华科传热学2009

1、2022-3-101主讲：刘志春主讲：刘志春能源与动力工程学院能源与动力工程学院华中科技大学华中科技大学2022-3-1029-1 数值计算的基本思想数值计算的基本思想*9-2 流动与传热的数值计算流动与传热的数值计算9-3 Saints2D软件简介软件简介2022-3-103首先，我们以导热问题为例，介绍计算区域离首先，我们以导热问题为例，介绍计算区域离散化的概念、内节点与边界节点方程式的建立散化的概念、内节点与边界节点方程式的建立方法、节点方程组的求解过程，以及非稳态导方法、节点方程组的求解过程，以及非稳态导热问题的显示与隐示差分格式。热问题的显示与隐示差分格式。然后，介绍在上述思想的基础

2、上开发的流动与然后，介绍在上述思想的基础上开发的流动与传热计算软件传热计算软件Saints2DSaints2D，并给出传热问题虚拟，并给出传热问题虚拟实验的计算示例。实验的计算示例。 2022-3-104数值求解数值求解通常是对微分方程直接进行数值积分通常是对微分方程直接进行数值积分或者把微分方程转化为一组代数方程组再进行或者把微分方程转化为一组代数方程组再进行求解。这里要介绍的是后一种方法。求解。这里要介绍的是后一种方法。 如何实现从微分方程到代数方程的转化又可以如何实现从微分方程到代数方程的转化又可以采用不同的数学方法，如采用不同的数学方法，如有限差分法有限差分法、有限元有限元法法和和边界

3、元法边界元法等。这里仅向读者简要地介绍用等。这里仅向读者简要地介绍用有限差分方法从微分方程确立代数方程的处理有限差分方法从微分方程确立代数方程的处理过程。过程。 2022-3-105有限差分法的基本思想有限差分法的基本思想是把原来在时间和空间是把原来在时间和空间坐标中连续变化的物理量（如温度、压力、速坐标中连续变化的物理量（如温度、压力、速度和热流等），用有限数目的离散点上的数值度和热流等），用有限数目的离散点上的数值集合来近似表达。集合来近似表达。有限差分的数学基础是用有限差分的数学基础是用差商差商代替代替微商微商（导（导数）。数）。几何意义是用函数在某区域内的平均变化率代几何意义是用函数在

4、某区域内的平均变化率代替函数的真实变化率。替函数的真实变化率。 2022-3-106在图中可以看出有限差分表示的温度场与真实在图中可以看出有限差分表示的温度场与真实温度场的区别。图中用温度场的区别。图中用T0、T1、T2表示连续表示连续的温度场的温度场T；x为为步长步长，它将区域的，它将区域的x方向划分方向划分为有限个数的区域，为有限个数的区域，x0、x1、x2，它们，它们可以相等，也可以不相等。可以相等，也可以不相等。 当当x相等时，相等时，T1处的真实处的真实变化率变化率a可以用平均变化可以用平均变化率率b、c或或d来表示，其中来表示，其中b、c和和d分别表示三种不同差分别表示三种不同差分

5、格式下的温度随时间的分格式下的温度随时间的变化率变化率 x1x2x3x0 xTT3T2T1T0 b d c a2022-3-107b为为向后差分格式向后差分格式xxxTxTdxdT)()(111c为为向前差分格式向前差分格式 111()()dTT xxT xdxxd为为中心差分格式中心差分格式 x1x2x3x0 xTT3T2T1T0 b d c axxxTxxTdxdT2)()(1112022-3-108这种差分格式可推广到高阶微商（导数）。这种差分格式可推广到高阶微商（导数）。对于对于二阶导数的差分格式二阶导数的差分格式可以在一阶差分格式可以在一阶差分格式的基础上得出：的基础上得出： 211

6、1212)()()(2)(xxxTxTxxTdxTd这样处理后，反映温度场随时间、空间连续变这样处理后，反映温度场随时间、空间连续变化的微分方程就可以用反映离散点间温度线性化的微分方程就可以用反映离散点间温度线性变化规律的代数方程来表示。变化规律的代数方程来表示。当利用相应的数学办法求解这些代数方程组之当利用相应的数学办法求解这些代数方程组之后，我们就能获得离散点上的温度值。这些温后，我们就能获得离散点上的温度值。这些温度值就可以近似表示温度场的连续的温度分布。度值就可以近似表示温度场的连续的温度分布。 2022-3-109从上面的分析不难看出，当我们要对流动与传从上面的分析不难看出，当我们要

7、对流动与传热问题进行数值求解时一定要采取热问题进行数值求解时一定要采取三个大的步三个大的步骤骤，即：，即：a) 研究区域的离散化；研究区域的离散化；b) 散点（节点）差分方程的建立；散点（节点）差分方程的建立；c) 节点方程（代数方程）的求解。节点方程（代数方程）的求解。2022-3-10101 时间与空间的离散化时间与空间的离散化 进行数值求解时进行数值求解时,首先是首先是在所研究的时间和空间在所研究的时间和空间区域内把时间和空间分区域内把时间和空间分割成为有限大小的小区割成为有限大小的小区域。域。图表示了长柱体矩形截图表示了长柱体矩形截面上区域离散化的情况。面上区域离散化的情况。K-1K-

8、1时刻时刻SWENPSWENPSWENPxyK K时刻时刻K+1K+1时刻时刻xy2022-3-1011对于给定的空间区域，对于给定的空间区域，在在x方向上的步长为方向上的步长为x，在在y方向上的步长为方向上的步长为y，用它们作为空间尺度可用它们作为空间尺度可以将矩形区域划分成纵以将矩形区域划分成纵横交错的网格系统。横交错的网格系统。计算区域就被这些网格线分隔成一系列的小的计算区域就被这些网格线分隔成一系列的小的区域，称为区域，称为控制面积控制面积，对于三维情况则为，对于三维情况则为控制控制体积体积或或控制容积控制容积，因而在一般意义上称之为，因而在一般意义上称之为控控制体制体；控制体的中心点

9、称为；控制体的中心点称为节点节点。K-1K-1时刻时刻SWENPSWENPSWENPxyK K时刻时刻K+1K+1时刻时刻xy2022-3-1012控制体的形状是随着坐标系的不同而改变的，控制体的形状是随着坐标系的不同而改变的，这里的控制体是一个个的矩形面积。网格的步这里的控制体是一个个的矩形面积。网格的步长在每一个方向上可以均匀划分，也可以不均长在每一个方向上可以均匀划分，也可以不均匀的划分；所得到网格，相应地被称为匀的划分；所得到网格，相应地被称为均匀网均匀网格格或者或者非均匀网格非均匀网格。选用不同的步长和不同的划分方法，可以将同选用不同的步长和不同的划分方法，可以将同一区域划分出不同大

10、小、不同数目的控制区域，一区域划分出不同大小、不同数目的控制区域，以及不同数目的以及不同数目的节点数节点数。 2022-3-1013获得每个节点上的温度值，就是导热数值计算获得每个节点上的温度值，就是导热数值计算的目的。，随着步长的不断减小，节点数目的的目的。，随着步长的不断减小，节点数目的不断增加，由节点温度表示的离散温度场就更不断增加，由节点温度表示的离散温度场就更接近连续温度场，但计算工作量也会增加。接近连续温度场，但计算工作量也会增加。 在时间方向上离散化的步长常用在时间方向上离散化的步长常用来表示，来表示，的选取也是可大可小的，也可以随时间的进程的选取也是可大可小的，也可以随时间的进

11、程而变化。显然，无限小的时间步长而变化。显然，无限小的时间步长亦会使得亦会使得离散温度变化接近连续的温度改变，但随之而离散温度变化接近连续的温度改变，但随之而来的是相应的计算工作量将会增加。来的是相应的计算工作量将会增加。2022-3-10142 节点方程的建立节点方程的建立建立节点差分方程可采用不同的方法，主要分建立节点差分方程可采用不同的方法，主要分为两大类：为两大类： 第一类包括第一类包括泰勒泰勒(Taylor)级数展开法级数展开法和和多项式多项式拟合法拟合法，它偏重于从数学的角度进行推导，其，它偏重于从数学的角度进行推导，其优点是便于对离散方程进行数学特性分析，但优点是便于对离散方程进

12、行数学特性分析，但缺点是变步长网格的离散方程形式复杂、导出缺点是变步长网格的离散方程形式复杂、导出过程的物理概念不清晰、不能保证差分方程具过程的物理概念不清晰、不能保证差分方程具有守恒特性。有守恒特性。 2022-3-1015第二类包括第二类包括控制体热平衡法控制体热平衡法和和控制容积积分法控制容积积分法，其优点是推导过程的物理概念清晰、离散方程其优点是推导过程的物理概念清晰、离散方程系数具有一定物理意义、保证差分方程具有守系数具有一定物理意义、保证差分方程具有守恒特性，但缺点是不便于对离散方程进行数学恒特性，但缺点是不便于对离散方程进行数学特性分析。特性分析。 下面我们采用控制体热平衡法来建

13、立节点方下面我们采用控制体热平衡法来建立节点方程。程。 2022-3-1016 内节点方程内节点方程控制体热平衡法控制体热平衡法建立节点方程的建立节点方程的过程过程是将能量守是将能量守恒方程应用于控恒方程应用于控制体，建立该节制体，建立该节点与周围节点之点与周围节点之间的能量平衡关间的能量平衡关系式。系式。(i+1,j)E(i, j)oyx(i-1,j)W(i,j-1)S (i,j+1)NP x x y y再利用傅里叶导热定律，最后获得控制体节点再利用傅里叶导热定律，最后获得控制体节点温度与周围节点温度之间的关系式。温度与周围节点温度之间的关系式。 2022-3-1017考察图中的节点考察图中

14、的节点P及其控制体，由能量平衡关及其控制体，由能量平衡关系应有系应有 EVNSEW式中，式中，W、E、S和和N分别为邻分别为邻近节点近节点W、E、S和和N通过传导方式传通过传导方式传给节点给节点P的热流量；的热流量；V为单位时间控制为单位时间控制体内热源发热量；体内热源发热量；(i+1,j)E(i, j)oyx(i-1,j)W(i,j-1)S (i,j+1)NP x x y yE为控制体单位时间内热能的增加量。为控制体单位时间内热能的增加量。2022-3-1018由导热傅里叶定律，由导热傅里叶定律，在线性温度分布的在线性温度分布的假设下，时刻假设下，时刻K周周围节点传给节点围节点传给节点P的热

15、流量分别为：的热流量分别为： 1 .)(1 .)(1 .)(1 .)(xTTyxTTyyTTxyTTxKPKNNKPKSSKPKEEKPKWW(i+1,j)E(i, j)oyx(i-1,j)W(i,j-1)S (i,j+1)NP x x y y2022-3-1019控制体的发热流量控制体的发热流量其中其中qV为为内热源强度内热源强度，即单位时间单位体积的，即单位时间单位体积的内热源发热量。内热源发热量。 1yxqVV控制体单位时间的内能增加量为控制体单位时间的内能增加量为 1111yxTTcEyxTTcEKPKPKPKP或前者为时间上的向前差分，而后者为时间上向前者为时间上的向前差分，而后者为

16、时间上向后差分。以上关系式中温度后差分。以上关系式中温度T的上标为所在时的上标为所在时刻，下标为所在空间位置。刻，下标为所在空间位置。 2022-3-1020假设假设x=y，经整理可以得出二维非稳态导热，经整理可以得出二维非稳态导热问题的内节点的两种差分格式的差分方程，即问题的内节点的两种差分格式的差分方程，即 a) 显式差分格式显式差分格式 cqTxaTTTTxaTVKPKNKSKEKWKP)41 ()(221定义定义网格傅里叶数网格傅里叶数 2aFox其物理意义是表征控制体的导热性能与热储蓄其物理意义是表征控制体的导热性能与热储蓄性能之间的对比关系，反映控制体温度随时间性能之间的对比关系，

17、反映控制体温度随时间变化的动态特性。显式差分格式简化为变化的动态特性。显式差分格式简化为 2022-3-1021b) 隐式差分格式隐式差分格式 cqTFoTTTTFoTVKPKNKSKEKWKP)41 ()(1cqTTTTTFoFoTVKPKNKSKEKWKP1)(411或改写为或改写为 111111()1 4KKKKKKVPWESNPqTFo TTTTTcFo2022-3-1022显示差分格式最突出的优点是节点温度表达式显示差分格式最突出的优点是节点温度表达式的右边只涉及的右边只涉及K K时刻时刻( (前一时刻前一时刻) )的节点温度值，的节点温度值，只要知道前一时刻周围节点的温度值就可以求

18、只要知道前一时刻周围节点的温度值就可以求出该节点的出该节点的K K +1+1时刻时刻( (当前时刻当前时刻) )的温度值；的温度值；cqTFoTTTTFoTVKPKNKSKEKWKP)41 ()(12022-3-1023隐示差分格式温度表达式的右端除了隐示差分格式温度表达式的右端除了K K时刻时刻( (前前一时刻一时刻) )的节点温度值以外，还含有的节点温度值以外，还含有K K +1+1时刻时刻( (当前时刻当前时刻) )的温度值，这就意味着必须同时计的温度值，这就意味着必须同时计算当前时刻所有节点的温度值，即必须联立求算当前时刻所有节点的温度值，即必须联立求解解K K +1+1时刻所有节点的

19、差分方程组，计算工作时刻所有节点的差分方程组，计算工作量增大也就是显而易见的了。量增大也就是显而易见的了。 111111()1 4KKKKKKVPWESNPqTFo TTTTTcFo2022-3-1024虽然显示差分格式计算比较方便，但它存在一虽然显示差分格式计算比较方便，但它存在一个缺点，即计算式中个缺点，即计算式中F FO O值必须满足一定的条件值必须满足一定的条件才不至于引起数值计算出现不收敛的问题，这才不至于引起数值计算出现不收敛的问题，这在数值计算中称为在数值计算中称为差分格式的不稳定性差分格式的不稳定性。这里差分方程稳定性的条件是式这里差分方程稳定性的条件是式(9-2)(9-2)中

20、的变量中的变量T T前面的系数必须大于或等于零，分析一下差分前面的系数必须大于或等于零，分析一下差分方程中的各项系数，有方程中的各项系数，有cqTFoTTTTFoTVKPKNKSKEKWKP)41 ()(141Fo2022-3-102541Fo此式称为此式称为显示差分格式的稳定性判据显示差分格式的稳定性判据，从中看，从中看出时间步长和空间步长是相互制约的。为了获出时间步长和空间步长是相互制约的。为了获得较为精确的节点温度值，空间步长得较为精确的节点温度值，空间步长xx的选的选择不能太小，按照稳定性判据的要求势必会使择不能太小，按照稳定性判据的要求势必会使时间步长时间步长也要相应地不能太大，因而

21、必须也要相应地不能太大，因而必须在增加节点数目的同时增多时间间隔，从而使在增加节点数目的同时增多时间间隔，从而使计算工作量加大。计算工作量加大。 隐示差分格式是无条件稳定的隐示差分格式是无条件稳定的。2022-3-1026这里指出，以上的讨论及结果适用于对非稳态这里指出，以上的讨论及结果适用于对非稳态导热问题。对于稳态导热问题，其实应更简单导热问题。对于稳态导热问题，其实应更简单, ,只需要在式只需要在式(9-2)(9-2)或者式或者式(9-3)(9-3)中中 Fo令令或或均可得到二维稳态导热问题的内节点方程式均可得到二维稳态导热问题的内节点方程式 4)(412xqTTTTTVNSEWP202

22、2-3-1027 边界节点方程边界节点方程以对流换热边界为例，从流体以对流换热边界为例，从流体侧来看，应用牛顿冷却公式，侧来看，应用牛顿冷却公式，流体与壁面之间的对流换热量流体与壁面之间的对流换热量()1Eeh TTy 再者，从控制体侧来看，应用再者，从控制体侧来看，应用傅里叶定律，假定壁面处的温傅里叶定律，假定壁面处的温度梯度取向后差分格式，则应度梯度取向后差分格式，则应有有 1/2KePETTyx SNWPTyyxxh h2022-3-1028只要网格步长只要网格步长x足够小，两者的结果应该是一足够小，两者的结果应该是一致的。从而可消去未知量致的。从而可消去未知量Te，得到，得到)/(21

23、2)(xhTTKPE代入式代入式(9-1)，得到对流换热边界节点的两种差，得到对流换热边界节点的两种差分格式，即分格式，即 a) 显式差分格式显式差分格式 122()(13)22KKKKKVPWSNPBiBiqTFoTTTFoFoTFo TcBiBi显式差分格式的稳定性判据为显式差分格式的稳定性判据为 )2/(231BiBiFo2022-3-1029b) 隐式差分格式隐式差分格式 112()1(56)/(2)2KKKKKVPWSNPBiqTFo TTTTFo TcFoBiBiBi隐式差分格式仍然是无条件稳定的隐式差分格式仍然是无条件稳定的 定义定义网格毕渥数网格毕渥数 xhBi其物理意义是体现

24、控制体和环境间的换热性能其物理意义是体现控制体和环境间的换热性能与其导热性能之间的对比关系。与其导热性能之间的对比关系。 2022-3-1030这里指出，上面的结果是在对流换热边界的情这里指出，上面的结果是在对流换热边界的情况下得到的，但经过简单处理，可直接用于绝况下得到的，但经过简单处理，可直接用于绝热边界条件与恒壁温边界条件。热边界条件与恒壁温边界条件。 BiconstTTe令令，则有，则有即可简化为恒壁温边界条件下对应的差分格式。即可简化为恒壁温边界条件下对应的差分格式。 0Bi0h而令而令，即，即即可简化为绝热边界条件下对应的差分格式。即可简化为绝热边界条件下对应的差分格式。 2022

25、-3-1031由上面的讨论可以看出，对应于离散温度场的由上面的讨论可以看出，对应于离散温度场的每一个节点均可以列出相应的差分方程，这样每一个节点均可以列出相应的差分方程，这样就可以得出就可以得出与节点数目相同的一组代数方程组与节点数目相同的一组代数方程组。当联立求解这个代数方程组时，最后就可以得当联立求解这个代数方程组时，最后就可以得出每一个节点的温度值。出每一个节点的温度值。 2022-3-10323 节点方程的求解节点方程的求解 一般情况下，差分方程组是线性代数方程组，一般情况下，差分方程组是线性代数方程组，而线性代数方程组是可以用而线性代数方程组是可以用直接法直接法和和迭代法迭代法求求解

26、的。常用的直接法有解的。常用的直接法有高斯消元法高斯消元法、列主元素列主元素消去法消去法和和矩阵求逆法矩阵求逆法，而迭代法常用的有，而迭代法常用的有高斯高斯赛德尔迭代赛德尔迭代和和超（欠）松弛迭代超（欠）松弛迭代。 2022-3-1033nnnnjnjnnininjijiinnjjnnjjbTaTaTaTabTaTaTaTabTaTaTaTabTaTaTaTa.22112211222222121111212111迭代求解该方程组的思路迭代求解该方程组的思路为，寻找一个由（为，寻找一个由（T1，T2，Tn）组成的列向量，使其收敛于某一）组成的列向量，使其收敛于某一个极限向量（个极限向量（T1*,

27、 T2*,，Tn*）,且该极限向量且该极限向量就是该就是该方程的精确解方程的精确解。 2022-3-1034当 这 个 线 性 代 数 方 程 组 的 系 数 项当 这 个 线 性 代 数 方 程 组 的 系 数 项 ai i 0（i=1,2,n）时，可将其改写成迭代形式，有：）时，可将其改写成迭代形式，有： nnnnnjnjnnnniininjijiiiinnjjnnjjaTaTaTaTabTaTaTaTaTabTaTaTaTaTabTaTaTaTaTabT/ ).(./ ).(./ ).(/ ).(11221122112222323121221111313212112022-3-1035步

28、骤是，合理选择（假设）各节点的初始温度，步骤是，合理选择（假设）各节点的初始温度，将其作为第零次迭代的近似温度值，记为将其作为第零次迭代的近似温度值，记为Ti(0) (i=1,2,n);将将Ti(0)代入上式的右端，得到第一次迭代的近代入上式的右端，得到第一次迭代的近似值似值Ti(1) ；之后将之后将Ti(1)再代入上式的右端，则得出第二次再代入上式的右端，则得出第二次的近似值的近似值Ti(2) ；如此反复进行下去，直至进行到如此反复进行下去，直至进行到K次，使相邻次，使相邻的两次近似解的两次近似解Ti(K+1)和和Ti(K) (i=1,2,n)之间的偏之间的偏差小于预先设定的小量差小于预先设

29、定的小量时，即满足时，即满足 Ti(K+1) Ti(K) 或或 (Ti(K+1) Ti(K) )/ Ti(K) 2022-3-10361 交错网格系统交错网格系统 2 通用输运方程及离散化通用输运方程及离散化 3 压力修正方程：压力修正方程：SIMPLE算法算法 4 紊流壁面法则紊流壁面法则 2022-3-1037通用输运方程通用输运方程 *SyryxurxtrSoftware for Arbitrary Integration of Navier-Stokes Equation with a Turbulence and Porous Media Simulator Saints2D中的坐标

30、系统与控制方程中的坐标系统与控制方程2022-3-1038Saints2D中的非独立变量中的非独立变量 Dependent variables 2022-3-1039选用适当的参数选用适当的参数refrefrefTuL和,带带*的特征值表示无量纲量的特征值表示无量纲量 )/(,/,/*refrefrefrefuLttLyyLxxrefrefrefuwwuuuu/,/,/*2*/ )(refrefuppprefrefTTTT/ )(*2*/refukk )/(3*refrefLu)/(*refrefuL2*/refLKK 平面问题轴对称问题1*yr 2022-3-1040Saints2D中的标准

31、源项中的标准源项S* Normalized source terms 2022-3-10412*2*2*2*2*2*)(222rwrrxwrxyuyxuPt其中其中*2*2*1ReyTcxTcGrGxxTt 特征长度特征长度refL 选取任一有代表性的长度如平板长度、管子选取任一有代表性的长度如平板长度、管子直径直径 Reference length 有量纲和无量纲计算有量纲和无量纲计算2022-3-1042 特征速度特征速度refu强迫对流强迫对流: 选取一个确定的速度尺度，如进口选取一个确定的速度尺度，如进口处平均流速处平均流速 自然对流：考虑浮力与惯性力之间的平衡关自然对流：考虑浮力与惯

32、性力之间的平衡关系，可取系，可取混合对流：选取两者中的较大者，或是有助于混合对流：选取两者中的较大者，或是有助于得到较好结果的那一个得到较好结果的那一个 refrefrefLTguReference velocity2022-3-1043 特征温差特征温差refT已知壁面温度时：选取进口处平均温度与壁面温已知壁面温度时：选取进口处平均温度与壁面温度之差度之差 已知壁面热流密度时：考虑壁面热流密度与对流项已知壁面热流密度时：考虑壁面热流密度与对流项的平衡的平衡 wBrefTTTrefwTTTT/ )(*refpwrefucqT/注意：无量纲形式的热流密度为注意：无量纲形式的热流密度为PrRe0*

33、NuStuTcqnTqrefrefpwnwReference temperature difference2022-3-1044 其它无量纲数其它无量纲数 /RerefrefLu23/refrefLTgGr2*/refLKDaK雷诺数雷诺数 格拉晓夫数格拉晓夫数 Darcy数数 Dimensionless numbers2022-3-1045 有量纲计算有量纲计算 简单地将所有特征数设为它们各自的单位量简单地将所有特征数设为它们各自的单位量 suLtCTsmumLrefrefreforefrefref1/,1,/1,1前述无量纲数则为前述无量纲数则为 )/(1Re2smLurefref2222

34、3)/()/1 ()/(smKsmgLTgGrrefrefDimensional solutions 2022-3-1046基本思想基本思想：根据边界处速度矢量是已知还是未知来：根据边界处速度矢量是已知还是未知来划分边界条件的类型。划分边界条件的类型。1 速度已知边界、速度未知边界的概念速度已知边界、速度未知边界的概念2022-3-1047 速度已知边界速度已知边界： 在边界处，除压力已知外，速度矢量及所有变量的在边界处，除压力已知外，速度矢量及所有变量的值或密度都是已知的值或密度都是已知的 。固体壁面固体壁面x边界：边界： 固体壁面固体壁面y边界：边界： constconstxTconstTu, 0或constconstTconstu,流动进口边界：流动进口边界： constconstyTconstTu, 0或Known-velocity boundary 2022-3-1048 速度未知边界速度未知边界： 在