计算流体力学的“百家争鸣”

来源：LBM与流体力学（ID：LBMCFD），作者：卢比与钢蛋。

在春秋战国时期，华夏大地上涌现出了很多著名的学派并流传至今，而后人也把这个思想文化迸发、群星闪耀的局面称为“百家争鸣”。回到我们最熟悉的CFD领域，作为RANS方法能够生存的最重要基石，无数的学者提出了各种各样的湍流模型，可谓百花齐放。可是迄今为止，仍然没有一款湍流模型能够一统天下，于是，在各种CFD软件中下拉湍流模型选项的时候也总是能看到很多种选择，这让有选择困难症的小伙伴们非常无奈。今天，我们就和大家一起聊一聊湍流模型的百家争鸣。

幸好对于大部分工程问题，我们并不需要完全求解湍流。比如工程上更关心流动的压力损失和平均速度分布，而非湍流的细节。实验做完五年以后，雷诺幡然醒悟：既然流动未可知，不妨使用统计学的思想——对N-S方程进行平均，把瞬时速度u分解为时均速度ū 和脉动速度u’，代入N-S方程即可得到雷诺平均的N-S方程，也就是RANS。

然而雷诺平均的N-S方程似乎更复杂了，除了平均速度的应力，上式中还多了脉动应力项，称之为雷诺应力，成为新的拦路虎。不过冥冥之中自有天意，在雷诺提出对N-S方程进行平均的十年前， Boussinesq 已经将湍流脉动引起的切应力类比成了牛顿内摩擦定律，即用粘度乘以速度梯度来表示湍流脉动引起的切应力，这就是大名鼎鼎的涡粘性假设：雷诺应力=μt*(əū/əy)，其中，μt 为涡粘性系数。至此，通往计算流体力学的大门终于打开了一条细缝。

流体的粘性来自于分子自由运动产生的掺混，与分子自由程密切相关；对于涡粘性，也可以定义类似的参数。于是，普朗特定义混合长度l 为湍流脉动掺混的长度，其物理意义为流体微团耗散前所经历的距离。

于是，脉动速度可以表示为混合长度与壁面法向速度梯度的乘积，而涡粘性系数则可以相应的表述为μt =ρl ² *(əū/əy)。只要知道了混合长度，便可以明确涡粘性系数，进而求解RANS方程。

然而，混合长度的准确值也很难得知，于是普朗特继续发扬“跟着感觉走，天下在我手”的科学精神，大胆的认为混合长度与到壁面的距离成正比，从而得到了CFD领域第一个实用的涡粘模型。尽管普朗特的理论有东拼西凑的嫌疑，但是它使得湍流模型和计算流体力学实现了真正意义上的破局，并引领了其后近百年时光的CFD大发展。

1978年， NASA Ames 研究中心的Baldwin 和Lowmax 在AIAA 第16届航空航天科学会议上发表经典的论文《Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows》，提出了一种更合理的混合长度模型。

Baldwin 和Lowmax 基于湍流边界层内外层的流动差异，将湍流粘性系数也相应的定义成了两个公式。内层延续普朗特的理论，而外层使用了更加符合实际的Clauser 公式。当然在实际使用中，两位学者也做了一些更新。后来很多学者使用B-L 模型计算了翼型的绕流、平板流动、台阶流，甚至是超声速的流动，仿真结果和实验都呈现比较好的一致性。

作为一种典型的代数模型，B-L 模型不需要引入额外的方程即可求解RANS方程，其计算量较小，因此很多商用CFD软件仍然会在湍流模型库中保留B-L 模型。不过需要特别说明的是，B-L 模型只适用于小曲率、无分离的流动，不适合分离流、喷流或逆压梯度较大的流动。

前人栽树，后人乘凉。于是一种携带着湍流演化思想的湍流模型便应运而生了，那便是两位CFD领域的顶级大神B.E. Launder 和D.B. Spalding 以及他们的博士生们在20世纪70年代初陆续提出来的k-epsilon 模型。

k-epsilon 模型也属于涡粘性模型，不过它变换了另外一种更加湍流的思想来求解混合长度，进而求解涡粘性系数：Kolmogorov 的各向同性假设告诉我们湍流脉动速度u'=v'=w'，因此湍动能k 则可简化为3/2*u'²，而另外一个重要的参数epsilon 表示流体微团的耗散率，即湍动能/耗散时间。

而定义为流体微团耗散前所经过距离的混合长度l，则可以顺势定义为耗散时间和流体微团速度的乘积。当然流体微团的速度和脉动速度是相当的，通过经验系数来表示，于是我们可以用湍动能和湍流耗散率来表示湍流粘性系数。

面对多出来的两个未知量k 和epsilon，则分别建立其输运方程进行封闭，这便是我们熟悉的k 方程和epsilon 方程。不过需要特别指出的是，补充的epsilon 方程中包含了经验系数，这些系数一般由实验得出，因此限制了k-epsilon 模型的适用范围。实际工程应用中，也有很多工程师会根据特定的应用场景修改相应的经验参数来匹配实验的结果。

k-epsilon 模型是对混合长度理论的进一步发展，因此相对于零方程模型适用于更复杂的流动。不过通过上述的推导，大家也可以看出来k-epsilon 模型仍然有拼凑的嫌疑。它建立于湍流各向同性假设基础之上，因此适合于完全发展的湍流，但不适用于逆压梯度较高的流动，以及喷流和二次流。因此，还有学者提出了Realizable 以及RNG k-epsilon 模型，在壁面边界层的处理以及分离流动的求解有了明显的增强。

该模型定义了比耗散率omega，并通过k 和omega 来描述湍流粘性系数。从方程结构来看，omega 方程和epsilon 方程非常相似，但是两者之间有一个重要的区别——epsilon 方程中的经验系数来源于粘性底层中的阻尼函数，而omega 方程中的经验系数却不需要。这个阻尼函数的精度在强逆压梯度的流动中存在一定的精度问题，因此k-epsilon 模型不适合于处理逆压梯度较大的流动，而k-omega 模型却能够很好的驾驭。

此处需要特别指出的是，如果把omega的表达式直接代入epsilon方程，直接推导出来的方程比omega方程多了一个交叉扩散项，而这一项正是两者之间的差异来源。

k-omega 模型的确能够很好的改善强逆压梯度附近的流动求解，但是它也存在一个比较明显的缺陷：当自由来流中湍动能出现微小的变化时，湍流粘度系数会出现剧烈的变化，导致流动的计算出现明显的不合理；同时表面摩擦系数也会出现一定的变化，影响分离位置的判断。反观k-epsilon 湍流模型对于自由来流的湍动能变化，则没有那么敏感。

既然k-omega 和k-epsilon 各有所长，那就把它们结合起来吧，于是一个更加适用的湍流模型就诞生了。1992年，NASA Ames 研究中心的Menter 博士提出了著名的k-omega SST 湍流模型，在omega 方程后面添加了一个交叉扩散项，这个交叉扩散项刚好就是k-omega 和k-epsilon 模型之间的差异，而Menter 博士则聪明的在这个交叉扩散项上面乘了一个混合函数，通过混合函数则可以轻易的控制湍流模型的变身：当混合函数为1时，omega 方程还是omega 方程，用于边界层附近的求解，而当混合函数为0时，omega 方程则变身成了epsilon 方程，用于主流的求解，当然两者之间通过一个混合区域来进行更加光顺的过渡。

凭借着更好的适用性和计算精度，k-omega SST 模型成为了越来越多的CFD 工程师的选择，尤其在航空航天和叶轮机械领域的出镜率极高。

一方程模型当然是有的，只是它并非是零和二之间的过渡，而是另外一种更加特立独行的湍流模型。波音公司的大神Spalart 基于翼型计算的丰富经验，于1994年和Allmaras 一起提出了著名的S-A 模型。该模型不再使用湍动能和湍流耗散率计算涡粘性系数，而是直接导出涡粘性系数的输运方程。

这种模式更适合于平均流场中有剧烈变化的湍流，比如几何曲率明显变化、存在激波等工况，因此广泛应用于航空航天和叶轮机械领域。不过标准的S-A 模型并没有针对一般工业流域的流动进行校准，尤其是某些自由剪切流动比如平面射流，可能会产生较大的误差。因此，它是一种更专用的湍流模型。

不过笔者倒是更愿意相信，随着硬件和CFD方法的提升，LES 甚至是DNS 会慢慢融入并替代现在的RANS。到那时，或许我们会时常怀念起曾经和这些湍流模型一起挣扎的岁月，而这些湍流模型也会因为我们的怀念而变成胸口的朱砂痣。