Solving two PDEs in parallel (linked boundary conditions) with two xmesh parameters

I am trying to solve two 1D heat equations plus convection, in parallel, that are joined by a flux condition at their right and left boundaries respectively. I believe I've achieved this using the PDEPE solver with relative ease.

However, I would like the PDEs to have separate xmesh parameters, as the physical domains they model are different lengths. This doesn't seem to be supported by PDEPE in its current form. I have looked into the PDEPE function (and the paper the function is based on) and cannot fathom how to implement it. I am sure this only requires simple adaptations to the orginal function, but I can't work it out!

I've included a skeleton version of the functions for PDEPE in my code below. I removed a lot of the other supporting code to focus in on the 'nuts and bolts'. I suspect it won't help to solve the actual problem, but it might help to give you some context.

Could anyone offer any help?

_____________________________________________________________________________________

m=0;

x = linspace(0,10,50);
t = linspace(0,30,75); 

sol = pdepe(m,governing_pdes,initial_cond,boundary_cond,x,t);

function [c,f,s] = governing_pdes(x,t,u,dudx)
c=[1;1];
f = [(alpha_rod*10^6);(alpha_leg*10^6)].*dudx;
s = [-(h_rod*A_rod/(k_nick/1000))*u(1);-(h_leg*A_leg/(k_cop/1000))*u(2)];
end

function u0 = initial_cond(x)
u0 = [0;0];
end

function [pl,ql,pr,qr] = boundary_cond(xl,ul,xr,ur,t)
if t >= 0 & t <= Pulse_Span
    pl = [(theta*(V_p^2/R)/A_rod);-(k_cop/1000)*ur(1)];
    ql = [1;1];
    pr = [-(k_nick/1000)*ul(2);-(1-theta)*(V_p^2/R)/A_leg];
    qr = [1;1];
else
    pl = [0;-(k_cop/1000)*ur(1)];
    ql = [1;1];
    pr = [-(k_nick/1000)*ul(2);0];
    qr = [1;1];
end
end

[Torsten]

If you define two pdes in one problem code for pdepe, the spatial domains in which the two pdes are solved must be the same. Or did I misunderstand your question ?

[Neil Parke]

You did not misunderstand. It isn't possible to create a custom function (by adapting the PDEPE function) to solve these two equations with separate spatial domains then?

If this isn't possible, I suspect I will have to do the spatial descretisation myself and then use ODE15 solver (as PDEPE does). Do you have any suggestions for other approaches please do say...? I am open to suggestions...

[Neil Parke]

Do you think I could solve the two spatial domains using two separate PDEPE calls, but link the boundaries somehow? I imagine that would be even more complex than implementing two xmesh parameters in a custom function...

[Torsten]

What is/are the interface condition(s) between the two regions ? Could you supply the two PDEs with inital and boundary conditions + transmission condition(s) in a mathematical notation ? It's easier than translating your MATLAB code in mathematical formulae.

[Bill Greene]

I also do not understand the problem you are trying to solve.

I assume it is not as simple as doing something like this?

x=[linspace(0,L1,25) linspace(L1+del,L2,50)];

[Neil Parke]

Of course. Please see the photos of my (attempt at) rigorous notation attached. This should give you some physical intuition. If you would like clarification on any points please ask.

As a brief description, I want to model two non-adiabatic rods of different materials, diameters, and crucially lengths that are joined at one boundary through conduction, lose heat via convection, plus have heat flux entering at their other boundary. I am aware of the perils of accurate meshing over the join of the two rods, but right now I am ignoring it for the sake of getting a (not so) basic working model.

@Bill Greene Unfortunatly it is not that simple. I have attempted the solution you suggested; however, then I lose the ability to differentiate between the two rods through linking boundary conditions.

[Torsten]

Since the two rods have different diameters, this is not an 1d, but at least a 2d problem in cylindrical coordinates.

[Neil Parke]

I am happy to neglect radial differences for now, as the effect is negligible compared to that of the length. My data actually lines up fairly well with experimental data from the system, it is just being able to alter the lengths of the two domains that is an issue for me right now.

[Torsten]

If differences in radius can be neglected (especially at the contact point), you can just proceed as @Bill Greene suggested. Use one unknown function to solve for, include the contact point of the two regions in your mesh and define

function [c,f,s] = governing_pdes(x,t,u,dudx)
  c = 1;
  if x < length_first_rod
    f = (alpha_rod*10^6)*dudx;
    s = -(h_rod*A_rod/(k_nick/1000))*u(1);
  else if x > length_first_rod
    f = (alpha_leg*10^6)*dudx;
    s = -(h_leg*A_leg/(k_cop/1000))*u(1);
  else
    f = 0.5*((alpha_rod*10^6)+(alpha_leg*10^6))*dudx;
    s = 0.5*( -(h_rod*A_rod/(k_nick/1000))+(-(h_leg*A_leg/(k_cop/1000))))
  end
end  

I think you know how to proceed for the boundary part.

[Neil Parke]

If this does work, then @Bill Greene I stand corrected and apologies. I had it in my head that I required boundary conditions to couple the equations, but if it's one equation you don't need to couple them!!

@Torsten thank you very much, I really appreciate it. Regarding boundary conditions, I would have to remove the conduction conditions between the two models I had in my original post (leaving just the heat fluxes divided up by theta at their respective sides), correct?

[Torsten]

It's necessary to set two transmission conditions at the interface. I think the finite element method used in "pdepe" will automatically ensure continuity of temperature and heat flux if you set up the code as suggested above. If you want to set different transmission conditions, you cannot use "pdepe" or some other software (like the PDE toolbox), but you will have to discretize the equations and the transmission conditions explicity and use ODE15S to solve the resulting system of ordinary differential equations (method-of-lines).

Maybe of help: Instationary heat conduction in two regions with different material properties - solved by pdepe and the method-of-lines:

code = 1;
test(code)

code = 2;
test(code)

function [] = test(code)
  % Set parameters
  x1start = 0.0;
  x1end = 0.25;
  x2start = x1end;
  x2end = 1.0;
  nx1 = 50;
  nx2 = 150;
  x1 = linspace(x1start,x1end,nx1).';
  x2 = linspace(x2start,x2end,nx2).';
  x = [x1;x2(2:end)];
  rho1 = 100;
  cp1 = 10;
  lambda1 = 0.4;
  rho2 = 20;
  cp2 = 25;
  lambda2 = 10;

  % Set initial conditions
  Tstart = 200;
  Tend = 400;
  T0 = 0;
  
  % Set interval of integration
  tspan = 0:10:100;
  
  if code == 1
    m = 0;
    sol = pdepe(m,@(x,t,u,dudx)governing_pdes(x,t,u,dudx,x1end,rho1,cp1,lambda1,rho2,cp2,lambda2),@(x)initial_cond(x,x1start,x1end,x2start,x2end,Tstart,Tend,T0),@(xl,ul,xr,ur,t)boundary_cond(xl,ul,xr,ur,t,Tstart,Tend),x,tspan);
    u = sol(:,:,1);
    % Plot results
    figure(1)
    plot(x,[u(1,:);u(2,:);u(3,:);u(4,:);u(5,:);u(6,:);u(7,:);u(8,:);u(9,:);u(10,:);u(11,:)])
  elseif code==2
    nodes = nx1+nx2-1;
    y0 = zeros(nodes,1);
    y0(1) = Tstart;
    y0(2:nx1+nx2-2) = T0;
    y0(nx1+nx2-1) = Tend;
    [T,Y] = ode15s(@(t,y)fun(t,y,x1,nx1,x2,nx2,lambda1,rho1,cp1,lambda2,rho2,cp2),tspan,y0,odeset("RelTol",1e-3,"AbsTol",1e-3));
    % Plot results
    figure(2)
    plot(x,Y.')  
  end
end

function [c,f,s] = governing_pdes(x,t,u,dudx,x1end,rho1,cp1,lambda1,rho2,cp2,lambda2)
  if x < x1end
    c = rho1*cp1;
    f = lambda1*dudx;
    s = 0.0;
  elseif x > x1end
    c = rho2*cp2;
    f = lambda2*dudx;
    s = 0.0; 
  else
    c = 0.5*(rho1*cp1+rho2*cp2);
    f = 0.5*(lambda1+lambda2)*dudx;
    s = 0.0;    
  end
end

function u0 = initial_cond(x,x1start,x1end,x2start,x2end,Tstart,Tend,T0)
   if x==x1start
    u0 = Tstart;
   elseif x==x2end
    u0 = Tend;
   else
    u0 = T0;
   end
end

function [pl,ql,pr,qr] = boundary_cond(xl,ul,xr,ur,t,Tstart,Tend)
  pl = ul - Tstart;
  ql = 0.0;
  pr = ur - Tend;
  qr = 0.0;
end

function dy = fun(t,y,x1,nx1,x2,nx2,lambda1,rho1,cp1,lambda2,rho2,cp2)
  T1 = y(1:nx1);          % Temperature zone 1
  T2 = y(nx1:nx1+nx2-1);  % Temperature zone 2

  % Compute temperature in ghost points
  xg1 = x1(end)+(x1(end)-x1(end-1));
  xg2 = x2(1)-(x2(2)-x2(1));
  % Set up linear system of equations for [Tg1, Tg2]
  a11 = lambda1/(xg1-x1(nx1-1));
  a12 = lambda2/(x2(2)-xg2);
  b1 = lambda1*T1(nx1-1)/(xg1-x1(nx1-1))+lambda2*T2(2)/(x2(2)-xg2);

  a21 = (lambda1/(xg1-x1(nx1))/((x1(nx1)+xg1)/2 - (x1(nx1)+x1(nx1-1))/2))*rho2*cp2;
  a22 = (-lambda2/(x2(1)-xg2)/((x2(2)+x2(1))/2 - (x2(1)+xg2)/2))*rho1*cp1;
  b2 = (lambda1*T1(nx1)/(xg1-x1(nx1))+lambda1*(T1(nx1)-T1(nx1-1))/(x1(nx1)-x1(nx1-1)))/...
        ((x1(nx1)+xg1)/2-(x1(nx1)+x1(nx1-1))/2)*rho2*cp2 +...
       (lambda2*(T2(2)-T2(1))/(x2(2)-x2(1))-lambda2*T2(1)/(x2(1)-xg2))/...
        ((x2(2)+x2(1))/2 - (x2(1)+xg2)/2)*rho1*cp1;

  A = [a11,a12;a21,a22];
  b = [b1;b2];
  % Solve linear system for [Tg1, Tg2]
  sol = A\b;

  Tg1 = sol(1);
  Tg2 = sol(2);

  % Solve heat equation in the two zones
  % Zone 1
  T1 = [T1;Tg1];
  x1 = [x1;xg1];
  dT1dt(1) = 0.0;
  for i=2:numel(x1)-1
    dT1dt(i) = (lambda1*(T1(i+1)-T1(i))/(x1(i+1)-x1(i)) -...
                lambda1*(T1(i)-T1(i-1))/(x1(i)-x1(i-1)))/...
                ((x1(i+1)+x1(i))/2-(x1(i)+x1(i-1))/2)/(rho1*cp1);
  end
  % Zone 2
  for i=2:numel(x2)-1
    dT2dt(i) = (lambda2*(T2(i+1)-T2(i))/(x2(i+1)-x2(i)) -...
                lambda2*(T2(i)-T2(i-1))/(x2(i)-x2(i-1)))/...
                 ((x2(i+1)+x2(i))/2-(x2(i)+x2(i-1))/2)/(rho2*cp2);
  end
  dT2dt(end+1) = 0.0;

  % Return time derivatives
  dy = [dT1dt(1:end),dT2dt(2:end)];
  dy = dy.';
end

[Neil Parke]

@Torsten, This works perfectly. I really appreciated it, thank you so much.

[Bill Greene]

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.

[Neil Parke]

@Bill Greene Did you mean DO NOT place grid points at the discontinuity?

[Torsten]

No, place a grid point exactly at the discontinuity.

@Bill Greene meant that the "else" case in the if-statement I wrote is superfluous (because the pde function is not called for x = x1end).

[Neil Parke]

@Torsten @Bill Greene Thanks again to you both. I have combined your answers below.

Answer 1

[Neil Parke]

Combined answer from @Torsten & @Bill Greene:

It's necessary to set two transmission conditions at the interface. I think the finite element method used in "pdepe" will automatically ensure continuity of temperature and heat flux if you set up the code as suggested above. If you want to set different transmission conditions, you cannot use "pdepe" or some other software (like the PDE toolbox), but you will have to discretize the equations and the transmission conditions explicity and use ODE15S to solve the resulting system of ordinary differential equations (method-of-lines).

The pdepe docs recommend you place an xmesh point at the discontinuity. 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 in the IF statement.

Maybe of help: Instationary heat conduction in two regions with different material properties - solved by pdepe and the method-of-lines:

code = 1;
test(code)

code = 2;
test(code)

function [] = test(code)
% Set parameters
x1start = 0.0;
x1end = 0.25;
x2start = x1end;
x2end = 1.0;
nx1 = 50;
nx2 = 150;
x1 = linspace(x1start,x1end,nx1).';
x2 = linspace(x2start,x2end,nx2).';
x = [x1;x2(2:end)];
rho1 = 100;
cp1 = 10;
lambda1 = 0.4;
rho2 = 20;
cp2 = 25;
lambda2 = 10;
% Set initial conditions
Tstart = 200;
Tend = 400;
T0 = 0;

% Set interval of integration
tspan = 0:10:100;

if code == 1
    m = 0;
    sol = pdepe(m,@(x,t,u,dudx)governing_pdes(x,t,u,dudx,x1end,rho1,cp1,lambda1,rho2,cp2,lambda2),@(x)initial_cond(x,x1start,x1end,x2start,x2end,Tstart,Tend,T0),@(xl,ul,xr,ur,t)boundary_cond(xl,ul,xr,ur,t,Tstart,Tend),x,tspan);
    u = sol(:,:,1);
    % Plot results
    figure(1)
    plot(x,[u(1,:);u(2,:);u(3,:);u(4,:);u(5,:);u(6,:);u(7,:);u(8,:);u(9,:);u(10,:);u(11,:)])
elseif code==2
    nodes = nx1+nx2-1;
    y0 = zeros(nodes,1);
    y0(1) = Tstart;
    y0(2:nx1+nx2-2) = T0;
    y0(nx1+nx2-1) = Tend;
    [T,Y] = ode15s(@(t,y)fun(t,y,x1,nx1,x2,nx2,lambda1,rho1,cp1,lambda2,rho2,cp2),tspan,y0,odeset("RelTol",1e-3,"AbsTol",1e-3));
    % Plot results
    figure(2)
    plot(x,Y.')  
end
end

function [c,f,s] = governing_pdes(x,t,u,dudx,x1end,rho1,cp1,lambda1,rho2,cp2,lambda2)
if x < x1end
    c = rho1*cp1;
    f = lambda1*dudx;
    s = 0.0;
else
    c = rho2*cp2;
    f = lambda2*dudx;
    s = 0.0;   
end
end

function u0 = initial_cond(x,x1start,x1end,x2start,x2end,Tstart,Tend,T0)
if x==x1start
    u0 = Tstart;
elseif x==x2end
    u0 = Tend;
else
    u0 = T0;
end
end

function [pl,ql,pr,qr] = boundary_cond(xl,ul,xr,ur,t,Tstart,Tend)
pl = ul - Tstart;
ql = 0.0;
pr = ur - Tend;
qr = 0.0;
end

function dy = fun(t,y,x1,nx1,x2,nx2,lambda1,rho1,cp1,lambda2,rho2,cp2)
T1 = y(1:nx1);          % Temperature zone 1
T2 = y(nx1:nx1+nx2-1);  % Temperature zone 2
% Compute temperature in ghost points
xg1 = x1(end)+(x1(end)-x1(end-1));
xg2 = x2(1)-(x2(2)-x2(1));
% Set up linear system of equations for [Tg1, Tg2]
a11 = lambda1/(xg1-x1(nx1-1));
a12 = lambda2/(x2(2)-xg2);
b1 = lambda1*T1(nx1-1)/(xg1-x1(nx1-1))+lambda2*T2(2)/(x2(2)-xg2);
a21 = (lambda1/(xg1-x1(nx1))/((x1(nx1)+xg1)/2 - (x1(nx1)+x1(nx1-1))/2))*rho2*cp2;
a22 = (-lambda2/(x2(1)-xg2)/((x2(2)+x2(1))/2 - (x2(1)+xg2)/2))*rho1*cp1;
b2 = (lambda1*T1(nx1)/(xg1-x1(nx1))+lambda1*(T1(nx1)-T1(nx1-1))/(x1(nx1)-x1(nx1-1)))/...
    ((x1(nx1)+xg1)/2-(x1(nx1)+x1(nx1-1))/2)*rho2*cp2 +...
    (lambda2*(T2(2)-T2(1))/(x2(2)-x2(1))-lambda2*T2(1)/(x2(1)-xg2))/...
    ((x2(2)+x2(1))/2 - (x2(1)+xg2)/2)*rho1*cp1;
A = [a11,a12;a21,a22];
b = [b1;b2];
% Solve linear system for [Tg1, Tg2]
sol = A\b;
Tg1 = sol(1);
Tg2 = sol(2);
% Solve heat equation in the two zones
% Zone 1
T1 = [T1;Tg1];
x1 = [x1;xg1];
dT1dt(1) = 0.0;
for i=2:numel(x1)-1
    dT1dt(i) = (lambda1*(T1(i+1)-T1(i))/(x1(i+1)-x1(i)) -...
        lambda1*(T1(i)-T1(i-1))/(x1(i)-x1(i-1)))/...
        ((x1(i+1)+x1(i))/2-(x1(i)+x1(i-1))/2)/(rho1*cp1);
end
% Zone 2
for i=2:numel(x2)-1
    dT2dt(i) = (lambda2*(T2(i+1)-T2(i))/(x2(i+1)-x2(i)) -...
        lambda2*(T2(i)-T2(i-1))/(x2(i)-x2(i-1)))/...
        ((x2(i+1)+x2(i))/2-(x2(i)+x2(i-1))/2)/(rho2*cp2);
end
dT2dt(end+1) = 0.0;
% Return time derivatives
dy = [dT1dt(1:end),dT2dt(2:end)];
dy = dy.';
end