Modeling diffusion through two layers

[robby991]

Hi, I am looking to model the diffusion of a particle through 2 layers using the finite difference method in matlab. I already successfully have modeled this in a mono-layer configuration, however I am looking to code it through two layers.

I have two layers, one is 1 micron (1E-6) and the other is 10 nanometers (10E-9), for a total of 1.01E-6.

I am going to set this up in one mesh that has a total length of 1.01E-6 (for simplicity I will exclude boundary conditions because I know how to implement those). So I will first set up a matrix of all zeros with those dimensions.

The spacing between nodes is defined as:

x = linspace(0,Length,numx);
S = zeros(numt,numx);
dx = S(2)-S(1);

For the region of the first layer, the diffusion I previously described as a normal second order diffusion represented in the central difference scheme ("j" is the time domain and "i" space domain):

 for j=1:numt-1  
        
    for i=2:numx-1
      S(j+1,i) = S(j,i) + (dts/dxs^2)*Ds*(S(j,i+1) - 2*S(j,i) + S(j,i-1)) 
         end
    
    S(j+1,1)=S(j,1)+dts*Ds*2*(S(j,2)-S(j,1))./dxs.^2-((Vmax*dts*S(j,i))/(Km+S(j,i)));      %Neumann Boundary Condition
    
end

For the interface of layer k and k+1, the spatial derivative gives:

dP/dt = [Dk+1*Pk+1 - (Dk+1 + Dk)Pj + Dk*Uk-1]/dx^2

I am unsure how to code this in for loops because I need to transistion diffusion equations after 1E-6, and I don't know exactly how many space steps (dx) this corresponds to. For example, the first for loop for the first diffusion equation for the first layer will be, where the question mark is the point at the interface, then the interface diffusion equation above will be implemented:

 for i=2:[B]?[/B]      
S(j+1,i) = S(j,i) + (dts/dxs^2)*Ds*(S(j,i+1) - 2*S(j,i) + S(j,i-1)) 
         end

I can calculate it by hand but I know there is a way to do it by programming. Is this the correct approach to do it all in one mesh? Any help on this would be greatly appreciated. Thank you!

 

[kuan]

I am confused! It seem for me this problem can be split in 2 diffusion problems with the coupling at the boundary. For the stability, the simplest one is dt / dx^2 <=1. I cannot remember what it is called. You find it in the Numerical Recipes.

 

[robby991]

That is my question. Can I set up two separate matrix's and couple inbetween them with the equation I posted? Or does it have to all be in the same matrix?

 

[kuan]

Hmmm, I am not sure I understand you properly? Do you ask if these 2 systems are ought be solved at the same time or not? 2 systems are coupled, thus generally they need to be solved at the same time.

 

[robby991]

Yes they are coupled. I am asking if I can use 2 separate matrix's. One for the first layer and one for the second layer.

 

[kuan]

hmmm, i might see your question now. Can you set up 2 pseudo-grids at the boundary? One for the diffusion in domain 1 and the other one for domain 2.

 

[yus310]

If you are going to calculate and analyze the concentration in two layers, one 2D or 3D matrix will do. Make sure to get your equations right, and apply the Gauss Siedel with an appriopriate convergence criteria, where you initially set an initial vector for C(t,x,y) at t=0 is {0... so on).. right out your equations in a loop for the chosen number of iterations, then right out your criteria at the end of this loop, and then finally... wait until it converges. Make sure to take into account whether or not the interface between the two boundaries is perfect, meaning does D*dC1/dt=D*Dc2/dt at the interface or does C1=C2 at the interface. Include boundary conditions in the loop as well. You should get a good answer in no time using one matrix and the gauss siedel.

best

 

[robby991]

Thank you for the comments. I think I am going to try it forst with finite difference because I already have working code for the diffusion equation. I set up a matrix, space and time, and I have the following questions:

numx = 30;                      %number of grid points in space
numx_tot = 2*numx;
numt = 1000;                    %number of time steps to be iterated over 
tmax = .0045;                   %number of time steps to be iterated over 

Length1 = 1E-6;                 %length of grid
Length2 = 10E-9;
Length_tot = Length1+Length2;

S1 = zeros(numt,numx_tot);      %Substrate matrix, initialize to zero
xS = [0:Length1/numx:Length1 (Length1+Length2/numx):Length2/numx:Length1+Length2]; %Substrate vector of x values, to be used for plotting
zS = xS/Length_tot;
ts = linspace(0,tmax,numt)';    %vector of t values, to be used for plotting

Should the spacing on the space domain (xS) be the same for both layers? As you can see from my code the number of intervals is the same but the spacing is different. In addition, each layer has a different diffusion coefficient, which I scaled to its respective layer (Length1 or 2). I am confused how to implement this. Should the grid spacing be the same in two layers so the scaling of the diffusion coefficients will be with respect to the whole layer? I am concerned about doing this because as you can see the second layer is about 3 orders of magnitude smaller than the first.

 

[yus310]

It is much easier if the mesh is even.. hence delta x is the same. It doesn't matter how long or short the layers are, the respective diffusion coefficients, treated as known variables e.g. D1=20, D2=30, will take care of that.. you job is to set up the mesh, make the differential equation easier to solve via finite differences, input the known variables, and then apply a numerical technique to solve it.

Hope this helps.

 

[robby991]

Thanks. I am having difficulty generating a vector in space that has evenly spaced points for the 2 different lengths. For example, if I have length 1 as 1E-6 and length 2 as 10E-9 which are fixed constants, I would set up the vector as:

numx=30;
Length1 = 1E-6;
Length2 = 10E-9;
Length_tot=Length1+Length2;
x = [0:Length_tot/numx:Length1 (Length1+Length_tot/numx):Length_tot/numx:Length_tot];

This will not generate equally spaced points from 0 to 1E-6, then 1E-6 to 1.01E-6 UNLESS the spacing between the points is a common integer to both. However, when doing finite difference schemes the number of points in space will vary until you get a stable result, and the probability of the stable result having a "numx" value that gives correct equal spacing is very unlikely. How to alleviate the situation?

 

[yus310]

Take the d1+d2=d3, or add the two layers together set equal to d3 (total distance), then choose an odd integer of points, meaning n=101, then set dx=d3/(n-1), hence you get dx to be uniform! This makes meshing easy, e.g. if the d3=1, and n=3, then dx=1/(3-1)=1/2, which makes sense as from pt. 1 to pt. 2 d=1/2, from pt.2 to pt.3 d=1/2, hence d3=1/2+1/2=1!

 

[robby991]

Ok got it. I've implemented the for loops for iterations through layer1, across the interface, and through the second layer. The diffusion equation through the 1st layer is slightly different than the 2nd layer, so I was wondering exactly how to code the loops. The diffusion through the first layer is as follows:

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

Explicit central difference in Matlab:

S(j+1,i) = S(j,i) + (dts/dxS^2)*Ds1*(S(j,i+1) - 2*S(j,i) + S(j,i-1))-((Vmax*dts*S(j,i))/(Km+S(j,i)));

at the boundary is:

in Matlab:

S(j+1,i) = S(j,1) + (dts/dxS^2)*Ds2*(S(j,i+1) - (Ds2+Ds1)*S(j,i) + Ds1*S(j,i-1))

Through the 2nd layer:

S(j+1,i) = S(j,i) + (dts/dxS^2)*Ds1*(S(j,i+1) - 2*S(j,i) + S(j,i-1));

My specific question is how do I incorporate the term at the end of the first diffusion layer (VmaxS/Km+S) at the bounday. Right now, the finite difference at the boundary (node before and after) does not include that term.

Can someone take a look at my for loops and verify if they are correct? Specifically the loop for at the boundary.

S(1,:) = S0;                        
P(:,numx)=0;

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

%iterate central difference equation% 

for j=1:numt-1  
    
      %2nd Derivative Central Difference Substrate Iteration% 
     
    for i=2:k-1
      S(j+1,i) = S(j,i) + (dts/dxS^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/dxP^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         %where k equals node of Length1, boundary of 2 layers
        
        S(j+1,i) = S(j,1) + (dts/dxS^2)*Ds2*(S(j,i+1) - (Ds2+Ds1)*S(j,i) + Ds1*S(j,i-1))-((Vmax*dts*S(j,i))/(Km+S(j,i))); 
    end
    
    for i=k+1:numx-1
        S(j+1,i) = S(j,i) + (dts/dxS^2)*Ds1*(S(j,i+1) - 2*S(j,i) + S(j,i-1));
    end
    
    S(j+1,1)=S(j,1)+dts*Ds1*2*(S(j,2)-S(j,1))./dxS.^2-((Vmax*dts*S(j,i))/(Km+S(j,i)));      %Neumann Boundary Condition
    
end

 

[yus310]

you need to state the mathematical equations, variables, and BC's before asking this. I really cannot tell anything from the computer program.

 

[robby991]

The system is 2 layers, from x = 0 to x = Length1, from x = Length1 to Length_tot, where Length_tot is Length1 + Length2.

In the region x = 0 to x = Length1, the diffusion of species S is defined as the diffusion equation evaluated with second order finite difference scheme with diffusion coefficient Ds1 with a non-linear term subtracted from the end (the non-linear term can be evaluated explicitly like in the code, as confirmed by literature and working single layer code):

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

From the region Length1 to Length_tot, the species S follows the classic diffusion equation, with diffusion coefficient Ds2:

dS/dt = Ds2*d^2S/dx^2

At the boundary,k, between the 2, the finite difference scheme is defined as (see attached):

dS/dt = Ds2*S(k+1)-(Ds2+Ds1)*S(k)+Ds1*S(k-1)/dx^2

The boundary conditions are as follows:

Neumann boundary condition dS/dx=0 at x = 0, which is in the body of the loop.
Dirchelet boundary condition at x = Length_tot is some constant "S0".

I just need help in merging the 2 diffusion equations. The central difference scheme at the boundary k describes merging of the two layers, I am just unsure if I am implementing it correctly in the code.

 

[robby991]

Hi everyone. I am still having issues implementing diffusion through 2 layers in Matlab. I attached my Matlab code using the information in my previous post. If anyone has any suggestions I'd appreciate the feedback. Somehow, I do not think I am implementing the equations right in my code. There is two species, S and P that I am modeling. This code runs and gives me some output, however I was wondering if someone can let me know if my loop at the k boundary are ok.

%{ 
 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 = 51;                         %number of grid points in space
numt = 20000;                      %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 = 1E-6;                    %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

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

%******** 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-2
      S(j+1,i) = S(j,i) + (dts/zSx^2)*Ds1*(S(j,i+1) - 2*S(j,i) + S(j,i-1)); 
      P(j+1,i) = P(j,i) + (dtp/zPx^2)*Dp1*(P(j,i+1) - 2*P(j,i) + P(j,i-1));
    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+2: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]

Your program is confusing... I'm going to share with you my program on something very similar, but heat transfer (same equations though)

n=29; (Defining nodes)
dx=5/(n-1); (defining dx, which is Total length/(n-1)
Y=0:dx:5; Defining the length of Y OR Z
alpha=0.5; (constants)
stab=(0.5*dx*dx)/(alpha); (need stability criterion for 0.5*dx*dx/alpha>dt
dt=0.03; (Hence from stability we choose dt)
g=0:dt:250;
[o,d]=size(g);
th=zeros(d,n); (setting up our matrix of interest, t,x)
tau=(0.5*dt)/(dx*dx);
for k=1
for i=1:1:n
th(k,i)=1;
end
end

for k=2:1:d
th(k,1)=0;
th(k,n)=2;
end
(THESE ALL ARE BOUNDARY CONDITIONS & INITIAL CONDITiONS

for t=1:1:d-1
for i=2:1:n-1
th(t+1,i)=((th(t,i+1)+th(t,i-1)-(2*th(t,i)))*(tau))+th(t,i);
end
end

(EVALUATING INSIDE THE WIRE! A LOOP!)

tc=10; %Choose t for Viewing, Steady at tc=1E5 or 5 s
for i=1:1:n
thetha(i)=th(tc,i);
end
hold on
figure (1)
plot(Y,thetha)
xlabel ('x (cm)')
ylabel('thetha(t,x)')

Point is
a) set up a mesh (in your case it is like a wire, two layers)
b) set up a matrix of interest, C(t,x)
c) determine the stability for your formulation to work in this case (your case)
dt<(coefficient for difussion*0.5/(dx*dx)
D) hence determine dt, determine your t length (0:dt:t)
E) Set up your initial conditions C(0,X)=?
F) Set up Boundary Conditions if Known (Keep outside the loop)
G) Input Unknown BC's & Inner Node Calculations inside mesh using finite differences, thus calculate these in a loop that changes with increasing t
H) Finally plot your results

for C,D,E,F,G you may have different results... so point is for example you have a perfect interface, at the two layers so D1*dc1/dx=D2*dc2/dx, so you can use finite differences to set this up in the loop to calculate concurrently as your other inner nodes change, etc.

Get back to me if you need more help.

 

[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.

 

[yus310]

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:

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:

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

And at the boundary of the 2 layers:

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:

% 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:

%{ 
 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.

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

Yes, the interface node is:

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:

 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