modeling diffusion through two layers


by robby991
Tags: diffusion, layers, modeling
robby991
robby991 is offline
#19
Mar29-11, 09:00 AM
P: 41
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
yus310 is offline
#20
Mar29-11, 09:45 AM
P: 81
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
robby991 is offline
#21
Mar29-11, 10:45 AM
P: 41
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
yus310 is offline
#22
Mar29-11, 11:06 AM
P: 81
send me the code, repost again with equations, trust me I will solve it using my method in matlab
robby991
robby991 is offline
#23
Mar29-11, 02:17 PM
P: 41
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
robby991 is offline
#24
Mar30-11, 09:31 PM
P: 41
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
yus310 is offline
#25
Apr19-11, 10:46 AM
P: 81
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
robby991 is offline
#26
Apr19-11, 12:13 PM
P: 41
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


Register to reply

Related Discussions
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