Using PDEPE for Ficks second law ?

Hello together,

i tried to use the PDEPE-Solver for Fick's second law of diffusion:

But I realy got confused with the initial and boundary conditions.

Is it possible to set the initial condition as c(x>0,0)=0

and the two boundaries as c(x=0,t)=10 und dc(x=l,t)/dt=0 ?

clc
clear
global dcoeff flux
L = 0.01; % m length of slab
tend = 36000; % s total of 10 hours
dcoeff = 10^-10; % m^2/s
flux = 10; % surface concentration wt%
m = 0;
x = linspace(0,L,20);
t = linspace(0,tend,20);
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);
u = sol(:,:,1);
surf(x,t,u)
xlabel('distance')
ylabel('time')
figure
ptot(x,u(:,:))
xlabel('distance')
ylabel('wt%')
figure
plot(t,u(:,:))
xlabel('time')
ylabel('wt%')
function [c,f,s] = pdefun(x,t,u,DuDx)
global dcoeff
c = 1/dcoeff;
f = DuDx;
s = 0;
end
function u0 = icfun(x)
if x > 0 
u0 = 0;
else 
   u0 = 10, 
end
end
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
pl = ul-10;
ql = 0;
pr = 0;
qr = 1;
end

Thank you for every hint.

[Bill Greene]

You say you want a BC . Is this a typo? This is a strange BC that would have no other effect than to set the value of c there equal to the initial condition.

[Stefan]

Hello Bill,

thank you for your comment! To be honest, I tried to reproduce following example (Fick's 2. Law) as an easy exercise for the beginning. I thought the BC: would set a constant c if x >= l ?

The biggest issue at the moment is that the If-Loop as IC isn't working right, it always sets u0=10 for all x at t=0 but i would like to set it as described in the example described above.

(To avoid confusio u in the code corresponds to c above)

function u0 = icfun(x)
if x > 0 
u0 = 0;
else 
   u0 = 10, 
end

Thank you very much in advance

Answer 1

[Stefan]

Ok I could solve the problem by myself!

I just mixed the IC and BC up! The If-Loop wasn't necessary.

I added the Code if someone is face the same topic!

Regards, Stefan

clc
clear
global dcoeff flux
L = 10; % m length of slab
tend = 50; % s total of 10 hours
dcoeff = 1; % m^2/s
flux = 10; % surface concentration wt%
m = 0;
x = linspace(0,L,100);
t = linspace(0,tend,100);
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);
u = sol(:,:,1);
surf(x,t,u)
xlabel('distance')
ylabel('time')

figure
contourf(x,t,u,[1:9],'ShowText','on')
ax = gca;
ax.YDir = 'reverse';
xlabel('x')
ylabel('time')


function [c,f,s] = pdefun(x,t,u,DuDx)
global dcoeff
c = 1/dcoeff;
f = DuDx;
s = 0;
end
function u0 = icfun(x)

   u0 = 0; 

end
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
pl = ul-10;
ql = 0;
pr = 0;
qr = 1;
end