Pdepe cylindrical coordinates with Neumann BC
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1632701/image.png)
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 Kommentare
Akzeptierte Antwort
Torsten
am 1 Mär. 2024
Bearbeitet: Torsten
am 1 Mär. 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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu PDE Solvers finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!