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

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

Answer 1

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

[Phil. U.]

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?

Answer 2

[Bill Greene]

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