# modeling diffusion through two layers

by robby991
Tags: diffusion, layers, modeling
 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?
 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
 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.
 P: 81 send me the code, repost again with equations, trust me I will solve it using my method in matlab
 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, %}
 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');
 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
 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

 Related Discussions Nuclear Engineering 15 Advanced Physics Homework 2 Introductory Physics Homework 1 General Astronomy 8 General Engineering 0