radial diffusion pde boundary conditions

27 views (last 30 days)
I am writing a spherical diffusion pde and looking at examples for help. I am wondering what some terms mean.
For example, i think icfun is supposed to identitify initial conditions (concentrations at time t=0)? is this correct?
and bcfun i think are the spatial boundary conditions. but pdepe solves in 1-d so what are the four terms pl, ql, pr, and qr?? do these terms change for spherical problems, when m=2 for pdepe?
what is pdefun on the pdepe input values?
any help is appreciated. thank you

Accepted Answer

Torsten
Torsten on 13 Jul 2022
Edited: Torsten on 13 Jul 2022
Here is an example for the 1d spherical diffusion equation
dT/dt = 1/r^2 * d/dr (r^2*D*dT/dr)
with boundary conditions
dT/dr = 0 at r = 0
T = 275 at r = 1
and initial condition
T = 0 for 0 <= r < 1
T = 275 for r = 1.
integrated for 1000 time units.
r = linspace(0,1,25);
t = linspace(0,1000,25);
m = 2;
sol = pdepe(m,@heatsphere,@heatic,@heatbc,r,t);
u = sol(:,:,1);
surf(r,t,u)
xlabel('r')
ylabel('t')
zlabel('u(r,t)')
view([150 25])
plot(t,sol(:,1))
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at center of sphere')
function [c,f,s] = heatsphere(r,t,u,dudr)
D = 1e-4;
c = 1;
f = D*dudr;
s = 0;
end
%----------------------------------------------
function u0 = heatic(r)
u0 = 0;
if r==1
u0 = 275;
end
end
%----------------------------------------------
function [pl,ql,pr,qr] = heatbc(rl,ul,rr,ur,t)
pl = 0.0;
ql = 1.0;
pr = ur - 275;
qr = 0;
end
%----------------------------------------------
  3 Comments
Torsten
Torsten on 13 Jul 2022
Edited: Torsten on 13 Jul 2022
"l" refers to the left boundary point, "r" refers to the right boundary point.
ul is T at r=0, ur is T at r=1.
The general boundary condition is
p + q*f = 0.
In the above case, f = D*dT/dr, thus
p + q*D*dT/dr = 0.
Left boundary point (r=0) with boundary condition dT/dr = 0 gives
pl = 0, ql = 1, e. g., because this leads to 0 + 1*D*dT/dr = 0 or dT/dr = 0.
Right boundary point (r=1) with boundary condition T = 275 gives
pr = ur - 275, qr = 0 because this leads to (T - 275) + 0*D*dT/dr = 0 or T = 275.
Juliana
Juliana on 13 Jul 2022
wow okay that makes sense. i cannot thank you enough!!

Sign in to comment.

More Answers (0)

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!