aliafrouzian / implement-boundary-condition-and-diffusion

Hello Guys, I was wondering if you could take a look at my code and give me your feedback if it is correct or not. I have some doubts about the results. Also, I don't know how to implement boundary condition and diffusion equations in my code. The algorithm is attached. Thank you so much.
0 stars 0 forks source link

How to implement boundary condition and diffusion? #1

Open aliafrouzian opened 4 years ago

aliafrouzian commented 4 years ago

Hello Guys, I was wondering if you could take a look at my code and give me your feedback if it is correct or not. I have some doubts about the results. Also, I don't know how to implement boundary condition and diffusion equations in my code. The algorithm is attached. trying to upload initial draft of code as well. Thank you so much. Algorithm.pdf

aliafrouzian commented 4 years ago

n=100; % number of iterations Dp=30e-6; % Pillar diameter [micrometer] h=Dp; % Thickness is equal to diameter of pillar dx=h/n; % Discritization (delta x) alpha=1.6563e-4; % Thermal diffusity of silver [m^2/s] Phi_o=0.2; % Initial porosity at time t=0 Phi_st=0.12; % Critical porosity (should be lower than Phi_o) t_final=22560; % total time (sec) t1=0.1t_final; % time that there is no temperature gradient dt=10; % time increment [sec] R=8.314; % Universal gas constant [J/mol.K] G=8300; % Activation energy [J/mol] B= 118e9; % Effective Bulk Modulus [Pa] To=298.15; % Initial temperature [K] T_f=523.15; % Final temperature (heating process) mu_ost=0.047; % mu_o = mu_ostexp(-G/RTo) etha=0.1; % Flow mobility delta_T=1; % Temperature change [celecius/min] DPDF=0.01; % In Eq. 13 AA=0; % Considering phi(y) %---------------------------------------------------------------------- % OUTPUT %---------------------------------------------------------------------- % Mu_o & Mu_pr mu_o = mu_ost exp(-G/(RTo)); mu_pr = mu_o G/(R*To^2);

% Initializing Phi & Theta_s

t=0:dt:t_final; y=dx/2:dx:h-dx/2; for jj=1:length(y) for ii=1:length(t) phi=ones(jj,ii); theta_s=zeros(jj,ii); T=zeros(jj,ii); end end phi(:,1)=Phi_o; phi(:,2:length(t))=0; % Entering Loops % Defining Temperature

for i=1:length(y) for j=1:length(t) if t(j)<t1 T(i,j)=To-(y(i)/h)delta_T(1-t(j)/t1); else T(i,j)=To; end end end

% Defining Porosity

for i=1:length(y) for j=1:length(t) phi(i,j+1)=dt(-mu_oexp(-G/(RT(i,j))DPDF))+ phi(i,j); end end

% Defining Theta_f Kapa=-ethaalphaB*h^2;

% Without considering second derivative of phi(y) beta=-Kapaalpha; theta_f=zeros(length(y)+1,length(t)); for i=2:length(y) for j=1:length(t)-1 if phi(i-1,j)>Phi_st theta_f(i-1,j+1)=theta_f(i-1,j)+ beta(theta_f(i+1,j)-... 2theta_f(i,j)+ theta_f(i-1,j))dt/(dx^2); else theta_f(i,j+1)=theta_f(i,j); end end end