Pdepe cylindrical coordinates with Neumann BC
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
Accepted Answer
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
on 3 Mar 2024
Great! thank you.
Does cl refers to the value at the left boundary? It fixed the error, but there is one thing I am not sure. I attached the solution for 20 hrs. After some point since the BC value at the left decreases, I expect negative slope at the interface. I believe when the slope is negative, the BC value on the left end supposed to be smaller than the value on the right of the interface (referring to the green line at t=20hrs I would expect a smaller value for the horizontal line). Do you have any idea on that?
%Sample Input: AdvDiff_Cyl(4E-14,4.2E-8,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-7:R_O]; % create spatial and temporal solution mesh
t = [0:1: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;
% time = [1 2 5 10 20]*3600/5
for time = [1 2 5 10 20]*3600
h1(i) = plot(r/R_C,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2);
hold on
Conc_vessel= interp1(NP1(:,1), NP1(:,2), time, 'spline', 'extrap');
plot([0 1], [Conc_vessel Conc_vessel],'-', 'color', colors(i,:),'LineWidth',2);
hold on
i = i + 1;
end
rectangle('Position',[0.99 0 0.01 5])
xlabel('Distance [r/R_N]', 'Interpreter','Tex');
ylabel('Concentration [\mug/ml]','Interpreter','Tex');
title('Numerical Solution of ADE using pdepe','Interpreter','Tex');
legend([h1(1), h1(2), h1(3), h1(4), h1(5)], '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: Transvascular Transport at R = R_c
pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (cl/Kav));
ql = 1; %ignored by solver since m=1
% BC2: Concentration gradient is zero at the oter edge
pr = 0;
qr = 1;
end
end
Does cl refers to the value at the left boundary?
Yes.
You should plot
Kav*interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')
, not
interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')
at the left interface to compare with Conc.
@Torsten Thank you!
Sorry for bothering you but I have one more problem that I couldn't adress.
On the surface plot below the shorter surface is the BC at r=R_c and the longer surface is the values inside the domain. It seems the values inside the domain are starting to be greater than the BC values after some time, which normally should be imposible. Do you think I am formulating sth wrong? For some other input parameters this error/difference between BC and domain even greater.
AdvDiff_Cyl(1E-11,4.2E-8,1.0,0,[0 1; 1*3600 1; 6*3600 1; 20*3600 1; 48*3600 1; 95*3600 1]);
function Conc = AdvDiff_Cyl(Deff,P,Kav,Phi,NP1)
% Define the parameters
R_C = 19.2*10^-6/2/pi;
R_O = 50*R_C;
tsim = 24 * 3 * 3600; % simulation time in [s]
r = [R_C:1E-7: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 ; %2-D Cylindrical coordinates with azimuthal angular symmetry
[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() %2D Surface Kymograph Plot
% ZZ = Kav*interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap'); %Constant or Specified BC
ZZ = Kav *(1+sin(2*pi/24/3600*(t+Phi*3600)))/2; %Sin BC
surf(r(1:10:numel(r))/R_C,t(1:10:numel(t))/3600,Conc((1:10:numel(t)),(1:10:numel(r))),'edgecolor','none')
hold on
surf([0 1],t/3600, ZZ'.*[1 1],'edgecolor','none')
xlabel('Distance [r/R_N]'); xlim([0 25])
ylabel('Time [hr]'); ylim([0 tsim/3600]); yticks([0:24:tsim/3600])
colorbar; caxis([0, 1]);
view(2)
figure() %1D plot
i = 1;
% time = [1 2 5 10 20]*3600/5
for time = [1 2 4 8 12 24]*3600
h1(i) = plot(r/R_C,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2);
hold on
Conc_vessel = interp1(NP1(:,1), NP1(:,2), time, 'spline', 'extrap'); %Constant or Specified BC
% Conc_vessel = (1+sin(2*pi/24/3600*time))/2; %Sin BC
plot([0 1], Kav*[Conc_vessel Conc_vessel],'-', 'color', colors(i,:),'LineWidth',2);
hold on
i = i + 1;
end
rectangle('Position',[0.8 0 0.4 5],'EdgeColor',"none",'FaceColor',[1 0 0 0.3])
ylim([0 1])
xlim([0 20])
xlabel('Distance [r/R_N]', 'Interpreter','Tex');
ylabel('Concentration [\mug/ml]','Interpreter','Tex');
title('Numerical Solution of ADE using pdepe','Interpreter','Tex');
legend([h1(1), h1(2), h1(3), h1(4), h1(5) h1(6)],'1 hr','2 hr','4 hr' ,'8 hr','12 hr','24 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.5;
end
% Define the boundary condition
function [pl,ql,pr,qr] = advdiff_bc(xl,cl,xr,cr,t)
% BC1: Transvascular Transport at R = R_c
pl = P * (1+sin(2*pi/24/3600*t)/2 - (cl/Kav)); %Sin BC
ql = 1;
% BC2: Concentration gradient is zero at the outer edge
pr = 0;
qr = 1;
end
end
Torsten
on 26 Mar 2024
Please correct your code (see above).
Ali
on 27 Mar 2024
If you add the following lines at the end of your code
figure(2)
for i = 1:10:numel(t)
[~,dudx(i)] = pdeval(m,r,sol(i,:),r(1));
bc_cond(i) = Deff*dudx(i) + P * (1+sin(2*pi/24/3600*t(i))/2 - sol(i,1)/Kav);
end
plot(t(1:10:numel(t)),bc_cond(1:10:numel(t)))
you will see that bc_cond is approximately 0.
Thus your boundary condition at r = R_C that you prescribe as
Deff*dc/dr + P * (1+sin(2*pi/24/3600*t)/2 - c/Kav);
is satisfied.
I don't know how you conclude from this boundary condition that c should not exceed 1.
Maybe you mean
pl = P * ((1+sin(2*pi/24/3600*t))/2 - (cl/Kav)); %Sin BC
ql = 1;
instead of
pl = P * (1+sin(2*pi/24/3600*t)/2 - (cl/Kav)); %Sin BC
ql = 1;
Ali
on 28 Mar 2024
(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;
More Answers (0)
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
