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

23 views (last 30 days)
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

Answers (2)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 3 Nov 2021
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 Comment
Phil. U.
Phil. U. on 3 Nov 2021
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?

Sign in to comment.


Bill Greene
Bill Greene on 4 Nov 2021
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

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!