Pdepe cylindrical coordinates with Neumann BC

34 views (last 30 days)
Ali
Ali on 1 Mar 2024
Edited: Torsten on 28 Mar 2024 at 11:45
I try to solve 1D diffusion eqn with Neumann(flux) BC ar r=R_c. I couldn't use variable Conc in the BC function. It gives 'Unrecognized function or variable' error. I appreciate for any ideas on how to resolve.
Example input: AdvDiff_Cyl(4E-11,4.2E-10,0.4,[0 6; 1*3600 5; 6*3600 3.75; 20*3600 1.85; 48*3600 0.7; 95*3600 0.3]);
function Conc = AdvDiff_Cyl(Deff,P,Kav,NP1)
% Define the parameters
R_C = 19.2*10^-6/2/pi; % radius of microvessel [m]
R_O = 10*R_C; % radius of surrounding outer tissue [m]
tsim = 24 * 3600 % simulation time in [s]
r = [R_C:1E-6:R_O]; % create spatial and temporal solution mesh
t = [0:5:tsim];
% Solve the PDE numerically using pdepe which solves the pde using ode15s
m = 1 ; %1D cylindrical model
[sol] = pdepe(m, @advdiff_pde, @advdiff_initial, @advdiff_bc,r,t);
Conc = sol(:,:,1);
colors = [ 0 0.4470 0.7410;
0.8500 0.3250 0.0980;
0.9290 0.6940 0.1250;
0.4940 0.1840 0.5560;
0.4660 0.6740 0.1880;
0.3010 0.7450 0.9330;
0.6350 0.0780 0.1840];
figure()
i = 1;
for time = [5/3600 1 2 5 10 20]*3600/5
plot(r/R_C,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2)
hold on
i = i + 1;
end
xlabel('Distance [r/R_N]', 'Interpreter','Tex');
ylabel('Concentration [\mug/ml]','Interpreter','Tex');
title('Numerical Solution of ADE using pdepe','Interpreter','Tex');
legend('Initial Condition', '1 hr', '2 hr','5 hr','10 hr','20 hr','Location','northeast')
lgd.FontSize = 6;
legend('boxoff')
% Define the DE
function [c,f,s] = advdiff_pde(r,t,Conc,dcdr)
c = 1;
f = Deff*dcdr;
s = 0
end
% Define the initial condition
function c0 = advdiff_initial(r)
c0 = 0;
end
% Define the boundary condition
function [pl,ql,pr,qr] = advdiff_bc(xl,cl,xr,cr,t)
% BC1:
pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (Conc/Kav));
ql = Deff; %ignored by solver since m=1
% BC2:
pr = 0;
qr = 1;
end
end

Accepted Answer

Torsten
Torsten on 1 Mar 2024
Edited: Torsten on 1 Mar 2024
You will have to use
pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (cl/Kav));
ql = 1; %ignored by solver since m=1
instead of
pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (Conc/Kav));
ql = Deff; %ignored by solver since m=1
  8 Comments
Ali
Ali on 28 Mar 2024 at 2:50
I should clarify my point furher. What I try to model is cyclic concentration BC (1+sin(2*pi/24/3600*t)/2 at r=R_c which is in between 0 and 1, diffusing into the domain r>R_c.
That's why physically I expect sinusoidal concentration BC in between 0 and 1 for (0<r<R_c) to be greater or equal to the concentration values inside the domain (r>R_c) after being diffused according to the pde in cylindrical coordinates.
Thank you for your help
Torsten
Torsten on 28 Mar 2024 at 11:45
Edited: Torsten on 28 Mar 2024 at 11:45
(1+sin(2*pi/24/3600*t)/2
is between 1/2 and 3/2.
That's why I corrected
pl = P * (1+sin(2*pi/24/3600*t)/2 - (cl/Kav)); %Sin BC
ql = 1;
to
pl = P * ((1+sin(2*pi/24/3600*t))/2 - (cl/Kav)); %Sin BC
ql = 1;

Sign in to comment.

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!