MATLAB: Solve 1-D interfacial mass transfer using pdepe

mass transferpdepe

Hi Everyone,
I am trying to use pdepe to solve a diffusion problem and Im having issues trying to set my left side boundary condition.
 dP=0.04; %thickness polymer layer [cm]
 d=100; %times in days [d]
 tt=d*86400; %time in seconds [s]
 x = linspace(0,dP,5); 
 t = linspace(0,tt,100); 
 m = 0;
 sol = pdepe(m,@equa,@IC,@BC,x,t);
 u =sol(:,:,1) 
 function [c,f,s,algo] = equa(x,t,u,DuDx)
 ∂Cp/∂t=D*∂/∂x(∂Cp/∂x)
 c = 1;
 f = D*DuDx;
 s = 0;
 end  
 function u0 = IC(x)
 u0 = Cpo;  %Initial  concentration[microg/cm^3] 
 end
 function [pl,ql,pr,qr] = BC(xl,ul,xr,ur,t)
 pl =0;
 ql =1;
 pr =h*((ur/K)-Cinf);
 qr =D;
 end 
  the first boundary condition (0,t)
  ∂Cp/∂x=0 
  the second boundary condition (x,t)
  h*(Cp/K-Cinf)+D*∂Cp/∂x=0
Rigth now looks like it is working using Cinf as a constant but actually Cinf should change and increase with time (accumulation).
Cinf=A/V*integral(∂Cp/∂x (t) dt, 0,t)
I really dont know how to solve it this way
 Could someone please guide me?

Best Answer

First of all one error in your actual code:
qr=1 instead of qr=D.
If you set qr=D, your boundary condition reads
h*(Cp/K-Cinf)+D^2*∂Cp/∂x=0
Concerning your question about Cinf there are two ways to solve this problem:
1. Discretize your PDE equation in space, add the equation
dCinf/dt = A/V*dCp/dx(@x=dP)
to the system and solve all the equations using ODE15S. Look up "method-of-lines" for more details.
2. Iterative procedure:
Fix a list of constant output times t_i.
a) Solve your equation for a constant Cinf.
b) From the solution, evaluate dCp/dx(@x=dP) at the output times t_i and use "cumtrapz" to build an approximation to Cinf(t_i)=A/V*integral_{t=0}^{t=t_i} dCp/dx(@x=dP) dt.
c) Now hand these two vectors (t_i and C_inf(t_i)) to BC and use "interp1" to interpolate Cinf at the time t requested by the solver.
d) Compare with the first solution. If error<eps, stop, else go to step b)
Best wishes
Torsten.
Related Question