传热学计算机实习指导书本指导书是为配合本科生传热学课中计算机应用方面的教学而编写的。

应用计算机解决工程实际问题，是现代工程技术人员所必备的技能。

在传热学课程中引入计算机实习的目的，是使学生初步掌握用计算机求解传热问题的技能，从而提高学生应用计算机解决工程实际问题的能力。

大量的传热问题能够用计算机求解。

研究如何用计算机求解传热问题的专门知识数值传热学（或称计算传热学）已经发展成了传热学的一个分支学科。

传热学课中所涉及的只是数值传热学的初步知识。

因此，本次计算机实习也仅仅是作为数值传热学的入门。

本指导书给出了三个练习题及相应的算法。

这三个练习题分别涉及了一维稳态导热、二维稳态导热和一维非稳态导热。

要求学生在掌握问题的数值计算方法的基础上，独立编写计算机程序并用所编的程序计算出这三个练习题的数值结果。

1 练习题一：一维稳态导热的数值计算1．1 物理问题图1示出了一个等截面直肋，处于温度t ∞=80℃的流体中。

肋表面与流休之间的对流换热系数为h=45W/m 2.℃，肋基处温度t w =300℃，肋端绝热。

肋片由铝合金制成，其导热系数为λ=110W/m ℃，肋片厚度为δ=0.01m ，高度为H=0.1m 。

试计算肋内的温度分布及肋的总换热量。

1． 2数学描述及其解析解引入无量纲过余温度∞∞--=t t t t w θ，则以无量纲温度θ描述的肋片导热微分方程及其边界条件为:0,1,00222=∂∂====-xH x x m dxd θθθθ (1-2)(1-1)(1-3)其中Ahpm λ=(其中符号含义与教科书杨世铭陶文铨编著《传热学》相同，以下同)。

上述数学模型的解析解为:()()[]()()()m H th t t mhpm H ch H x m ch t t t t w w ∞∞∞-=-⋅-=-φ (1-4)按式(1-4)计算得到的在肋内各点的温度由表1给出。

1． 3数值离散1．3．1区域离散在对方程(1-1)~(1-3)进行数值离散之前,应首先进行计算区域的离散。

计算区域的离散如图1所示，总节点数取N 。

1．3．2微分方程的离散由于方程(1-1)在计算区域内部处处成立，因而对图1所示的各离散点亦成立。

对任一节点i 有：0222=-⎪⎪⎭⎫ ⎝⎛i im dx d θθ 用θ在节点i 的二阶差分代替θ在节点i 的二阶导数,得:022211=-∆+--+i i i i m xθθθθ整理上式成迭代形式:()112221-++∆+=i i i x m θθθ (i=2,3,…,N-1) (1-5)1．3．3边界条件离散上面得到的离散方程式(1-5)，对所有内部节点都成立，因此每个内部节点都可得出一个类似的方程。

事实上，式(1-5)表达的是一个代数方程组。

但这个方程组的个数少于未知数i θ (i=1,2, ……,N)的个数。

因此，还需要根据边界条件补充进两个方程后代数方程组才封闭。

左边界(x=0)为第一类边界条件，温度为已知，因此可以根据式(1-2)直接补充一个方程为：11==-=∞w w t t θθ 右边界为第三类边界条件，由图1中边界节点N 的向后差分来代替式(1-3)中的导数，得：01=∆--xN N θθ将此式整理为迭代形式，得：1-=N N θθ1．3．4最终的离散格式()111221211--+=+∆+===N N i i i w xm θθθθθθθ (i=2,3,…,N-1) （1-6）1．3．5代数方程组的求解及其程序代数方程组有各种求解方法，较为有效而简便的方法是高斯-赛德尔迭代方法。

式（1-6）已给出了代数方程组的迭代形式。

在实际计算中，应首先假定一个温度场的初始分布，即给出各节点的温度初值：)1(,...,,0100201==w Nθθθθθ将这些初值代入方程组（1-6）中进行迭代计算，直至收敛。

假如第K 步迭代已完成，即KNK K θθθ,...,,21为已知，则K+1次迭代的计算式为： ()1111112211121+-++-+++=+∆+==K N K N K i K i K iwK xm θθθθθθθ (i=2,3,…,N-1) （1-7） 根据式（1-7）编写程序的工作由学生自行完成。

计算结果可与解析解比较。

2 练习题二：二维稳态导热的数值计算2． 1 物理问题图2示出了一矩形区域，其边长L=W=1，假设区域内无内热源，导热系数为常数，三个边温度为T 1=0，一个边温度为T 2=1，求该矩形区域内的温度分布。

2． 2 数学描述对上述问题的微分方程及其边界条件为：1,10,00,10,0021112222=============∂∂+∂∂T T y T T y T T x T T x y Tx T (2-2)作为参考，以下给出该问题的解析解:⎪⎭⎫⎝⎛⋅⎪⎭⎫ ⎝⎛⋅⎪⎭⎫ ⎝⎛⋅--=--∑∞=W L n sh y L n sh x L n n T T T T n n ππππsin )1(121121 (2-3) 表2列出了由式(2-3)计算得到的在平面区域内各不同位置的温度值。

2．3 数值离散2.3.1区域离散区域离散如图2所示，x 方向总节点数为N ，y 方向总节点数为M ，区域内任一节点用i,j 表示。

2.3.2方程的离散对于图2中所有的内部节点方程(2-1)都适用，因此可写为：0,22,22=⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂ji j i y t x t用i,j 节点的二阶中心差分代替上式中的二阶导数，得：(2-1)02221,,1,2,1,,1=∆+-+∆+--+-+yT T T xT T T j i j i j i ji j i j i上式整理成迭代形式：()()()()1,1,222,1,1222,22-+-++∆+∆∆++∆+∆∆=j i j i j i j i j i T T yx x T T y x y T (i=2,3,…,N-1),(j=2,3,…,M-1)补充四个边界上的第一类边界条件得：1,1T T j =(j=1,2,…,M) 1,T T j N = (j=1,2,…,M) 11,T T i = (i=1,2,…,N) 2,T T M i =(i=1,2,…,N)2.4计算程序计算程序由学生自行完成。

3 练习题三：一维非稳态导热的数值计算非稳态导热问题由于有时间变量，其数值计算出现了一些新的特点。

在非稳态导热微分方程中，与时间因素相关的非稳态项是温度对时间的一阶导数，这给差分离散带来了新的特点。

由于这个特点，可以采用不同的方法构造差分方程，从而得到几种不同的差分格式，即所谓的显式、隐式和半隐式。

我们仍然从一个具体问题出发来研究非稳态导热问题的数值计算。

3．1 问题一块无限大平板（如图3所示），其一半厚度为L=0.1m ，初始温度T 0=1000℃，突然将其插入温度T ∞=20℃的流体介质中。

平板的导热系数λ=34.89W/m ℃，密度ρ=7800kg/m 3，比热c=0.712310⨯J/kg ℃，平板与介质的对流换热系数为h=233W/m 2.℃，求平板内各点的温度分布。

3．2 数学描述由于平板换热关于中心线是对称的，仅对平板一半区域进行计算即可。

坐标x 的原点选在平板中心线上，因而一半区域的非稳态导热的数学描述为：()∞-=∂∂-==∂∂===∂∂=∂∂T T h xTL x x Tx T T x Ta T λττ,0,0,0022 该数学模型的解析解为:()02cos cos sin sin 210F n n nn n n n e L x T T T T μμμμμμ-∞=∞∞⎪⎭⎫ ⎝⎛+-+=∑（3-5）其中20La F τ=,n μ为方程i B ctg /μμ=的根，λhLB i =。

表3给出了在平板表面(x=L)处由式(3-5)计算得到的不同时刻的温度值。

3．3数值离散3．3．1计算区域的离散一维非稳态导热指的是空间坐标是一维的。

若考虑时间坐标，则所谓的一维非稳态导热实(3-1) (3-2) (3-3) (3-4)际上是二维问题(见图4)，即：有时间坐标τ和空间坐标x 两个变量。

但要注意，时间坐标是单向的，就是说，前一时刻的状态会对后一时刻的状态有影响，但后一时刻的状态却影响不到前一时刻，图4示出了以x 和τ为坐标的计算区域的离散，时间从τ=0开始，经过一个个时层增加到K 时层和K+1时层。

3．3．2微分方程的离散对于i 节点，在K 和K+1时刻可将微分方程(3-1)写成下面式子：122122++⎪⎪⎭⎫⎝⎛∂∂=⎪⎭⎫ ⎝⎛∂∂⎪⎪⎭⎫⎝⎛∂∂=⎪⎭⎫⎝⎛∂∂K iK iKi Ki x T a T x T a T ττ将式(3-6)~(3-8)的左端温度对时间的偏导数进行差分离散为：ττττ∆-=⎪⎭⎫⎝⎛∂∂∆-=⎪⎭⎫⎝⎛∂∂+++Ki K iK iKi K iK i T T T T T T 111观察式(3-8)和(3-9)，这两个式子的右端差分式完全相同，但在两个式子中却有不同含义。

对式(3-8)，右端项相对i 点在K 时刻的导数Ki T ⎪⎭⎫⎝⎛∂∂τ是向前差分。

而在式(3-9)中，右端项是I 点在K+1时刻的导数1+⎪⎭⎫⎝⎛∂∂K iT τ的向后差分。

将式(3-8)和(3-9)分别代入式(3-6)和(3-7)，并将式(3-6)和(3-7)右端关于x 的二阶导数用相应的差分代替，则可得到下列显式和隐式两种不同的差分格式： 显式:Ki K i K i K i fT T f fT T 111)21(-+++-+= (3-10)(K=0,1,2, ………, i=2,3,…,N-1)全隐式：()KiK i K i K i T fT fT fT +++=+-+++11111211 (3-11) (K=0,1,2, …………… i=2,3,…,N-1)以上两式中的2x a f ∆∆=τ。

从式(3-10)可见，其右端只涉及K 时刻的温度，当从K=0（即τ=0时刻）开始计算时，在K=0时等号右端都是已知值，因而直接可计算出K=1时刻各点的温度。

由K=1时刻的各点的温度值，又可以直接利用式(3-10)计算K=2时刻的各点的温度，这样一个时层一个时层的往下推，各时层的温度都能用式(3-10)直接计算出来，不要求解代数方程组。

而对于式(3-11) 等号右端包含了与等号左端同一时刻但不同节点的温度，因而必须通过求解代数方程组才能求得这些节点的温度值。

3．3．3边界条件的离散(3-6)(3-7)(3-8)(3-9)对于式(3-3)和(3-4)所给出的边界条件，可以直接用差分代替微分，也可以用元体平衡法给出相应的边界条件，亦有显式和隐式之分。

通常，当内部节点采用显式时，边界节点也用显式离散；内部节点用隐式时，边界节点亦用隐式。