pdepe function maximum value

Hello,

I want to solve Fick's second law of diffusion using the pdepe function, where w is the concentration of a species in a solid and D is the diffusion coefficient

• ∂w/∂t=∂/∂x(D ∂w/∂x)

I defined the following parameters

• c=1 m=0 s=0 and f=D ∂w(x,t)/∂x

I made a function to define D

My boundary conditions are (at x=0 and x=L)

• -D ∂w(0,t)/∂x=beta
• -D ∂w(L,t)/∂x=0

Since boundaries are defined as p+qf=0 I've set

• pl=beta
• ql=1
• pr=0
• qr=1

My code seems to work just fine. The concentration w at point x=0 grows with time and numerically, the way I’ve written the code, it has no limit but physically it cannot surpass a saturation concentration Cp. Is there a way to set a maximum to the function? I do not know how to translate this to the code or even where I can add this condition.I would like a code that sets w(0,t)=Cp whenever, w(0,t) ≥Cp

Thanks

[Torsten]

Usual boundary condition at x=0:

D*dw/dx = k*(w-Cp)

Best wishes

Torsten.

[FM]

I've seen the boundary condition written that way but my problem is that I have no way of knowing parameter k, that's why I am using the boundary condition where the diffusion flow D*dw/dx equals the flow (beta) that arrives at the solid at x=0,

Thank you

[Torsten]

Then you could make several computations for different k-values ...

[FM]

So with my current boundary condition there is no way of setting a maximum? The only way is to use the boundary condition you suggested?

[Torsten]

ql = 1.0;
if ul <= Cp
  pl = beta;
else
  pl = 0.0;
end

[FM]

I did try that but the concentration at x=0 does not stay constant once it reaches Cp, instead it starts decreasing with time and even becomes negative

[Torsten]

if ul <= Cp
  ql = 1.0;
  pl = beta;
else
  ql = 0.0;
  pl = ul - Cp;
end

If this does not work, stop the solver when ul reaches Cp and restart it with the new boundary condition.

Best wishes

Torsten.

[FM]

Thank you so much Torsten! I stopped the solver and then restarted it with a new boundary condition as you suggested and it works.

However I had to make a whole other function for the new boundary condition. Is there a way to put a condition on the boundary condition? Like if a chioce=1 then use the first boundary condition and if choice=2 then use the second boundary condition? If I try to put choice as an input for my function containing the boundary condition it says there are too many input arguments.

[Torsten]

sol = pdepe(m,pdefun,icfun,@(xl,ul,xr,ur,t)bcfun(xl,ul,xr,ur,t,ichoice),xmesh,tspan)
...

[pl, ql, pr, qr] = bcfun(xl,ul,xr,ur,t,ichoice)
...

[FM]

Thank you for your help! i really appreciate it

[FM]

Sorry to bother you again Torsten but I have a problem when I try to launch pdepe twice in a same matlab file. I launch it once with the first boundary conditions and then when the concentration reaches Cp I launch it again but with the new boundary conditions. However it stops at the first launch. If I run my code directly on the command window it works fine, (it runs pdepe twice) so I don't think there is a problem with my code. Do you know why this happens?

Thank you in advance

[Torsten]

Please show your code.