系统总结1.7.2湍流模型的优缺点及适用范围

目前，湍流的数值模拟方法可分为^[31]：①直接数值模拟（Direct Numerical Simulation，DNS）；②非直接数值模拟。

直接数值模拟是指直接用瞬时的N-S方程对湍流进行计算，无需任何简化或近似，理论上可得到比较准确的计算结果，但由于要模拟最小尺度涡，计算网格的分辨率要足够高，而计算区域的尺寸应足以容纳最大尺度涡，整个计算域的网格节点数的量级非常大，基于目前计算机技术的发展水平，对于水泵这样复杂的三维高雷诺数湍流的数值计算还不现实。

非直接数值模拟包括：①大涡模拟（Large Eddy Simulation，LES）；②雷诺平均法（Reynolds Averaged Navier-Stokes，RANS）；③统计平均法。大涡模拟是近年来研究比较热门的一种模拟方法，其基本思路^[31]就是通过过滤函数，将涡分为大涡和小涡，大涡通过直接数值模拟进行计算，而小涡对大涡的影响通过近似的模型来考虑。大涡模拟对计算机的要求也比较高，但低于DNS方法。近些年，大涡模拟得到了很快速的发展，Fluent等商业软件也提供了大涡模拟的方法，但由于其对计算机性能的要求较高，因此应用于实际工程中还不是很成熟。统计平均法^[31]是基于湍流相关函数的统计理论，采用相关函数及谱分析的方法来研究湍流结构；统计理论主要涉及小尺度涡的运动，这种方法在工程上应用较少。

现代工程数值计算中更多的是通过求解时均化的N-S方程（RANS方法）来进行计算的。由于时均化的N-S方程组不封闭，需要引入额外的方程来封闭该方程组，这就提出了湍流模型的概念。所谓的湍流模型^[36-39]就是依据湍流的理论知识、实验数据以及经验，对雷诺平均运动方程中的雷诺应力项建立表达式或方程，并对雷诺应力方程中的一些项进行合理的模化，从而使湍流的平均雷诺方程组封闭的理论。根据对雷诺应力做出的假定或处理方式的不同，可将湍流模型分为雷诺应力模型和涡黏模型。涡黏模型是在Boussinesq涡黏假设的基础上提出的，从给出的方程个数来看，可将涡黏模型分为零方程模型（代数模型）、一方程模型和两方程模型。下面从湍流模型的优、缺点及适用范围等方面进行系统地归纳和总结。

1.Boussinesq涡黏假设

鉴于分子黏性系数与切应力之间的关系，Boussinesq于1877年提出了引用一个涡黏性系数μ_t，使湍流的雷诺应力与湍流场中的平均速度梯度之间建立下述关系：

涡黏性系数μ_t与分子动力黏性系数有相同的量纲，但与分子黏度不同的是它还依赖于流动情况，这是涡黏性假设的弱点^[32]。

2.零方程模型

零方程模型又称为代数方程模型，它通过直接建立雷诺应力与平均速度之间的代数关系，不涉及高阶统计量微分方程。这类模型比较简单、直观，无须附加湍流特性的微分方程；但它属于当地平衡型模型，不能反映上游历史的影响，且没有考虑对流和扩散的影响，在一些比较简单的流动，如射流、边界层、管流及喷管流动等^[40-42]有较为理想的计算结果，对于表面曲率或压力梯度很大的流动情况，其计算结果不太理想。常见的零方程模型有Prandtl混合长度理论，基于Prandtal混合长度理论发展出了C-S模型、B-L模型及J-K模型等。

（1）Prandtl混合长度理论

Prandtl于1925年提出了混合长度理论^[43，44]。Prandtl假设湍动黏性系数ν_t正比于时均速度梯度和混合长度l_m，且ν_t等于时均速度梯度和混合长度l_m的乘积。

Prandtl混合长度理论在预测一些简单流动时取得了较为成功的结果，如管流、槽流、边界层和自由剪切流动等，但它不适用于计算有回流的流动。

（2）C-S模型

C-S模型^[45]是由Cebeci和Smith于1967年提出的双层代数湍流模型。该湍流模型的涡黏性系数在内、外层分别采用不同的表达式进行表示。内层涡黏性系数的计算表达式可表示为：

外层涡黏性系数的计算表达式为：

式中 l_D——衰减长度系数；

δ*——边界层位移厚度；

u——外部自由流速度；

γ——间隙因子。

（3）B-L模型

B-L模型^[46]是由Baldwin和Lomax两位学者于1978年提出的平衡型双层代数湍流模型，它的提出是为了解决边界层属性如δ，δ_ν*，U_e等很难被确定的计算。该模型具有形式简单，计算量小，不需要考虑边界层厚度，健壮易收敛，应用方便等优点^[47]；适用于模拟附体流动^[48]，对于较小的局部分离流动，该模型也有一定的模拟能力，尤其对亚临界流动的计算比较有效。但在比较明显的分离流、激波附面层干扰、多物面等流动情况中的计算结果较差，且在多块网格流场计算中，计算物面距离比较困难；另外，当采用B-L模型计算跨声速分离流动时，其激波总是较实际位置靠后^[49，50]。

B-L模型的计算形式类似于C-S模型，两者之间的主要区别在于外层涡黏性系数的表达式不同。B-L模型的内层涡黏性系数μ_t，i可表示为：

外层涡黏性系数可表示为：

式中 y——距壁面的法向距离；

Ω_V——速度旋度（涡量）的模；

w——物面值；

F_max——函数F（y）的最大值；

y_max—F（y）——取得最大值所对应的物面法向距离；

q——当地速度的模；

Г——间歇因子。

（4）J-K模型

1985年，Johnson D A和King L S提出了非平衡型代数湍流模型J-K85模型[51]。该模型采用涡黏性假设，将涡黏性分布同最大雷诺切应力的分布联系起来。J-K模型虽然仍采用涡黏假设，却具有雷诺应力模型的特点，求解常微分方程比一方程、两方程模型中的偏微分方程要简单且节省计算时间，其计算量稍多于B-L模型。J-K85模型能较好地计算二维翼型分离流动，但该模型在模拟附体平衡流动时激波位置过于靠前，且对边界层厚度的确定非常困难。

Abid R等几位研究者于1990年提出了J-K90A模型[52，53]。该模型对分离流动有较好的模拟能力。Johnson D A和Coakley T J于1990年提出了J-K90J模型来改进J-K85模型模拟附体平衡流动激波位置靠前的缺点，该模型考虑到压缩性影响，采用双曲正切函数形式的混合内层湍流黏性系数。1992年，Johnson D A在J-K90J模型的基础上提出了J-K92模型[54]，该模型采用相似的内层混合黏性系数，将平衡流内层黏性系数修改为B-L模型内层黏性系数，同时对附面层厚度进行重新定义，同样采用双曲正切函数对流场湍流黏性系数进行合成。J-K92模型具有更好的模拟大攻角、大分离流动的能力，它能较准确地模拟背风面涡的结构与强度，同时能更好地模拟壁面压力分布。由于J-K模型考虑了最大黏性切应力非平衡性质的影响，因此，它能有效地模拟流动的上游历史效应，与此同时，J-K模型在模拟强激波诱导分离、大攻角分离流动时要优于B-L模型^[55]。

3.一方程模型

在零方程模型中，湍动黏度和混合长度l_m都把雷诺应力和当地平均速度梯度相联系，是一种局部平衡的概念，忽略了对流和扩散的影响。一方程湍流模型突破了混合长度假定的局限性。比较常见的一方程湍流模型有k方程模型、S-A模型及B-B模型等。

（1）k方程模型

只保留k方程为微分方程，而其他二阶脉动相关量均由代数方程来表示。这种模型对复杂湍流的揭示能力有所降低。用普朗特混合长度来确定湍动能耗散率ε，进而把ε方程简化成代数方程，湍流长度l、湍动能k和耗散率ε三者之间的关系式为：

涡黏性系数ν_t和湍动能k的输运方程为：

输运方程中的经验常数为：C_μ=0.09，C_D=0.09，σ_k=1。

（2）S-A模型

S-A模型^[56，57]是由Spalart和Allmaras于1992年提出的涡黏性一方程湍流模型。该模型的计算量较小，对物面处的网格要求不是很高^[58]，具有良好的鲁棒性、稳定性和收敛性，对一定范围内的分离流动模拟能力要好于B-L代数方程模型；但对于较大的分离流动，S-A模型显得比较保守，对流动分离会有抑制作用，计算精度也会有所降低。最初的S-A模型应用于如飞机的翼剖面和机翼流动等^[59，60]外流的数值计算中，近年来不少学者通过对其进行改进而将其应用到叶轮机械的内流数值计算中^[61]。S-A模型在壁面限制流、附着流、薄层自由剪切流以及受到逆压梯度作用的边界层流动中得到了较多的验证，取得了比较理想的计算结果。

S-A模型忽略了包含湍动能k的雷诺应力项，而是通过对变量（与涡黏性有关的变量）构造输运方程^[62]来获得涡黏性系数ν_t。S-A模型的涡黏性系数ν_t和变量的输运方程可表示为：

变量ν～的输运方程中的右端分别是梯度耗散项、生成项和湍动能的壁面破坏项。

式中，d为到物面的距离，输运方程中的经验常数为：C_b1=0.1355，C_b2=0.622，σ=2/3，κ=0.41，C_ω1=3.239，C_ω2=0.3，C_ω3=2.0，C_ν1=7.1。

（3）B-B湍流模型

B-B模型^[63-65]是由Barrett S.Baldwin和Timothy J.Barth于1990年提出的一方程湍流模型。其输运方程可表示为：

D₁和D₂为阻尼函数。输运方程中的经验常数为：C_ε1=1.2，C_ε2=2.0，C_μ=0.09，A0⁺=26，A2⁺=10，κ=0.41，σ_ε2=0.7。

4.两方程模型

较为完善的湍流模型是求解两个变量（速度尺度和长度尺度）的两方程湍流模型。随着计算流体动力学的快速发展，出现了各种不同的两方程湍流模型。下面对几种常见的两方程湍流模型进行介绍。

（1）标准k-ε模型

1972年，Launder和Spalding提出了标准k-ε模型^[66]，该模型是在涡黏理论基础上发展起来的，是从实验现象中总结出来的半经验公式模型。标准k-ε模型的涡黏系数包含了流动的部分历史效应，具有简单、稳定、经济、合理的计算精度、适用范围广等优点，自从被提出之后，它就成为工程流场计算的主要工具之一，且得到了广泛的检验，能够满足较大范围内的工程精度要求。但标准k-ε模型属于高雷诺数湍流模型，只适合于充分发展的完全湍流流动的数值计算。它也具有比较明显的缺点^[67，68]：

1）标准k-ε模型假定雷诺应力和当时当地的平均切变率成正比，所以不能反映雷诺应力沿流向的松弛效应；

2）标准k-ε模型具有明显的各向同性特性，不能反映雷诺应力的各向异性特征；

3）平均涡量会影响雷诺应力的分布，这种影响在分离流动中十分明显，而标准k-ε模型不能反映平均涡量对雷诺应力分布的影响。

这些缺点使它在计算复杂剪切湍流场时有明显的不足，且不能很好地预测有回流、大曲率流线及强旋流等流动情况的湍流特性。标准k-ε模型对湍动能k和耗散率ε构造输运方程。湍动能k的输运方程是通过精确的方程推导得到的；耗散率ε的输运方程则是通过物理推理，数学上模拟相似原形方程得到的。湍动能k和耗散率ε的定义公式可表示为：(www.xing528.com)

标准k-ε模型的湍流涡黏性系数ν_t、湍动能k和耗散率ε的输运方程可表示为：

G_k是由于平均速度梯度引起的湍动能k的生成项，σ_k和σ_ε分别为湍动能k和耗散率ε的湍流普朗特数，输运方程中的常数是通过对典型流动的实验结果和算例结果进行最佳拟合而得到的：C_ε1=1.44，C_ε2=1.92，C_μ=0.09，σ_k=1.0，σ_ε=1.3。

（2）RNG k-ε模型

由于标准k-ε模型具有明显的各向同性特征，在计算强旋流或带有弯曲壁面的流动时会出现失真等缺点，Yakhot和Orzag于1986年提出了RNG k-ε模型[69，70]，随后，Yakhot和Smith于1992年修正了RNG k-ε模型。RNG k-ε模型是对瞬时的N-S方程使用重整化群的数学方法推导出来的模型，在形式上与标准k-ε模型非常相似，但在以下方面进行了改进[71]：

1）标准k-ε模型中的系数是由试验结果和典型算例拟合得到的经验常数，而RNG k-ε模型中的系数则由重整化群理论计算所得；

2）RNG k-ε模型对耗散率ε的输运方程进行了修正，增加了R项来考虑高湍流应变率的影响，即在一定程度上考虑了湍流的各向异性效应和湍流的漩涡特性，从而改善了对较为复杂的湍流如旋转流和大曲率流等流动的预报精度。

这些特点使得RNG k-ε模型比标准k-ε模型在更广泛的流动中具有更高的可信度和精度。RNG k-ε模型的缺点是该模型是为高张力引起的湍流黏度降低而设计的。RNG k-ε模型的涡黏性系数ν_t、湍动能k和耗散率ε的输运方程可表示为：

输运方程中的常数为：C_ε1=1.42，C_ε2=1.68，C_μ=0.0845，σ_k=0.7179，σ_ε=0.7179，η₀=4.38，β=0.012。

（3）Realizable k-ε模型

标准k-ε模型在对湍流耗散率ε的输运方程模化时过于粗糙，在某些情况下不能给出合理的湍流尺度，所以标准k-ε模型对于具有较高主流剪切率、较大曲率壁面或有流动分离情况的流场计算往往得不到理想的结果。Shih和Liou等针对这些缺陷进行修正，于1995年提出了Realizable k-ε模型^[72]。该模型是一种新型的高雷诺数湍流模型，它在如下方面进行了改进：

1）Realizable k-ε模型的湍流黏度表达式中的系数C_μ表示为湍流时间尺度与应变张量和旋转张量的函数，从而引入了与旋转和曲率相关的内容；

2）Realizable k-ε模型对耗散率ε的输运方程采用了新的模化方法，使雷诺主应力和切应力满足可实现性条件；

3）对湍动能k和湍流耗散率ε输运方程中的模型系数进行了修改。

Realizablek-ε模型在定义湍流黏度时考虑了平均旋度的影响，这样使得在计算旋转和静态流动区域时不能提供自然的湍流黏度。Realizable k-ε模型的这些特性使得该模型在对平板射流和圆柱射流的发散比率预测上得到了更为精确的结果，且对于旋转流动、强逆压梯度的边界层流动、流动分离和复杂的二次流等有很好的表现。Realizable k-ε模型已被有效地应用于旋转均匀剪切流、平面混合流、平面射流、圆形射流、管道内充分发展流动、边界层流动以及带分离的流动等^[73]。Realizable k-ε模型的涡黏性系数ν_t、湍动能k和湍流耗散率ε的输运方程可表示为：

σ_k，σ_ε分别是湍动能k和湍流耗散率ε的湍流普朗特数，Ω_ij为旋转坐标系中的时均转动速率张量。输运方程中的常数为：A₀=4.0，C₂=1.9，σ_k=1.0，σ_ε=1.2。

（4）低雷诺数k-ε模型

上述几种k-ε模型都是适用于高雷诺数流动计算的两方程模型，对于近壁面低雷诺数流动区域的计算需要加入壁面函数，而壁面函数的表达式主要是根据比较简单的平行流动边界层的实测资料归纳出来的。Jones和Launder^[74，75]认为，低雷诺数的流动主要体现在黏性底层中，流体的分子黏性起着绝对支配地位，故通过对高雷诺数k-ε模型进行修改，使其可应用于各种雷诺数流动区域的计算：

1）通过控制方程中的扩散系数项同时包括湍流扩散系数和分子扩散系数两部分来体现分子黏性的影响；

2）在系数计算公式中引入湍流雷诺数Re_t=ρk²/（ηε），使控制方程的有关系数考虑不同流态的影响；

3）在k方程中应考虑壁面附近湍动能耗散的各向异性特性。

低雷诺数k-ε模型的涡黏性系数ν_t、湍动能k和耗散率ε的输运方程可表示为：

当Re_t很大时，f_μ、f₁和f₂均趋近于1；当局部湍流的Re_t数小于150时，使用低雷诺数k-ε模型，此时，充分发展的湍流核心区及黏性底层采用同一套公式进行计算，但由于黏性底层的速度梯度较大，所以黏性底层的网格要求更密^[31]。

（5）k-ω模型

k-ε模型中的湍流耗散率ε缺乏自然边界条件，需要引入额外的边界条件，这些边界条件要么渐进不一致，要么数值计算困难，且k-ε模型对逆压梯度缺乏敏感，常常过多地模拟切应力，推迟分离。为克服k-ε模型的这些固有缺陷，Wilcox提出了形式更为简单的k-ω模型[76，77]，该模型是积分到壁面的两方程线性涡黏性模型，它不需要确定物面法向距离，提高了对逆压梯度的灵敏度，降低了近壁区的计算难度，且不包含k-ε模型中需要的复杂非线性阻尼函数，因此更加精确和稳定。

k-ω模型不会阻止近壁区的耗散减弱，在处理靠近物面区域的湍动能时，可能会得到渐进不一致的解；对于自由剪切流动，k-ω模型的计算结果过分地依赖于自由流的ω值[78]。

k-ω模型在尾流、混合流动、平板绕流、圆柱绕流、放射状喷射流、墙壁束缚流等流动情况中得到了较为满意的检验结果^[79]。

k-ω模型采用湍动能k和比耗散率ω的输运方程来确定涡黏性系数，其涡黏性系数ν_t、湍动能k和比耗散率ω的输运方程可表示为：

输运方程中的常数为：γ^*=1，γ=5/9，β_k=0.09，β_ω=0.075，σ_k=1/2，σ_ω=1/2。

（6）BSL k-ω模型

为了保持Wilcox k-ω模型在近壁面区域的健壮性和计算的精确性，且在边界层外边缘和自由流区充分利用标准k-ε模型的优势^[56]，提出了BSL k-ω模型^[80]。BSL k-ω模型的基本思想就是：在靠近壁面的区域采用比较精确的k-ω模型，而在其余边界层和自由流区域采用标准k-ε模型，为了得到统一的计算公式，需要将标准k-ε模型转换成k-ω模型的形式，然后通过加权函数来处理两套公式的系数。整合后的表达式与k-ω模型的主要差别就在于比耗散率ω的输运方程中增加了一项附加的交叉扩散项，且输运方程的模型常数发生了改变，变换后的标准k-ε模型具有类似于Wilcox k-ω模型的性能。BSL k-ω模型的输运方程可表示为：

在外边界层和自由流区域，Wilcox k-ω模型对自由流的ω值有着较强的敏感性，而标准k-ε模型则可避免这一不足，所以在加权函数F₁的构造中，将加权函数F₁设计成在子层或边界层的对数区域等于1，在尾流区域逐渐变换为零，即意味着在靠近壁面的区域激活k-ω模型，而在远离壁面或尾流的外部区域和自由剪切流区域中激活标准k-ε模型。

（7）SST k-ω模型

SST k-ω模型^[81-83]是在BSL k-ω模型的基础上改进了涡黏性的定义，以考虑湍流主切应力输运的影响，这样，SST k-ω模型使得对逆压梯度流动的预测（如分离流）得到了重要的改进^[84]。SST k-ω模型的湍动能k和比耗散率ω的输运方程表达式为：

加权函数F₁的表达式为：

同时，由于原始k-ω模型没有考虑湍流切应力的输运。Menter认为这会导致对于涡黏性的过分估计。因此他提出应该使用以下公式对涡黏性进行限制：

F₂是类似F₁的混合函数，用来修正F₁在自由剪切流中的误差。

（8）q-ω模型

Coakley于1983年提出了两方程线性涡黏性模型q-ω模型^[85]，该模型采用q和ω的输运方程来确定湍流尺度，与传统的两方程线性涡黏性模型相比，q-ω模型采用了新的近壁项以及更多地考虑了应力松弛影响。采用新的近壁项增强了计算的起始阶段和尖前缘、后缘处的计算稳定性；增加对应力松弛影响的考虑能够更准确地预测激波/边界层相互作用等复杂流动。

q-ω模型的涡黏性系数ν_t、q和比耗散率ω的输运方程可表示为：

式中 y——到最近壁面的距离。

输运方程中的常数为：C_μ=0.09，C₁=0.405D+0.045，C₂=0.92，Pr_q=1.0，Prω=1.3，α=0.0065。

5.雷诺应力模型

（1）雷诺应力模型

k-ε模型是建立在Boussinesq的涡黏性假设基础上的，该模型无法描述湍流能量转变为平均运动能量的传递过程和涡黏性的各向异性张量特性。为此，通过直接建立雷诺应力输运方程来封闭湍流流动方程组，发展出了雷诺应力模型（Reynolds Stress Equation Model，RSM）。RSM模型是指直接求解雷诺应力控制微分方程的二阶湍流模型，与k-ε模型相比，RSM模型计算得到的湍流特性更为合理一些，普遍性和预测能力也均优于其他二阶模型；但它过于复杂，需要求解的微分方程数目较多，计算机耗时较长。经模化后的雷诺应力方程可表示为：

湍动能k及耗散率ε的输运方程为：

式中 D——对应下标输运量的扩散项；

P——对应下标的源项；

PS——雷诺应力方程中的压力-应变项。

输运方程中的经验常数为：C_k=0.09～0.11，C₁=1.5～2.2，C₂=0.4～0.5，C_ε=0.07～0.09，C_ε1=1.41～1.45，C_ε2=1.9～1.92。

（2）代数应力模型

为了便于工程应用，在保存RSM模型湍流各向异性特点的基础上将其雷诺应力的微分方程简化为代数表达式，提出了代数应力模型（Algebraic Stress Equation Model，ASM）^[86]，该模型的计算量比雷诺应力模型要少，但因为比k-ε模型多解6个代数方程组，其计算量还是要远大于k-ε模型。代数应力模型在带有扭曲和出现二次流的方形管道、三角形管道以及带有强旋的管道流动中得到了较为满意的计算结果。

不考虑浮力作用、系统无旋转、忽略固体壁面的反射影响时，代数应力方程可表示为：

代数应力模型中的湍动能k和耗散率ε的输运方程分别为：

输运方程中的常数为：C_μ=0.09，C_ε1=1.44，C_ε2=1.92，σ_k=0.82，σ_ε=1.0。