System of PDE-ODE with multi point boundary conditions- How to can i solve it?

I have attached a pdf which includes the system i am trying to solve. To summerise, it is a system of PDEs and ODEs. The system consists of several regions that are joined together (see the figure in the pdf) and each region has different physical properties however the governing equation remains the same. The boundary conditions are provided at the intersection of each region. I am not sure how to solve this problem. I am struggling to implement this in matlab. I kindly ask your support and i thank you in advance.

[Torsten]

I don't understand which kind of answer you expect. You give us 5 PDEs and 2 ODEs and ask for our support. We won't code this problem for you. So ask a specific question for a specific problem you encounter (preferably with your code so far) and we might be able to help you.

And yes - the method of lines is applicable here (also for mulitpoint boundary value problems). But you need not only one, but two transmission conditions at the internal boundaries for each solution variable.

[Sachin Hegde]

Hi Torsten, I am sorry if my question was misleading. i am not asking here for the code. I have basic understanding of matlab but i am not so well versed with all the methods, syntaxes... that are available here.

My question was, which method can solve such problems (you mentioned method of lines and this was the answer i was looking for. Thank you). I have tried to understand MOL through simple implementations of PDE as an example. However the internal boundary conditions is what i am struggling at. You mentioned now, that there needs to be two transmission conditions at internal boundaries. I am not quite clear by what you mean. Could you please explain it in detail?

I genuinly thank you for the solutions you have been providing.

[Torsten]

You mentioned now, that there needs to be two transmission conditions at internal boundaries. I am not quite clear by what you mean. Could you please explain it in detail?

You have a second-order PDEs - thus you must provide two conditions at the interfaces to fix the solution. These conditions usually are continuity of the field variable and its flux:

uleft = uright
Dleft * du/dx(x->xleft) = Dright * du/dx(x->xright)

You may have other conditions (which might even cause a discontinuity at the interface), but you need two of them - at least if the field variable is solved for in both regions.

Maybe the description here is more detailed:

https://www.rpi.edu/dept/chem-eng/WWW/faculty/plawsky/Comsol%20Modules/TwoRegionDiffusion/TwoRegionDiffusion.html

[Sachin Hegde]

Thank you Torsten. This clears things up. In my problem i want continuity between the regions. I will figure things out and write the complete code. If there are still problems, then i will get back to you. Thank you and have a nice day!

[Sachin Hegde]

I found this sample code from one of your answers Interface boundary condition in transient heat problem - MATLAB Answers - MATLAB Central (mathworks.com). I am guessing that you are suggesting me this approach for interface BCs.

[Torsten]

Yes of course - I recommend my answer :-)

[Sachin Hegde]

Hi Torsten, in your sample code, i am not able to understand what you did for a11,a12,b1,a21,a22 and b2. Could you please elaborate it to me? How did you get those equations?

Your 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).';

rho1 = 100;

cp1 = 10;

lambda1 = 0.4;

rho2 = 20;

cp2 = 25;

lambda2 = 10;

% Set initial conditions

Tstart = 200;

Tend = 400;

T0 = 0;

nodes = nx1+nx2-1;

y0 = zeros(nodes,1);

y0(1) = Tstart;

y0(2:nx1+nx2-2) = T0;

y0(nx1+nx2-1) = Tend;

% Set time span of integration

tspan = 0:10:100;

[T,Y] = ode15s(@(t,y)fun(t,y,x1,nx1,x2,nx2,lambda1,rho1,cp1,lambda2,rho2,cp2),tspan,y0);

% Plot results

x = [x1;x2(2:end)];

plot(x,Y.')

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

[Torsten]

For this, read the comment before my answer.

[Sachin Hegde]

Hi Torsten, yes i am looking at it as well. If i am correct, equations 3 and 4 are implemented for dT2dt. The computation of ghost points are represented equations 1 and 2. Am i correct?

[Torsten]

The system of equations to determine T_l_ghost and T_r_ghost (rhs of equation (1) = rhs of equation (2) and rhs of equation (3) = rhs of equation (4)) is written as

a11*T_l_ghost + a12*T_r_ghost = b1  (1)
a21*T_l_ghost + a22*T_r_ghost = b2  (2)

and is then solved for (T_l_ghost,T_r_ghost) (which are named (Tg1,Tg2) in the code) as

[T_l_ghost,T_r_ghost] = [a11 a12;a21 a22]\[b1;b2]

[Sachin Hegde]

I understood that part. What i didnt get is, how did you decompose and got the equation for b

[Torsten]

I understood that part. What i didnt get is, how did you decompose and got the equation for b

By collecting all terms of

rhs of equation (1) = rhs of equation (2)

and

rhs of equation (3) = rhs of equation (4)

that don't contain T_ghost_l or T_ghost_r on the right-hand side.

[Sachin Hegde]

Thank you, now i understand. I have one last question on your solution. How did you descretize the PDE? I usually use FDM for second order DEs. However yours seem to be somewhat different method. How did you arrive at these two equations?

rho_l*cp_l * dT/dt (xi-) = lambda_l*((T_l_ghost - Ti)/(x_l_ghost - xi) - (Ti - T_l(n-1))/(xi - x_l(n-1)))/((x_l_ghost+xi)/2 - (xi+x_l(n-1))/2) (3)

rho_r*cp_r * dT/dt (xi+) = lambda_r*((T_r2) - Ti)/(x_r2 - xi) - (Ti - T_r_ghost)/(xi - x_r_ghost))/((x_r2+xi)/2 - (xi+x_r_ghost)/2) (4)

[Torsten]

It's d^2/dx^2(i) ~ (u(i+1)-2*u(i)+u(i-1))/h^2 for a possibly non-uniform gridding. Assume (x_l_ghost - xi) = (xi - x_l(n-1)) = h, and you'll arrive at the usual approximation for the second derivative.

[Sachin Hegde]

Perfect, i tried to solve it on paper and it checks out. Thank you. So this makes it clear for continuity, however what if there is no continuity at the interfaces. Lets say that in the right zone the PDE applies but not on the left. In such cases how do we define the condition at the interface?

Consider equation: di/dx = -aj (i = -sigma*d(phi)/dx). a is a constnat, j is current density, sigma a constant, phi is electron potential

Boundary condition (---- Continuity; xxxx No calculation)

| Zone 1 | Zone 2 | Zone 3 | Zone 4 | Zone 5 |

phi=0 ------------------ i=0 xxxxx i=0 --------------------i=I

For Zone 3 there is no calculation. However there is a BC at the interfaces. How can i implement the BC at the interface in this case?.

[Torsten]

In this case, since phi in zones 1 and 2 and phi in zones 4 and 5 are not coupled, these should be two independent boundary value problems to be solved. There are no transmission conditions here, but only 4 x 1 boundary conditions (one at the start of zone 1, one at the end of zone 2, one at the start of zone 4 and one at the end of zone 5), just as you wrote above.

[Sachin Hegde]

Thank you. that clears everything! I am now trying to model so that i understand it by practice. However i am getting an error (Warning: Failure at t=3.571866e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.938894e-18) at time t. ). If you could have a look it would be very helpful. The PDE i chose is very similar to your example code.

My PDE is: rho*Cp*dT/dt = kt*d2T/dx2 + Qtot; where kt is thermal conductivity, Qtot is heat source

At the moment i chose Qtot to be constants, but in the future i will calculate it when i have other PDEs solved.

Boundary condition (---- Continuity)

| Zone 1 | Zone 2 | Zone 3 | Zone 4 | Zone 5 |

T=343 -------------------------------------------------------T=343

My code:

clc
clear all
close all

L = [160 10 25 10 160]*1e-6; % [m] MEA layer thicknesses
kt = [1.6 0.27 0.3 0.27 1.6]; % [W/mK] Thermal conductivity 
rcp = [1.58 1.56 1.9 1.56 1.58]*1e6;    % [J/m3K] Volumetric heat capacity 
Qtot = [0 10 100 120 15];     % [W] Heat Source

% Initial mesh
Nd = numel(L); % number of domains
nm = 11;
% x = interp1(0:Nd, Lsum, linspace(0, Nd, Nd*10+1));
% x = sort([x Lsum(2:end-1)]); % duplicate interface nodes
x1 = linspace(0,L(1),11)';
x2 = linspace(x1(end),x1(end)+L(2),nm).';
x3 = linspace(x2(end),x2(end)+L(3),nm).';
x4 = linspace(x3(end),x3(end)+L(4),nm).';
x5 = linspace(x4(end),x4(end)+L(5),nm).';
x = [x1 x2 x3 x4 x5];

% Set initial conditions
T_start = 343;
T_end = T_start;
T0 = 0;

nodes = Nd*nm-1;
y0 = zeros(nodes,1);
y0(1) = T_start;
y0(2:nodes-1) = T0;
y0(nodes) = T_end;

% Set time span of integration
tspan = 0:10:100;

[T,Y] = ode15s(@(t,y)fun(t,y,x1,x2,x3,x4,x5,x,kt,rcp,nm,L,Qtot),tspan,y0);

Warning: Failure at t=3.610610e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.938894e-18) at time t.

%Plot results


% Function
function dy = fun(t,y,x1,x2,x3,x4,x5,x,kt,rcp,nm,L,Qtot)

T1 = y(1:nm);
T2 = y(nm:2*nm-1);
T3 = y(2*nm:3*nm-1);
T4 = y(3*nm:4*nm-1);
T5 = y(4*nm:5*nm-1);

% Compute temperature at ghost points
for i = 1:numel(L)-1
    xg(i,1) = x(end,i)+(x(end,i)-x(end-1,i));
    xg(i,2) = x(1,i+1)-(x(2,i+1)-x(1,i+1));
    
    %Setup linear system of equations
    a11 = kt(i)/(xg(i,1)-x(nm-1,i));
    a12 = kt(i+1)/(x(2,i+1)-x(1,i+1));
    b1 = kt(i)*T1(nm-1)/(xg(i,1)-x(nm-1,i))+kt(i+1)*T2(2)/(x(2,i+1)-xg(i,2));
    
    a21 = (kt(i)/(xg(i,1)-x(nm,i))/((x(nm,i)+xg(i,1))/2 - (x(nm,1)+x(nm-1,i))/2)+Qtot(i))*rcp(i+1);
    a22 = (-kt(i+1)/(x(1,i+1)-xg(i,2))/((x(2,i+1)+x(1,i+1))/2 - (x(1,i+1)+xg(i,2))/2)+Qtot(i+1))*rcp(i);
    b2  = ((kt(i)*T1(nm)/(xg(i,1)-x(nm,i))+kt(i)*(T1(nm)-T1(nm-1))/(x(nm,i)-x(nm-1,i)))/...
        ((x(nm,i)+xg(i,1))/2-(x(nm,i)+x(nm-1,i))/2)+Qtot(i))*rcp(i+1) +...
        ((kt(i+1)*(T2(2)-T2(1))/(x(2,i+1)-x(1,i+1))-kt(i+1)*T2(1)/(x(1,i+1)-xg(i,2)))/...
        ((x(2,i+1)+x(1,i+1))/2 - (x(1,i+1)+xg(i,2))/2)+Qtot(i+1))*rcp(i);
    
    A = [a11,a12;a21,a22];
    b = [b1;b2];
    
    % Solve linear system for [Tg1, Tg2]
    sol = A\b;
    Tg(i,1) = sol(1);
    Tg(i,2) = sol(2);
    
end

% Governing Equation
% AGDL
T1 = [T1;Tg(1,1)];
x1 = [x1;xg(1,1)];
dT1dt(1) = 0;
for i = 2:numel(x1)-1
    dT1dt(i) = ((kt(1)*(T1(i+1)-T1(i))/(x1(i+1)-x1(i)) -...
        kt(1)*(T1(i)-T1(i-1))/(x1(i)-x1(i-1)))/...
        ((x1(i+1)+x1(i))/2-(x1(i)+x1(i-1))/2)+Qtot(1))/(rcp(1));
end

% ACL
T2 = [T2;Tg(2,1)];
x2 = [x2;xg(2,1)];
for i = 2:numel(x2)-1
    dT2dt(i) = ((kt(2)*(T2(i+1)-T2(i))/(x2(i+1)-x2(i)) -...
        kt(2)*(T2(i)-T2(i-1))/(x2(i)-x2(i-1)))/...
        ((x2(i+1)+x2(i))/2-(x2(i)+x2(i-1))/2)+Qtot(2))/(rcp(2));
end
% PEM
T3 = [T3;Tg(3,1)];
x3 = [x3;xg(3,1)];
for i = 2:numel(x3)-1
    dT3dt(i) = ((kt(3)*(T3(i+1)-T3(i))/(x3(i+1)-x3(i)) -...
        kt(3)*(T3(i)-T3(i-1))/(x3(i)-x3(i-1)))/...
        ((x3(i+1)+x3(i))/2-(x3(i)+x3(i-1))/2)+Qtot(3))/(rcp(3));
end
% CCL
T4 = [T4;Tg(4,1)];
x4 = [x4;xg(4,1)];
for i = 2:numel(x4)-1
    dT4dt(i) = ((kt(4)*(T4(i+1)-T4(i))/(x4(i+1)-x4(i)) -...
        kt(4)*(T4(i)-T4(i-1))/(x4(i)-x4(i-1)))/...
        ((x4(i+1)+x4(i))/2-(x4(i)+x4(i-1))/2)+Qtot(4))/(rcp(4));
end
% CGDL
for i = 2:numel(x5)-1
    dT5dt(i) = ((kt(5)*(T5(i+1)-T5(i))/(x5(i+1)-x5(i)) -...
        kt(5)*(T5(i)-T5(i-1))/(x5(i)-x5(i-1)))/...
        ((x5(i+1)+x5(i))/2-(x5(i)+x5(i-1))/2)+Qtot(5))/(rcp(5));
end
dT5dt(end+1) = 0;

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

[Torsten]

You compute dT2dt(2:10), dT3(2:10),dT4(2:10), but you return dy = [...,dT2dt(1:10),dT3dt(1:10),dT4dt(1:10),...].

So something in your numbering must be wrong.

[Sachin Hegde]

Thank you, I was able to solve it yesterday (i had made mistake with the number of nodes). Now the code runs however, the results are not correct. Compared to your example code, lambda is replaced by kt and i have additional heat source (Qtot). When i turn heat source to zero then it should be exactly like your example. But In my result, there seems to be no effect of time. Right from the second time step the values for T remain constant. Also i am getting negative values which doesnt make sense either (see figure). I have went through my code several times from yesterday and i am unable to figure out the source of the problem.

Current Code:

clc
clear all
close all

L = [160 10 25 10 160]*1e-6; % [m] MEA layer thicknesses
kt = [1.6 0.27 0.3 0.27 1.6]; % [W/mK] Thermal conductivity 
rcp = [1.58 1.56 1.9 1.56 1.58]*1e6;    % [J/m3K] Volumetric heat capacity 
Qtot = [0 10 100 120 15]*1e3;     % [W] Heat Source

% Initial mesh
Nd = numel(L); % number of domains
nm = 11; % 
x1 = linspace(0,L(1),11)';
x2 = linspace(x1(end),x1(end)+L(2),11).';
x3 = linspace(x2(end),x2(end)+L(3),11).';
x4 = linspace(x3(end),x3(end)+L(4),11).';
x5 = linspace(x4(end),x4(end)+L(5),11).';
x = [x1 x2 x3 x4 x5];

% Set initial conditions
T_start = 343;
T_end = T_start;
T0 = 343;

nodes = (Nd*nm)-(Nd-1);
y0 = zeros(nodes,1);
y0(1) = T_start;
y0(2:nodes-1) = T0;
y0(nodes) = T_end;

% Set time span of integration
tspan = 0:1:10;

[T,Y] = ode15s(@(t,y)fun(t,y,x1,x2,x3,x4,x5,x,kt,rcp,nm,L,Qtot,Nd),tspan,y0);

%Plot results
xplot = [x1;x2(2:end);x3(2:end);x4(2:end);x5(2:end)];
plot(xplot,Y.');

% Function
function dy = fun(t,y,x1,x2,x3,x4,x5,x,kt,rcp,nm,L,Qtot,Nd)

T1 = y(1:nm);
T2 = y(nm:2*nm-(Nd-4));
T3 = y(2*nm-(Nd-4):3*nm-(Nd-3));
T4 = y(3*nm-(Nd-3):4*nm-(Nd-2));
T5 = y(4*nm-(Nd-2):5*nm-(Nd-1));

% Compute temperature at ghost points
for i = 1:numel(L)-1
    xg(i,1) = x(end,i)+(x(end,i)-x(end-1,i));
    xg(i,2) = x(1,i+1)-(x(2,i+1)-x(1,i+1));
    
    %Setup linear system of equations
    a11 = kt(i)/(xg(i,1)-x(nm-1,i));
    a12 = kt(i+1)/(x(2,i+1)-x(1,i+1));
    b1 = kt(i)*T1(nm-1)/(xg(i,1)-x(nm-1,i))+kt(i+1)*T2(2)/(x(2,i+1)-xg(i,2))+Qtot(i+1)-Qtot(i);
    
    a21 = (kt(i)/(xg(i,1)-x(nm,i))/((x(nm,i)+xg(i,1))/2 - (x(nm,1)+x(nm-1,i))/2)+Qtot(i))*rcp(i+1);
    a22 = (-kt(i+1)/(x(1,i+1)-xg(i,2))/((x(2,i+1)+x(1,i+1))/2 - (x(1,i+1)+xg(i,2))/2)+Qtot(i+1))*rcp(i);
    b2  = ((kt(i)*T1(nm)/(xg(i,1)-x(nm,i))+kt(i)*(T1(nm)-T1(nm-1))/(x(nm,i)-x(nm-1,i)))/...
        ((x(nm,i)+xg(i,1))/2-(x(nm,i)+x(nm-1,i))/2)+Qtot(i))*rcp(i+1) +...
        ((kt(i+1)*(T2(2)-T2(1))/(x(2,i+1)-x(1,i+1))-kt(i+1)*T2(1)/(x(1,i+1)-xg(i,2)))/...
        ((x(2,i+1)+x(1,i+1))/2 - (x(1,i+1)+xg(i,2))/2)+Qtot(i+1))*rcp(i);
    
    A = [a11,a12;a21,a22];
    b = [b1;b2];
    
    % Solve linear system for [Tg1, Tg2]
    sol = A\b;
    Tg(i,1) = sol(1);
    Tg(i,2) = sol(2);
    
end

% Governing Equation
% AGDL
T1 = [T1;Tg(1,1)];
x1 = [x1;xg(1,1)];
dT1dt(1) = 0;
for i = 2:numel(x1)-1
    dT1dt(i) = ((kt(1)*(T1(i+1)-T1(i))/(x1(i+1)-x1(i)) -...
        kt(1)*(T1(i)-T1(i-1))/(x1(i)-x1(i-1)))/...
        ((x1(i+1)+x1(i))/2-(x1(i)+x1(i-1))/2)+Qtot(1))/(rcp(1));
end

% ACL
T2 = [T2;Tg(2,1)];
x2 = [x2;xg(2,1)];
for i = 2:numel(x2)-1
    dT2dt(i) = ((kt(2)*(T2(i+1)-T2(i))/(x2(i+1)-x2(i)) -...
        kt(2)*(T2(i)-T2(i-1))/(x2(i)-x2(i-1)))/...
        ((x2(i+1)+x2(i))/2-(x2(i)+x2(i-1))/2)+Qtot(2))/(rcp(2));
end
% PEM
T3 = [T3;Tg(3,1)];
x3 = [x3;xg(3,1)];
for i = 2:numel(x3)-1
    dT3dt(i) = ((kt(3)*(T3(i+1)-T3(i))/(x3(i+1)-x3(i)) -...
        kt(3)*(T3(i)-T3(i-1))/(x3(i)-x3(i-1)))/...
        ((x3(i+1)+x3(i))/2-(x3(i)+x3(i-1))/2)+Qtot(3))/(rcp(3));
end
% CCL
T4 = [T4;Tg(4,1)];
x4 = [x4;xg(4,1)];
for i = 2:numel(x4)-1
    dT4dt(i) = ((kt(4)*(T4(i+1)-T4(i))/(x4(i+1)-x4(i)) -...
        kt(4)*(T4(i)-T4(i-1))/(x4(i)-x4(i-1)))/...
        ((x4(i+1)+x4(i))/2-(x4(i)+x4(i-1))/2)+Qtot(4))/(rcp(4));
end
% CGDL
for i = 2:numel(x5)-1
    dT5dt(i) = ((kt(5)*(T5(i+1)-T5(i))/(x5(i+1)-x5(i)) -...
        kt(5)*(T5(i)-T5(i-1))/(x5(i)-x5(i-1)))/...
        ((x5(i+1)+x5(i))/2-(x5(i)+x5(i-1))/2)+Qtot(5))/(rcp(5));
end
dT5dt(end+1) = 0;

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

[Torsten]

You must either have a differential equation for dT1dt(x1(end)) or dT2dt(x2(1)) to compute T at the first interface.

You must either have a differential equation for dT2dt(x2(end)) or dT3dt(x3(1)) to compute T at the second interface.

You must either have a differential equation for dT3dt(x3(end)) or dT4dt(x4(1)) to compute T at the third interface.

You must either have a differential equation for dT4dt(x4(end)) or dT5dt(x5(1)) to compute T at the fourth interface.

But you don't have these 4 differential equations in your code.

[Sachin Hegde]

x1 = [x1;xg(1,1)]; This line includes actual interface point x1(end) when i say numel(x1)-1 right?

Also in your code, i do not see dT1dt(x1(end)) or dT2dt(x2(1)) to compute at the interface. I formulated in the similar way. So i am not understanding what you mean.

Your sample code:

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;

[Torsten]

Yes, you are correct. I overlooked that the xi vectors are enlarged by the ghost point. Thus the loops run up to the original length of the xi vectors.

[Sachin Hegde]

Yes. Initially i thought there is some mistake in the heat source values, but when i make Qtot to zero, then it should be a normal heat flow problem without any heat source within the zones. Even with that, the results are not plausible.

[Sachin Hegde]

Oh i found a mistake. In the line for b1 and b2, i took always T1 and T2 instead of i within for loop. I have corrected it, but now again i have error.

Updated code:

clc
clear all
close all

L = [160 10 25 10 160]*1e-6; % [m] MEA layer thicknesses
kt = [1.6 0.27 0.3 0.27 1.6]; % [W/mK] Thermal conductivity 
rcp = [1.58 1.56 1.9 1.56 1.58]*1e6;    % [J/m3K] Volumetric heat capacity 
Qtot = [0 10 100 120 15]*1e3;     % [W] Heat Source

% Initial mesh
Nd = numel(L); % number of domains
nm = 11;
% x = interp1(0:Nd, Lsum, linspace(0, Nd, Nd*10+1));
% x = sort([x Lsum(2:end-1)]); % duplicate interface nodes
x1 = linspace(0,L(1),11)';
x2 = linspace(x1(end),x1(end)+L(2),nm).';
x3 = linspace(x2(end),x2(end)+L(3),nm).';
x4 = linspace(x3(end),x3(end)+L(4),nm).';
x5 = linspace(x4(end),x4(end)+L(5),nm).';
x = [x1 x2 x3 x4 x5];

% Set initial conditions
T_start = 343;
T_end = T_start;
T0 = 0;

nodes = (Nd*nm)-(Nd-1);
y0 = zeros(nodes,1);
y0(1) = T_start;
y0(2:nodes-1) = T0;
y0(nodes) = T_end;

% Set time span of integration
tspan = 0:1:10;

[T,Y] = ode15s(@(t,y)fun(t,y,x1,x2,x3,x4,x5,x,kt,rcp,nm,L,Qtot,Nd),tspan,y0);

Warning: Failure at t=1.606731e-04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-19) at time t.

%Plot results
xplot = [x1;x2(2:end);x3(2:end);x4(2:end);x5(2:end)];
plot(xplot,Y.');

% Function
function dy = fun(t,y,x1,x2,x3,x4,x5,x,kt,rcp,nm,L,Qtot,Nd)

T1 = y(1:nm);
T2 = y(nm:2*nm-(Nd-4));
T3 = y(2*nm-(Nd-4):3*nm-(Nd-3));
T4 = y(3*nm-(Nd-3):4*nm-(Nd-2));
T5 = y(4*nm-(Nd-2):5*nm-(Nd-1));
T_coll = [T1 T2 T3 T4 T5];

% Compute temperature at ghost points
for i = 1:numel(L)-1
    xg(i,1) = x(end,i)+(x(end,i)-x(end-1,i));
    xg(i,2) = x(1,i+1)-(x(2,i+1)-x(1,i+1));
    
    %Setup linear system of equations
    a11 = kt(i)/(xg(i,1)-x(nm-1,i));
    a12 = kt(i+1)/(x(2,i+1)-x(1,i+1));
    b1 = kt(i)*T_coll(nm-1,i)/(xg(i,1)-x(nm-1,i))+kt(i+1)*T_coll(2,i+1)/(x(2,i+1)-xg(i,2))+Qtot(i+1)-Qtot(i);
    
    a21 = (kt(i)/(xg(i,1)-x(nm,i))/((x(nm,i)+xg(i,1))/2 - (x(nm,1)+x(nm-1,i))/2)+Qtot(i))*rcp(i+1);
    a22 = (-kt(i+1)/(x(1,i+1)-xg(i,2))/((x(2,i+1)+x(1,i+1))/2 - (x(1,i+1)+xg(i,2))/2)+Qtot(i+1))*rcp(i);
    b2  = ((kt(i)*T_coll(nm,i)/(xg(i,1)-x(nm,i))+kt(i)*(T_coll(nm,i)-T_coll(nm-1,i))/(x(nm,i)-x(nm-1,i)))/...
        ((x(nm,i)+xg(i,1))/2-(x(nm,i)+x(nm-1,i))/2)+Qtot(i))*rcp(i+1) +...
        ((kt(i+1)*(T_coll(2,i+1)-T_coll(1,i+1))/(x(2,i+1)-x(1,i+1))-kt(i+1)*T_coll(1,i+1)/(x(1,i+1)-xg(i,2)))/...
        ((x(2,i+1)+x(1,i+1))/2 - (x(1,i+1)+xg(i,2))/2)+Qtot(i+1))*rcp(i);
    
    A = [a11,a12;a21,a22];
    b = [b1;b2];
    
    % Solve linear system for [Tg1, Tg2]
    sol = A\b;
    Tg(i,1) = sol(1);
    Tg(i,2) = sol(2);
    
end

% Governing Equation
% AGDL
T1 = [T1;Tg(1,1)];
x1 = [x1;xg(1,1)];
dT1dt(1) = 0;
for i = 2:numel(x1)-1
    dT1dt(i) = ((kt(1)*(T1(i+1)-T1(i))/(x1(i+1)-x1(i)) -...
        kt(1)*(T1(i)-T1(i-1))/(x1(i)-x1(i-1)))/...
        ((x1(i+1)+x1(i))/2-(x1(i)+x1(i-1))/2)+Qtot(1))/(rcp(1));
end

% ACL
T2 = [T2;Tg(2,1)];
x2 = [x2;xg(2,1)];
for i = 2:numel(x2)-1
    dT2dt(i) = ((kt(2)*(T2(i+1)-T2(i))/(x2(i+1)-x2(i)) -...
        kt(2)*(T2(i)-T2(i-1))/(x2(i)-x2(i-1)))/...
        ((x2(i+1)+x2(i))/2-(x2(i)+x2(i-1))/2)+Qtot(2))/(rcp(2));
end
% PEM
T3 = [T3;Tg(3,1)];
x3 = [x3;xg(3,1)];
for i = 2:numel(x3)-1
    dT3dt(i) = ((kt(3)*(T3(i+1)-T3(i))/(x3(i+1)-x3(i)) -...
        kt(3)*(T3(i)-T3(i-1))/(x3(i)-x3(i-1)))/...
        ((x3(i+1)+x3(i))/2-(x3(i)+x3(i-1))/2)+Qtot(3))/(rcp(3));
end
% CCL
T4 = [T4;Tg(4,1)];
x4 = [x4;xg(4,1)];
for i = 2:numel(x4)-1
    dT4dt(i) = ((kt(4)*(T4(i+1)-T4(i))/(x4(i+1)-x4(i)) -...
        kt(4)*(T4(i)-T4(i-1))/(x4(i)-x4(i-1)))/...
        ((x4(i+1)+x4(i))/2-(x4(i)+x4(i-1))/2)+Qtot(4))/(rcp(4));
end
% CGDL
for i = 2:numel(x5)-1
    dT5dt(i) = ((kt(5)*(T5(i+1)-T5(i))/(x5(i+1)-x5(i)) -...
        kt(5)*(T5(i)-T5(i-1))/(x5(i)-x5(i-1)))/...
        ((x5(i+1)+x5(i))/2-(x5(i)+x5(i-1))/2)+Qtot(5))/(rcp(5));
end
dT5dt(end+1) = 0;

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

Answer 1

[Torsten]

% Set physical parameters
L = [160 10 25 10 160]*1e-6; % [m] MEA layer thicknesses
kt = [1.6 0.27 0.3 0.27 1.6]; % [W/mK] Thermal conductivity
rcp = [1.58 1.56 1.9 1.56 1.58]*1e6;    % [J/m3K] Volumetric heat capacity
Qtot = [0 10 100 120 15]*1e3;     % [W] Heat Source
%Qtot = zeros(1,5);
% Number of mesh points in zones
nx = [11 11 11 11 11];
% Set boundary conditions
T_start = 343;
T_end = T_start;
% Set initial conditions
T0 = 0;
% Set time span of integration
tspan = 0:0.001:0.025;

% Create grid in domains
% Number of domains
Nd = numel(L); 
% Grid vectors
xstart = 0;
xend = L(1);
x{1} = linspace(xstart,xend,nx(1)).';
for i = 2:Nd
  xstart = xend;
  xend = xstart + L(i);
  x{i} = linspace(xstart,xend,nx(i)).';
end

% Number of ODEs to be solved
nodes = sum(nx) - (Nd-1); 

% Create initial solution vector 
y0 = zeros(nodes,1);
y0(1) = T_start;
y0(2:nodes-1) = T0;
y0(nodes) = T_end;

% Call solver
[T,Y] = ode15s(@(t,y)fun(t,y,Nd,nx,x,kt,rcp,Qtot),tspan,y0);

% Plot results
xplot = [x{1};x{2}(2:end);x{3}(2:end);x{4}(2:end);x{5}(2:end)];
plot(xplot,Y.');

% Function
function dy = fun(t,y,Nd,nx,x,kt,rcp,Qtot)

  % Write y in local variables
  ileft = 1;
  iright = nx(1);
  T{1} = y(ileft:iright);
  for j = 2:Nd
    ileft = iright ;
    iright = ileft + nx(j) - 1;
    T{j} = y(ileft:iright);
  end

  % Compute values in ghost points
  for j = 1:Nd-1
    T1 = T{j};
    T2 = T{j+1};
    x1 = x{j};
    x2 = x{j+1};
    lambda1 = kt(j);
    lambda2 = kt(j+1);
    rhocp1 = rcp(j);
    rhocp2 = rcp(j+1);
    nx1 = nx(j);
    nx2 = nx(j+1);
    Q1 = Qtot(j);
    Q2 = Qtot(j+1);

    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))*rhocp2;
    a22 = (-lambda2/(x2(1)-xg2)/((x2(2)+x2(1))/2 - (x2(1)+xg2)/2))*rhocp1;
    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)+Q1)*rhocp2 +...
         ((lambda2*(T2(2)-T2(1))/(x2(2)-x2(1))-lambda2*T2(1)/(x2(1)-xg2))/...
          ((x2(2)+x2(1))/2 - (x2(1)+xg2)/2)+Q2)*rhocp1;
    A = [a11,a12;a21,a22];
    b = [b1;b2];
    % Solve linear system for [Tg1, Tg2]
    sol = A\b;
    Tg1(j) = sol(1);
    Tg2(j) = sol(2);
    Xg1(j) = xg1;
    Xg2(j) = xg2;
  end

  % Compute time derivatives
  dTdt{1}(1) = 0.0;
  for j = 1:Nd
    if j < Nd
      Tj = [T{j};Tg1(j)];
      xj = [x{j};Xg1(j)];
    else
      Tj = T{j};
      xj = x{j};
    end
    lambdaj = kt(j);
    rhocpj = rcp(j);
    Qj = Qtot(j);
    for i=2:numel(xj)-1
      dTdt{j}(i) = ((lambdaj*(Tj(i+1)-Tj(i))/(xj(i+1)-xj(i)) -...
                     lambdaj*(Tj(i)-Tj(i-1))/(xj(i)-xj(i-1)))/...
                     ((xj(i+1)+xj(i))/2-(xj(i)+xj(i-1))/2)+Qj)/(rhocpj);
    end
  end
  dTdt{Nd}(end+1) = 0.0;

  % Return time derivatives
  dy = dTdt{1}(1:end);
  for j = 2:Nd
    dy = [dy,dTdt{j}(2:end)];
  end
  dy = dy.';
 end

[Sachin Hegde]

Wow! Wonderful! I need to go through the code once to understand where i went wrong. Thank you very much for your support Torsten. Have a nice weekend!

[Torsten]

"Thank you" here means accepting the answer.