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

3 次查看（过去 30 天）
Neil Parke
Neil Parke 2023-4-7
编辑： Neil Parke 2023-4-12
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
  16 个评论
显示 14更早的评论
Torsten
Torsten 2023-4-12
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
Neil Parke 2023-4-12
编辑：Neil Parke 2023-4-12
Ran in:
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

更多回答（0 个）

类别

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

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by