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

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

Answer 1

[Torsten]

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.

[Ali Aykut]

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]

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.