一维非稳态导热的数值解法
一、导热问题数值解法的认识
（一）背景
所谓求解导热问题，就是对导热微分方程在规定的定解条件下的积分求解。这样获得的解称为分析解。近100年来，对大量几何形状及边界条件比较简单的问题获得了分析解。但是，对于工程技术中遇到的许多几何形状或边界条件复杂的导热问题，由于数学上的困难目前还无法得出其分析解。另一方面，在近几十年中，随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展十分迅速，并得到日益广泛的应用。这些数值方法包括有限差分法、有限元法及边界元法等。其中，有限差分法物理概念明确，实施方法简便，本次大作业即采用有限差分法。
（二）基本思想
把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场，用有限个离散点上的值的集合来代替，将连续物理量场的求解问题转化为各离散点物理量的求解问题，将微分方程的求解问题转化为离散点被求物理量的代数方程的求解问题。
（三）基本步骤
（1）建立控制方程及定解条件。根据具体的物理模型，建立符合条件的导热微分方程和边界条件。
（2）区域离散化。用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。每一个节点都可以看成是以它为中心的一个小区域的代表，将小区域称之为元体。
（3）建立节点物理量的代数方程。建立方法主要包括泰勒级数展开法和热平衡法。 （4）设立迭代初场。 （5）求解代数方程组。
（6）解的分析。对于数值计算所获得的温度场及所需的一些其他物理量应作仔细分析，以获得定性或定量上的一些结论。对于不符合实际情况的应作修正。
二、问题及求解
（一）题目
一厚度为0.1m的无限大平壁，两侧均为对流换热边界条件，初始时两侧流体温度与壁内温度一致，已知两侧对流换热系数分别为h1=11 W/m2K、h2=23W/m2K，tf1?tf2?t0?5℃；壁的导热系数?=0.43W/mK，导温系数a=0.3437×10-6 m2/s。如果一侧的环境温度tf1突然升高为50℃并维持不变，计算在其它参数不变的条件下，平壁内温度分布及两侧壁面热流密度随时间的变化规律（用图形表示）。
1
（二）分析
为无限大平壁，可认为是一维导热问题。在某瞬时边界条件发生突变，因此是非稳态问题。因此该题模型为第三类边界条件下常物性、无内热源大平壁的一维非稳态导热问题。
（三）求解过程
1、求解域的离散
建立二维坐标，以X为横坐标，时间?为纵坐标。令空间步长?x=0.001m，时间步长
??=0.5 s。于是将计算区域划分为100等份，得到101个空间节点，编号为n=1101（由
于Matlab语句要求，编号n必须为正数，因此不能为习惯上的0两个空间节点为边界点。
本题采用温度方程的显示差分格式，为保证节点温度方程求解的稳定性，必须要求
100）。n=1和n=101这
Fo??11且Fo??。 22Bi??2a???0.3437?10?6?0.51Fo??==0.17< 22?x20.001Bi?1?Bi?2h1??x11?0.001==0.0256 ?0.43h??x23?0.001?2==0.0535
?0.4311且Fo??，能保证其稳定性。 22Bi??2100）
易知，满足Fo??2、建立节点温度差分方程（显式差分格式） 内部节点：（n?2i?1iii tn?Fo?tn?1?tn?1??1?2Fo??tn
??边界节点1：（n=1）
i?1iii t1?2Fo?t2?Bi?1tf1??1?2Bi?1Fo??2Fo??t1
??边界节点101：（n=101）
i?1iii t101?2Fo?t100?Bi?2tf2??1?2Bi?2Fo??2Fo??t101
??3、壁面温度的计算
在编写程序时，需要用到tw1在稳态情况下的值，计算如下（设在稳态下）：
q?tf1?tf21?1??h1?h2?50?5?122.634 w/m2
10.11??110.43232
又因为q?h1(tf1?tw1)， 故tw1?tf1?4、计算框图
5、程序段 ①计算温度分布
dx=0.001,dt=0.5; （设置空间步长和时间步长） B1=11*dx/0.43; B2=23*dx/0.43; F=3.437e-007*dt/dx^2;
for n=1:101; （赋初值，即初始时刻温度） t(1,n)=5; end;
m=1; （设置时间节点的初值） while 38.85-t(m,1)>0.01 （误差小于等于0.01） t(m+1,1)=2*F*t(m,2)+2*F*B1*50+(1-2*B1*F-2*F)*t(m,1); （壁面1处的边界节点方程） for n=2:100;
t(m+1,n)=(1-2*F)*t(m,n)+F*(t(m,n-1)+t(m,n+1)); （内部节点方程） end;
t(m+1,101)=2*F*t(m,100)+2*F*B2*5+(1-2*B2*F-2*F)*t(m,101); （壁面2处的边界节点方程） m=m+1;
t(m,1) （显示壁面1处的节点温度） m （显示循环的时间步长数） end;
3
q122.634?50??38.85℃ h111设定m=1时各节点的温度，即赋初值 计算m=2时n=1～101各节点的温度 计算m=3时各节点的温度 循环直至计算结果达到允许的误差范围 ②曲线生成 x=0:100;
hold on （将不同时间的温度分布曲线画在同一 plot(x,t(96000,:)) 坐标图中） plot(x,t(76800,:)) plot(x,t(64800,:)) plot(x,t(52800,:)) plot(x,t(44400,:)) plot(x,t(34800,:)) plot(x,t(26400,:)) plot(x,t(18000,:)) plot(x,t(9600,:)) plot(x,t(7200,:)) plot(x,t(6000,:)) plot(x,t(4800,:)) plot(x,t(4200,:)) plot(x,t(3600,:)) plot(x,t(3000,:)) plot(x,t(2400,:)) plot(x,t(1800,:)) plot(x,t(1200,:)) plot(x,t(720,:)) plot(x,t(480,:)) plot(x,t(240,:)) plot(x,t(120,:)) hold off
i=1:96115; q1(i)=11*(50-t(i+1,1)); i=1:96115; plot(i,q1)
i=1:96115; q2(i)=23*(t(i+1,101)-5); i=1:96115; plot(i,q2)
（四）结果及分析
用Matlab进行计算，总运行次数m=96115。①温度分布图如下：
4
1处的热流密度随时间的变化图） 2处的热流密度随时间的变化图） （画出壁面 （画出壁面 温度
节点数
各条曲线为不同时刻的温度分布曲线，分别为1min，2min，4min，6min，10min，15min，20min，25min，30min，35min，40min，50min，60min，80min，150min，220min，290min，370min，440min，540min，640min，800min，左下方为1min对应的曲线，右上方为800min对应的曲线。壁面1（高温一侧）的温度上升很快，从左向右温升速率依次降低。当导热进行足够长的时间后，温度基本呈线性分布，即与稳态导热的情况相同。由图知，tw1与计算结果相同。
tw2?q122.634?tf2??5?10.33℃，计算结果也与图中值相同。 h223 ②
q15