计算流体力学从RANS到DNS

01 初识湍流

提起湍流，大家想到的是什么呢，是奔涌的河流？是袅袅炊烟？还是飞机划破天际时的两道白杠？仔细想想，湍流似乎真的无处不在，伴随着我们生活，又影响着我们的生活，那么关于湍流的科学研究是从哪里开始的呢？

相信很多小伙伴们第一次正式接受到湍流较为“学术”的描述，来自于流体力学课程中的雷诺实验——当圆管中流动的雷诺数逐渐增加大时，平稳的流动便逐渐演化为杂乱无章的流动，即为湍流。而教科书上2300这个临界雷诺数便成为我们遇见湍流的开始。

02 雷诺数的困惑

正如上面所述，雷诺数的大小通常被认为是衡量流态的标准。雷诺数的定义大家都不陌生：Re=ρVL/μ，其中ρ为密度，V为特征速度，L为特征长度，μ为动力粘度系数；其中ρ，V和μ都很好理解，特征长度L呢？书上说，对于平板边界层流动，特征长度为平板起点到某位置的长度；而对于机翼或者圆柱绕流，则选取弦长或者直径作为特征长度：

可以看到特征长度从某种意义上表达了固体几何对流动的扰动尺度。更大的特征长度意味着对流场有更多的扰动和能量注入，流动越趋向于湍流，而此时雷诺数也越大。因此我们可以得出结论，湍流属于高雷诺数的流动。

不过需要说明的是，对于某个特定的工况，尽管可以用临界雷诺数来近似判断流动是否趋向于湍流，但此临界雷诺数的值并不是放之四海而皆准的。比如对于我们熟知的圆管实验，自雷诺之后，有很多学者通过增加管壁光滑程度和减小前方来流扰动的方式，成功将临界雷诺数大大延后，甚至达到107这个量级。

03 湍流即是漩涡？

“流体不经搓，一搓就成涡”，这句广为流传的话被很多人拿来形容湍流产生的原因。直观的经验告诉我们，剪切会形成涡，比如用勺子搅拌牛奶。只是漩涡就代表着湍流吗？

上图为某风洞喷口的剪切层涡量分布。在剪切层最开始发展的阶段，流动呈现明显的规律性，涡的大小和强度都比较单一，在涡的内部仍可以认为是层流的状态。而随着流体向下游发展，掺混作用增强，大涡不断破碎为小涡，小涡则进一步破碎并逐渐消耗，化为流体的内能，这个阶段才能称之为完全发展的湍流。而在湍流的内部，也并非完全是杂乱无章的，而是存在一定的相似结构，称之为拟序结构。因此，复杂的、多尺度的、混沌变化的流体运动状态才是湍流的特征。

04 湍流中涡的尺度

湍流中涡的尺度通常被分为如下四个层级：宏观尺度（与边界相关的主流尺度），积分尺度（大涡的典型尺度，一般与宏观尺度同量级），泰勒尺度（无损的进行能量传递的惯性尺度），耗散尺度（最小涡的尺度，小于该尺度的涡将被耗散）。

对于常见的湍流流动来说，耗散尺度和宏观尺度可能相差很多个数量级。为了界定不同的尺度，学者们又定义了?EI和?DI两个参数作为泰勒尺度的上限和下限。对于高雷诺数的流动，?EI约为积分尺度的1/6，而?DI约为耗散尺度的60倍。

那么这些尺度和流场有什么相关性呢？

简而言之，大涡的尺度取决于边界条件（比如对于整车的尾流，其大涡的尺度和车高为同一量级），而小涡的尺度取决于湍流雷诺数。积分尺度和耗散尺度大致符合关系式：?0/η ~ ReL^(3/4)（湍流雷诺数ReL=ρVL/μ，其中的L和V分别为积分尺度和湍流脉动速度）。可见湍流雷诺数越高，耗散尺度越小（这也是LES和DNS在高雷诺数条件下需要更小的 格才能准确的捕捉湍流结构的原因）。

既然耗散尺度取决于湍流雷诺数，因此在给定特征长度（大涡尺度）和特征速度（湍流脉动速度）的条件下，流体的粘性几乎决定了耗散涡的大小，即湍动能在什么样的尺度上耗散。至于湍动能的耗散量，对于稳定的湍流来说，应等于大涡的湍动能生成量。

05 湍流问题的数学困局

尽管我们可以通过观察和实验来了解湍流的物理现象，可是唯有准确的数学表达才能够实现数值的计算。提到流体的数学描述，当然逃不开“虐我千百遍”的N-S方程。

为了省点心力，我们以不可压流体的N-S方程为例，其张量形式如下：

作为非线性方程界的扛把子，N-S方程只有在非常苛刻的条件下，比如极低雷诺数的定常流动等少数情况下，才可能得到稳定的解析解。而随着雷诺数的增大，湍流问题的复杂性使得N-S方程解的性质发生了变化，目前还看不到求出复杂湍流问题解析解的苗头。

不过很多工程问题，我们并不需要完全求解湍流，比如对于大部分的管路流动，工程上更关心管道内的压力损失和平均速度分布，而非湍流的细节。因此学者们转换思路——“对于未知的自然规律，不妨使用统计学的方法”——对N-S方程进行平均，把瞬时速度u分解为时均速度ū和脉动速度u’。这个点子是雷诺最先想出来的，因此平均后的N-S方程被称为雷诺平均的N-S方程，也就是RANS：

然而雷诺平均的N-S方程似乎更复杂了，除了平均速度的应力，上式中还多了脉动应力项，称之为雷诺应力，成为新的拦路虎。不过“他大舅他二舅都是他舅”，既然平均速度的应力可以描述为粘度乘以速度梯度的方式，那么雷诺应力应该也可以，于是1877年Boussinesq提出了影响深远的涡粘性假设：雷诺应力=μt*(?ū/?y)，其中μt表征了湍流脉动引起的切应力效应，称为涡粘性系数。至此，雷诺应力和平均流动之间通过这种举(生)一( 搬)反( 硬)三( 套)的方式建立起来了，为湍流的数值计算指明了新的方向。

06 计算流体力学之RANS时代

有了涡粘性假设，这似乎意味着我们可以通过数值方法求解复杂的湍流问题了。不过涡粘性系数主导了湍动能的耗散，如果无法准确给出，后果将非常严重（试想一个系统中，能量在不断生成，却无法耗散掉——只吃不拉是会出毛病的）。

然而对于绝大部分情况，μt是未知的，而且涡粘性在边界层附近变化很大。尽管Boussinesq指明了湍流数值计算的方向，然而道路上却充满了沼泽和泥泞，直到咱们广义的祖师爷——普朗特于1924年提出了混合长度理论，湍流的计算从数学表达到工程应用这座桥梁才逐渐变得清晰。

01零方程模型

我们知道流体的粘性来自于分子自由运动产生的掺混，与分子运动自由程密切相关；而对于涡粘性，也可以类似的定义湍流脉动掺混的长度，称之为混合长度l，其物理意义为流体微团耗散前所经历的距离，因此脉动速度可以表示为混合长度与法向速度梯度的乘积，而涡粘性系数则可以相应的表述为μt=ρl^2 *(?ū/?y)。如果知道了混合长度，便可以明确涡粘性系数，进而求解RANS方程。

然而混合长度的准确值也很难得知，于是普朗特继续发扬了“跟着感觉走，天下在我手”的科学精神，大胆的认为混合长度与到壁面的距离成正比，从而得到了CFD领域第一种实用的涡粘模型，即著名的“混合长度模型”。1978年， Baldwin和Lowmax基于湍流边界层内外层的流动差异提出了更合理的B-L模型，即针对湍流边界层的内层和外层分别定义混合长度。

混合长度模型是代数模型，相当于直接用代数公式定义了涡粘性系数，被称为零方程模型，即不引入额外的方程即可求解RANS方程。虽然B-L模型只适用于小曲率、无分离的流动，但其计算量很小，于是很多商用CFD软件仍然会在湍流模型库中保留B-L模型。

02两方程模型

而我们熟知的k-epsilon模型及其变种（如k-omega模型等），也属于涡粘性模型。该模型针对混合长度继续演化，将其表示为湍动能k、湍流耗散率epsilon和湍流脉动速度u’的函数，而涡粘性系数可由k和epsilon导出。

前苏联数学家Kolmogorov在1941年提出了湍流在惯性子区之下的各向同性假设（又称K41理论），即湍流脉动速度u’=v’=w’，而湍动能k则可简化为3/2*u’^2。在此基础上学者们分别建立了k和epsilon两个输运方程，完成了模型的封闭。既然有输运，因此基于两方程模型算出的涡粘性，包含了湍流运动的部分历史效应，是对混合长度模型的改进。

标准的k-epsilon模型由Launder 和Spalding于1974年提出，其k方程和epsilon方程的表达式分别如下：

两方程模型本质上是在湍流各向同性假设基础上，对混合长度理论的进一步发展，因此适用于更复杂的流动。不过熟悉k-epsilon模型的CFDer都知道，其补充的epsilon方程中包含了经验系数，这些系数一般由特定的实验得出，因此限制了k-epsilon模型的适用范围。当然，k-epsilon模型也有诸多进化，主要是对难以精确求解的湍流耗散率epsilon的解法进行改进。

既然有零方程模型和两方程模型了，那中间的一方程呢？一方程模型当然是有的，只是它并非是零和二之间的过渡。

03一方程模型

波音公司的大神Spalart基于翼型计算的丰富经验，于1994年和Allmaras一起提出了著名的S-A模型。该模型不再使用湍动能和湍流耗散率计算涡粘性系数，而是直接导出涡粘性系数的输运方程：

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

07 LES逐步走向实用化

严格来讲，RANS其实是为了工程化应用而拼凑出来的数学模型，放弃了对所有非定常湍流信息的模拟，转而寻求平均意义下的流动结果。其涡粘性假设也意味着简化处理由湍流脉动引起的雷诺应力。因此无论 格画的多细，RANS也只能求解出平均流场的流动结构。

为了更好的解析小尺度的脉动对流动的影响，学者们又发展了大涡模拟方法（LES）。其核心思想是对N-S方程进行某种过滤，把小于过滤尺度的脉动使用模型表达，而大于过滤尺度的涡则直接求解，从而分辨出更多的流动细节。

LES的理念最早由Smagorinsky于1963年提出，但由于硬件资源消耗过大，直到近年来才得到一些工程应用。LES对小尺度脉动的模型化处理借鉴了Boussinesq的涡粘性假设——使用小于 格尺度的应力，即亚格子应力来代表小尺度的耗散。

08 终极利器DNS

有人说，统计方法是人类面对未知的无奈之举——如果我们有无限的计算资源，当然就不用考虑RANS和LES了。我们可以在极小的空间和时间尺度上求解流动从而获得湍流全部的流场细节。这种对所有尺度的湍流都进行直接求解的办法称之为直接数值模拟（DNS），最小的 格尺度要小于耗散尺度。

不同于RANS和LES在时空分辨率上的缩水，DNS直接捕捉了湍流的全部细节。不过，大家可以使用前述的耗散尺度和积分尺度之间的关系式，大致估算一下最小的 格尺度，再估算一下计算量，估计就认识到什么是天文数字了。

目前DNS只能应用于简单几何的流动机理研究，比如上图的准二维叶型的边界层流动。而Spalart更是预测直到2080年，DNS才会进入工程实用领域。若没有诸如量子计算机等Bug级突破，估计小编有生之年也不指望能用上DNS。当然，随着LBM方法在CFD领域的大放异彩，由于其极高的并行效率和天然瞬态的算法，DNS的进程很有可能会大大的加快。

09 DNS是否是湍流计算的终结？

尽管针对具体的流动问题，DNS能够得出N-S方程的特解，并得到所有尺度的流动细节。然而就湍流计算而言， DNS仍然任重而道远，因为数值求解不仅依赖于算法，还依赖于准确的初场和边界条件。而湍流问题的高度非线性，又导致其结果对初始条件非常敏感。而真实湍流的初始条件是很难准确给出的，比如想要完美预测下一时刻的天气变化，就要知道当前时刻详细的大气条件——这几乎是不可能的。

即便是使用DNS得出了湍流的精确流场，也只是针对某一特例的数值计算结果。就湍流理论的研究而言，DNS只是一个新的起点。

近百年前，湍流理论的创始人之一Richardson曾用一首诗来描述湍流：

Big whirls have little whirls

That feed on their velocity

And little whirls have lesser whirls

And so on to viscosity

谨借此诗献给每一位“湍流虐我千百遍，我待湍流如初恋”的CFDer！