Pdepe cylindrical coordinates with Neumann BC
10 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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
More Answers (0)
See Also
Categories
Find more on PDE Solvers in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!