传热学数值计算

导热问题的数值解法1 、重点内容：① 掌握导热问题数值解法的基本思路；② 利用热平衡法和泰勒级数展开法建立节点的离散方程。

2 、掌握内容：数值解法的实质。

3 、了解内容：了解非稳态导热问题的两种差分格式及其稳定性。

由前述3 可知，求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。

但是，对于工程中几何形状及定解条件比较复杂的导热问题，从数学上目前无法得出其分析解。

随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展得十分迅速，并得到广泛应用，并形成为传热学的一个分支——计算传热学（数值传热学），这些数值解法主要有以下几种：（１）有限差分法（ 2 ）有限元方法（ 3 ）边界元方法数值解法能解决的问题原则上是一切导热问题，特别是分析解方法无法解决的问题。

如：几何形状、边界条件复杂、物性不均、多维导热问题。

分析解法与数值解法的异同点：相同点：根本目的是相同的，即确定① t=f(x ，y ，z) ；②。

不同点：数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散点的数值。

§4-1 导热问题数值求解的基本思想及内节点离散方程的建立实质对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。

该方法称为数值解法。

这些离散点上被求物理量值的集合称为该物理量的数值解。

2 、基本思路：数值解法的求解过程可用框图4-1 表示。

由此可见：1 ）物理模型简化成数学模型是基础；2 ）建立节点离散方程是关键；3 ）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。

一数值求解的步骤如图4-2 （a ），二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤如下：1 建立控制方程及定解条件控制方程：是指描写物理问题的微分方程针对图示的导热问题，它的控制方程（即导热微分方程）为：（a ）边界条件：x=0 时，x=H 时，当y=0 时，当y=W 时，区域离散化（确立节点）用一系列与坐标轴平行的网格线把求解区域划分成若干个子区域，用网格线的交点作为需要确定温度值的空间位置，称为节点( 结点) ，节点的位置用该节点在两个方向上的标号m ，n 表示。

Nu = 2+0.6(Re^1/2)(Pr^1/3) 。

F=Q/kK*△tm F 是换热器的有效换热面积。

Q 是总的换热量。

k 是污垢系数一般取0.8-0.9K。

是传热系数。

△tm 是对数平均温差。

传热学三种传热方式可以分开学。

传热学相较于理论力学，工程热力学，流体力学而言还是比较简单的，一般大学生掌握了高等数学完全可以自学的。

学习传热学必须有耐心，了解几种换热方式和常见的几个常数公式（努谢尔特数、格拉晓夫数、伯努利常数，傅里叶常数，而且常常推导下几个常用常数公式间的关系，你会惊奇地发现他们其实不少是远亲的），其实解决传热学问题绝大多数都是在和导热系数较劲，有时候是直接涉及。

扩展资料：
在热对流方面，英国科学家牛顿于1701年在估算烧红铁棒的温度时，提出了被后人称为牛顿冷却定律的数学表达式，不过它并没有揭示出对流换热的机理。

传热学作为学科形成于19世纪。

1804年，法国物理学家毕奥在热传导方面得出的平壁导热实验结果是导热定律的最早表述。

稍后，法国的傅里叶运用数理方法，更准确地把它表述为后来称为傅里叶定律的微分形式。

1860年，基尔霍夫通过人造空腔模拟绝对黑体，论证了在相同温度下以黑体的辐射率(黑度)为最大，并指出物体的辐射率与同温度下该物体的吸收率相等，被后人称为基尔霍夫定律。

Fe aE De 2
Thermal
Fw Fe aP Dw De aW aE Fe Fw 2 2
2、对方程的几点说明 由于连续性，Fe=Fw， aP aW aE（只是在流 场满足连续性条件时才具有这一性质）； 方程 aP P aE E aW W 隐含着分段线性分布的含
讨论只有对流项和扩散项存在时的一维稳态问题，控 制方程为：
u j ( ) S x j x j x j
d d d u ( ) dx dx dx
d 连续方程： u 0 dx
u const
任务：导出相应方程的离散化形式
义，也是熟知的中心差分格式（用左右节点值表示
界面上的值以及界面上的导数值）； 方程必须遵守四项基本法则，否则会产生灾难性的 结果。
2018-11-24
太 原 理 工 大 学
7 /70
Thermal
例如：设 De Dw 1， Fe Fw 4 若E、W给定，即可由离散方程求得P 。
即， F 2D时，有可能使 aE 或aW 为负 产生不切实际的结果
2018-11-24
太 原 理 工 大 学
8 /70
Thermal
这就是中心差分格式求解对流换热问题时仅限于低
Fw Fe Re(低的F/D)的原因 . aP Dw De aW aE Fe Fw 2 2
原通用方程可改写为
u j ( )S x j x j x j
对于已知的ρ、uj、Γ及S（常量）的分布，任何解及
+c 将同时满足方程，故系数和的法则仍然适用。
2018-11-24
太 原 理 工 大 学

导热问题的数值求解方法数值解法的基本思想是用空间和时间区域内有限个离散点（称为节点）上温度的近似值，代替物体内实际的连续温度分布，然后由导热方程和边界条件推导出各节点温度间的相互关系的代数方程组（称为离散方程），求解此方程组，得到节点上的温度值，此即物体中温度场的解。

只要节点分布的足够稠密，数值解就有足够的精度。

求解导热问题的数值方法有有限差分法及有限元法，近几年又发展了边界元法和有限分析法。

数值方法适用于求解各种导热问题，不管物体的几何形状有多复杂，不管线性或非线性问题，都能使用。

由于计算机的飞速发展，计算技术软件发展也很快，数值方法的的地位越来越重要。

1 数值求解的基本思路及稳态导热内节点离散方程的建立一、 解法的基本思路1、基本思路：数值解法的求解过程可用框图4-1表示。

由此可见：1）物理模型简化成数学模型是基础；2）建立节点离散方程是关键；3）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。

二、稳态导热中位于计算区域内部的节点离散方程的建立方法1、基本方法方法：①泰勒级数展开法；②热平衡法。

1）泰勒级数展开法如图4-3所示，以节点(m,n)处的二阶偏导数为例，对节点(m+1,n)及(m-1,n)分别写出函数t 对(m,n)点的泰勒级数展开式：对(m+1,n)：+∂∂∆+∂∂∆+∂∂∆+∂∂∆+=+444333,222,,,12462x t x x t x x t x x t x t t n m n m n m n m （a ）对（m-1,n ）：+∂∂∆+∂∂∆-∂∂∆+∂∂∆-=-444333,222,,,12462x t x x t x x t x x t xt t n m n m n m n m （b ）（a ）+（b ）得： +∂∂∆+∂∂∆+=+-+444,222,,1,1122x t x x t x t t t n m n m n m n m 变形为n m x t,22∂∂的表示式得：n m x t,22∂∂)(0222,1,,1x x t t t nm n m n m ∆+∆+-=-+ 上式是用三个离散点上的值计算二阶导数n m x t ,22∂∂的严格表达式，其中：)(02x ∆―― 称截断误差，误差量级为2x ∆在数值计算时，用三个相邻节点上的值近似表示二阶导数的表达式即可，则相应的略去)(02x ∆。

二、主要的数值解法
有限差分法（FDM），主要介绍方法 有限元法（FEM） 边界元法（BEM）
传热学 Heat Transfer
三、数值解法的基本步骤
传热学 Heat Transfer
以一个二维稳态导热问题为例，介绍数值求解导热问题的具体 过程，重点是节点离散方程的建立和代数方程组的迭代求解。
（1）物理问题 二维矩形域稳态无内热源，常物性的导热问题 y
传热学 Heat Transfer
4-1 导热问题数值求解的基本思想
一、数值解法的本质
数值解法是用物理问题所 涉及的空间和时间区域内有 限个离散点（称为节点）的 物理量近似值来代替物体内 实际连续的物理量分布，将 连续物理量分布函数的求解 问题转化为各节点物理量值 的求解问题。
传热学 Heat Transfer
传热学 Heat Transfer
传热学 Heat Transfer
主讲老师：王舫 适用专业：能源与动力工程专业
传热学 Heat Transfer
第四章 热传导问题的数值解法
§4-0 引言
1 求解导热问题的三种基本方法： (1) 理论分析 (2)实验 (3)数值计算
2 三种方法的特点 3 三种方法的基本求解过程
传热学heattransfer?本章教学内容41导热问题数值求解基本思想42内节点离散方程的建立43边界节点离散方程的建立及代数方程的求解44非稳态导热问题的数值解法传热学heattransfer41导热问题数值求解的基本思想一数值解法的本质数值解法是用物理问题所涉及的空间和时间区域内有限个离散点称为节点的物理量近似值来代替物体内实际连续的物理量分布将连续物理量分布函数的求解问题转化为各节点物理量值的求解问题
传热学 Heat Transfer

一维非稳态导热的数值解法一、导热问题数值解法的认识（一）背景所谓求解导热问题，就是对导热微分方程在规定的定解条件下的积分求解。

这样获得的解称为分析解。

近100年来，对大量几何形状及边界条件比较简单的问题获得了分析解。

但是，对于工程技术中遇到的许多几何形状或边界条件复杂的导热问题，由于数学上的困难目前还无法得出其分析解。

另一方面，在近几十年中，随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展十分迅速，并得到日益广泛的应用。

这些数值方法包括有限差分法、有限元法及边界元法等。

其中，有限差分法物理概念明确，实施方法简便，本次大作业即采用有限差分法。

（二）基本思想把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场，用有限个离散点上的值的集合来代替，将连续物理量场的求解问题转化为各离散点物理量的求解问题，将微分方程的求解问题转化为离散点被求物理量的代数方程的求解问题。

（三）基本步骤（1）建立控制方程及定解条件。

根据具体的物理模型，建立符合条件的导热微分方程和边界条件。

（2）区域离散化。

用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。

每一个节点都可以看成是以它为中心的一个小区域的代表，将小区域称之为元体。

（3）建立节点物理量的代数方程。

建立方法主要包括泰勒级数展开法和热平衡法。

（4）设立迭代初场。

（5）求解代数方程组。

（6）解的分析。

对于数值计算所获得的温度场及所需的一些其他物理量应作仔细分析，以获得定性或定量上的一些结论。

对于不符合实际情况的应作修正。

二、问题及求解（一）题目一厚度为0.1m 的无限大平壁，两侧均为对流换热边界条件，初始时两侧流体温度与壁内温度一致，1205f f t t t ===℃；已知两侧对流换热系数分别为h 1=11 W/m 2K 、h 2=23W/m 2K ，壁的导热系数λ=0.43W/mK ，导温系数a=0.3437×10-6 m 2/s 。

2 h∆y
=
h∆y 2 2 + λ t m ,n+1 + t m ,n−1 + 2(t m−1,n + Bi ⋅ t f ) 2(2 + Bi )
9
λ
tf
3）方程组的求解 ） 对应每个未知量（一个节点温度） 一条方程 一条方程（ 对应每个未知量（一个节点温度）→一条方程（一 个节点方程） 方程组有唯一解 个节点方程）→方程组有唯一解 解法：（ ：（1） 解法：（ ）矩阵法 （2）迭代法：高斯 赛德尔迭代 ）迭代法：高斯-赛德尔迭代
λ∆y
t m−1,n − t m ,n ∆x ∆x t m ,n−1 − t m ,n +λ + ∆y ⋅ h(t f − t m ,n ) = 0 2 ∆y ∆y ∆x t m ,n+1 − t m ,n +λ ∆y 2
如取正方形网络
∆x = ∆y
上式简化为： 上式简化为：
t m ,n =
2t m−1,n + t m ,n+1 + t m ,n−1 +
对于非稳态导热问题，除了空间上进行网格划分外， 对于非稳态导热问题，除了空间上进行网格划分外， 还要把时间分割成许多间隔。 还要把时间分割成许多间隔。
4
2）有限元法 ）
把整个求解域离散成为有限个子域， 把整个求解域离散成为有限个子域，每一子域内运 用变分法， 用变分法，即利用与原问题中微分方程相等价的变 分原理来进行推导，从而使原问题的微分方程组退 分原理来进行推导， 化到代数联立方程组，得到数值解。 化到代数联立方程组，得到数值解。 有限元法和差分法都是常用的数值计算方法， 有限元法和差分法都是常用的数值计算方法，差分 法计算模型对于不规则的几何形状难以应用。 法计算模型对于不规则的几何形状难以应用。有限 元法能够很好地适应复杂的几何形状、 元法能够很好地适应复杂的几何形状、复杂的材料 特性和复杂的边界条件。 特性和复杂的边界条件。

例题 2
某逆流套管式换热器，刚投入工作时的运行参数为：
t1 360C,t1 300C,t2 30C,t2 200C 已知 qm1cp1=2500 W/K, k = 800 W/(m2.K)。运行一年后发现， 在 qm1cp1，qm2cp2，及入口温度不变的情况下，由于积垢使 得冷流体只能加热到162℃. 确定此情况的
(d)由式 kAtm 求出换热量 ;
(e)比较 与 ，如果相差较大，再重新假设流体出口温度， 重复上述计算，直到满意为止。
10.5 传热的强化与削弱（自学）
传热工程技术是根据现代工业生产和科学实践的需要而发展 起来的科学与工程技术，其主要任务是按照工业生产和科学 实践的要求来控制和优化热量传递过程。
和换热量 。
计算步骤：
(a) 先假设一个流体的出口温度，热平 衡方程式求出换热量 和 另一个流体的 出口温度；
kAtm
qm1cp1 t1 t1
qm2cp2 t2 t2
(b) 根据流体的进、出口4个温度求平均温差 tm ；
(c) 计算换热面两侧的表面传热系数 h1, h，2 进而求得总传热系数k；
tm (tm )ctf
教程中图10-23~10-26分别给出了管壳式换热器和交叉流式 换热器的 。
值取决于无量纲参数 P和 R： P tc tc , th tc
R th th tc tc
式中：下标h、c分别表示两种流体，上角标 ` 表示进口，`` 表示出口，图表中均以P为横坐标，R为参量。
1. 通过平壁的传热
K
1
1
1
h1 h2
KAt
说明： (1) h1和h2的计算；（2）如果计及辐射，换热系 数应该采用等效换热系数(总表面传热系数)

传热学是研究热量如何通过传导、对流和辐射进行传递的学科。

在传热学中，有一些常用的表达式，如Nu数、Re数、Pr数和Gr数，它们分别表示不同的传热特性。

本文将对这些表达式的含义进行详细的介绍。

一、 Nu数的含义Nu数是Nusselt数的缩写，它表示流体中的对流传热能力。

Nu数的计算公式为：Nu = hL/k其中，h是对流传热系数，L是特征长度，k是流体的导热系数。

Nu 数是对流传热与导热的比值，它越大表示对流传热能力越强，反之则表示导热能力较强。

Nu数的大小与流体的性质、流动状态和流体与固体界面的情况有关。

二、 Re数的含义Re数是Reynolds数的缩写，它表示流体的流动状态。

Re数的计算公式为：Re = ρVD/μ其中，ρ是流体密度，V是流体流速，D是特征长度，μ是流体的动力黏度。

Re数反映了流体的惯性力与黏性力之间的比值，它的大小决定了流体的流动状态，当Re数较小时，流体呈现层流状态，当Re数较大时，流体呈现湍流状态。

Re数对流体的流动特性以及传热和传质过程都有重要影响。

三、 Pr数的含义Pr数是Prandtl数的缩写，它表示流体的热传导能力与动力黏度之间的比值。

Pr数的计算公式为：Pr = μCp/κ其中，μ是动力黏度，Cp是定压比热，κ是流体的导热系数。

Pr数越大，流体的热传导能力越强，而动力黏度的影响越小，反之则动力黏度的影响越大。

Pr数的大小对对流传热和边界层的发展都有重要影响。

四、 Gr数的含义Gr数是Grashof数的缩写，它表示自然对流传热的能力。

Gr数的计算公式为：Gr = gβΔTL^3/ν^2其中，g是重力加速度，β是体积膨胀系数，ΔT是温度差，L是特征长度，ν是运动黏度。

Gr数的大小决定了自然对流传热的强弱，当Gr数较大时，自然对流传热能力越强，当Gr数较小时，传热能力较弱。

总结在传热学中，Nu数、Re数、Pr数和Gr数是常用的表达式，它们分别代表了对流传热能力、流体流动状态、热传导能力与动力黏度之间的比值以及自然对流传热的能力。

传热学
数值计算
4、一无限大平壁厚度为 0.3m，其导热系数为λ =36.4 W/ (m·K)。平壁两侧表面均给定 为第三类边界条件，即 h 1 =60W/ (m 2 ·K)，t f 1 =25℃；h 2 =300W/ (m 2 ·K)，t f 2 =215℃。 当平壁中具有均匀的内热源 q v =2 10 W / m 时，试计算沿平壁厚度的稳态温度分布。
k
0
k k 1
11
k 2Fo (t 10 t f Bi ) (1 2Fo 2Bi Fo) t11 t11
0
四、
计算过程
⑴ 设定初值：
t (1~11) 35 ℃;
36W / m k ;
Bi h / ;
根据不同的 Fo 计算Δ τ ：
qv x 2

)/2
h2 x
h2 x TRB


qv x 2

) /(1

)
|T[i]-t[i]|<=EPS
NO
IT=IT+1
YES
打印“t[i]” ， “IT”
YES
IT>K
NO
停止
⑷ 程序与计算结果 #include"iostream.h" #include"iomanip.h" #include"math.h" #include"stdio.h" #define N 16 void main() { int M,i,IT,flag; //定义节点个数
计算结果：
各节点温度： 节点 1 2 411.24 10 422.65 3 420.33 11 414.23 4 427.22 12 403.62 5

传热系数计算公式传热系数是指单位时间内，单位面积的热量与温度差之间的比值。

它描述了物体传热的快慢程度，是传热过程的重要参数。

根据传热形式的不同，传热系数有不同的计算公式。

当传热方式是传导传热时，我们可以使用傅立叶定律计算传热系数。

傅立叶定律表示，通过单位面积传导的热量与温度梯度之间成正比，可以表示为：q = -kA(dT/dx)其中，q表示单位时间内传导的热量，k表示传导热系数，A表示传热面积，(dT/dx)表示温度梯度。

传导热系数k可以通过实验测量得到，也可以通过材料的性质计算得到。

当传热方式是对流传热时，我们可以使用庙卡定律计算传热系数。

庙卡定律表示，对流传热的热流密度与温度差之间成正比，可以表示为：q=hAΔT其中，q表示单位时间内传导的热量，h表示对流传热系数，A表示传热面积，ΔT表示温度差。

对流传热系数h可以通过实验测量得到，也可以通过流体的性质和流动情况计算得到。

对于辐射传热方式，我们可以使用斯特藩-玻尔兹曼定律计算传热系数。

斯特藩-玻尔兹曼定律表示，辐射传热的热流密度与温度之间成正比，可以表示为：q=εσA(T1^4-T2^4)其中，q表示单位时间内传导的热量，ε表示表面发射率，σ表示斯特藩-玻尔兹曼常数，A表示传热面积，T1和T2分别表示辐射体和接受体的温度。

表面发射率ε可以通过表面的材料性质计算得到。

总的来说，传热系数的计算公式和传热方式有关。

一般情况下，物体传热的方式是由传导、对流和辐射三种方式共同作用，因此传热系数是这三种传热系数的总和：h总=h传导+h对流+h辐射其中h传导、h对流和h辐射分别表示传导、对流和辐射传热系数。

在实际应用中，为了保持传热系数的连续性，可以通过换热系数来表示总的传热能力。

传热系数的计算是热力学和传热学中的重要内容，它影响着热工设备和系统的设计和运行。

通过合理地计算传热系数，可以提高热工设备的传热效率，减少能源损失，提高能源利用率。

因此，准确计算传热系数对于工程实际具有重要意义。

1.2 传热与流动问题数值计算的基本思想及应用举例
1.2.1 数值解基本思想（基于连续介质假设）
把原来在空间与时间坐标中连续的物理量的场 （如速度场、温度场、浓度场等），用一系列有限 个离散点（称为节点，node)上的值的集合来代替； 通过一定的原则建立起这些离散点上变量值之间关 系的代数方程(称为离散方程，discretization equation)；求解所建立起来的代数方程以获得所求 解变量的近似值。
u v w 0 x y z
div( U ) 0 t
称为流动无散（度）条件 (Zero divergence)。
2. 动量守恒方程
对上图所示的微元体分别在三个坐标方向上应用 Newton第2定律（F=ma)在流体中的表现形式： [微元体内动量的增加率]=[作用在微元体上各种力之和] 假设流体中切应力与正应力满足Stokes假定：应 力与应变成线性关系，可得u-动量方程如下：
为流体的动力粘度 ， 称为流体的第2分子粘度。
上式右端部分可进一步转化：
v u p u u w (divU 2 ) [ ( )] [ ( )] Fx x x y x y z z x x
u u u u v w ( ) ( ) ( ) ( ) ( ) ( ) (divU ) x x y y z z x x y x z x x p Fx u u u x div( gradu ) Su grad (u ) i j k x y z
Elliptic
的函数。 椭圆型 （回流型） 抛物型 （边界层）
0,
b 4ac
2
0, 0,
Parabolic

传热学数值计算作业数值解程序：tw1=40 %三边温度tw2=100 %一边温度正弦变化幅度l1=40 %板长L1:40厘米l2=20 %板宽L2:20厘米m=41 %分划成40*20的网格n=21k=2dx=l1/(m-1)c=ones(n,m)for i=1:ma2(i)=tw1+tw2*sin(pi*dx*(i-1)/l1)c(1,i)=tw1 ,c(n,i)=a2(i)endfor j=1:nc(j,1)=tw1c(j,m)=tw1endwhile (abs(c(j,i)-k)>0.0001)k=c(j,i)for i=2:m-1for j=2:n-1c(j,i)=0.25*(c(j,i-1)+c(j,i+1)+c(j-1,i)+c(j+1,i)) endendend数值解中各网格点的温度值：数值二维温度分布图像：解析解程序： tw1=40 tw2=100 l1=40 l2=20 p=40 q=20 x(1)=0 for i=1:px(i+1)=x(i)+1 end y(1)=0 for j=1:qy(j+1)=y(j)+1 endfor i=1:p+1 for j=1:q+1n(j,i)=tw1+tw2*sinh(pi*y(j)/l1)*sin(pi*x(i)/l1)/sinh(pi*l2/l1) end end各网格点用解析式得到的温度值：50L1/cmnumerical calculation 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e解析二维温度分布图像：误差分析：取x=21，即位于板长一半处，温度随y （宽度）的变化曲线。

c1(:,1) 取自于数值解， c1(:,2) 取自于解析解 c1(:,1) c1(:,2) 40.0000 40.0000 43.3106 43.4164 46.6465 46.8538 50.0313 50.3335 53.4889 53.8771 57.0430 57.5062 60.7178 61.2434 64.5376 65.1117 68.5273 69.1350 72.7122 73.3381 77.1187 77.7470 81.7736 82.3888 86.7050 87.2922 91.9423 92.4875 97.5162 98.0068 103.4592 103.8840 109.8058 110.1555 116.5925 116.8600 123.8586 124.0388 131.6461 131.7363 140.0000 140.000050L1/cmanalytical method 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e误差曲线：由相对误差公式：d1= (c1(:,2) -c1(:,1))./ c1(:,2) 可得： d1 = 0 0.0024 0.0044 0.00600.00720.0081 0.0086 0.0088 0.0088 0.0085 0.0081 0.0075 0.0067 0.0059 0.0050 0.0041 0.0032 0.0023 0.0015 0.0007 0结论：数值解与解析解吻合很好。

1
1求解导热问题的三种基本方法：
(1) 理论分析法；(2) 数值计算法；(3) 实验法
2三种方法的基本求解过程
(1)所谓理论分析方法，就是在理论分析的基础上，直接对微分方
程在给定的定解条件下进行积分，这样获得的解称之为分析解，或叫理论解；
(2)数值计算法，把原来在时间和空间连续的物理量的场，用有限
个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，从而获得离散点上被求物理量的值；并称之为数值解；
(3) 实验法就是在传热学基本理论的指导下，采用对所研究对象的
传热过程所求量的方法
3 三种方法的特点
(1) 分析法
a 能获得所研究问题的精确解，可以为实验和数值计算提供比较依据；
b 局限性很大，对复杂的问题无法求解；
c 分析解具有普遍性，各种情况的影响清晰可见
(2) 数值法：
在很大程度上弥补了分析法的缺点，适应性强，特别对于复杂问题更显其优越性；与实验法相比成本低
(3) 实验法: 是传热学的基本研究方法，a 适应性不好；b 费用昂贵。

hx
1 ti, j
ti1, j
ti, j1
x2 i, j
2
2hx
t
f
（c）内部角点
g
2
hx
3
ti,
j
2
ti1, j ti, j1
ti1, j
ti, j1
3x2 i, j
2
2hx
tf
三 节点差分方程的求解
1） 直接解法：通过有限次运算获得精确解的方 法，如：矩阵求解，高斯消元法。 2） 迭代法：先对要计算的场作出假设（设定初 场），在迭代计算中不断予以改进，直到计算前 的假定值与计算结果相差小于允许值为止的方法， 称迭代计算收敛。
4-2 稳态导热问题的数值计算
（6） 解的分析
通过求解代数方程，获得物体中的温度分布， 根据温度场应进一步计算通过的热流量，热应力及 热变形等。因此，对于数值分析计算所得的温度场 及其它物理量应作详细分析，以获得定性或定量上 的结论。
4-2 稳态导热问题的数值计算
建立离散方程的常用方法：
(1) Taylor（泰勒）级数展开法； (2) 多项式拟合法； (3) 控制容积积分法； (4) 控制容积平衡法(也称为热平衡法)
j
y
y
x
x
i
I
除 i=1 的左边界上各节点的温度已知外，其余 (i-1)j 个 节点均需建立离散方程，共有 (i-1)j 个方程，则构成一 个封闭的代数方程组。
4-2 稳态导热问题的数值计算
1 ）线性代数方程组：代数方程一经建立，其中各 项系数在整个求解过程中不再变化； 2 ）非线性代数方程组：代数方程一经建立，其中 各项系数在整个求解过程中不断更新； 3 ）是否收敛判断：是指用迭代法求解代数方程是 否收敛，即本次迭代计算所得之解与上一次迭代计 算所得之解的偏差是否小于允许值。

对物理问题进行数值解法的基本思路可以概括 为：把原来在时间、空间坐标系中连续的物理量的 场，如导热物体的温度场等，用有限个离散点上的 值的集合来代替，通过求解按一定方法建立起来的 关于这些值的代数方程，来获得离散点上被求物理 量的值，该方法称为数值解法。 这些离散点上被求物理量值的集合称为该物理 量的数值解。
1 (tm 1,n tm 1,n tm ,n 1 tm,n 1 ) 4
一阶
4.2.2 控制容积平衡法(热平衡法)
基本思想：对每个有限大小的控制容积应用能量守
恒，从而获得温度场的代数方程组，它从基本物理
现象和基本定律出发，不必事先建立控制方程，依 据能量守恒和Fourier导热定律即可。 能量守恒： 流入控制体的总热流量＋控制体内热源生成热 ＝ 流出控制体的总热流量＋控制体内能的增量
x y
tm , n
y
2x 2 qw x 1 3x 2 (2tm 1,n 2tm ,n 1 tm ,n 1 tm 1,n ) 6 2
讨论关于边界热流密度的三种情况： （1）绝热边界
即令上式 qw 0 即可。
（2）qw 值不为零
流入元体，qw 取正，流出元体，qw 取负 （3）对流边界 此时 qw h(t f tm,n ) ，将此表达式代入上述方程，并 将此项中的
tm , n 1 (tm 1,n tm 1,n tm ,n 1 tm ,n 1 ) 4
（4） 设立迭代初场
代数方程组的求解方法有直接解法与迭
代解法，传热问题的有限差分法中主要采用
迭代法。采用迭代法求解时，需对被求的温
度场预先设定一个解，这个解称为初场，并
在求解过程中不断改进。
4.1.2 物理问题的数值求解过程

传热学作业数值计算1 程序内容： 数值计算matlab程序内容：>> tw1=10; % 赋初值赋初值 tw2=20; c=1.5; p2=20; p1=c*p2; L2=40; L1=c*L2; deltaX=L2/p2; a=p2+1; b=p1+1; =ones(a,b)*5; m1=ones(a,b); m1(a,2:b-1)=zeros(1,b-2); m1(2:a,1)=zeros(a-1,1); m1(2:a,b)=zeros(a-1,1); m1(1,:)=ones(1,b)*2; k=0; max1=1.0; tn= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw1+tw2*sin(pi*(j-1)/(b-1)); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 t2=ones(a,b); %求解析温度场求解析温度场for i=a:-1:1 for j=1:1:b y=deltaX*(a-i); x=deltaX*(j-1); t2(i,j)=tw1+tw2*sin(pi*x/L1)*(sinh(pi*y/L1))/(sinh(pi*L2/L1)); end end t2 迭代次数k =706 数值解温度场 数值解每次迭代的最大误差max1 =9.8531e-07 解析温度场 t2 解析温度场取第11行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差数值计算matlab 程序内容：程序内容： >> tw1=10; tw2=20; c=1.5; p2=20; p1=c*p2; L2=20; deltaX=L2/p2; L1=c*L2; a=p2+1; b=p1+1; =ones(a,b)*5; m1=ones(a,b); m1(a,2:b-1)=zeros(1,b-2); m1(2:a,1)=zeros(a-1,1); m1(2:a,b)=zeros(a-1,1); m1(1,:)=ones(1,b)*2; k=0; max1=1.0; tn= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw2; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 tx=ones(a,b); for i=1:1:a for j=1:1:b y=(a-i)*deltaX; x=(j-1)*deltaX; m=sym('m'); g=(((-1)^(m+1)+1)/m)*sin(m*pi*x/L1)*sinh(m*pi*y/L1)/sinh(m*pi*L2/L1); h=symsum(g,m,1,100); tx(i,j)=2*h*(tw2-tw1)/pi+tw1; end end tx 迭代次数k = 695 数值解温度场 数值解每次迭代的最大误差max1 =9.8243e-07 解析温度场 tx = 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像： 数值温度场 图像图像 ：解析温度场tx图像：数值解与解析解的误差数值解与解析解的误差程序内容： 数值计算matlab程序内容：>> t0=90; =10; L=10; c=0.25; p2=20; p1=p2/c; B=c*L; d=0.5*B; h=10; a=p2+1; b=p1+1; deltaX=B/p2; lambda=160; Bi=h*deltaX/lambda; =ones(a,b)*10; m1=ones(a,b)*3; m1(2:a-1,1)=zeros(a-2,1); m1(a,2:b-1)=ones(1,b-2); m1(1,2:b-1)=ones(1,b-2)*6; m1(2:a-1,b)=ones(a-2,1)*2; m1(1,b)=ones(1,1)*4; m1(a,b)=ones(1,1)*5; m1(1,1)=7; m1(a,1)=8; tn= ; max1=1.0; k=0; while ( max1>1e-6) k=k+1; max1=0; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=t0; case 1 tn(i,j)=(2*tn(i-1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 2 tn(i,j)=(2*tn(i,j-1)+tn(i-1,j)+tn(i+1,j)-4* )/(4+2*Bi)+ ; case 3 tn(i,j)=0.25*(tn(i,j-1)+tn(i,j+1)+tn(i-1,j)+tn(i+1,j)); case 4 tn(i,j)=(tn(i,j-1)+tn(i+1,j)-2* )/(2*Bi+2)+ ; case 5 tn(i,j)=(tn(i,j-1)+tn(i-1,j)-2* )/(2*Bi+2)+ ; case 6 tn(i,j)=(2*tn(i+1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 7 tn(i,j)=t0; case 8 tn(i,j)=t0; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k ta=ones(a,b); Bi1=h*d/lambda; sbi=sqrt(Bi1); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end x=deltaX*(j-1); ta(i,j)=(cosh(sbi*(L-x)/d)+sbi*sinh(sbi*(L-x)/d))*(t0- )/(cosh(sbi*L/d)+sbi*sinh(sbi*L/d))+ ; end end ta 迭代次数k =1461 数值解温度场 解析温度场 ta 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点曲线为第11行的解析解的直线，散点为其数值解的点行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像如下图像如下数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差程序内容：数值计算matlab程序内容：>> tw=10; L2=15; c=0.75; L1=L2/c; p2=24 ; p1=p2/c; deltaX=2*L2/p2; a=p2+1; b=p1+1; lambda=16; qv0=24; =ones(a,b)*5; m1=ones(a,b); m1(1,:)=zeros(1,b); m1(2:a,b)=zeros(a-1,1); m1(2:a,1)=zeros(a-1,1); m1(a,2:b-1)=zeros(1,b-2); tn= ; max1=1.0; k=0; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=tw; case 1 tn(i,j)=0.25*(tn(i-1,j)+tn(i+1,j)+tn(i,j-1)+tn(i,j+1)+qv0*(deltaX^2)/lambda); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k; tx=ones(a,b); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end if j>(b+1)/2 x=(j-(b+1)/2)*deltaX; else x=-((b+1)/2-j)*deltaX; end m=sym('m'); xi=(2*m-1)*pi/2; g=((-1)^m)/(xi^3)*(cosh(xi*y/L1)/cosh(xi*L2/L1))*cos(xi*x/L1); h=symsum(g,m,1,100); tx(i,j)=2*qv0*L1^2/lambda*h+qv0*(L1^2-x^2)/(2*lambda)+tw; end end tx 数值温度场 解析温度场tx 取第13行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点曲线为第13行的解析解的直线，散点为其数值解的点解析解 第13行的误差=[数值解(13行) –解析解(13行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差。