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

6 次查看（过去 30 天）
Sachin Hegde
Sachin Hegde 2023-7-25
编辑： Torsten 2023-7-28
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.
  24 个评论
显示 22更早的评论
Sachin Hegde
Sachin Hegde 2023-7-28
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
Sachin Hegde 2023-7-28
编辑：Torsten 2023-7-28
Ran in:
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

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2023-7-28
编辑：Torsten 2023-7-28
Ran in:
% 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

更多回答（0 个）

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by