1.2 CFD基本模型

流体流动所遵循的物理定律，是建立流体运动基本方程组的依据。这些定律主要包括质量守恒、动量守恒、动量矩守恒、能量守恒、热力学第二定律，加上状态方程、本构方程。在实际计算时，还要考虑不同的流态，如层流与湍流。

1.2.1 基本控制方程

1．系统、控制体与常用运算符

在流体力学中，系统是指某一确定流体质点集合的总体。系统以外的环境称为外界，分隔系统与外界的界面，称为系统的边界。系统通常是研究的对象，外界则用来区别于系统。系统将随系统内质点一起运动，系统内的质点始终包含在系统内，系统边界的形状和所围空间的大小，则可随运动而变化。系统与外界无质量交换，但可以有力的相互作用及能量（热和功）交换。

控制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不动形状任意的空间体积。包围这个空间体积的边界面，称为控制面。控制体的形状与大小不变，并相对于某坐标系固定不动。控制体内的流体质点组成并非不变的。控制体既可通过控制面与外界有质量和能量交换，也可与控制体外的环境有力的相互作用。

本书用到一些数学运算符如下所示

上述公式中，称为哈密顿算子。

式中，Δ＝▽•▽＝▽2称为拉普拉斯算子。

2．质量守恒方程（连续性方程）

在流场中，流体通过控制面A1流入控制体，同时也会通过另一部分控制面A2流出控制体，在这期间控制体内部的流体质量也会发生变化。按照质量守恒定律，流入的质量与流出的质量之差，应该等于控制体内部流体质量的增量，由此可导出流体流动连续性方程的积分形式为

式中，V表示控制体；A表示控制面。等式左边第一项表示控制体V内部质量的增量；第二项表示通过控制表面流入控制体的净通量。

根据数学中的奥-高公式，在直角坐标系下可将其化为微分形式如下

对于不可压缩均质流体，密度为常数，则有

对于圆柱坐标系，其形式为

对于不可压缩均质流体，密度为常数，则有

3．动量守恒方程（运动方程）

动量守恒是流体运动时应遵循的另一个普遍定律，描述为：在一给定的流体系统，其动量的时间变化率等于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动方程，或N-S方程，其微分形式表达如下

式中，Fbx，Fby，Fbz分别是单位质量流体上的质量力在三个方向上的分量；pyx，pzx，pxy等是流体内应力张量的分量。

动量守恒方程在实际应用中有许多表达形式，其中比较常见有如下几种：

（1）可压缩黏性流体的动量守恒方程

（2）常黏性流体的动量守恒方程

（3）常密度常黏性流体的动量守恒方程

（4）无黏性流体的动量守恒方程（欧拉方程）

（5）静力学方程

（6）相对运动方程

在非惯性参考系中的相对运动方程是研究像大气、海洋及旋转系统中流体运动所必须考虑的。由理论力学得知，绝对速度va为相对速度vr及牵连速度ve之和：

式中，ve＝v0＋Ω×r，v0为运动系中的平动速度，Ω是其转动角速度，r为质点矢径。

而绝对加速度aa为相对加速度ar、牵连加速度ae及科氏加速度ac之和：

其中。

将绝对加速度代入运动方程（1-49），即得到流体的相对运动方程

4．能量守恒方程

将热力学第一定律应用于流体运动，把上式各项用有关的流体物理量表示出来，即是能量方程。

式中，；keff是有效热传导系数，keff＝k＋kt，其中kt是湍流热传导系数，根据所使用的湍流模型来定义；Jj′是组分j的扩散流量；Sh包括了化学反应热以及其他用户定义的体积热源项；上面方程右边的前三项分别描述了热传导、组分扩散和黏性耗散带来的能量输运。

1.2.2 湍流模型

湍流是自然界广泛存在的流动现象。大气、海洋环境的流动；飞行器和船舰的绕流；叶轮机械、化学反应器、核反应器中的流体运动都是湍流。湍流流动的核心特征是其在物理上近乎无穷多的尺度和数学上强烈的非线性，这使得人们无论是通过理论分析、实验研究还是计算机模拟来彻底认识湍流都非常困难。回顾计算流体力学的发展，特别是活跃的20世纪80年代，人们提出和发展了一大批高精度、高分辨率的计算格式，相当成功地解决了Euler方程的数值模拟，可以说Euler方程数值模拟方法的精度已接近于它有效使用范围的极限；同时，研究者还发展了一大批有效的网格生成技术及相应的软件，具体实现了工程计算所需要的复杂外形的计算网格；且随着计算机的发展，无论从计算时间还是从计算费用考虑，Euler方程都已能适用于各种实践所需。在此基础上，研究者于20世纪80年代还进行了求解可压缩雷诺平均方程及其三维定态黏流流动的模拟。90年代业界又开始了一个非定常黏流流场模拟的新局面，这里所说的黏流流场具有高雷诺数、非定常、不稳定、剧烈分离流动的特点，显然需要继续探求更高精度的计算方法和更实用可靠的网格生成技术。但研究湍流机理，建立相应的模式，并进行适当的模拟仍是解决湍流问题的重要途径。

1．湍流模型分类

湍流流动模型很多，但大致可以归纳为以下三类。

第一类是湍流输运系数模型，即将速度脉动的二阶关联量表示成平均速度梯度与湍流黏性系数的乘积，用笛卡儿张量表示为

模型的任务就是给出计算湍流黏性系数μt的方法。根据建立模型所需要的微分方程的数目，可以分为零方程模型（代数方程模型），单方程模型和双方程模型。

第二类是抛弃了湍流输运系数的概念，直接建立湍流应力和其他二阶关联量的输运方程。

第三类是大涡模拟。前两类是以湍流的统计结构为基础，对所有涡旋进行统计平均。大涡模拟把湍流分成大尺度湍流和小尺度湍流，通过求解三维修正的Navier-Stokes方程，得到大涡旋的运动特性，而对小涡旋运动还采用上述的模型。

实际求解中，选用什么模型要根据具体问题的特点来决定。选择的一般原则是精度高，应用简单，节省计算时间，具有通用性。

Fluent提供的湍流模型包括单方程（Spalart-Allmaras）模型、双方程模型（标准k-ε模型、重整化群（RNG）k-ε模型、可实现（Realizable）k-ε模型）、k-ω（包括Standard k-ω和SST k-ω）模型及雷诺应力模型和大涡模拟（见图1-1）。

2．平均量输运方程

雷诺平均就是把Navier-Stokes方程中的瞬时变量分解成平均量和脉动量两部分。对于速度，有

式中，和分别是平均速度和脉动速度（i＝1，2，3）。

类似地，对于压力等其他标量，也有

式中，ϕ表示标量，如压力、能量、组分浓度等。

把上面的表达式代入瞬时的连续与动量方程，并取平均（去掉平均速度上的横线），可以把连续与动量方程写成如下的笛卡儿坐标系下的张量形式：

上面两个方程称为雷诺平均的Navier-Stokes（RANS）方程。它们和瞬时Navier-Stokes方程有相同的形式，只是速度或其他求解变量变成了时间平均量。额外多出来的项是雷诺应力，表示湍流的影响。

对于密度变化的流动过程，如燃烧问题，需要采用法夫雷（Favre）平均才可以求解。法夫雷平均就是除了压力和密度本身以外，所有变量都用密度加权平均。变量的密度加权平均定义如下

式中，符号～表示密度加权平均，对应于密度加权平均值的脉动值用表示，有。显然，这种脉动值的简单平均值不为零，但它的密度加权平均值等于零，即：。

为了求解方程（1-62），必须模拟雷诺应力项以使方程封闭。通常的方法是应用Boussinesq假设，认为雷诺应力与平均速度梯度成正比，表达如下

Boussinesq假设被用于Spalart-Allmaras单方程模型和k-ε双方程模型。这种近似方法好处是与求解湍流黏性系数有关的计算时间比较少。例如，在Spalart-Allmaras单方程模型中只多求解一个表示湍流黏性的输运方程；在k-ε双方程模型中只需多求解湍动能k和耗散率ε两个方程，湍流黏性系数用湍动能k和耗散率ε的函数来描述。Boussinesq假设的不足之处是假设μt是个各向同性标量，对于一些复杂流动，该条件并不严格成立，所以具有其应用局限性。

另外的近似方法是求解雷诺应力各分量的输运方程。这也需要额外再求解一个标量方程，通常是耗散率ε方程。这就意味着对于二维湍流流动问题，需要多求解4个输运方程，而三维湍流问题需要多求解7个方程，需要较多的计算时间，要求更高的计算机内存。

在很多情况下基于Boussinesq假设的模型很好用，而且计算量并不是很大。但是，如果湍流场各向异性很明显，如强旋流动以及应力取得的二次流等流动中，求解雷诺应力RSM模型可以得到更好的结果。

3．常用湍流模型简介

1）单方程（Spalart-Allmaras）模型

单方程模型求解变量是，表征出了近壁（黏性影响）区域以外的湍流运动黏性系数。的输运方程为

式中，Gv是湍流黏性产生项；Yv是由于壁面阻挡与黏性阻尼引起的湍流黏性的减少；和Cb2是常数；ν是分子运动黏性系数。

湍流黏性系数，其中，fv1是黏性阻尼函数，定义为：。而湍流黏性产生项Gv模拟如下：，其中，Cb1和k是常数，d是计算点到壁面的距离；。在Fluent软件中，考虑到平均应变率对湍流产生也起到很大作用，，其中，，平均应变率。

在涡量超过应变率的计算区域计算出来的涡旋黏性系数变小。这适合涡流靠近涡旋中心的区域，那里只有“单纯”的旋转，湍流受到抑止。包含应变张量的影响更能体现旋转对湍流的影响。忽略了平均应变，估计的涡旋黏性系数产生项偏高。

湍流黏性系数减少项Yv为，其中，是常数，在计算r时用到的受平均应变率的影响。

上面的模型常数在Fluent软件中默认值为：。

2）标准k-ε模型（Standard k-ε）

标准k-ε模型需要求解湍动能及其耗散率方程。湍动能输运方程是通过精确的方程推导得到，但耗散率方程是通过物理推理，数学上模拟相似原形方程得到的。该模型假设流动为完全湍流，分子黏性的影响可以忽略。因此，标准k-ε模型只适合完全湍流的流动过程模拟。标准k-ε模型的湍动能k和耗散率ε方程为如下形式：

式中，Gk表示由于平均速度梯度引起的湍动能产生；Gb是由于浮力影响引起的湍动能产生；YM为可压速湍流脉动膨胀对总的耗散率的影响。湍流黏性系数。

在Fluent中，作为默认值常数，C1ε＝1.44，C2ε＝1.92，C3ε＝0.09，湍动能k与耗散率ε的湍流普朗特数分别为σk＝1.0，σε＝1.3。

3）重整化群k-ε模型（RNG k-ε）

重整化群k-ε模型是对瞬时的Navier-Stokes方程用重整化群的数学方法推导出来的模型。模型中的常数与标准k-ε模型不同，而且方程中也出现了新的函数或新项。其湍动能与耗散率方程与标准k-ε模型有相似的形式：

式中Gk、Gb、YM参数与标准k-ε模型中相同。αk和αε分别是湍动能k和耗散率ε的有效湍流普朗特数的倒数。湍流黏性系数计算公式为：，其中，。对于前面方程的积分，可以精确到有效雷诺数（涡旋尺度）对湍流输运的影响，这有助于处理低雷诺数和近壁流动问题的模拟。对于高雷诺数，上面方程可以给出：。这个结果和标准k-ε模型的半经验推导给出的常数Cμ＝0.09非常近似。在Fluent中，如果是默认设置，用重整化群k-ε模型时候是针对的高雷诺数流动问题。如果对低雷诺数问题进行数值模拟，必须进行相应的设置。

4）可实现k-ε模型（Realizable k-ε）

可实现k-ε模型的湍动能及其耗散率输运方程为

式中，，η＝Sk/ε；C2和C1ε是常数；σk，σε分别是湍动能及其耗散率的湍流普朗特数。在Fluent中，作为默认值常数，C1ε＝1.44，C2＝1.9，σk＝1.0，σε＝1.2。

在上述方程中，Gk表示由于平均速度梯度引起的湍动能产生，Gb表示用于浮力影响引起的湍动能产生，YM表示可压速湍流脉动膨胀对总的耗散率的影响。

该模型的湍流黏性系数与标准k-ε模型相同。不同的是，黏性系数中的Cμ不是常数，而是通过公式计算得到：，其中，，式中，为在角速度ωk旋转参考系下的平均旋转张量率；模型常数，。从这些式中发现，Cμ是平均应变率与旋度的函数。在平衡边界层惯性底层，可以得到Cμ＝0.09，与标准k-ε模型中采用的常数一样。

该模型适合的流动类型比较广泛，包括有旋均匀剪切流，自由流（射流和混合层），腔道流动和边界层流动。对以上流动过程模拟结果都比标准k-ε模型的结果好，特别是可实现k-ε模型对圆口射流和平板射流模拟中，能给出较好的射流扩张角。

双方程模型中，无论是标准k-ε模型、重整化群k-ε模型还是可实现k-ε模型，三个模型有类似的形式，即都有k和ε的输运方程，它们的区别在于：①计算湍流黏性的方法不同；②控制湍流扩散的湍流Prandtl数不同；③ε方程中的产生项和Gk关系不同。但都包含了相同的表示由于平均速度梯度引起的湍动能产生Gk，用于浮力影响引起的湍动能产生Gb；可压速湍流脉动膨胀对总的耗散率的影响YM。

湍动能产生项

式中，Prt是能量的湍流普特朗数，对于可实现k-ε模型，默认设置值为0.85，T为温度。对于重整化群k-ε模型，Prt＝1/α，α＝1/Prt＝k/μCp。热膨胀系数，对于理想气体，浮力引起的湍动能产生项变为

5）雷诺应力模型（RSM）

雷诺应力模型是求解雷诺应力张量的各个分量的输运方程。具体形式为

式中，为雷诺应力；u为速度，左边的第二项是对流项Cij；右边第一项是湍流扩散项；第二项是分子扩散项；第三项是应力产生项Pij；第四项是浮力产生项Gij；第五项是压力应变项ϕij；第六项是耗散项εij；第七项是系统旋转产生项Fij。

在式（1-76）中，不需要模拟，而需要模拟以封闭方程。下面简单对几个需要模拟项进行说明。

可以用Delay and Harlow的梯度扩散模型来模拟，但这个模型会导致数值不稳定，在Fluent中是采用标量湍流扩散模型：

式中，湍流黏性系数用来计算，根据Lien and Leschziner，σk＝0.82，这和标准k-ε模型中选取1.0有所不同。

压力应变项ϕij可以分解为三项，即

ϕij，1，ϕij，2和分别是慢速项，快速项和壁面反射项，具体表述可以见参考文献[2]。

浮力引起的产生项Gij模拟为

耗散张量εij模拟为

式中，是马赫数；标量耗散率ε用标准k-ε模型中的采用的耗散率输运方程求解。

6）大涡模拟（LES）

湍流中包含了不同时间与长度尺度的涡旋。最大长度尺度通常为平均流动的特征长度尺度。最小尺度为Komogrov尺度。LES的基本假设是：①动量、能量、质量及其他标量主要由大涡输运；②流动的几何和边界条件决定了大涡的特性，而流动特性主要在大涡中体现；③小尺度涡旋受几何和边界条件影响较小，并且各向同性；大涡模拟过程中，直接求解大涡，小尺度涡旋模拟，从而使得网格要求比DNS低。

LES的控制方程是对Navier-Stokes方程在波数空间或者物理空间进行过滤得到的。过滤的过程是去掉比过滤宽度或者给定物理宽度小的涡旋，从而得到大涡旋的控制方程。

其中，τij为亚网格应力，。

很明显，上述方程与雷诺平均方程很相似，只不过大涡模拟中的变量是过滤过的量，而非时间平均量，并且湍流应力也不同。

7）湍流模型选择策略

湍流模型选取的需要考虑的因素：流体是否可压、建立特殊的可行的问题、精度的要求、计算机的能力、时间的限制。

Fluent软件中各湍流模型的应用范围及特点如下，供读者选择时参考。

（1）Spalart-Allmaras模型

应用范围：这个模型被Spalart-Allmaras提出，用来解决因湍流动粘滞率而修改的数量方程，主要是壁面约束的流动，而且已经显示出很好的效果。该模型多用于航空领域、透平机械中，对于有壁面约束的空气流动问题应优先选用此模型。

特点：①Spalart-Allmaras模型是相对简单的单方程模型，只需求解湍流黏性的输运方程，不需要求解当地剪切层厚度的长度尺度；由于没有考虑长度尺度的变化，对一些流动尺度变换比较大的流动问题不太适合；比如平板射流问题，从有壁面影响流动突然变化到自由剪切流，流场尺度变化明显等问题。②该模型的输运变量在近壁处的梯度要比模型中的小，表现出对网格粗糙带来数值误差不太敏感。③该模型不能去预测均匀衰退，各向同性湍流等复杂的工程流动问题。

（2）k-ε模型

①标准的k-ε模型（Standard k-ε）

在Fluent中，标准k-ε模型系数由经验公式给出，计算比较稳定，自从被Launder和Spalding提出之后，就变成了工程分析计算中的主要工具。现有研究表明，对于简单的充分发展湍流问题，该模型完全适用，而且比其他湍流更快收敛。

应用范围：适合完全湍流（充分发展的湍流）的流动过程模拟，特别是对高Re的湍流有效，包含黏性热、浮力、压缩性等选项，但模拟旋流和绕流时有缺陷。

②重整化群k-ε模型（RNG k-ε）

重整化群k-ε模型来源于严格的统计技术，与标准模型相比，有一些改进：A. RNG k-ε模型在ε方程中加了一个条件，有效地改善了精度。B.考虑到了湍流漩涡，提高了在这方面的精度。C. RNG理论为湍流Prandtl数提供了一个解析公式，而标准k-ε模型是用户提供的常数。D.标准k-ε模型是一种高雷诺数的模型，RNG理论提供了一个考虑低雷诺数流动黏性的解析公式。这些公式的作用取决于正确地对待近壁区域。

应用范围：除强旋流过程无法精确预测外，其他流动都可以使用此模型，如模拟射流撞击、分离流、二次流和旋流等，能取得较标准模型更高的精度，但中等复杂流动受到各向同性涡旋黏度假设限制。

③可实现的k-ε模型（Realizable k-ε）

可实现的k-ε模型与标准k-ε模型相比，有两个的主要不同点：A.可实现的k-ε模型为湍流黏性增加了一个公式。B.为耗散率增加了新的传输方程，该方程来源于为层流速度波动而作的精确方程。

应用范围：除强旋流过程无法精确预测外，其他流动都可以使用此模型，包括有旋均匀剪切流，自由流（射流和混合层），腔道流动和边界层流动。

（3）k-ω模型

本书末列出k-ω模型表达式，具体见参考文献[10]。

①标准的k-ω模型

两个输运方程求解k与ω。对于有界壁面和低雷诺数的可压缩性和剪切流动，能取得较好的模拟效果，尤其是圆柱绕流问题、放射状喷射、混合流动等，包含转捩、自由剪切和压缩性选项。

②剪切压力传输k-ω（SST k-ω）模型

剪切压力传输k-ω模型合并了来源于ω方程中的交叉扩散，湍流黏度考虑到了湍流剪应力的传播，而且模型常量不同，这些特点使得该模型在近壁区域的自由流、逆压梯度流动、翼型、跨音速激波等方面有着更高的精度。

（4）雷诺应力模型（RSM）

雷诺应力模型同时考虑了连续性方程、动量方程、输运方程和各向异性湍流剪切力方程，计算收敛较难，计算耗时较大。主要应用于复杂3D流动（如弯曲管道、旋转、旋流燃烧、旋风等流动），尤其是强旋流运动，若要考虑雷诺应力各向异性时，必须应用该模型。

（5）大涡模拟（LES）

大涡模拟将湍流过程分为大尺度脉动和小尺度（小于惯性区尺度）脉动，把小尺度脉动认为是湍流耗散（作为额外的应力项（亚网格应力）），通过动量方程进行传递，将计算的重点放在大尺度脉动上，可以表现出高精度的湍流细节，但需要更为精细的网格，及较大的计算耗时。

（6）近壁面处理方法的选择

壁面对湍流流动影响较大，壁面的不光滑及其他因素均会对流动产生影响，在离壁面很近的地方，黏性力将抑制流体切线方向速度变化，同时壁面也阻碍正常的波动。RNG k-ε在近壁面的外部区域，湍流动能受平均流速的影响而增大，湍流运动加剧。而k-ε模型、RSM模型、LES模型仅适用于湍流核心区域（远离壁面处），故需要考虑如何使这些模型能够用于壁面边界层处的流动。Fluent软件的湍流模型提供了三种近壁处理方法：标准壁面函数法、非平衡壁面函数法和加强壁面函数法，具体选择可以参考下面的说明。

标准壁面函数法：应用广泛，计算精度与经济性合理，但属于高Re数的经验公式，不适用于低Re数区的流动，也没有考虑质量蒸发、▽p和体积力的影响，三维计算也较差。

非平衡壁面函数法：考虑了▽p和非平衡方程及分离、再附着、撞击现象。但在低Re区的流动、伴随质量蒸发的流动、▽p过大的流动、强体积力流动及三维流动模拟方面存在不足。

加强壁面函数法：不再依赖于经验壁面法则，较好地应用于复杂流动及低Re区的流动，但要求精细的网格、较大的CPU与内存。

1.2.3 初始条件和边界条件

流体动力学（CFD）计算中，初始条件和边界条件的正确设置是关键的一步。现有的CFD软件都提供了现成的各种类型的边界条件，这里对有关的初始条件和边界条件作简单介绍。

1．初始条件

初始条件是计算初始给定的参数，即t＝t0时给出各未知量的函数分布，如：

当流体运动定常时，无初始条件问题。

2．边界条件

边界条件是流体力学方程组在求解域的边界上，流体物理量应满足的条件。例如，流体被固壁所限，流体就不应有穿过固壁的速度分量；在水面这个边界上，大气压强认为是常数（一般在距离不大的范围内可如此）；在流体与外界无热传导的边界上，流体与边界之间无温差等。由于各种具体问题不同，边界条件一般要保持恰当：①保持在物理上是正确的；②要在数学上不多不少，刚好能用来确定积分微分方程中的积分常数，而不是矛盾的或有随意性。

通常流体边界分为流固分界面和流流（液液、液气）分界面。

1）流固分界面边界条件

飞机、船舶在空气或水中运动时的流固分界面，水在岸边及底部的流固分界面，均属这一类。一般而言，流体在固体边界上的速度依流体有无黏性而定，对于黏性流体，流体将粘附于固体表面（无滑移）：

式中，v|F是流体速度，v|S是固壁面相应点的速度。上式表明，在流固边界面上，流体在一点的速度等于固体在该点的速度。对于无黏性流体，流体可沿界面滑移，即有速度的切向分量，但不能离开界面，也就是流体的法向速度分量vn|F等于固体的法向速度分量vn|S：

另外，也可视所给条件，给出无温差条件：

式中，v|F是流体温度，v|S是固壁面相应点的温度。

2）液液分界面边界条件

如果密度不同的两种液体的分界面就属于这一类。一般而言，对分界面两侧的液体情况经常给出的条件是

对应力及传导热情况给出的条件是

3）液气分界面边界条件

液气边界面最典型的是水与大气的分界面，即自由面，由于自由面本身是运动和变形的，而且它的形状常常也是一个需要求解的未知函数。因此就有一个自由面的运动学条件问题。设自由面方程为

并假定在自由面上的流体质点始终保持在自由面上，则流体质点在自由面上一点的法向速度，应该等于自由面本身在这一点的法向速度。经过一系列推导（见文献[2]），得到自由液面运动学条件：

如果要考虑液气边界上的表面张力，则在界面两侧，两种介质的压强差与表面张力有如下关系：

这就是自由面上的动力学条件，当不考虑表面张力时，有

式中，pa为大气压强。

4）无限远的条件

流体力学中的很多问题，流体域是无限远的。例如，飞机在空中飞行时，流体是无界的。如果将坐标系取在运动物体上，这时无限远处的边界条件为

其中下标∞表示无穷远处的值。