Run pdepe solver in adjacent domains: introduce flux BC at the interface

3 次查看（过去 30 天）
Ali Aykut
Ali Aykut 2023-11-18
编辑： Torsten 2023-11-20
Ran in:
I have heat diffusion equation in spherical coordinates. I try to model two neighbouring domains (interface at r=r_T). For that I added an if condition for introducing the pde. However, I also need to introduce additional flux BC at the interface of materials. At origin (r=0) and (r=R) I have 0 flux BC. At the interface I need to intoduce an additional BC as:
I believe I need to modify advdiff_bc function but I not sure how to. Thank you in advance for your help!
clc
close all
rho = 1000; % fluid density (kg/m^3)
mu = 1.002E-3; % fluid dynamic viscosity
k_B = 1.38E-23; % Boltzman constant (J/K)
kin_v = mu/rho; % kinematic viscosity
T = 293; % temperature of the medium (K)
d_p = 38-9; % diameter of the particle (m)
Dtherotical_SE = k_B * T / (3 * pi * mu * d_p) %calculated using stokes-einstein eqn
Dtherotical_SE = 1.4764e-20
NP1 = [0 6/6; 1*3600 5/6; 6*3600 3.75/6; 20*3600 1.85/6; 48*3600 0.7/6; 95*3600 0.3/6];
[r,t,c] = AdvDiff_sph(4E-14, 1.65E4, 1.74E4, 0.4, NP1, Dtherotical_SE);
plot(r,[c(2,:);c(20,:);c(200,:);c(2000,:);c(20000,:);c(end,:)])
%% 1D Spherical Drug Transport Model
% c = AdvDiff_sph(4E-14,40*133.2,1.65E4,0.4,[0 6; 1*3600 5; 6*3600 3.75; 20*3600 1.85; 48*3600 0.7; 95*3600 0.3])
function [r,t,Conc] = AdvDiff_sph(Deff,SVt,SVn,Kav,NP1,D_se)
% Define the parameters
R_T = 0.005; % radius [m]
R_N = 0.25; % radius of surrounding normal [m]
tsim = 12 * 3600; % simulation time in [s]
r = [0:1E-3:R_N]; % create spatial and temporal solution mesh
t = [0:1:tsim];
% Solve the PDE numerically using pdepe which solves the pde using ode15s
m = 2; %the 3D spherical model with azimuthal and zenith angular symmetries
[sol] = pdepe(m, @advdiff_pde, @advdiff_initial, @advdiff_bc,r,t);
Conc = sol(:,:,1);
function [c,f,s] = advdiff_pde(r,t,Conc,dcdr)
c = 1;
f = Deff*dcdr;
if r < 0.005
s = 20* D_se* SVt*(interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')-(Conc/Kav));
else
s = 4.2E-10* SVn*(interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')-(Conc/Kav));
end
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)
% take zero slope/flux BC on both origin and very far from interest
pl = 0;
ql = 1;
pr = 0;
qr = 1;
end
end

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2023-11-18
移动：Torsten 2023-11-18
If you make r = R_T a point of the r-vector, the finite element solver pdepe will automatically respect continuity of c and dc/dr at the interface at r = R_T. So no need to change your code from above.
The pdepe docs recommend you place grids point exactly at the locations of any discontinuities in the coefficients. If you do this, your pde function will never be called with x=discontinuity_location so you don't need to explicitly handle this case.
  2 个评论
Ali Aykut
Ali Aykut 2023-11-20
编辑：Torsten 2023-11-20
Thank you for the information, I didn't know pdepe was already checking the contiunity of profile and its slope.
However, the profile I found is significantly planar, however I believe it should be dependent on r since its in spherical domain with source term. Do you think it seems reasonable?
clc
% close all
% clear all
rho = 1000; % fluid density (kg/m^3)
mu = 1.002E-3; % fluid dynamic viscosity
k_B = 1.38E-23; % Boltzman constant (J/K)
kin_v = mu/rho; % kinematic viscosity
T = 310; % temperature of the tissue 37C (K)
d_p = 38-9; % diameter of the particle (m)
Dtherotical_SE = k_B * T / (3 * pi * mu * d_p) %calculated using stokes-einstein eqn
NP1 = [0 6/6; 1*3600 5/6; 6*3600 3.75/6; 20*3600 1.85/6; 48*3600 0.7/6; 95*3600 0.3/6];
c = AdvDiff_sph(4E-14, 1.65E4, 1.74E4, 0.4, NP1, Dtherotical_SE);
%% 1D Spherical Drug Transport Model
% c = AdvDiff_sph(4E-14,40*133.2,1.65E4,0.4,[0 6; 1*3600 5; 6*3600 3.75; 20*3600 1.85; 48*3600 0.7; 95*3600 0.3])
function c = AdvDiff_sph(Deff,SVt,SVn,Kav,NP1,D_se)
% ic, timedata, concBC, cutoff, a, t0)
% Inputs:
% fxu : Solute velocity: retardation times solution velocity [m/s]
% t0 : Time [mins]
% Define the parameters
R_T = 0.005 % radius of tumor tissue [m]
R_N = 0.25 % radius of surrounding normal tissue [m]
tsim = 24 * 3600 % simulation time in [s]
% Deff = 4E-14 % Effective diffusivity [m2/s]
% Pi = 40 * 133.32; % Capillary fluid pressure mmHg to [Pa] taken from Lecture Notes
% SVi = 1.65E4; %Transvascular transport surface to volume ratio (Tumor) [1/m]
% Kav = 0.4; % Available volume fraction (Tumor)
r = [0:1E-5:R_N]; % create spatial and temporal solution mesh
t = [0:5:tsim];
% Provide the drug plasma concentration to be used as the BC [time conc];
% NP1 = [0 6; 1*3600 5; 6*3600 3.75; 20*3600 1.85; 48*3600 0.7; 95*3600 0.3];
% Solve the PDE numerically using pdepe which solves the pde using ode15s
m = 2; %the 3D spherical model with azimuthal and zenith angular symmetries
[sol] = pdepe(m, @advdiff_pde, @advdiff_initial, @advdiff_bc,r,t);
Conc = sol(:,:,1);
% x0 = [0.3:0.01:2.5]*1E3;
% x0 = [cutoff*1.6347*1E-3:0.01:3]*1E3;
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];
% % Plot the numerical solution
% h = surf(x,t./60,c);
% set(h,'edgecolor','none')
% xlabel('Distance x [\mum]', 'Interpreter','Tex');
% ylabel('Time t [min]', 'Interpreter','Tex');
% zlabel('Normalized Concentration [C/C_{max}]','Interpreter','Tex');
% title('Numerical Solution of ADE Using pdepe','Interpreter','Tex');
%
figure()
% x = x + cutoff*1.6347E-6;
i = 1;
for time = [5/3600 1 2 5 10 20]*3600/5
% c_pde(i,:) = interp1(x*1E6,c(time,:),x0,'linear','extrap');
plot(r/R_T,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2)
hold on
i = i + 1;
end
% title('Solution at 1 hr intervals')
% legend('Numerical, 20 mesh points','Analytical','Location','South')
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')
xlim([0 4])
lgd.FontSize = 6;
legend('boxoff')
% xlabel('Distance x')
% ylabel('c(x,2)')
% delC/delt = D_eff * del^2 C/delx^2 + (uR) * delC/delx
% Define the wave equation
function [c,f,s] = advdiff_pde(r,t,Conc,dcdr)
c = 1;
f = Deff*dcdr;
% Cv = interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap');
if r < 0.005
s = 4.2E-10* SVt*(interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')-(Conc/Kav));
else
% s = 0;
s = 4.2E-10* SVn*(interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')-(Conc/Kav));
% 4.2E-10
% 20* D_se
end
%this advective velocity might be a linear function such that the input
%ADE has 3 inputs, D_eff, Vinlet_1, V_inlet2
% s = -(a(2)*1E-7-((a(2)*1E-7)/(a(3)*1E3)*t))*dcdx;
% s = -(0*1E-7-((0*1E-7)/(a(3)*1E3)*t))*dcdx;
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)
% Constant BC on the left (=1)
% pl = cl - 1;
% take zero slope/flux BC on both origin and very far from interest
pl = 0;
ql = 1;
pr = cr;
qr = 0;
% pr = 0;
% qr = 1;
end
end
Torsten
Torsten 2023-11-20
编辑：Torsten 2023-11-20
However, the profile I found is significantly planar, however I believe it should be dependent on r since its in spherical domain with source term. Do you think it seems reasonable?
Your source term(s) don't depend on r - so why do you expect a profile in r-direction ?
E.g. if you have the equation
dc/dt = 1/r^2*d/dr(r^2*D*dc/dr) + 1
and you start with a concentration of 0, the c-profile will remain flat because a constant rate of 1 mol/(m^3*s) of c is formed throughout the sphere.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by