Analytical solution of 1-D transient heat transfer

Question

Nour 2024-4-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2106286-analytical-solution-of-1-d-transient-heat-transfer

编辑： Torsten 2024-4-14

Hello,

I am trying to solve a 1D transient heat equation using the analytical solution for radii from 0 to 5 cm, with a convective bounday conditions as shown in the picture. I have attached my discretization method as well to give a better insight into my problem. I have already coded the solution with out including convenction boundary condition but with a constant T_inf= 298 K.

Can you please help me edit my code to include the convection part?

-- Governing equation used

-- Convection boundary condition to be included

clc 
clear all; 
%-----------------------------Input---------------------------------------- 
% Inputs 
r0=0.05;                         % wall thickness(m)
k=0.1;                           % conductivity(W/m-K)
rho=821.42;                      % density(kg/m^3)                   
cp=1399;                         % specific heat capacity(J/kg-K)
T_ini=353;                       % initial temperature(K) 
T_inf=298;                       % external temperature(K) 
h=5;                             % heat transfer coefficient(w/m^2-K) 
t_f=1000
lampda1=1.694                   % Eigen root
C1=1.3787                       % constant
R=0.05;
T_i=353;
T_inf=298;
r=linspace(0,0.05,50);
t=linspace(0,100,t_f);
%-----------------------------Main parameters -----------------------------
% Parameters
Bi=(h*r0)/k                    % Biot number
alpha=k/(rho*cp)               % Thermal diffusivity
F0=(alpha*t)/(r0*r0)           % Fourier number
r_div=50;
epsilon=1.694;
C_1=1.3787;
                         
T=zeros(1,r_div);
u=zeros(1,r_div);
t=480;
%-----------------------------Temperature calculation ---------------------
for j=1:5
At least one END is missing. The statement beginning here does not have a matching end.
for i=1:50
u(i)=epsilon*r(i) /R;   
T(i)=T_inf+(T_i-T_inf)*C_1*exp(-epsilon^2*alpha*t/R^2)*besselj(0,u(i));
end

6 个评论
显示 4更早的评论隐藏 4更早的评论

Nour 2024-4-13

在 MATLAB Online 中打开

Thank you for the reply. Really sorry, I coppeid the approximate solution code that I wrote earlier by mistake.

The following is my exact solution that is related to the graph.

 
clear all
clc
% constant properties by demand 
density = 821.42;     % density(kg/m^3)
rho=density;
cp=1399;              % specific heat capacity(J/kg-K)
k=0.1;                % conductivity(W/m-K)
alpha=k/(rho*cp);      % Thermal diffusivity (m^2/s)
% %constant epsilon Bi for all other analysis 
% infinte cyclinder 
nd=20;
Fo=zeros(1,nd); % Fourier number
BT=zeros(1,nd);
theta0_st=zeros(1,nd);
%***********************************************************************
%solution of the trancendental equation x(Ji(x)/Jo(x))=Bi
%x=epsilon
nt=20; % number of roots required 
n=20; % number of iterations 
x=zeros(1,n);
f_x=zeros(1,n);
fp_x=zeros(1,n);
Bi=2.5;
%initial guess value
x(1)=1;
% finding x_0
xg_1=1;
xg_2=2;
x_0=(xg_1+xg_2)/2;
% 
x(1)=x_0;
roots=zeros(1,nt);
for k=1:nt
    for i=1:n-1
        
f_x(i)=x(i)*besselj(1,x(i))/besselj(0,x(i))-Bi;
fp_x(i)=x(i)+x(i)*besselj(1,x(i))^2/besselj(0,x(i))^2;
x(i+1)=x(i)-f_x(i)/fp_x(i);
roots(k)=x(i+1);
    end
    
xg_1=roots(k);
xg_2=roots(k)+5;
x_0=(xg_1+xg_2)/2;
% 
x(1)=x_0;
  
  
end
epsilon=roots;
%*************end of solution of transecndental equation******************
C=zeros(1,nt);
PFT=zeros(1,nt);
SMT=zeros(1,nd);
h=5;
R_i=0.05;
k=0.1;
Bi=h*R_i/k;
% ft_p=zeros(1,nd);
r=linspace(0,R_i,nd);
t=50000;
T_i=400;
T=zeros(1,nd);
for i=1:nd
   
Tw=298; % working temperature in K
for n=1:nt
  
   C(n)=(2/epsilon(n))*(besselj(1,epsilon(n))/((besselj(0,epsilon(n)))^2 + (besselj(1,epsilon(n)))^2));
end
%% 
%% 
for n=1:nt
     
    PFT(n)=C(n)*exp(-epsilon(n)^2*alpha*t/R_i^2)*besselj(0,epsilon(n)*r(i)/R_i);
    %PFT(n)=(C(n)/epsilon(n))*(1-exp(-epsilon(n)^2*k*t_contact(i)/(rho*Cp*R_i^2)))*besselj(1,epsilon(n));
    
end
SMT(i)=sum(PFT);
T(i)=Tw+(T_i-Tw)*SMT(i);
end

Torsten 2024-4-13

编辑：Torsten 2024-4-13

I don't understand your code from above.

T is a vector - so you only compute the r-profile for one give time t (in your case t = 50000) ?

And if you plot the temperature profile for smaller times (e.g. t = 100) , the temperature at r = R_i does not equal T_w = 298 K for all times t as you claim.

Or is the code you posted an attempt to implement the convective boundary condition ?

Nour 2024-4-13

编辑：Torsten 2024-4-14

在 MATLAB Online 中打开

yes, I am now struggling to add the convection boundary condition to the same code.

Actually I am just a beginner and these are my first trials. I was changing the time every time, save the value, and then plot them as follows.

Thank you a lot for your patience!

T_0=[399.348105708207 400.220997362980 399.823899267308 400.157085055238 399.853772767351 400.138907622977 399.866656443988 400.128678977623 399.875572112861 400.120262237576 399.884092056048 400.111071455955 399.894645024123 400.098108624175 399.911899650103 400.072636544409 399.955270636652 399.981215858787 400.234856351442 397.386534761319]

T_0 = 1x20

399.3481 400.2210 399.8239 400.1571 399.8538 400.1389 399.8667 400.1287 399.8756 400.1203 399.8841 400.1111 399.8946 400.0981 399.9119 400.0726 399.9553 399.9812 400.2349 397.3865

<mw-icon class=""></mw-icon>

T_100=[399.999999230084 400.000000268964 399.999999783514 400.000000191480 399.999999826333 400.000000157780 399.999999858242 400.000000124585 399.999999894301 400.000000083422 399.999999873225 399.999997772070 399.999945739678 399.999098778850 399.989556948428 399.914856521539 399.505742260982 397.923154547764 393.537920248508 384.648819048777]

T_100 = 1x20

400.0000 400.0000 400.0000 400.0000 400.0000 400.0000 400.0000 400.0000 400.0000 400.0000 400.0000 400.0000 399.9999 399.9991 399.9896 399.9149 399.5057 397.9232 393.5379 384.6488

<mw-icon class=""></mw-icon>

T_600=[399.999883833492 399.999829774393 399.999604832226 399.998958507437 399.997222201791 399.992765755409 399.981860358590 399.956501011340 399.900598983991 399.783965784248 399.553876193179 399.124915245478 398.369306722075 397.111634660331 395.132978110469 392.188893420349 388.042627631852 382.509713544762 375.504394130675 367.074828552094]

T_600 = 1x20

399.9999 399.9998 399.9996 399.9990 399.9972 399.9928 399.9819 399.9565 399.9006 399.7840 399.5539 399.1249 398.3693 397.1116 395.1330 392.1889 388.0426 382.5097 375.5044 367.0748

<mw-icon class=""></mw-icon>

T_1800=[399.190501334017 399.148212305607 399.017444481274 398.786457266990 398.435582809434 397.937169878682 397.255683223726 396.348111646569 395.164846259286 393.651172628276 391.749473812938 389.402166358795 386.555293647069 383.162591487097 379.189734947228 374.618391275875 369.449659021043 363.706481627073 357.434690829772 350.702457235095]

T_1800 = 1x20

399.1905 399.1482 399.0174 398.7865 398.4356 397.9372 397.2557 396.3481 395.1648 393.6512 391.7495 389.4022 386.5553 383.1626 379.1897 374.6184 369.4497 363.7065 357.4347 350.7025

<mw-icon class=""></mw-icon>

T_3600=[390.805227064895 390.675779000604 390.286017493707 389.631786740209 388.706477511719 387.501491077936 386.006865145261 384.212038201660 382.106722046482 379.681846109104 376.930531814971 373.849051217311 370.437721882964 366.701690070457 362.651556927952 358.303807952486 353.681014246139 348.811784878551 343.730462361549 338.476567080242]

T_3600 = 1x20

390.8052 390.6758 390.2860 389.6318 388.7065 387.5015 386.0069 384.2120 382.1067 379.6818 376.9305 373.8491 370.4377 366.7017 362.6516 358.3038 353.6810 348.8118 343.7305 338.4766

<mw-icon class=""></mw-icon>

T_18000=[320.784154829704 320.738258431128 320.600846434660 320.372748583693 320.055341615340 319.650540052191 319.160783399654 318.589019837082 317.938686515191 317.213686595838 316.418363192915 315.557470394871 314.636141569840 313.659855173590 312.634398298245 311.565828215802 310.460432184841 309.324685801294 308.165210184555 306.988728298620]

T_18000 = 1x20

320.7842 320.7383 320.6008 320.3727 320.0553 319.6505 319.1608 318.5890 317.9387 317.2137 316.4184 315.5575 314.6361 313.6599 312.6344 311.5658 310.4604 309.3247 308.1652 306.9887

<mw-icon class=""></mw-icon>

T_50000=[298.890663675783 298.888869371201 298.883497300502 298.874579920011 298.862171081446 298.846345670453 298.827199104151 298.804846691191 298.779422858871 298.751080252716 298.719988714898 298.686334148649 298.650317276681 298.612152302323 298.572065482773 298.530293624494 298.487082511280 298.442685276008 298.397360727463 298.351371643918]

T_50000 = 1x20

298.8907 298.8889 298.8835 298.8746 298.8622 298.8463 298.8272 298.8048 298.7794 298.7511 298.7200 298.6863 298.6503 298.6122 298.5721 298.5303 298.4871 298.4427 298.3974 298.3514

<mw-icon class=""></mw-icon>

figure

r = linspace(0,0.05,20)

r = 1x20

0 0.0026 0.0053 0.0079 0.0105 0.0132 0.0158 0.0184 0.0211 0.0237 0.0263 0.0289 0.0316 0.0342 0.0368 0.0395 0.0421 0.0447 0.0474 0.0500

<mw-icon class=""></mw-icon>

plot(r,T_0,'LineWidth',2)

hold on

plot(r,T_100,'LineWidth',2)

hold on

plot(r,T_600,'LineWidth',2)

hold on

plot(r,T_1800,'LineWidth',2)

hold on

plot(r,T_3600,'LineWidth',2)

hold on

plot(r,T_50000,'LineWidth',2)

%legend('0 sec','100 sec','600 sec','1800 sec','3600 sec','50000')

xlabel('r (m)')

ylabel('T (K)')

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-4-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2106286-analytical-solution-of-1-d-transient-heat-transfer#answer_1440821

编辑：Torsten 2024-4-14

在 MATLAB Online 中打开

I'd compare with the solution from "pdepe":

m = 1;

Ri = 0.05;

nd = 20;

r = linspace(0,Ri,nd);

t = [0 100 600 1800 3600 50000];

sol = pdepe(m,@heatcyl,@heatic,@heatbc,r,t);

u = sol(:,:,1);

plot(r,u(1:6,:),'LineWidth',2)

grid on

function [c,f,s] = heatcyl(x,t,u,dudx)

rho = 821.42;

cp = 1399;

k = 0.1;

c = rho*cp;

f = k*dudx;

s = 0;

end

%----------------------------------------------

function u0 = heatic(x)

u0 = 400;

end

%----------------------------------------------

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)

pl = 0; %ignored by solver since m=1

ql = 0; %ignored by solver since m=1

alpha = 5;

Ta = 298.0;

pr = alpha*(ur-Ta);

qr = 1;

end

%----------------------------------------------

But your results look ok at a quick glance.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Analytical solution of 1-D transient heat transfer

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Analytical solution of 1-D transient heat transfer

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

6 个评论
显示 4更早的评论隐藏 4更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论