Using PDEPE for Ficks second law ?

6 次查看(过去 30 天)
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.
  2 个评论
Bill Greene
Bill Greene 2021-6-11
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
Stefan 2021-6-14
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

请先登录,再进行评论。

采纳的回答

Stefan
Stefan 2021-6-14
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

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Manage Products 的更多信息

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by