Why does solving the heat equation with MATLAB (pdepe) yield a completely different result than the Heisler chart (analytical solution)?

8 次查看（过去 30 天）
Phil. U.
Phil. U. 2021-11-2
回答： Bill Greene 2021-11-4
Hi,
I want to solve the heat equation with the pdepe solver for a slab / wall. In contrast to the example provided by MATLAB, a heat transfer coeffcient α is present at both sides of the wall. Therefore, I adopted the answer given in this thread, the resulting MATLAB-code looks as follows (code is working):
x=linspace(0,0.1,5);
t=linspace(0,3600,12);
m = 0;
sol=pdepe(m,@PDE,@IC,@BC,x,t); %solution
surf(sol); % plot
function u0=IC(x)
u0=0; % Initial temperature = 0°C
end
function [pl,ql,pr,qr]=BC(xl,ul,xr,ur,t)
T_ext = 50; % Ambient temperature = 50 °C
alpha = 10; % Heat transfer coefficient in W/m^2K
ql = 1; pl = -alpha*(ul-T_ext); % Left boundary
qr = 1; pr = alpha*(ur-T_ext); % Right boundary
end
function [c,f,s]=PDE(x,t,u,DuDx)
rho = 1000; % Density in kg/m^3
cp = 1000; % Specific heat in J/kgK
k = 1; % Heat transfer coefficient in W/mK
c=rho*cp;
f=k*DuDx;
s=0;
end
The resulting plot (x-axis = width of the wall (1 = left side, 3 = mid, 5 = right side), y-axis = time (1 = 5 min, 2 = 10 min, ... 12 = 60 min), z-axis = temperature) looks as follows and makes sense to me:
The derived solution tells me that a slab with
  • a width of 10 cm / 0.1 m,
  • an initial temperature of 0 °C,
  • an ambient temperature of 50 °C,
  • a heat transfer coefficient (at both sides of the wall) of 10 W/m²K,
  • a density of 1,000 kg/m³,
  • a specific heat capacity of 1,000 J/kgK,
  • a heat transfer coefficient of 1 W/mK,
will have a midplane temperature of 20,96 °C after 1 h / 3,600 s:
sol(size(sol,1),3)
I want to double-check this solution with the Heisler chart for a slab (source):
For the parameters given above, the Fourier-number is , the Biot-number is , and the reciprocal Biot-number is (note that only one half of the slab´s width (i.e. 0.05 m) must be used as a characteristic dimension for a slab).
Plugging these parameters into the Heisler chart (see blue arrow in the upper left diagram where I derived the value of 0.55) gives me
This value (27.5 °C) significantly deviates from the MATLAB solution (20,96 °C) derived above (however, the Heisler chart solution should be valid since ). Can you help me to spot my mistake and / or give me a reasonable explanation for this huge deviation?
Thanks a lot in advance,
Phil

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

回答（2 个）

Sulaymon Eshkabilov
In fact, in your matlab code solution, you have obtained the solution in the two corners T = 27.019 C (about) that is approx equal to your calcs using Heisler chart.
  1 个评论
Phil. U.
Phil. U. 2021-11-3
Sorry, I don´t understand this comment. The corner temperatures is (according to my solution)
sol(size(sol,1),1)
or
sol(size(sol,1),5)
27.027 °C (which is close to the value that you are mentioning (27.019 °C) but not identical). Anyhow, the textbook / Heisler chart approach provides a solution for the midplane temperature (and not the corner temperature). Can you therefore please elaborate why you think that the corner temeratures are of importance here?

请先登录，再进行评论。

Bill Greene
Bill Greene 2021-11-4
I referred to section 5.5 of Bergman (the source of your pdf file) and wrote a function that computes the solution the Heisler chart is supposedly based on.
I then modified your posted m-file to compare the pdepe solution with this approximate analytical solution. Here are a couple of comparison plots
And here is my modified version of your m-file
function matlabAnswers_11_3_2021
rho = 1000; % Density in kg/m^3
cp = 1000; % Specific heat in J/kgK
k = 1; % Heat transfer coefficient in W/mK
T_ext = 50; % Ambient temperature = 50 °C
hc = 10; % Heat transfer coefficient in W/m^2K
T0=0; % initial temperature
nx=11;
nx2=ceil(nx/2);
x=linspace(0,0.1,nx);
t=linspace(0,3600,20);
solAnal=heatTransferConvectionBCAnal(k, rho, cp, hc, T0, T_ext, x, t);
m = 0;
pde=@(x,t,u,DuDx) PDE(x,t,u,DuDx,rho,cp,k);
bc=@(xl,ul,xr,ur,t) BC(xl,ul,xr,ur,t,hc,T_ext);
ic = @(x) IC(x, T0);
sol=pdepe(m,pde,ic,bc,x,t); %solution
%surf(sol); % plot
figure; plot(t,sol(:,nx2),t,solAnal(:,nx2));
title 'Center temperature as a function of time'
legend('numerical', 'anal');
ylabel 'Temperature'
xlabel 'Time'
figure; plot(x,sol(end,:),x,solAnal(end,:)); grid
legend('numerical', 'anal');
ylabel 'Temperature'
xlabel 'x'
title 'Temperature at final time'
end
function u0=IC(x, T0)
u0=T0; % Initial temperature = 0°C
end
function [pl,ql,pr,qr]=BC(xl,ul,xr,ur,t,hc,T_ext)
ql = 1; pl = -hc*(ul-T_ext); % Left boundary
qr = 1; pr = hc*(ur-T_ext); % Right boundary
end
function [c,f,s]=PDE(x,t,u,DuDx,rho,cp,k)
c=rho*cp;
f=k*DuDx;
s=0;
end
function T=heatTransferConvectionBCAnal(k, rho, cp, hc, T0, T_ext, x, t)
x=x(:)';
t=t(:);
L=(x(end)-x(1))/2;
xStar=(x-L)/L;
Bi=hc*L/k;
alpha=k/(rho*cp);
f= @(x) x*tan(x)-Bi;
z=fzero(f, 1);
F0=alpha*t/L^2;
C1=4*sin(z)/(2*z + sin(2*z));
theta0=C1*exp(-z^2*F0)*cos(z*xStar);
T=(T0-T_ext)*theta0 + T_ext;
end

类别

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

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by