Solve 1-D interfacial mass transfer using pdepe

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?

[Bill Greene]

I would be interested in trying to help you with this. Could you reformat your code so that it is more readable. And supply a runnable version with Cinf=constant; the version you have posted is missing values D, Cpo, etc.

Also, please send me a private message with a little background on this example.

Answer 1

[Torsten]

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.

[Edgar Ortiz]

I ve been thinkig about the second aproach for days but is not clear to me how to evaluate dCp/dx(@x=dP) at t_i

Thank you so much Torsten

[Torsten]

https://de.mathworks.com/help/matlab/ref/pdeval.html

Best wishes

Torsten.