CFD基础课程系列-第五章（3）

CAE大拿

CFD基础课程系列

第五章  热流体仿真基础知识（3）

     在这个《CFD基础课程系列》里，针对刚刚开始，或者将要开始进行热流体仿真的工程师，我们尽量通过通俗易懂的语言和直观的现象来阐述CFD的概念。在第五章的第一部分，介绍了热流体的运动方程式，有限体积法的概念以及仿真区域选定的思想方法。在第二部分，我们介绍了热流体仿真中必不可少的仿真空间的网格划分，边界条件和初始条件设置的思想方法。在第五章的第三部分我们将介绍热流体仿真的计算方法。着重介绍具有代表性的湍流的计算模型以及它们的适用范围。

5.7矩阵的计算方法

     如前所说，各个单元的物理量是由邻接单元的值计算而来。如图5.20所示，单元3的物理量的计算将由以下算式得出。

                        

图5.20 邻接单元的位置关系和方程式

     上图是2维问题，所以上下左右的邻接单元有四个。3维问题的话，前后方向也有邻接单元，所以有六个邻接单元。如果是流动方程，常数部分为外力，如果是能量方程，常数部分为发热量的值。

     对每个单元都建立了这个关系。所以如果有100万个单元的话，需要求解100万个联立方程式。这个方程式可以用矩阵形式表示，矩阵的英语是“matrix”，所以方程式的求解也被称为“矩阵求解”。矩阵的求解方法有直接法和迭代法两大类。热流体仿真一般采用迭代法。

     这里，为了简单说明，我们采用以下2维的联立方程式，介绍迭代法中比较早且简单的雅克比迭代法。

           2x – y= 2

             x – y= –1                                (1)

这个联立方程式的解是x = 3, y = 4。

以下，为了方便理解，式（1）变形为

           x = 0.5y  + 1  

           y = x + 1                                   (2)

            

在用迭代法求解时，首先适当地假设x和y的初始值。这里我们把x和y的初始值定为0并代入式（2）。

          x = 0.5 × 0 + 1 = 1

          y = 0 + 1 = 1                                (3)

由此，第一次迭代得到x = 1, y = 1 。再次将x和y值代入式（2）。

   

          x = 0.5 × 1 + 1 = 1.5

          y = 1 + 1 = 2                                  (4)

                                                                       

第二次迭代得到的解是x = 1.5, y = 2 ，比第一次迭代的结果x =１, y = 1 更接近正确解。这样的操作重复进行的话，x和y的解会如下表所示。

表5.1 雅克比迭代法的计算结果

从表5.1可以看到，迭代过程重复进行的话，x和y的值将逼近联立方程的解。

     这里为了解说的简单，我们利用雅克比法来说明迭代法的概念。一般来说会采用更先进，收敛更高效的迭代方法。迭代次数越多，结果与正解（严密解）就越接近。但是要得到与正解（严密解）完全一致的结果并不容易。而且对于许多问题，我们并不知道什么是正解。实际计算时，当第k-1次的迭代结果与第k次的计算结果的差小于标准值时，迭代计算就终止了。这个标准值被称为收敛判定值。需要注意的是，在收敛不充分的情况下得到的结果很有可能与实际结果相差很大。

5.8时间历程和时间步长

     热流体仿真有定常分析和非定常分析两种。

     在定常分析时，只计算随时间的变化结果不再变化为最终状态的结果。因此，定常分析中没有时间的概念，达到最终状态之前的中间结果没有实际的物理意义。只有达到了最终状态的结果才有意义。

     而在非定常分析时，整个时间过程被切分成许多时间段，现象从某个时间瞬态到下一个瞬态一步一步地计算。因为需要计算现象的整个时间历程，中间结果（每个时间瞬态的结果）是有意义的。可以得到从初始条件到最终状态的变化过程。

     从一个瞬间到下一个瞬间的时间间隔被称为时间步长。在仿真一定时间段的现象变化时，时间步长越大，时间瞬间的计算次数就越少，可以减少一定時間段里现象变化的计算次数，但是代价是计算精度会下降。因此，在非定常分析时，设定合理的时间步长是很重要的。

        众所周知，一天后的天气预报会比一周以后天气预报的精度要高。从这个事实，就很容易理解时间步长越大精度越差的道理了。

      在设定时间步长时可以参考由下式定义的参数，这个参数被称为克朗数（Courant数），也被称为CFL数。

          C=uΔt/ΔL                                      (5)

C是克朗数，u是流速，Δt是时间步长、ΔL是单元长度。

图5.21 克朗数相关参数

克朗数的物理意义是，在一个设定的计算步长中流体经过了多少个单元。

     比如，如果克朗数为1，如图5.22所示，在一个时间步长里，流体正好从单元的左边流到单元的右边，在下一个步长，流体会流入邻接的单元。这样的网格密度可以得到比较高的计算精度。

    如果克朗数为10，如图5.22所示，流体在一个步长要通过10个单元，其中通过9个单元的过程没有计算。这与克朗数为1，网格密度扩大10倍所能得到的计算精度差不多。

图5.22 时间步长（克朗数）与流动距离

     克朗数很大程度取决于单元的尺寸和通过这个单元的流动速度。因此，即便是同样的流速单元小的地方，或者同样的单元尺寸流速较大的地方，克朗数都可能增大。

     很多场合，仿真对象的单元尺寸和流速在不同的位置不一样，克朗数也随位置的不同而变化。在设定时间步长时，在考虑节省计算时间的同时，应该尽量使可能出现的最大克朗数不要过大。

5.9湍流解析

     在基础课程（3.3）讨论层流和湍流时讲述到现实中大部分的流动属于湍流。图5.23是湍流的大致模样。一旦发生的涡流会不断地分裂产生更小的涡流，最终由于粘性转换成热能而消失。在这个过程中流场里存在大大小小的涡流，使得流动和温度分布受到影响。

图5.23 湍流的示意图

湍流的描述方程根据湍流处理尺寸不同，可以分成3大类计算模型。

・雷诺平均法（Reynolds Average Navier-Stokes RANS）

・大涡模拟法（Large Eddy Simulation LES）

・直接数值模拟法（Direct Numerical Simulation DNS）

以下，对这些方法做简单的说明。

5.9.1雷诺平均法（Reynolds Averaged Navier-Stokes, RANS）

     热流体分析的时候，工程意义上大部分场合是关心平均的流动和温度分布，或者对物体加载的力的平均值。如图5.24所示，本来非定常的湍流现象可以用时间平均的方法来考察。这个平均化的操作过程被称为雷诺平均。

图5.24 雷诺平均

     基于经过雷诺平均后的Navier-Stokes方程式求解平均流动的方法被称为雷诺平均法（Reynolds Average Navier-Stokes），取英语的第一个字母，简称RANS法。

     在基于经过雷诺平均后的方程式里，包含了表示湍流影响的雷诺应力项。这一项表达了平均流动的湍流效果，但是仅仅靠雷诺平均后的方程式还不能求得雷诺应力，还需要有表现雷诺应力的近似方法。这种近似方法一般称为湍流模型，具有代表性的是标准k-ε模型。

     标准k-ε湍流模型无法模拟随时间变化的微小脉动或湍流涡（比如瞬态流动释放出来的涡，像卡门涡），上述现象不能用雷诺平均方法来模拟。RANS最主要的优势是计算负荷低，它是上述3种方法中计算负荷最低的计算方法，在实际工程中应用广泛。

5.9.2 大涡模拟法（Large Eddy Simulation LES）

     相对小漩涡，大漩涡携带更多的能量，对流场的影响更大。这个方法的基本操作就是低通滤波，对基础方程进行空间过滤，对过滤尺寸以上的大涡流不加近似，用直接法模拟，而对过滤尺寸以下的小涡流使用简化模型。简单地说用单元尺度进行筛选，大于单元尺寸的漩涡用直接法求解。这个方法被称为大涡模拟法（Large Eddy Simulation），取英语第一个字母，简称LES法。

      一般来说，以单元的尺寸作为过滤尺度，用直接法计算的大尺寸称为Grid Scale，GS，用平均模型计算的小尺寸成为subgrid scale, SGS。

图5.25 大涡模拟法的过滤尺度：Grid Scale和 Subscale Grid

网格分割的越细，利用平均法计算的比例就越小，虽然计算负荷会增大，可以得到与现实更接近的结果。

5.9.3 直接数值模拟法（Direct Numerical Simulation, DNS）

     湍流分析最单纯且严密的方法是对基础方程式不做任何简化直接进行求解。这样的解法称为直接数值模拟法，取英语的第一个字母，DNS法。

     虽然原理上最简单，但是用DNS进行正确的流场分析，需要大大提高网格的分别率，使其能够准确地捕捉不断分裂派生的小尺度涡流。这个最小涡流尺度取决于粘性大小的，被称为柯尔莫哥洛夫尺度（Kolmogorov scale）。雷诺数越大的流动，粘性的影响就越小，进而最小涡流的尺度就越小，仿真的网格密度就需要越大。3维DNS法所需的网格数可以由下式来预估。

                      (6)

N是网格数，Re是瑞诺数。如图5.26所示，用式（6）来预估几种典型仿真所需要的网格数。

图5.26 雷诺数与网格数的关系

    从图5.26可知，即便是雷诺数较低的流动，需要数千亿以上的网格规模。

     根据目前计算机的计算速度和内存，用DNS求解规模如此庞大的计算问题并不容易。因此DNS大多数用在以研究目的的情况。工程问题的实用计算主要还是采用前述的RANS和LES方法。

编后记

     这个CFD基础课程连载系列分9次将热流体仿真相关的概念，思想方法做了介绍。感谢读者一直支持和阅读这个连载。

     近年来，随着计算方法和计算机性能的日新月异，热流体仿真软件已经成为设计/开发/制造工艺不可或缺的工具之一。而且，在日常业务中需要进行热流体仿真的设计人员，工程师也在不断增加。学好用好热流体仿真软件并不容易，希望这个连载能帮助那些刚刚开始，或者将要开始进行热流体仿真的工程师加深对热流体仿真的概念和思想方法的理解，并从中获益。