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

8 views (last 30 days)
Phil. U. on 2 Nov 2021
Answered: Bill Greene on 4 Nov 2021
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?
Phil

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.
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?

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
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

### Categories

Find more on Heat Transfer in Help Center and File Exchange

R2021b

### Community Treasure Hunt

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

Start Hunting!