Why does solving the heat equation with MATLAB (pdepe) yield a completely different result than the Heisler chart (analytical solution)?
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788105/image.png)
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/787705/image.png)
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).
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788110/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788115/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788120/image.png)
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 ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788125/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788125/image.png)
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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/787730/image.png)
Thanks a lot in advance,
Phil
0 Commenti
Risposte (2)
Sulaymon Eshkabilov
il 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.
Bill Greene
il 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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/789900/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/789904/image.png)
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 Commenti
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!