modeling diffusion through two layers

yus310

Quote by robby991

Thank you for the reply, I went through your code and it gave me different ideas on how to call initial and boundary conditions. This was only considered as a 1 layer system correct? One problem I cannot figure out, maybe you can help me, is how to vectorize this code (mine and yours) Right now we are both initializing a matrix where the computations are performed, while all the data is being stored inside that matrix. I was wondering how to change the code so as to use vectors, so only the current and previous vectors are saved, then at the end plot the final vector. I am asking this because I am running out memory in my simulation because the matrix I have to initialize is so large for the algorithm to be stable. I feel it is a simple code adjustment but I cannot seem to get it working.

1) Set up your mesh, and number of nodes (n), find dx via total length/(n-1) and make a X matrix from 0:dx:total length
2) Set up your dt, via stability criterion which says dt<(dx*dx*0.5)/(coefficient for diffusion whatever it is)
3) set up a time iteration 0:dt:final time of interest
4) FInally make a matrix of C(x,t) where C is a matrix of zeroes that has dimensions D X N or where D is the number of time slots in the time iteration and N is the number of nodes or the number of x slots in the x-iteration'
5) Input your known BC's (constants, not dervivatives), and initial conditon, C(x,0)=?
6) Set up a time-marching loop, which says from t=2:1:D-1 (going from the time at 2 for the time iteration to time of D-1... to calculate along the layer)
7) Input your finite difference equations for the inner nodes (not the beginning and end.. two different inner node equations because of two layers)
8) At the new interface where the two layers meet set up the equation of D1*dc1/dx=D2*dc2/dx, where 1 is layer 1 and D2 is layer 2... this means that concentration diffusion through the common node of 1 and 2 is perfect,
9) Set up your other bc's if you beginning is insulated or no concentration coming from right (then use backwards difference)
10) Each time you calculate you need to calculate C(t+1,x)=C(t,x).... whatever finite difference, which means you store the next concentration at the next time using the previous times and previous values
11) Done! It should be show appropriate behavior then....

I'm attaching my heat transfer program it is exactly the same as yours except similar again just compare the two I sent you

%Spencer Nelle
%1D Transient Conduction - Simulation of Platinum Hot Wire in Molten Salt

clc
clear

%Constants:
k_plat=71.6; %W/mK Good Enough Approximation k(T)? But Then keeps in 300 Range!
k_salt=0.511; %W/mK

R=0.002699;
I=10; %Need to Be Measured
V=pi*(0.001^2)*(0.25)*(0.02);
L=0.02;
Qg=(I*I*R)/V;

Rho_plat=21450; %kg/m^3
Rho_salt=1928.1; %kg/m^3

Cp_plat=130; %J/kg K
Cp_salt=1290; %J/kg K

K=k_plat/k_salt;

alpha_plat=(k_plat/(Rho_plat*Cp_plat));
alpha_salt=(k_salt/(Rho_salt*Cp_salt));

%Establish Matrix/Nodes/Intervals
n=11; %(Number of Nodes)
m=(n-1); %(Number of Intervals)
I=(n+1)/2;

r1=.1275*10^-3; %m (Radius of Platinum Wire)
r2=20*r1; %m (radius of Salt)
dr=(r2+r1)/m; %(Length of interval along 1D radius)
r=0:dr:(r2+r1);
T0=300; %C

%Time Step
dt=0.00005;
time_set=0:dt:10; %seconds
[o,d]=size(time_set);

%Check Fourier Number for Stability
Fo1 = (alpha_plat*dt)/(dr*dr);
Fo2= (alpha_salt*dt)/(dr*dr);
%Initial Conditions
T=zeros(d,n);

%Initial Condition
for i=1:1:n
T(1,i)=T0;
end

%Calculation of Fourier Number to Find Stability Criterion is observed to see that it is met.
%Inputting the initial conditions (boundary + initial)
if Fo1 <= 1/2 && Fo2<= 1/2
for t=1:1:d-1
if t<=4001 %Time Steps corresponding to 2 seconds of pulse
for i=2:1:I-1 %Region I
T(t+1,i)=(((T(t,i)*((1/(dr*(dr*(i-1))))-(2/(dr*dr)))) + (T(t,i-1)*((-1/((dr*(dr*(i-1))))) + (1/(dr*dr)))) + (T(t,i+1)*(1/(dr*dr))) + (Qg/k_plat))*((alpha_plat*dt)))+T(t,i);
end

for i=I+1:1:n-1 %Region II
T(t+1,i)=(((T(t,i)*((1/(dr*(dr*(i-1))))-(2/(dr*dr)))) + (T(t,i-1)*((-1/((dr*(dr*(i-1))))) + (1/(dr*dr)))) + (T(t,i+1)*(1/(dr*dr))))*(alpha_plat*dt))+T(t,i);
end

%Left Node Boundary Condition Changes with Time with Pulse
T(t+1,1) = ((-1*Qg*2*dr*L)/(3*k_plat)) + (T(t+1,2)*(4/3)) - (T(t+1,3)*(1/3));

else %Time Steps corresponding to after pulse
for i=2:1:I-1 %Region I
T(t+1,i)=(((T(t,i)*((1/(dr*(dr*(i-1))))-(2/(dr*dr)))) + (T(t,i-1)*((-1/((dr*(dr*(i-1))))) + (1/(dr*dr)))) + (T(t,i+1)*(1/(dr*dr))))*((alpha_plat*dt)))+T(t,i);
end

for i=I+1:1:n-1 %Region II
T(t+1,i)=(((T(t,i)*((1/(dr*(dr*(i-1))))-(2/(dr*dr)))) + (T(t,i-1)*((-1/((dr*(dr*(i-1))))) + (1/(dr*dr)))) + (T(t,i+1)*(1/(dr*dr))))*((alpha_plat*dt)))+T(t,i);
end

%Left Node Boundary Condition Changes with Time NO Pulse
T(t+1,1) = (T(t+1,2)*(4/3)) - (T(t+1,3)*(1/3));
end

%Right Node Boundary Condition
T(t+1,n)= (1/3)*(4*T(t+1,n-1) - (T(t+1,n-2)));

%Interface
T(t+1,I)= (-1/((3*K)+3))*((K*T(t+1,I-2)) - 4*K*T(t+1,I-1) -4*T(t+1,I+1) + T(t+1,I+2));
end
else
('Please Input a New Time Step for this Program ')
end

tc=5E4; %Choose t for Viewing, Steady at tc=1E5 or 5 s
for i=1:1:n
TT(i)=T(tc,i);
end

plot (r,TT)
timeofinterest=time_set(tc)

robby991

Yus: Thanks this is great. I was going through the code and there are some things that are a little confusing to me. I understand the general heat equations and how you are calling the boundary conditions, however could you clearly state the problem at hand (purpose of pulse?) and boundary conditions?

Also, what I am trying to say is that we both initialize a matrix where we perform computations, in your case T=zeros(d,n); The data from the time and space stepping computations gets saved in this matrix after the program runs. This is bad for computations since the memory could get full (happened in my case). A better way to code this would be to use vectors, using the previous vector to calculate the new vector, then the new vector to calculate the next vector, and not saving the used vectors. So at the end of the stepping only 2 vectors remain. Do you understand?

yus310

you can do as you please... gauss siedel, time marching, either way one matrix is fine enough.. you should not run out of memory.. it means you are doing the problem wrong

robby991

Unfortunately there is no other option other than vectorizing. To keep the stability in my code, the dt has to be extremely small (total time is fixed), which makes the number of time steps in my matrix extremely large (on the order of 500000). Matlab either gives me a memory error or crashes when it runs. It runs when I reduce the number of time steps and increase x, however I cannot do this since my x has to be a fixed value. Using single vectors we can save the same amount of memory as a matrix in a vector, then keep overwriting the vector with new values.

yus310

send me the code, repost again with equations, trust me I will solve it using my method in matlab

robby991

The system is as follows.

Diffusion of 2 species, S and P, through 2 layers. First layer is 1E-6 second is 10E-9. The diffusion through layer 1 for S and P is:

Code:

dS/dt = Ds1*d^2S/dx^2 - (vmaxS/km +S)
dP/dt = Dp1*d^2P/dx^2 + (vmaxS/km +S)

Diffusion through the second layer is the standard diffusion equation:

Code:

dS/dt = Ds2*d^2S/dx^2 
dP/dt = Dp2*d^2P/dx^2

And at the boundary of the 2 layers:

Code:

Ds1*dS1/dx = Ds2*dS2/dx
Dp1*dP1/dx = Dp2*dP2/dx

Where the nonlinear part of layer 1 needs to be integrated into the interface diffusion which I have not done yet.

BC and ICs are as follows:

Code:

% Dirchelet Boundary Conditions, where Neumann BC dS/dx=0 is in body of for loop
S(:,numx) = S0;          %S(x=d,t)=S0
P(:,1) = 0;              %P(x=0,t)=0
P(:,numx) = 0;           %P(x=d,t)=0

%{
Initial conditions satisfied with martix initialization and BCs   
 
S(0,x) = 0, 0 <= x < d,
S(0,d) = S0,
P(0,x) = 0, 0 <=x < d,
%}

robby991

Here is my code:

Code:

%{ 
 Solution to the substrate and product diffusion equations using the 
central difference scheme with Neumann boundary conditions for S at x = 0, 
and dirchelet boundary conditions:
S(:,numx) = S0;          %S(x=d,t)=S0
P(:,1) = 0;              %P(x=0,t)=0
P(:,numx) = 0;           %P(x=d,t)=0

Initial conditions satisfied with martix initialization and BCs    
S(0,x) = 0, 0 <= x < d,
S(0,d) = S0,
P(0,x) = 0, 0 <=x < d,
%}

clear all;

numx = 201;                         %number of grid points in space
numt = 1000000;                      %number of time steps to be iterated over 

tmax = .1;
ts = linspace(0,tmax,numt)';       %vector of t values, to be used for plotting
dts = ts(2)-ts(1);                 %Define grid spacing in time
tp = linspace(0,tmax,numt)';       %vector of t values, to be used for plotting
dtp = tp(2)-tp(1);                 %Define grid spacing in time

Length1 = 1E-6;                    %length of layer 1
Length2 = 10E-9;                    %length of layer 2
Length_tot = Length1+Length2;      %total length of layer

S = zeros(numt,numx);              %Substrate matrix, initialize to zero     
dxS = Length_tot/(numx-1);         %Define grid spacing in space

xS = 0:dxS:Length1:dxS:Length_tot; %Substrate vector of x values
zS = xS/Length1;                   %scaled vector
zSx = zS(2)-zS(1);                 %scaled grid spacing
k = round(Length1/dxS+1);          %vector at the interface of the 2 layers

P = zeros(numt,numx);              %Product matrix, initialize everything to zero
dxP = Length_tot/(numx-1);         %Define grid spacing in space
xP = 0:dxP:Length1:dxP:Length_tot; %Substrate vector of x values
zP = xP/Length1;                   %Scaled vector
zPx = zP(2)-zP(1);                 %scaled grid spacing

S0=(1E-3);                         %mol/L initial substrate concentration

%********* 1st Layer ***********&
D1 = .019E-9;                      
Ds1 = D1/(Length1^2);              %scale the substrate diffustion coefficient

S_Stability1 = (Ds1*dts)/((zSx)^2) %check requirement Ds(dt)/dx^2 < .5, cm^2/sec
D2 = .1E-9;                         
Dp1 = D2/(Length1^2);              %scale the product diffusion coefficient                 

P_Stability1 = (Dp1*dtp)/((zPx)^2) %check requirement Dp(dt)/dx^2 < .5

V = 275E-6;                     %mol/cm^2sec
Vmax = V/Length_tot;                %Scale the velocity

K = 3E-3;                       %mol/cm^3
Km = K/Length_tot;                  %Scale the Michaelis Menten constant    
%**************************************%

%******** 2nd Layer********************%
D3 = .019E-9;                       %requirement Ds(dt)/dx^2 < .5, cm^2/sec
Ds2 = D3/(Length2^2);               %scale the substrate diffustion coefficient

D4 = .1E-9;                         %requirement Dp(dt)/dx^2 < .5
Dp2 = D4/(Length2^2);               %scale the product diffusion equation                 

%**************************************%

% Dirchelet Boundary Conditions, where Neumann BC dS/dx=0 is in body of for loop
S(:,numx) = S0;          %S(x=d,t)=S0
P(:,1) = 0;              %P(x=0,t)=0
P(:,numx) = 0;           %P(x=d,t)=0

%{
Initial conditions satisfied with martix initialization and BCs    
S(0,x) = 0, 0 <= x < d,
S(0,d) = S0,
P(0,x) = 0, 0 <=x < d,
%}

%2nd Derivative Central Difference Iteration% 
for j=1:numt-1  
        
    for i=2:k-1
      S(j+1,i) = S(j,i) + (dts/zSx^2)*Ds1*(S(j,i+1) - 2*S(j,i) + S(j,i-1))-((Vmax*dts*S(j,i))/(Km+S(j,i))); 
      P(j+1,i) = P(j,i) + (dtp/zPx^2)*Dp1*(P(j,i+1) - 2*P(j,i) + P(j,i-1))+((Vmax*dtp*S(j+1,i))/(Km+S(j+1,i)));
    end
    
    for i=k-1:k+1
      S(j+1,i) = S(j,i) + (dts/zSx^2)*(Ds2*S(j,i+1) - (Ds2+Ds1)*S(j,i) + Ds1*S(j,i-1)); 
      P(j+1,i) = P(j,i) + (dtp/zPx^2)*(Ds2*S(j,i+1) - (Dp2+Dp1)*P(j,i) + Dp1*P(j,i-1));  
    end
    
    for i=k+1:numx-1
      S(j+1,i) = S(j,i) + (dts/zSx^2)*Ds2*(S(j,i+1) - 2*S(j,i) + S(j,i-1));
      P(j+1,i) = P(j,i) + (dtp/zPx^2)*Dp2*(P(j,i+1) - 2*P(j,i) + P(j,i-1));
    end
    
      S(j+1,1)=S(j,1)+dts*Ds1*2*(S(j,2)-S(j,1))./zSx.^2;      %Neumann Boundary Condition
    
end

figure(1);
plot(zP,P(numt,:));
xlabel('Length');
ylabel('[P]');

figure(2);
plot(zS,S(numt,:));
xlabel('Length');
ylabel('[S]');

figure(3);   
i = 2*(96500)*Dp1*(P(:,2)/zPx);
plot(tp,i);
xlabel('Time');
ylabel('Current');

yus310

what is the equation for the last node? Backward difference, flux of concentration at the end?

Where is the interface node between the two layers.. the equation.. again I said assume

D1*dCA1/dx=D2*dCA2/dx

you need to add a nodal equation for the end and the interface

robby991

There is no flux at the end, it is a dirchelet boundary condition.

Code:

S(:,numx) = S0;          %S(x=d,t)=S0
P(:,numx) = 0;           %P(x=d,t)=0

Yes, the interface node is:

Code:

Ds1*dS1/dx = Ds2*dS2/dx
Dp1*dP1/dx = Dp2*dP2/dx

which I have evaluated in matlab as the following, where k is the interface node:

Code:

 for i=k-1:k+1
      S(j+1,i) = S(j,i) + (dts/zSx^2)*(Ds2*S(j,i+1) - (Ds2+Ds1)*S(j,i) + Ds1*S(j,i-1)); 
      P(j+1,i) = P(j,i) + (dtp/zPx^2)*(Ds2*S(j,i+1) - (Dp2+Dp1)*P(j,i) + Dp1*P(j,i-1));  
    end

Similar discussions for: modeling diffusion through two layers
Thread	Forum	Replies
Diffusion equation and neutron diffusion theory	Nuclear Engineering	15
Conductivity of a set of layers	Advanced Physics Homework	2
atmosphere layers	Introductory Physics Homework	1
Dynamics of the Sun's layers.	General Astronomy	8
COMSOL- Unsteady 1D Diffusion modeling	General Engineering	0

yus310	you can do as you please... gauss siedel, time marching, either way one matrix is fine enough.. you should not run out of memory.. it means you are doing the problem wrong