# Modeling diffusion through two layers

by robby991
Tags: diffusion, layers, modeling
 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: 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