华科传热学课件

2023/7/261传热学主讲：刘志春能源与动力工程学院华中科技大学2023/7/261传热学主讲：刘志春2023/7/262第九章流动与传热的数值计算§9-1

数值计算的基本思想*§9-2

流动与传热的数值计算§9-3Saints2D软件简介2023/7/262第九章流动与传热的数值计算§9-12023/7/263首先，我们以导热问题为例，介绍计算区域离散化的概念、内节点与边界节点方程式的建立方法、节点方程组的求解过程，以及非稳态导热问题的显示与隐示差分格式。然后，介绍在上述思想的基础上开发的流动与传热计算软件Saints2D，并给出传热问题虚拟实验的计算示例。2023/7/263首先，我们以导热问题为例，介绍计算区域离2023/7/264§9-1数值计算的基本思想数值求解通常是对微分方程直接进行数值积分或者把微分方程转化为一组代数方程组再进行求解。这里要介绍的是后一种方法。如何实现从微分方程到代数方程的转化又可以采用不同的数学方法，如有限差分法、有限元法和边界元法等。这里仅向读者简要地介绍用有限差分方法从微分方程确立代数方程的处理过程。2023/7/264§9-1数值计算的基本思想数值求解通常2023/7/265有限差分法的基本思想是把原来在时间和空间坐标中连续变化的物理量（如温度、压力、速度和热流等），用有限数目的离散点上的数值集合来近似表达。有限差分的数学基础是用差商代替微商（导数）。几何意义是用函数在某区域内的平均变化率代替函数的真实变化率。2023/7/265有限差分法的基本思想是把原来在时间和空间2023/7/266在图中可以看出有限差分表示的温度场与真实温度场的区别。图中用T0、T1、T2…表示连续的温度场T；Δx为步长，它将区域的x方向划分为有限个数的区域，Δx0、Δx1、Δx2…，它们可以相等，也可以不相等。当Δx相等时，T1处的真实变化率a可以用平均变化率b、c或d来表示，其中b、c和d分别表示三种不同差分格式下的温度随时间的变化率Δx1Δx2Δx3Δx0xTT3T2T1T0bdca2023/7/266在图中可以看出有限差分表示的温度场与真实2023/7/267b为向后差分格式c为向前差分格式

d为中心差分格式

Δx1Δx2Δx3Δx0xTT3T2T1T0bdca2023/7/267b为向后差分格式c为向前差分格式d为中2023/7/268这种差分格式可推广到高阶微商（导数）。对于二阶导数的差分格式可以在一阶差分格式的基础上得出：这样处理后，反映温度场随时间、空间连续变化的微分方程就可以用反映离散点间温度线性变化规律的代数方程来表示。当利用相应的数学办法求解这些代数方程组之后，我们就能获得离散点上的温度值。这些温度值就可以近似表示温度场的连续的温度分布。2023/7/268这种差分格式可推广到高阶微商（导数）。这2023/7/269从上面的分析不难看出，当我们要对流动与传热问题进行数值求解时一定要采取三个大的步骤，即：a)

研究区域的离散化；b)散点（节点）差分方程的建立；c)节点方程（代数方程）的求解。2023/7/269从上面的分析不难看出，当我们要对流动与传2023/7/26101时间与空间的离散化

进行数值求解时,首先是在所研究的时间和空间区域内把时间和空间分割成为有限大小的小区域。图表示了长柱体矩形截面上区域离散化的情况。K-1时刻SWENPSWENPSWENPτxyK时刻K+1时刻ΔxΔy2023/7/26101时间与空间的离散化进行数值求解时2023/7/2611对于给定的空间区域，在x方向上的步长为Δx，在y方向上的步长为Δy，用它们作为空间尺度可以将矩形区域划分成纵横交错的网格系统。计算区域就被这些网格线分隔成一系列的小的区域，称为控制面积，对于三维情况则为控制体积或控制容积，因而在一般意义上称之为控制体；控制体的中心点称为节点。K-1时刻SWENPSWENPSWENPτxyK时刻K+1时刻ΔxΔy2023/7/2611对于给定的空间区域，在x方向上的步长为2023/7/2612控制体的形状是随着坐标系的不同而改变的，这里的控制体是一个个的矩形面积。网格的步长在每一个方向上可以均匀划分，也可以不均匀的划分；所得到网格，相应地被称为均匀网格或者非均匀网格。选用不同的步长和不同的划分方法，可以将同一区域划分出不同大小、不同数目的控制区域，以及不同数目的节点数。2023/7/2612控制体的形状是随着坐标系的不同而改变的2023/7/2613获得每个节点上的温度值，就是导热数值计算的目的。，随着步长的不断减小，节点数目的不断增加，由节点温度表示的离散温度场就更接近连续温度场，但计算工作量也会增加。在时间方向上离散化的步长常用Δτ来表示，Δτ的选取也是可大可小的，也可以随时间的进程而变化。显然，无限小的时间步长Δτ亦会使得离散温度变化接近连续的温度改变，但随之而来的是相应的计算工作量将会增加。2023/7/2613获得每个节点上的温度值，就是导热数值计2023/7/26142节点方程的建立建立节点差分方程可采用不同的方法，主要分为两大类：第一类包括泰勒(Taylor)级数展开法和多项式拟合法，它偏重于从数学的角度进行推导，其优点是便于对离散方程进行数学特性分析，但缺点是变步长网格的离散方程形式复杂、导出过程的物理概念不清晰、不能保证差分方程具有守恒特性。2023/7/26142节点方程的建立建立节点差分方程可采2023/7/2615第二类包括控制体热平衡法和控制容积积分法，其优点是推导过程的物理概念清晰、离散方程系数具有一定物理意义、保证差分方程具有守恒特性，但缺点是不便于对离散方程进行数学特性分析。下面我们采用控制体热平衡法来建立节点方程。2023/7/2615第二类包括控制体热平衡法和控制容积积分2023/7/2616①

内节点方程控制体热平衡法建立节点方程的过程是将能量守恒方程应用于控制体，建立该节点与周围节点之间的能量平衡关系式。(i+1,j)E(i,j)oyx(i-1,j)W(i,j-1)S

(i,j+1)NPxxyy再利用傅里叶导热定律，最后获得控制体节点温度与周围节点温度之间的关系式。2023/7/2616①内节点方程控制体热平衡法建立节点方2023/7/2617考察图中的节点P及其控制体，由能量平衡关系应有式中，Ф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)NPxxyyΔE为控制体单位时间内热能的增加量。2023/7/2617考察图中的节点P及其控制体，由能量平衡2023/7/2618由导热傅里叶定律，在线性温度分布的假设下，时刻K周围节点传给节点P的热流量分别为：(i+1,j)E(i,j)oyx(i-1,j)W(i,j-1)S

(i,j+1)NPxxyy2023/7/2618由导热傅里叶定律，在线性温度分布的假设2023/7/2619控制体的发热流量其中qV为内热源强度，即单位时间单位体积的内热源发热量。控制体单位时间的内能增加量为前者为时间上的向前差分，而后者为时间上向后差分。以上关系式中温度T的上标为所在时刻，下标为所在空间位置。2023/7/2619控制体的发热流量控制体单位时间的内能增2023/7/2620假设Δx=Δy，经整理可以得出二维非稳态导热问题的内节点的两种差分格式的差分方程，即a)

显式差分格式

定义网格傅里叶数

其物理意义是表征控制体的导热性能与热储蓄性能之间的对比关系，反映控制体温度随时间变化的动态特性。显式差分格式简化为2023/7/2620假设Δx=Δy，经整理可以得出二维非稳2023/7/2621b)

隐式差分格式

或改写为2023/7/2621b)隐式差分格式或改写为2023/7/2622显示差分格式最突出的优点是节点温度表达式的右边只涉及K时刻(前一时刻)的节点温度值，只要知道前一时刻周围节点的温度值就可以求出该节点的K+1时刻(当前时刻)的温度值；2023/7/2622显示差分格式最突出的优点是节点温度表达2023/7/2623隐示差分格式温度表达式的右端除了K时刻(前一时刻)的节点温度值以外，还含有K+1时刻(当前时刻)的温度值，这就意味着必须同时计算当前时刻所有节点的温度值，即必须联立求解K+1时刻所有节点的差分方程组，计算工作量增大也就是显而易见的了。2023/7/2623隐示差分格式温度表达式的右端除了K时刻2023/7/2624虽然显示差分格式计算比较方便，但它存在一个缺点，即计算式中FOΔ值必须满足一定的条件才不至于引起数值计算出现不收敛的问题，这在数值计算中称为差分格式的不稳定性。这里差分方程稳定性的条件是式(9-2)中的变量T前面的系数必须大于或等于零，分析一下差分方程中的各项系数，有2023/7/2624虽然显示差分格式计算比较方便，但它存在2023/7/2625此式称为显示差分格式的稳定性判据，从中看出时间步长和空间步长是相互制约的。为了获得较为精确的节点温度值，空间步长Δx的选择不能太小，按照稳定性判据的要求势必会使时间步长Δτ也要相应地不能太大，因而必须在增加节点数目的同时增多时间间隔，从而使计算工作量加大。隐示差分格式是无条件稳定的。2023/7/2625此式称为显示差分格式的稳定性判据，从中2023/7/2626这里指出，以上的讨论及结果适用于对非稳态导热问题。对于稳态导热问题，其实应更简单,只需要在式(9-2)或者式(9-3)中令或均可得到二维稳态导热问题的内节点方程式2023/7/2626这里指出，以上的讨论及结果适用于对非稳2023/7/2627②

边界节点方程以对流换热边界为例，从流体侧来看，应用牛顿冷却公式，流体与壁面之间的对流换热量再者，从控制体侧来看，应用傅里叶定律，假定壁面处的温度梯度取向后差分格式，则应有SNWPT∞ΔyΔxh2023/7/2627②边界节点方程以对流换热边界为例，从2023/7/2628只要网格步长Δx足够小，两者的结果应该是一致的。从而可消去未知量Te，得到代入式(9-1)，得到对流换热边界节点的两种差分格式，即a)

显式差分格式显式差分格式的稳定性判据为2023/7/2628只要网格步长Δx足够小，两者的结果应该2023/7/2629b)

隐式差分格式隐式差分格式仍然是无条件稳定的定义网格毕渥数其物理意义是体现控制体和环境间的换热性能与其导热性能之间的对比关系。2023/7/2629b)隐式差分格式隐式差分格式仍然是2023/7/2630这里指出，上面的结果是在对流换热边界的情况下得到的，但经过简单处理，可直接用于绝热边界条件与恒壁温边界条件。令，则有即可简化为恒壁温边界条件下对应的差分格式。而令，即即可简化为绝热边界条件下对应的差分格式。2023/7/2630这里指出，上面的结果是在对流换热边界的2023/7/2631由上面的讨论可以看出，对应于离散温度场的每一个节点均可以列出相应的差分方程，这样就可以得出与节点数目相同的一组代数方程组。当联立求解这个代数方程组时，最后就可以得出每一个节点的温度值。2023/7/2631由上面的讨论可以看出，对应于离散温度场2023/7/26323节点方程的求解

一般情况下，差分方程组是线性代数方程组，而线性代数方程组是可以用直接法和迭代法求解的。常用的直接法有高斯消元法、列主元素消去法和矩阵求逆法，而迭代法常用的有高斯－赛德尔迭代和超（欠）松弛迭代。2023/7/26323节点方程的求解一般情况下，差分方2023/7/2633迭代求解该方程组的思路为，寻找一个由（T1，T2，…，Tn）组成的列向量，使其收敛于某一个极限向量（T1*,T2*,…，Tn*）,且该极限向量就是该方程的精确解。2023/7/2633迭代求解该方程组的思路为，寻找一个由（2023/7/2634当这个线性代数方程组的系数项aii≠0（i=1,2,…,n）时，可将其改写成迭代形式，有：2023/7/2634当这个线性代数方程组的系数项aii≠02023/7/2635步骤是，合理选择（假设）各节点的初始温度，将其作为第零次迭代的近似温度值，记为Ti(0)(i=1,2,…,n);将Ti(0)代入上式的右端，得到第一次迭代的近似值Ti(1)

；之后将Ti(1)再代入上式的右端，则得出第二次的近似值Ti(2)

；如此反复进行下去，直至进行到K次，使相邻的两次近似解Ti(K+1)和Ti(K)(i=1,2,…,n)之间的偏差小于预先设定的小量ε时，即满足∣Ti(K+1)

－Ti(K)

∣≤ε或∣(Ti(K+1)

－Ti(K))/Ti(K)

∣≤ε2023/7/2635步骤是，合理选择（假设）各节点的初始温2023/7/2636*§9-2流动与传热的数值计算

1

交错网格系统2

通用输运方程及离散化3压力修正方程：SIMPLE算法4

紊流壁面法则2023/7/2636*§9-2流动与传热的数值计算12023/7/2637通用输运方程SoftwareforArbitraryIntegrationofNavier-StokesEquationwithaTurbulenceandPorousMediaSimulator

§9-3Saints2D软件简介

Saints2D中的坐标系统与控制方程2023/7/2637通用输运方程Softwarefor2023/7/2638Saints2D中的非独立变量Dependentvariables2023/7/2638Saints2D中的非独立变量Dep2023/7/2639选用适当的参数带*的特征值表示无量纲量2023/7/2639选用适当的参数带*的特征值表示无量纲量2023/7/2640Saints2D中的标准源项S*

Normalizedsourceterms2023/7/2640Saints2D中的标准源项S*No2023/7/2641其中①

特征长度

选取任一有代表性的长度如平板长度、管子直径

Referencelength有量纲和无量纲计算2023/7/2641其中①特征长度选取任一有代表性的长2023/7/2642②

特征速度强迫对流:选取一个确定的速度尺度，如进口处平均流速自然对流：考虑浮力与惯性力之间的平衡关系，可取混合对流：选取两者中的较大者，或是有助于得到较好结果的那一个Referencevelocity2023/7/2642②特征速度强迫对流:选取一个确定的2023/7/2643③特征温差已知壁面温度时：选取进口处平均温度与壁面温度之差已知壁面热流密度时：考虑壁面热流密度与对流项的平衡注意：无量纲形式的热流密度为Referencetemperaturedifference2023/7/2643③特征温差已知壁面温度时：选取进口处2023/7/2644④其它无量纲数

雷诺数格拉晓夫数Darcy数Dimensionlessnumbers2023/7/2644④其它无量纲数雷诺数格拉晓夫数2023/7/2645⑤有量纲计算简单地将所有特征数设为它们各自的单位量前述无量纲数则为Dimensionalsolutions2023/7/2645⑤有量纲计算简单地将所有特征数设为2023/7/2646基本思想：根据边界处速度矢量是已知还是未知来划分边界条件的类型。1速度已知边界、速度未知边界的概念2023/7/2646基本思想：根据边界处速度矢量是已知还是2023/7/2647①速度已知边界：在边界处，除压力已知外，速度矢量及所有变量的值或密度都是已知的。固体壁面x边界：固体壁面y边界：流动进口边界：Known-velocityboundary2023/7/2647①速度已知边界：在边界处，除压力已2023/7/2648②速度未知边界：在边界处已知速度是正态分布的，但其绝对值未知（变量值或密度都是未知的）。流动出口x边界：流动出口y边界：Unknown-velocityboundary2023/7/264