No.08 1026 4 导热问题的数值解法

1、第四章第四章热传导问题的数值解法热传导问题的数值解法2导热问题一般为导热问题一般为:边界条件初始条件导热微分方程).,( 0)(zyxftttc1. 分析解分析解(analytical method): 通过对上述方程积分求得通过对上述方程积分求得(有限有限情况情况)，也称，也称理论解理论解。2. 数值解数值解(numerical method): 用某种方式把微分方程化为关用某种方式把微分方程化为关于各个离散点于各个离散点(节点节点)的代数方程的代数方程,通过解代数方程获得问题近通过解代数方程获得问题近似解的方法。似解的方法。 3 把原来在时间、空间坐标系中把原来在时间、空间坐标系中连续的连

2、续的物理量的场，如导物理量的场，如导热物体的温度场等，热物体的温度场等，用有限个离散点上的值的集合来代替，用有限个离散点上的值的集合来代替，通过求解通过求解按一定方法按一定方法建立起来的关于这些值的代数方程，来建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。获得离散点上被求物理量的值。 连续连续 离散离散(任意情况任意情况)数值解法的基本思想：数值解法的基本思想： (m,n)(m-1,n)(m+1,n)tm,ntm-1,ntm+1,n41. 数学描述数学描述2. 区域离散化区域离散化3. 建立节点物理量的代数方程建立节点物理量的代数方程4. 设立迭代初场（初始值设定）设立迭代初场

3、（初始值设定）5. 求解代数方程组求解代数方程组6. 解的分析解的分析数值求解的基本步骤：数值求解的基本步骤：5物理问题数值求解过程的框图：物理问题数值求解过程的框图：建立控制方程及定解条件建立控制方程及定解条件确定节点（区域离散化）确定节点（区域离散化）建立节点物理量的离散方程建立节点物理量的离散方程设立温度场的迭代初值设立温度场的迭代初值求解离散方程求解离散方程是否收敛是否收敛解的分析解的分析改进初场改进初场是是否否61. 控制方程的建立及定解条件控制方程的建立及定解条件无限长棱柱（如图）导热、沿高度各截面的温度分布相同，无限长棱柱（如图）导热、沿高度各截面的温度分布相同，可简化为二维问题

4、。可简化为二维问题。)(constx ht qbtt 0aby7无限长棱柱（如图）导热、沿高度各截面的无限长棱柱（如图）导热、沿高度各截面的温度分布相同，可简化为二维问题。温度分布相同，可简化为二维问题。)(constbttbytthytyqxtaxxtxytyt )( 0 0 002222x ht qbtt 0aby8实现过程：实现过程：以二维矩形区域为例，所谓区域离散化，就是将研究区域分解成有限数量的小区域（单元），单元的顶点（或中心点）称作每个节点都有自己的控制区域，称作，控制体内所有特性都是均匀的，节点的温度代表每个控制体的温度。节点之间的距离称为节点之间的连线称为控制容积的分界面称为

5、节点的表示方法：节点的表示方法： （m，n）横坐标节点编号纵坐标节点编号xyxynm(m,n)MN2. 区域离散化区域离散化9网格网格(grid )划分划分节点节点(node): 网格线交点网格线交点.控制容积控制容积(control volume): 节点代表的区域节点代表的区域 ,其边界位于两点之间其边界位于两点之间.界面界面(interface): 控制容积的边界控制容积的边界.均分网格均分网格:constconst yx节点编号节点编号: 从小往大排从小往大排网格划分方法网格划分方法:方法方法1： 先确定节点，后定界面先确定节点，后定界面 方法方法2： 先确定界面，后定节点先确定界面，

6、后定节点xyxynm(m,n)MN10(1) Taylor级数展开法；级数展开法；(2) 控制容积平衡法控制容积平衡法(也称为热平衡法也称为热平衡法)(3) 多项式拟合法；(4) 控制容积积分法；3. 离散方程的建立离散方程的建立11(1) 泰勒级数展开法泰勒级数展开法根据泰勒级数展开式，用节点(m,n)的温度tm,n来表示节点(m+1,n)的温度tm+1,n! 3! 23,332,22, 1xxtxxtxxtttnmnmnmnmnm! 3! 23,332,22, 1xxtxxtxxtttnmnmnmnmnm用节点(m,n)的温度tm,n来表示节点(m-1,n)的温度tm-1,n(m,n)(m

7、,n+1)(m+1,n)(m,n-1)(m-1,n)12若取上面式右边的前三项，并将两式相加移项整理即得二阶导数的中心差分：)(222, 1, 1,22xoxtttxtnmnmnmnm截断误差)(2221,1,22yoytttytnmnmnmnm同理可得：(m,n)(m,n+1)(m+1,n)(m,n-1)(m-1,n)13对于二维稳态导热问题，在直角坐标中，其导热微分方程为：22220gttxy)(222, 1, 1,22xoxtttxtnmnmnmnm)(2221,1,22yoytttytnmnmnmnm0ytt 2txtt 2tj , i , v21j , ij , i1j , i2j

8、, 1ij , ij , 1i其节点方程为：(m,n)(m,n+1)(m+1,n)(m,n-1)(m-1,n)14(2)热平衡法 (控制容积平衡法)基本思想基本思想：对每个有限大小的控制容积应用能量守恒，从而获得温度场的代数方程组。它从基本物理现象和基本定律出发，不必事先建立控制方程，依据能量守恒和Fourier导热定律即可。能量守恒：流入控制体的总热流量控制体内热源生成热流入控制体的总热流量控制体内热源生成热= 控制体内能的增量控制体内能的增量即： 单位单位：igst W注意： 上面的公式对内部节点和边界节点均适用 所谓流入控制体的总热流量，是指从所有方向流入的151,1,1,10imnmn

9、m nm n 稳态、无内热源时：从所有方向流入控制体的总热流量i0内部节点：0右左下上(m, n)oyx(m-1,n)(m+1,n)(m,n-1) x x y y (m,n+1) (m,n+1) (m,n-1) (m+1,n) (m-1,n)igst 16以二维、稳态、有内热源的导热问题为例，此时：以二维、稳态、有内热源的导热问题为例，此时：0ig 下左右上(m, n)oyx(m-1,n)(m+1,n)(m,n-1) x x y y(m,n+1)xtyxtAdddd左可见：当温度场还没有求出来之前，我们并不知道所以，必须假设相邻节点间的温度分布形式，这里，我们假定温度呈分段线性分布，如下图所示

10、xt dd17(m,n)(m-1,n)(m+1,n)tm,ntm-1,ntm+1,nxttyxxttyxtynmnmnmnmnmnm, 1, 1, 1dd左可见，节点越多，假设的分段线性分布越接近真实的温度分布。此时：xttynmnm, 1右yttxnmnm,1,上yttxnmnm,1,下内热源：g Vx y m,n+1) (m,n-1) (m+1,n) (m-1,n)180g下左右上0, 1, 1,1,1,yxxttyxttyyttxyttxnmnmnmnmnmnmnmnmyx042,1,1, 1, 1xtttttnmnmnmnmnmxtttttnmnmnmnmnm21,1, 1, 1,4(

11、m,n)(m,n+1)(m+1,n)(m,n-1)(m-1,n)19xtttttnmnmnmnmnm21,1, 1, 1,4无内热源时变为：无内热源时变为：1,1, 1, 1,4nmnmnmnmnmttttt待求节点的温度前的系数一定等于其它所有相邻节点温度前的系数之和。这一结论也适用于边界节点。但这里不包括热流、或热流密但这里不包括热流、或热流密度和内热源之前的系数。度和内热源之前的系数。(m,n)(m,n+1)(m+1,n)(m,n-1)(m-1,n)204-2 一维无限大平板、稳态、常物性、无内热源、左侧第一类边界，右侧第三类，如右图所示，将其均匀分成三个控制体，试建立离散方程13212

12、23433441:2:03:04:()0wttttttxxttttxxtth ttx边界节点1234twth内部节点内部节点边界节点21形成如下代数方程组：11 112 21121 122 2221 12 2.n nn nnnnn nna ta ta tba ta ta tba ta ta tb代数方程组的通用形式为：ATb23234342:23:204:()wttttttth thtxx23421012100wtttxxhhtt 22写出所有内节点和边界节点的温度差分方程n个未知节点温度，n个代数方程式：代数方程组的求解方法：直接解法、迭代解法直接解法、迭代解法11 112 21121 12

13、2 2221 12 2.n nn nnnnn nna ta ta tba ta ta tba ta ta tbATb4-2 23通过有限次运算获得代数方程精确解; 矩阵求逆、高斯消元法迭代解法：迭代解法：先对要计算的场作出假设、在迭代计算过程中不断予以改进、直到计算结果与假定值的结果相差小于允许值。称迭代计算已经收敛。缺点：缺点：所需内存较大、方程数目多时不便、不适用于非线性问题迭代解法有多种：高斯-赛德尔迭代的特点：24在计算后面的节点温度时应按下式（采用最新值）在计算后面的节点温度时应按下式（采用最新值）例如：根据第例如：根据第 k 次迭代的数值次迭代的数值(k)n(k)2(k)1.ttt

14、、可以求得节点温度：可以求得节点温度：(1)( )( )( )112 211.kkkkn nta ta tb(1)(1)( )( )221 122(1)(1)(1)( )( )331 132 233(1)(1)(1)(1 12 211.kkkkn nkkkkkn nkkkknnnnnnta ta tbta ta ta tbta ta tat1)( )knb25判断迭代是否收敛的准则：判断迭代是否收敛的准则：)(max)() 1()()() 1()() 1(maxmaxmaxkkikikikikikikitttttttt631010 取值一般相对偏差允许的偏差；k 及及 k+1 表示迭代次数；表

15、示迭代次数；第第k 次迭代得到的最大值次迭代得到的最大值(k)maxt当有接近于零的当有接近于零的t 时，第三个较好时，第三个较好26获得收敛解的条件：获得收敛解的条件：主对角占优主对角占优111333231222321111312aaaaaaaaa热平衡法建立的离散方程自动满足上述条件。热平衡法建立的离散方程自动满足上述条件。27主要内容：主要内容：4-1 导热问题数值求解的基本思想及内节点离散方程的建立4-2 边界节点离散方程的建立及代数方程的求解4-3 非稳态导热问题的数值解法非稳态导热问题的数值解法28 4-3 非稳态非稳态导热问题的数值解法导热问题的数值解法 增加了非稳态项增加了非稳

16、态项 扩散项的处理方法与前一样扩散项的处理方法与前一样数学描述数学描述区域离散化区域离散化建立节点物理量的代数方程建立节点物理量的代数方程设立迭代初场设立迭代初场求解代数方程组求解代数方程组解的分析解的分析29空间坐标空间坐标 x: 1 N x 空间步长空间步长时间坐标时间坐标 : 1 I 时间步长时间步长（n,i）代表了时间空间区）代表了时间空间区域中的一个接点位置域中的一个接点位置 t(i)n一维非稳态网格划分一维非稳态网格划分30将温度函数将温度函数 t 在节点（在节点（n,i+1）和（）和（n,i-1）对点（）对点（n,i）作泰勒级）作泰勒级数展开数展开inininintttt,222

17、,)()1(2inininintttt,222,)()1(2从第一式得出从第一式得出)()1()()1(,)(inininininttOttt向前差分向前差分 forward difference31inininintttt,222,)()1(2inininintttt,222,)()1(2从第二式得出从第二式得出)1()()1()(,)(inininininttOttt向后差分向后差分 back difference32inininintttt,222,)()1(2inininintttt,222,)()1(2二级数相减得到二级数相减得到:2)(2)1()1(2)1()1(,inininin

18、inttOttt中心差分中心差分 central finite difference33常物性一维非稳态问题常物性一维非稳态问题: 非稳态项非稳态项-向前差分向前差分 扩散项扩散项-中心差分中心差分22xtat2)(1)()(1)()1(2xtttattininininin显示格式显示格式 explicit finite difference scheme可直接求解可直接求解.)(2)(1)(12)1(21inininintxattxat34如果扩散项也用如果扩散项也用(i+1)时层上的值来表示时层上的值来表示,则则2)1(1)1()1(1)()1(2xtttattininininin不能写成

19、与显示格式类似的形式，不能写成与显示格式类似的形式， 称为称为隐示格式隐示格式 implicit finite difference scheme不能直接求解，不能直接求解，必需联立求解必需联立求解.与显示格式比与显示格式比,计算工作量大计算工作量大,但它是时间步长没有限制但它是时间步长没有限制. 35边界节点的处理边界节点的处理)()1()()()(12)(iNiNiNfiNiNttxctthxttfiNiNiNtxchtxaxaxchtt22221)(122)()1(式中式中: 网格傅里叶数网格傅里叶数BiFoxhxcxch2fiNiNiNtBiFotFoFoBiFott22221)(1)

20、()1()()(1)(1)1(21inininintFottFot一维无限大平板一维无限大平板边界边界:内节点内节点:2xaFoV计算机实习：计算机实习：1、补充材料中计算机实习指导书、补充材料中计算机实习指导书: 练习题一、二、三练习题一、二、三班级名称班级名称上机上机人数人数机机 房房周周 次次日日 次次单单 元元传热学（能动学院） 42 405机房机房1335传热学（能动学院） 46 406机房机房1335传热学（能动学院） 42 405机房机房1435传热学（能动学院） 46 406机房机房1435传热学（能动学院） 42 405机房机房1535传热学（能动学院） 46 406机房机房15352、作业上交：、作业上交： 请于第14周发到： 邮件名称以姓名学号姓名学号命名，附件采用压缩文件；37、等截面直肋如右图所示,碳钢=43.2W/(m.K), 写出1，2，3，4节点的离散方程。 1030KmW 120C 25 C 2002000httx10 x123438、在图中所示的有内热源的二维导热区域中，一个界面绝热，一个界面等温温度为t0（包括节点4），其余两个界面与温度