0 引言

高雷诺数壁湍流是国防和航空航天领域关键的问题之一。对高雷诺数壁湍流的精确模拟和控制是指导飞行器航天器的设计和改进的基础，是国家安全、经济建设和科学技术中的重点研究对象之一，也是制约我国航空航天事业发展的基础科学难题之一。近些年，随着计算机技术突飞猛进的发展以及数值计算方法的进一步完善，数值模拟方法正逐渐成为湍流研究的重要手段之一。

经典的湍流数值模拟方法有三种：直接数值模拟方法(Direct Numerical Simulation, DNS)，大涡模拟方法(Large-eddy Simulation，LES)和雷诺平均方法(Raynolds-Averaged Navier-Stokes，RANS)。由于网格需求大，格式精度要求高，现阶段DNS方法只停留在简单流动的物理机理研究上^[1]。现在工程上最常用的方法还是更易于实现的RANS^[2-5]。但是雷诺应力模型普适性较差，特别是在大分离、非定常流动时效果不好。LES方法的计算量介于RANS和DNS之间，它可以用来模拟大分离、非定常流动，一般模型参数也相对较少^[6]。但是对于高雷诺数壁湍流而言，在壁面解析的情况下，LES的网格需求仍然很大^[7-8]，因此在工程应用中，壁面解析的LES使用得并不多。

为了更好地模拟大分离、非定常流动，人们提出了各种改进的方法。例如，壁面模化的LES方法(Wall-modelled LES，WMLES)^[9-10]，分离涡模拟方法(Detached Eddy Simulation, DES)等RANS-LES混合方法(Hybrid RANS-LES Method)^[11-14]。RANS-LES混合方法结合了RANS方法和LES方法的优点，同时又避免了这两种方法各自的缺点，对分离流动的模拟效果大大优于RANS方法。RANS-LES混合方法被认为是未来数十年湍流模拟最有希望的方法。虽然以DES为代表的RANS-LES混合方法被越来越多地用来模拟复杂流动问题，但是此类方法有一个最致命的问题就是在模拟附着壁湍流时，存在对数区不匹配现象(Log-layer Mismatch, LLM)。例如，如果用1997年原始的DES方法^[15]来模拟槽道湍流，会出现两段不同的对数律区间，其中底层的RANS区域有一段对数律，而外层的LES区域也有一段对数律。两段对数律的斜率非常接近，但是截距差别明显^[16]。此外，如果DES方法使用的网格不是很好，使得DES方法在很粗的网格里切换到LES分支，这样会使得模型化部分的应力偏小，而总体雷诺应力又没有能力弥补这个缺失，从而出现所谓的模型应力亏损(Modeled Stress Depletion, MSD)^[17]，甚至有可能出现网格引起的虚假分离(Grid-induced Separation，GIS)^[17-18]。为了克服这些缺点，研究人员相继提出了DDES(Delayed DES)^[17]、IDDES(improved DDES)^{[14, 18]}等修正方法。虽然IDDES方法在槽道湍流的验证中消除了LLM现象，但是其引入了新的尺度定义、新的混合函数和数个模型参数，从而使得其适用性大为降低。Spalart悲观地认为IDDES是DES方法发展的一个不幸的趋势^[12]。

为了消除LLM，本文作者及其合作者提出了约束大涡模拟方法(Constrained Large Eddy Simulation, CLES)^[19-20]。CLES方法同之前的DES类方法不同，它全场作LES计算，但是在靠近壁面的区域采用约束了的SGS模型。CLES方法解决了RANS-LES混合方法最主要的缺陷——对近壁区的模拟缺少小尺度脉动，从而消除了LLM。本文将首先详细介绍CLES方法的基本原理，而后对CLES方法应用过程中的一些问题展开讨论分析，最后对CLES方法未来的一些发展方向做了一点概括。对于CLES方法的数值验证，我们将简略带过，读者有兴趣可以参考一些已发表的文献。

1 CLES方法的基本原理

CLES方法是湍流约束思想^[21-22]在LES建模中的应用，本质上还是一种LES方法。与传统的LES方法不同的是，它在靠近壁面的内区需要使用约束的亚网格模型(Constrained Subgrid-scale Model，CSGS)。接下来，我们将首先推导CSGS模型的基本约束等式。不失一般性，我们接下来以不可压缩问题为例。对于三维非定常不可压缩常密度流体，其控制方程如下：

$ \begin{array}{l} \frac{\partial u_{j}}{\partial x_{j}}=0 \\ \frac{\partial u_{j}}{\partial t}+\frac{\partial\left(u_{i} u_{j}\right)}{\partial x_{i}}=-\frac{\partial p}{\partial x_{j}}+\nu \frac{\partial^{2} u_{j}}{\partial x_{i} \partial x_{i}} \end{array} $

(1)

可压缩问题情况类似，具体推导可参考文献[20]。

1.1 RANS控制方程

对上述方程(1)做系综平均运算〈〉，便可得到RANS方程组：

$ \begin{array}{l} \frac{\partial\left\langle u_{j}\right\rangle}{\partial x_{j}}=0 \\ \frac{\partial\left\langle u_{j}\right\rangle}{\partial t}+\frac{\partial\left\langle u_{i}\right\rangle\left\langle u_{j}\right\rangle}{\partial x_{i}}=-\frac{\partial\langle p\rangle}{\partial x_{j}}+\nu \frac{\partial^{2}\left\langle u_{j}\right\rangle}{\partial x_{i} \partial x_{i}}-\frac{\partial R_{i j}}{\partial x_{i}} \end{array} $

(2)

这里$ {R_{ij}} = \left\langle {{u_i}{u_j}} \right\rangle - \left\langle {{u_i}} \right\rangle \left\langle {{u_j}} \right\rangle $为雷诺应力张量。

1.2 LES控制方程

类似地，如果我们对方程(1)做滤波运算~，就可以得到LES方程组：

$ \begin{array}{l} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}=0 \\ \frac{\partial \tilde{u}_{j}}{\partial t}+\frac{\partial\left(\tilde{u}_{i} \tilde{u}_{j}\right)}{\partial x_{i}}=-\frac{\partial \widetilde{p}}{\partial x_{j}}+\nu \frac{\partial^{2} \tilde{u}_{j}}{\partial x_{i} \partial x_{i}}-\frac{\partial \tau_{i j}^{\mathrm{SGS}}}{\partial x_{i}} \end{array} $

(3)

这里，$ \tau _{ij}^{{\rm{SGS}}} \equiv \widetilde {{u_i}{u_j}} - {{\tilde u}_i}{{\tilde u}_j}$为亚网格应力张量。

从方程(2)、(3)可见，LES和RANS方程在形式上完全一致。这也是RANS/LES混合方法的最基本出发点。

1.3 CLES的约束等式

如果我们将方程(3)做一个系综平均，那么其动量方程会变成：

$ \frac{\partial\left\langle\tilde{u}_{j}\right\rangle}{\partial t}+\frac{\partial\left\langle\tilde{u}_{i} \tilde{u}_{j}\right\rangle}{\partial x_{i}}=-\frac{\partial\langle\tilde{p}\rangle}{\partial x_{j}}+\nu \frac{\partial^{2}\left\langle\tilde{u}_{j}\right\rangle}{\partial x_{i} \partial x_{i}}-\frac{\partial\left\langle\tau_{i j}^{\mathrm{SGS}}\right\rangle}{\partial x_{i}} $

(4)

如果我们定义$ R_{i j}^{\mathrm{LES}}=\left\langle\tilde{u}_{i} \tilde{u}_{j}\right\rangle-\left\langle\tilde{u}_{i}\right\rangle\left\langle\tilde{u}_{j}\right\rangle$为可解雷诺应力，那么方程(4)可以改写成：

$ \begin{array}{l} \frac{\partial\left\langle\tilde{u}_{j}\right\rangle}{\partial t}+\frac{\partial\left\langle\tilde{u}_{i}\right\rangle\left\langle\tilde{u}_{j}\right\rangle}{\partial x_{i}} \\ =-\frac{\partial\langle\tilde{p}\rangle}{\partial x_{j}}+\nu \frac{\partial^{2}\left\langle\tilde{u}_{j}\right\rangle}{\partial x_{i} \partial x_{i}}-\frac{\partial\left(\left\langle\tau_{i j}^{\mathrm{SGS}}\right\rangle+R_{i j}^{\mathrm{LES}}\right)}{\partial x_{i}} \end{array} $

(5)

在各态遍历假设条件下，我们有$ \langle \cdot \rangle = \langle \tilde \cdot \rangle $，比较方程(2)和方程(5)，我们可以得到CLES中关于雷诺应力的约束等式：

$\left\langle\tau_{i j}^{\mathrm{SGS}}\right\rangle+R_{i j}^{\mathrm{LES}}=R_{i j}^{\mathrm{RANS}} $

(6)

即：总雷诺应力可以分解为亚网格应力的平均值和可解雷诺应力两部分。对于可压缩问题，我们同样推导出了类似的关于雷诺应力的约束等式以及关于热流的约束等式^[20]。

1.4 CLES方法

有了上面的雷诺应力约束等式，我们就可以构建CSGS模型。对于亚网格应力张量，我们首先将其做一个分解：

$ \tau_{i j}^{\text {model }}=\left\langle\tau_{i j}^{\text {model }}\right\rangle+\left(\tau_{i j}^{\text {model }}\right)^{\prime} $

(7)

它的两部分可以分别给出：其平均部分可以由上面的雷诺应力约束等式(6)给出:

$ \left\langle\tau_{i j}^{\text {model }}\right\rangle=R_{i j}^{\mathrm{RANS}}-R_{i j}^{\mathrm{LES}} $

而其脉动部分我们可以简单地由已有的SGS模型来构造，以Smagorinsky模型^[23]为例：

$ \left(\tau_{i j}^{\operatorname{model}}\right)^{\prime}=C_{s}^{\prime}\left(\Delta^{2}|\tilde{S}| \widetilde{S}_{i j}-\left\langle\Delta^{2}|\widetilde{S}| \widetilde{S}_{i j}\right\rangle\right) $

这里，其系数C′_s可以沿用常系数Smagorinsky模型来获取，也可以结合Germano等式^[24-25]来动态求解。于是，我们就可以构建一个以Smagorinsky模型为基础的CSGS模型：

$ \begin{array}{c} \tau_{i j}^{\text {model }}=R_{i j}^{\mathrm{RANS}}-R_{i j}^{\mathrm{LES}}+C_{s}^{\prime}\left(\Delta^{2}|\tilde{S}| \widetilde{S}_{i j}-\right. \\ \left.\left\langle\Delta^{2}|\tilde{S}| \widetilde{S}_{i j}\right\rangle\right) \end{array} $

(8)

在CLES方法的具体应用的时候，在外区，我们可以使用传统的SGS模型，而在靠近壁面的内区，我们则可以使用CSGS模型。通过上面的推导，我们也发现，虽然我们全场都在做LES计算，但是我们在内区实际上做了一个隐含的RANS计算。由于CLES方法在内区和外区都开展LES计算，因此小尺度脉动能够在全场存在，不会出现RANS/LES混合方法那样由于缺乏近壁小尺度脉动而导致的LLM现象。

CLES方法在不可压缩/可压缩槽道湍流^[19-20]、周期丘陵绕流^[26]、不可压缩/可压缩圆柱绕流^{[19, 27-28]}等算例中得到了验证，并被用来模拟了大型客机外型在迎角14°时的流动状态^[29-31]。图 1给出了在不可压缩槽道湍流中验证结果。由图可见，当CLES方法用来模拟不可压缩槽道湍流时，能得到与DNS结果基本吻合的平均速度型。相反，DES方法出现了LLM，而传统的LES-DSM则会出现平均速度在对数区的整体上移。而且我们还发现，在槽道湍流中，使用代数涡黏模型和SA模型来预估雷诺应力，对结果影响不大。

2 关于CLES方法的一些讨论

CLES方法中，CSGS模型至关重要。

2.1 关于CSGS模型的思考

CSGS建模的想法可追溯到湍流中的约束关系的提出^[21-22]。图 2给出了一些常见的SGS模型的物理/数学考量。由图我们可以看出，大多数SGS模型都是有或多或少的物理/数学考量的。当我们将更多的物理/数学考量加入到SGS模型中，我们就能得到更好的模型。例如当将Germano等式这个数学约束加入到常系数Smagorinsky模型之后，我们就可以得到比常系数模型更准确的动态Smagorinsky模型(Dynamic Smagorinsky Model, DSM)。类似地，我们如果在动态混合相似性模型(dynamic mixed similarity model)基础上加入能流约束，我们就可以得到精确度更高的能流约束模型^[32]；在壁湍流的模拟中，如果我们在DSM基础上加入雷诺应力约束，就能得到雷诺应力约束模型^[19-20]。

需要特别指出的是，约束等式(6)在数学的意义上是一个必要条件，而不是一个充要条件，这和Germano等式类似。事实上，Meneveau^[33]曾经推导过亚网格应力模型需要满足的一些必要条件。在他推导的必要条件里，可以允许两者之间相差一个无散的常张量C_ij。在这里，我们考虑到在壁面上的相容条件，总雷诺应力、可解雷诺应力和亚网格雷诺应力都需要为0，因此C_ij=0。对于其他的约束，我们可以允许差一个常数张量，但是它们在实际计算中并没有贡献，因为雷诺应力是以散度的形式出现在动量方程中的。

此外，我们还需要说明的是基于雷诺应力约束的CSGS模型更多的是一种方法论，采用不同基准的SGS模型，可以得到不同的CSGS模型。例如，在文献[34]中，我们就考虑了基于Vreman模型^[35]的CSGS模型。

2.2 总雷诺应力的影响

在CSGS模型中，总雷诺应力R_ij^RANS对结果影响很明显。如何获取该总雷诺应力是CLES模拟的一个关键。如果有DNS数据，那么我们可以得到准确的总雷诺应力。但是大多数问题的模拟中，我们并没有DNS数据。那么此时，我们获取总雷诺应力就需要考虑精度。在简单的附着流动中，例如槽道湍流，大多数的RANS模型都可以对雷诺应力有一个不错的估计，因此不同RANS模型得到的结果基本相同。对于复杂的分离流动而言，此时总雷诺应力的预估就非常重要了。如果预估的雷诺应力不准确，那么相当于给问题设置一个错误的目标，结果肯定也不好。在文献[26]中，我们讨论了在周期丘陵绕流中不同雷诺应力模型对结果的影响。结果发现，不同雷诺应力模型预估的雷诺应力对结果影响非常大。全场求解RANS模型得到的雷诺应力和实际雷诺应力相差较大，结果很不好。相对而言，使用代数涡黏模型或者基于DES类的长度尺度改进了的RANS模型求解获得的雷诺应力结果更好。在更为复杂的流动计算中，我们更推荐后者。

2.3 CLES方法与传统LES和DES的联系

如果我们以常系数Smagorinsky模型为基础，此时CSGS模型(8)可以被改写为：

$ \begin{array}{c} \tau_{i j}^{\operatorname{model}}=C_{s} \Delta^{2}|\tilde{S}| \tilde{S}_{i j}+R_{i j}^{\mathrm{RANS}}-R_{i j}^{\mathrm{LES}}- \\ C_{s}\left\langle\Delta^{2}|\tilde{S}| \tilde{S}_{i j}\right\rangle \end{array} $

(9)

由此式，我们可以看出CSGS模型可以看成是在传统的LES方法的Smagorinsky模型基础上，在内区外加了一个张量$ R_{i j}^{\mathrm{RANS}}-R_{i j}^{\mathrm{LES}}-C_{s}\left\langle\Delta^{2}|\widetilde{S}| \tilde{S}_{i j}\right\rangle$，该张量在动量方程中就体现为一个外加力。在传统的SGS模型中，我们发现平均速度型在对数区会向上平移，其原因在于内区的亚网格模型雷诺应力偏小。CSGS模型采用约束的方法，于内区在传统的SGS模型基础上增加了一个附加力，从而弥补了传统的SGS模型在内区偏小的亚网格雷诺应力。该结果在文献[19]的分析中也得到印证了。

类似地，我们可以将方程式(8)作另一种改写：

$ \begin{aligned} \tau_{i j}^{\text {model }}=& R_{i j}^{\mathrm{RANS}}-\left[R_{i j}^{\mathrm{LES}}-C_{s}^{\prime}\left(\Delta^{2}|\widetilde{S}| \widetilde{S}_{i j}-\right.\right.\\ &\left.\left.\left\langle\Delta^{2}|\widetilde{S}| \widetilde{S}_{i j}\right\rangle\right)\right] \end{aligned} $

(10)

如果内区的总雷诺应力是由求解一方程或者两方程的DES型(即原始的RANS模型方程中的特征长度尺度按照DES方法做了修改)RANS模型方程来预估的，而外区的SGS模型也是使用DES的LES分支(也可以继续使用Smagorinsky模型)，那么此时CSGS模型可以看成是在DES方法的基础上，在内区加了一个附加力，该力为张量

$ - R_{ij}^{{\rm{LES}}} + C_s^\prime \left( {{\Delta ^2}|\tilde S|{{\tilde S}_{ij}} - \left\langle {{\Delta ^2}|\tilde S|{{\tilde S}_{ij}}} \right\rangle } \right) $

的散度。值得提出的是，Keating和Piomelli^[36]曾经提出过在界面处加入随机里的方法来减弱RANS/LES混合方法中的LLM现象，但这种附加的随机力本身与流场没有相关性，需要根据实际需要调节随机力的强度。CSGS模型事实上在整个内区都增加一个与流场相关的附加力。该附加力使得小尺度脉动在内区一致存在，并且在界面附近得到充分发展，从而可以消除LLM。事实上，在文献[19]中，我们分析了DES方法中的模型应力和可解应力的法向变化规律。同传统的LES和CLES不同，DES方法在内区是RANS分支，此时模型应力占主导，而可解部分的应力其实是随着壁面距离的增加而缓慢增加的。这种外区长度尺度的切换，从外部也改变了内部的雷诺应力分布，从而使得总体雷诺应力偏小。而CLES方法也可以视为在内区增加了一个力，这个力可以使得内区的脉动迅速增长起来，从而改变了内区的可解雷诺应力。

2.4 约束界面的影响

在CLES方法的应用中，其内外区界面的位置决定了使用CSGS模型的区域，从而对结果会产生不同程度的影响。可以预见的是，如果约束区域很小，那么其效果就会不好，当界面的位置设置于固壁上，CLES中CSGS模型的使用区域就会随着消失，CLES方法也会退化为传统的LES方法。而如果CLES中CSGS区域太大，一个问题是我们很难保证远离壁面的外区我们还能得到准确的雷诺应力，特别是针对复杂的分离流动而言。此外，约束区域的增大也会相应增加计算量，毕竟CSGS模型的计算量比传统的SGS模型要大。在槽道湍流等简单流动中，我们研究过不同约束界面位置对结果的影响。如分析的那样，如果约束界面位于缓冲区(y⁺=5.05)，那么平均速度会上移；如果界面位于对数区(y⁺=20.1和44.9), 结果基本相同^[19]。如果用CLES方法模拟复杂壁湍流，如何选取界面位置的确是一个需要研究的问题。目前大多数计算中，CLES的界面位置会参考DES的界面位置，并比DES的位置略微靠近壁面一点。

此外，我们还需要考虑内外区之间的衔接问题。虽然CSGS模型和传统的SGS模型在数学上不存在连续性问题，但是在实际计算中，两种不同的模型在固定的界面处连接在一起，势必会在界面处出现一些统计量不连续的情况。为了避免这个问题，我们可以参照一些RANS/LES混合方法的做法，可以在界面处使用一个连接函数将两个区域连接起来。一个最简单的尝试就是界面处的亚网格应力为两种模型的平均值。

2.5 系综平均的实现

在CSGS模型中，我们看到需要进行系综平均的操作，这是一个耗费巨大计算量的工作。在简单的流动里，系统本身存在统计均匀方向，此时我们可以用统计均匀方向的空间平均来近似系综平均。例如在槽道湍流中，我们可以使用面平均来近似系综平均。然而，在一般的复杂流动中，可能不存在统计均匀方向，这时我们无法采用空间平均来近似系综平均。在这类问题中，我们可以仿照Meneveau等^[37]构造Lagrangian平均来近似，但是Lagrangian平均往往计算复杂，计算量也很大，在复杂流动的计算中很难实用。另一个选择是采用Euler时间平均来近似系综平均：

$ \begin{array}{l} \langle f\rangle^{n+1}=\varepsilon f^{n+1}+(1-\varepsilon)\langle f\rangle^{n} \\ \varepsilon=\Delta t /(\Delta t+T) \end{array} $

(11)

这里Δt是计算的时间步长，T是Euler时间平均的特征间隔时间。在槽道湍流里，我们验证了时间平均能够得到和空间平均相近的结果^[19]。

3 CLES方法的未来研究建议

CLES方法经过近十年的发展，已经在很多经典问题中得到验证，并已经被用于模拟复杂的工程流动问题。未来，我们希望CLES方法能够被更多地应用到工程湍流问题的模拟中去。对于CLES方法的理论研究，我们认为一下几个方向值得研究：

(1) CLES的转捩预测能力的研究与验证。在层流-湍流转捩过程中，我们需要模型能够跟踪捕获这一变化过程。在层流区域，通过分析，我们不难发现CSGS模型会退化为τ_ij^model=R_ij^RANS，因为层流区几乎没有脉动，因此可解应力等项会消失。因此CSGS模型模拟转捩的能力就很大程度上依赖于预估的雷诺应力了。在早期的槽道流时间转捩模拟的研究中，我们发现如果将雷诺应力在转捩过程中的行为考虑进去，CLES的确可以具备自适应地模拟转捩的能力^[38]。在空间转捩问题中的表现如何，这个还需要进一步研究和验证。另外，在复杂的湍流问题的模拟中，也牵涉到剪切层的转捩问题，例如高雷诺数钝体绕流问题中，分离后的剪切层的发展演化。此类问题的精确模拟，对分离、再附等问题准确模拟至关重要，也需要深入的研究和验证。

(2) 复杂问题模拟中，如何给出准确的总雷诺应力的预估值？总雷诺应力是CLES方法中的目标值，因此对CLES方法来说，总雷诺应力的预估值至关重要，会直接影响CLES最终的模拟结果。在简单的壁湍流问题中，大多数RANS模型都可以给出较好的预测结果；在一般的复杂流动中，如何在CSGS应用区域中给出合理的总雷诺应力的预估值是值得研究的问题。在这里，基于机器学习的RANS建模可望提供一种新思路。

(3) 约束界面位置的选取方法。在过去的研究中，我们通常采用的固定的界面位置。我们知道的是，固定的界面位置不能太靠近壁面，否则约束效果不好。如何选取合理的界面位置是一个值得研究的问题。

(4) 约束大涡模拟在气动声学和气动光学中的应用。CLES能够模拟出比传统RANS-LES混合方法更加丰富而精细的小尺度流动结构，这些小尺度流动结构对光学传播有重要的影响，同时也是气动噪声的主要的声源。因此CLES方法在气动声学和气动光学的模拟方面有广阔的应用前景。