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

12 views (last 30 days)

Show older comments

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)

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

##### 0 Comments

### Answers (2)

Sulaymon Eshkabilov
on 3 Nov 2021

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

##### 0 Comments

### See Also

### Community Treasure Hunt

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

Start Hunting!