大涡模拟

气体速度场的大涡模拟。]]
大涡模拟(Large eddy simulation, LES) 是用于计算流体动力学的湍流数学模型。它最初由约瑟夫·斯玛格林斯基于 1963 年提出，用于模拟大气气流，并由迪尔多夫(1970)推广。 LES 目前广泛应用于多个研究领域，包括燃烧、声学、和大气边界层模拟。

通过数值求解Navier-Stokes 方程来模拟湍流需要解决非常广泛的时间和长度尺度，所有这些都会影响流场。这样的分辨率可以通过直接数值模拟(DNS) 来实现，但 DNS 的计算成本很高，难以模拟具有复杂几何形状或流动配置的实际工程系统，例如湍流喷射、泵、车辆和起落架。

LES 背后的主要思想是通过 Navier-Stokes 方程的低通滤波忽略最小长度尺度来降低计算成本。这种可以被视为时间和空间平均的低通滤波，有效地从数值解中去除了小尺度信息。然而，这些信息并不是无关紧要的，它对流场的影响必须被建模，由此衍生小尺度可以发挥重要作用的问题的活跃研究领域，例如近壁流，反应流、

过滤器定义和属性
LES 滤波器可应用于空间和时间场\phi(\boldsymbol{x},t)并执行空间滤波操作、时间滤波操作或两者一并进行。过滤后的流场加了上划线，定义为：

:
\overline{\phi(\boldsymbol{x},t)} = \displaystyle{
\int_{-\infty}^{\infty}} \int_{-\infty}^{\infty} \phi(\boldsymbol{r},\tau) G(\boldsymbol{x}-\boldsymbol{r},t - \tau) d\tau d \boldsymbol{r}

其中G是滤波器卷积核。这也可以写成：

:
\overline{\phi} = G \star \phi .

过滤器内核G有一个相关的截止长度尺度\Delta和截止时间尺度\tau_{c} .小于这些的尺度将被消掉。使用上述过滤器定义，任何流场\phi可以分为过滤和子过滤（用素数表示）部分，如

:
\phi = \bar{\phi} + \phi^{\prime} .

需要注意的是，大涡模拟滤波操作不满足雷诺算子的性质。

过滤的控制方程
LES的控制方程是通过过滤控制流场的偏微分方程得到的\rho \boldsymbol{u}(\boldsymbol{x},t) .不可压缩和可压缩 LES 控制方程之间存在差异，这导致了新滤波操作的定义。

不可压缩流动
对于不可压缩流动，对连续性方程和 Navier-Stokes 方程进行滤波，得到滤波后的不可压缩连续性方程，

:
\frac{ \partial \bar{u_i} }{ \partial x_i } = 0

和过滤后的 Navier-Stokes 方程，

:
\frac{ \partial \bar{u_i} }{ \partial t }

\frac{ \partial }{ \partial x_j } \left( \overline{ u_i u_j } \right)

= - \frac{1}{\rho} \frac{ \partial \overline{p} }{ \partial x_i }

\nu \frac{\partial}{\partial x_j} \left( \frac{ \partial \bar{u_i} }{ \partial x_j } + \frac{ \partial \bar{u_j} }{ \partial x_i } \right)

= - \frac{1}{\rho} \frac{ \partial \overline{p} }{ \partial x_i }

2 \nu \frac{\partial}{\partial x_j} \bar{S}_{ij},

其中\bar{p}是过滤后的压力场和\bar{S}_{ij}是使用过滤速度评估的应变率张量。非线性滤波平流项\overline{u_i u_j}是LES建模困难的主要原因。它需要未过滤的速度场，这是未知的，因此必须对其进行建模。下面的分析说明了非线性带来的困难，即它导致大小尺度之间的相互作用，防止尺度分离。

过滤后的平流项可以按照莱昂纳德(1975)拆分为：

:
\overline{u_i u_j} = \tau_{ij} + \overline{u}_i \overline{u}_j

其中\tau_{ij}是残余应力张量，因此过滤后的 Navier-Stokes 方程变为

:
\frac{ \partial \bar{u_i} }{ \partial t }

\frac{ \partial }{ \partial x_j } \left( \overline{u}_i \overline{u}_j \right)

= - \frac{1}{\rho} \frac{ \partial \overline{p} }{ \partial x_i }

2 \nu \frac{\partial}{\partial x_j} \bar{S}_{ij}
\frac{ \partial \tau_{ij} }{ \partial x_j }

与残余应力张量\tau_{ij}对所有未封闭的项进行分组。 Leonard 将这个应力张量分解为\tau_{ij} = L_{ij} + C_{ij} + R_{ij}并为每个项提供物理解释。 L_{ij} = \overline{ \bar{u}_{i} \bar{u}_{j} } - \bar{u}_{i} \bar{u}_{j} 是伦纳德张量，代表大尺度之间的相互作用；R_{ij} = \overline{u^{\prime}_{i} u^{\prime}_{j}} 为类雷诺应力项，表示子滤波器尺度 (SFS) 之间的相互作用；C_{ij} = \overline{\bar{u}_{i} u^{\prime}_{j}} + \overline{\bar{u}_{j} u^{\prime}_{i}}
是克拉克张量，表示大尺度和小尺度之间的跨尺度相互作用。提出了一种密度加权滤波操作，称为 Favre 滤波，定义为任意量\phi为：

:
\tilde{\phi} = \frac{ \overline{\rho \phi} }{ \overline{\rho} }

在不可压缩的限制下，它变成了正常的过滤操作。这使得质量守恒方程化为：

:
\frac{\partial \overline{\rho}}{\partial t} + \frac{ \partial \overline{\rho} \tilde{u_i} }{ \partial x_i } = 0.

然后可以将此概念扩展到编写可压缩流的 Favre 过滤动量方程。根据弗雷曼：

:
\frac{ \partial \overline{\rho} \tilde{u_i} }{ \partial t }

\frac{ \partial \overline{\rho} \tilde{u_i} \tilde{u_j} }{ \partial x_j }
\frac{ \partial \overline{p} }{ \partial x_i }
\frac{ \partial \tilde{\sigma}_{ij} }{ \partial x_j }

= - \frac{ \partial \overline{\rho} \tau_{ij}^{r} }{ \partial x_j }

\frac{ \partial }{ \partial x_j } \left( \overline{\sigma}_{ij} - \tilde{\sigma}_{ij} \right)

其中\sigma_{ij}是剪应力张量，对于牛顿流体，由下式给出：

:
\sigma_{ij} = 2 \mu(T) S_{ij} - \frac{2}{3} \mu(T) \delta_{ij} S_{kk}

\frac{ \partial }{\partial x_j} \left( \overline{\sigma}_{ij} - \tilde{\sigma}_{ij} \right)这一项表示评估粘度的子过滤器粘性贡献\mu(T)使用 Favre 过滤温度\tilde{T} . Favre 滤波动量场的亚网格应力张量：
\tau_{ij}^{r} = \widetilde{ u_i \cdot u_j } - \tilde{u_i} \tilde{u_j}

以此类推，莱昂纳德分解也可以写成滤波三重积的残余应力张量\overline{\rho \phi \psi} 。三重乘积可以使用 Favre 过滤运算符重写为\overline{\rho} \widetilde{\phi \psi} ，这是一个未封闭的项（它需要场\phi和\psi的信息, 当只有\tilde{\phi}和\tilde{\psi}两项是已知的时候）。它可以以类似于上面提到的\overline{u_i u_j}的方式分解，这导致子滤波器应力张量\overline{\rho} \left( \widetilde{\phi \psi} - \tilde{\phi} \tilde{\psi} \right) 。这个子过滤器项可以分为三种类型的相互作用的贡献：伦达德张量L_{ij}, 代表解析尺度之间的相互作用；克拉克张量C_{ij} ，表示已解析和未解析的尺度之间的相互作用；和雷诺张量R_{ij} ，它表示未解析的尺度之间的相互作用。

滤波动能方程
除了过滤的质量和动量方程之外，过滤动能方程可以提供额外的见解。可以过滤动能场以产生总过滤动能：

:
\overline{E} = \frac{1}{2} \overline{ u_i u_i }

过滤后的总动能可以分解为两项：过滤后的速度场的动能E_f ,

:
E_f = \frac{1}{2} \overline{u_i} \, \overline{u_i}

和剩余动能k_r ,

:
k_r = \frac{1}{2} \overline{ u_i u_i } - \frac{1}{2} \overline{u_i} \, \overline{u_i} = \frac{1}{2} \tau_{ii}^{r}

这样\overline{E} = E_f + k_r .

守恒方程为E_f可以通过将滤波后的动量传输方程乘以\overline{u_i}产生：

:
\frac{\partial E_f}{\partial t}

\overline{u_j} \frac{\partial E_f}{\partial x_j}
\frac{1}{\rho} \frac{\partial \overline{u_i} \bar{p} }{ \partial x_i }
\frac{\partial \overline{u_i} \tau_{ij}^{r}}{\partial x_j}
2 \nu \frac{ \partial \overline{u_i} \bar{S_{ij}} }{ \partial x_j }

\epsilon_{f}
\Pi

其中\epsilon_{f} = 2 \nu \bar{S_{ij}} \bar{S_{ij}}是粘性应力对过滤速度场的动能的耗散，并且\Pi = -\tau_{ij}^{r} \bar{S_{ij}}表示动能的子过滤尺度 (SFS) 耗散。

左边的项代表运输，右边的项是耗散动能的汇项。

LES的数值方法
大涡模拟涉及使用计算流体动力学求解离散滤波控制方程。 LES 从域大小中解析尺度L过滤器尺寸\Delta ，因此必须解决很大一部分高波数湍流波动。这需要高阶数值方案，或者如果使用低阶数值方案，则需要精细的网格分辨率。波普的第 13 章解决了网格分辨率有多精细的问题\Delta x需要解析过滤的速度场\overline{u}(\boldsymbol{x}) . Ghosal 发现，对于低阶离散化方案，例如在有限体积方法中使用的那些，截断误差可以与子滤波器尺度贡献相同，除非滤波器宽度\Delta比网格间距\Delta x大得多。虽然偶数阶方案具有截断误差，但它们是非耗散的，并且由于子滤波器尺度模型是耗散的，偶数阶方案不会像耗散方案那样强烈地影响子滤波器尺度模型的贡献。

过滤器实现
大涡模拟中的滤波操作可以是隐式的，也可以是显式的。隐式过滤认识到子过滤器比例模型将以与许多数值方案相同的方式消散。通过这种方式，可以假设网格或数值离散化方案是 LES 低通滤波器。虽然这充分利用了网格分辨率，并消除了计算子滤波器比例模型项的计算成本，但很难确定与一些数值问题相关的 LES 滤波器的形状。此外，截断误差也可能成为问题。

在显式滤波中， LES 滤波器应用于离散的 Navier-Stokes 方程，提供明确定义的滤波器形状并减少截断误差。然而，显式过滤比隐式过滤需要更精细的网格，并且计算量与(\Delta x)^4成正比。Sagaut (2006) 的第 8 章更详细地介绍了 LES 的数值问题。

(1) 提供准确的流动特性信息，即速度和湍流；

(2) 满足 Navier-Stokes 方程和其他物理；

(3)易于实施和适应不同情况。

目前，为 LES 生成入口条件的方法大致分为 Tabor 等人分类的两类：

产生湍流入口的第一种方法是根据具体情况合成它们，例如傅里叶技术、原理正交分解（POD）和涡流方法。合成技术试图在入口处构建具有合适的类湍流特性的湍流场，并使其易于指定湍流参数，例如湍流动能和湍流耗散率。此外，使用随机数生成的入口条件在计算上并不昂贵。然而，该方法存在一个严重的缺陷。合成的湍流不满足由 Navier-Stokes 方程控制的流体流动的物理结构。存在两类 SGS 模型；第一类是功能模型，第二类是结构模型。一些模型可能被归类为两者。

功能（涡粘）模型
功能模型比结构模型更简单，只关注以物理上正确的速率耗散能量。这些基于人工涡流粘度方法，其中湍流的影响集中在湍流粘度中。该方法将亚网格尺度上的动能耗散视为类似于分子扩散。在这种情况下，偏斜部分\tau_{ij}被建模为：

:
\tau_{ij}^r - \frac{1}{3} \tau_{kk} \delta_{ij} = -2 \nu_\mathrm{t} \bar{S}_{ij}

其中\nu_\mathrm{t}是湍流涡流粘度和\bar{S}_{ij} = \frac{1}{2} \left( \frac{\partial \bar{u}_i }{\partial x_j} + \frac{\partial \bar{u}_j}{ \partial x_i} \right)是应变率张量。

根据量纲分析，涡流粘度的单位必须为\left[ \nu_\mathrm{t} \right] = \frac{\mathrm{m^2}}{\mathrm{s}} 。大多数涡流粘度 SGS 模型将涡流粘度建模为特征长度尺度和特征速度尺度的乘积。

Smagorinsky-Lilly 模型
第一个成功开发的 SGS 模型是 Smagorinsky-Lilly SGS 模型，它由Smagorinsky 使用 Smagorinsky 模型确定了许多研究，每个研究都发现了不同的 Smagorinsky 常数值C针对不同的流量配置。为了试图为 SGS 模型制定一种更通用的方法，Germano 等提出了一个动态 Smagorinsky 模型，它使用了两个过滤器：一个网格 LES 过滤器，表示为\overline{f} ，以及一个测试 LES 滤波器\hat{f}，用于任何湍流的场f 。测试过滤器的尺度大于网格过滤器，并在 LES 表示的已经平滑的场上增加了对湍流场的额外平滑。将测试滤波器应用于 LES 方程（通过将“网格”滤波器应用于 Navier-Stokes 方程获得）会产生一组新的方程，它们的形式相同但 SGS 应力项\tau_{ij} = \overline{u_{i} u_{j}} - \bar{u}_{i} \bar{u}_{j}被T_{ij} = \widehat{\overline{u_{i} u_{j}}} - \hat{\bar{u}}_{i} \hat{\bar{u}}_{j}所替代。 Germano等注意到即使两者都由于存在未解析的尺度而无法精确计算，这两个张量之间仍然存在一个精确的关系。这种关系被称为 Germano 恒等式：

L_{ij} = T_{ij} - \hat{\tau}_{ij}.

这里 L_{ij} = \widehat{\bar{u}_{i} \bar{u}_{j}} - \widehat{\bar{u}_{i}} \widehat{\bar{u}_{j}}可以直接计算得出，因为它只涉及过滤的速度和测试过滤的操作。同一性的意义在于，如果假设湍流是自相似的，那么网格和测试级别的 SGS 应力具有相同的形式\tau_{ij} - (\tau_{kk}/3)\delta_{ij} = - 2 C \Delta^{2} |\bar{S}_{ij}| \bar{S}_{ij}和T_{ij} - (T_{kk}/3)\delta_{ij} = - 2 C \hat{\Delta}^{2} |\hat{\bar{S}}_{ij}| \hat{\bar{S}}_{ij} ，则 Germano 恒等式提供了一个方程，从中可以得到 Smagorinsky 系数C （不再是“常数”）可能得以推出。

为了做到这一点，需要在原始推导中引入x两个额外的步骤。首先，假设即使C原则上是可变的，变化足够慢，可以从过滤操作中移出： \widehat{C (.)} = C \widehat{(.)} 。其次，既然C是一个标量，Germano 恒等式与一个二阶张量（选择应变张量的速率）联系起来，以将其转换为一个标量方程，从而推出C。 Lilly 发现了一种从张量恒等式中推出 C的更好的方法。他指出，Germano 恒等式需要在空间中的每个点上满足单个量的九个方程（其中只有五个是独立的） C值。因此对C的推算条件多余了。他提议C通过最小化残差使用最小二乘拟合来确定。由此得：
C = \frac{ L_{ij} m_{ij} }{ m_{kl} m_{kl} }.

其中
m_{ij} = \alpha_{ij} - \widehat{\beta}_{ij}
。为简洁起见 \alpha_{ij} = - 2 \hat{\Delta}^{2} | \hat{\bar{S}} | \hat{\bar{S}}_{ij} , \beta_{ij} = - 2 \Delta^2 | \bar{S} | \bar{S}_{ij} 在 LES 模拟中实施该模型的最初尝试被证明是不成功的。首先，计算出的系数根本不像假设的那样“缓慢变化”，而且变化与任何其他湍流场一样多。其次，计算C可以是积极的，也可以是消极的。后一个事实本身不应被视为缺点，因为使用过滤的 DNS 字段的先验测试表明，本地子网格耗散率 - \tau_{ij} \bar{S}_{ij}即使在流体域上的积分始终为正，表示大尺度上的能量净耗散，在湍流场中的积分几乎与正数一样可能为负数。与涡粘度的严格正相反，正值的轻微优势导致观察到的净耗散。这种从小到大的能量“反向散射”确实对应于 Smagorinsky 模型中的负 C 值。然而，发现 Germano-Lilly 公式不能产生稳定的计算。通过在同质方向上平均分子和分母（在流动中存在这样的方向），采用了一种特别的方法

:
C = \frac{
\left\langle L_{ij} m_{ij} \right\rangle
}{
\left\langle m_{kl} m_{kl} \right\rangle
}.

当平均涉及足够大的统计样本时，计算C是正的（或至少极少负值）稳定的计算是可能的。简单地将负值设置为零（一个称为“剪裁”的过程），无论有或没有平均，也会提高计算稳定性。 Meneveau 提出对具有指数衰减“记忆”的拉格朗日流体轨迹进行平均。这可以应用于缺乏均匀方向的问题，并且如果进行平均的有效时间足够长，但又不足以消除感兴趣的空间不均匀性，那么这可以是稳定的。

Lilly 对 Germano 方法的修改，然后是统计平均或综合去除负粘度区域似乎是临时的，即使它可以“工作”。 Ghosal 等人提出了一种称为“动态定位模型”（DLM）的最小二乘最小化过程的替代公式。在这种方法中，首先定义一个量

:
E_{ij} = L_{ij} - T_{ij} + \hat{\tau}_{ij}

张量\tau_{ij}和T_{ij}替换为适当的 SGS 模型。然后，该张量表示子网格模型未能在每个空间位置符 Germano 恒等式的数量。在Lilly的方法中，C可以脱离帽符算子：

:
\widehat{C (.)} = C \widehat{(.)}

使得E_{ij}成为一个C的代数函数，然后通过要求确定E_{ij} E_{ij}被认为是C的函数具有最小的可能值。然而，由于C如此获得的结果与湍流中的任何其他波动量一样可变，原始假设C不能后天证明。在 DLM 方法中，通过不调用从测试过滤操作中删除C的步骤来避免这种不一致。相反，人们通过数量定义了整个流的全局误差：

:
E [ C ] = \int E_{ij} E_{ij} dV

其中积分范围在整个流体体积上。这个全局误差E[C(x,y,z,t)]是空间变化函数C(x,y,z,t) 的函数（这里的瞬时t是固定的，只是作为一个参数出现），C的结果需要让误差最小化。这个变分问题的解决方法是C必须满足第二类 Fredholm 积分方程

:
C (\boldsymbol{x}) = f ( \boldsymbol{x} ) +
\int K(\boldsymbol{x}, \boldsymbol{y}) C ( \boldsymbol{y} ) d\boldsymbol{y}

其中函数K(\boldsymbol{x}, \boldsymbol{y})和f ( \boldsymbol{x} )根据解析的项L_{ij},\alpha_{ij},\beta_{ij}定义，因此在每个时间步长和整个流体域的积分范围都是已知的。积分方程通过迭代程序进行数值求解，如果与预处理方案一起使用，则发现收敛通常很快。尽管这种变分方法消除了 Lilly 方法中固有的不一致性，但C(x,y,z,t)从积分方程获得的仍然显示出与负粘度相关的不稳定性。这可以通过坚持来解决E[C]受约束最小化C(x,y,z,t) \geq 0 .这导致了一个C的非线性方程：

:
C (\boldsymbol{x}) = \left[ f ( \boldsymbol{x} ) +
\int K(\boldsymbol{x}, \boldsymbol{y}) C ( \boldsymbol{y} ) d\boldsymbol{y} \right]_{+}

这里后缀 + 表示“正数部分”，即， x_{+} = (x + |x|)/2 。尽管这表面上看起来像“剪裁”，但它不是一个临时方案，而是约束变分问题的真正解决方案。这种 DLM(+) 模型被发现是稳定的，并且对于强制和衰减的各向同性湍流、通道流动和各种其他更复杂的几何形状产生了出色的结果。如果流动恰好具有均匀的方向（让我们说方向 x 和 z），那么可以引入假定 C = C(y,t) 。然后，变分方法立即产生 Lilly 的结果，对均匀方向进行平均，无需对先前结果进行临时修改。

DLM(+) 模型的一个缺点是它没有描述反向散射，而反向散射在分析 DNS 数据时被认为是真实的“事物”。开发了两种方法来解决这个问题。由于 Carati 等人的一种方法。类似于朗道的脉动流体力学理论，增加了一个由涨落耗散定理确定的振幅的脉动力。在第二种方法中，有人注意到任何“反向散射”能量出现在解析尺度中，只是以亚网格尺度中的能量为代价。 DLM 可以通过一种简单的方式进行修改，以考虑到这一物理事实，从而在本质上稳定的同时允许反向散射。在DLM 的这个 k 方程版本中，DLM(k) 替换在 Smagorinsky 涡流粘度模型中的 \Delta | \bar{S} | 一项，以 \sqrt{k} 作为适当的速度尺度。确定C的流程保持与“无约束”版本相同，除了张量 \alpha_{ij} = - 2 \hat{\Delta} \sqrt{K} \hat{\bar{S}}_{ij} , \beta_{ij} = - 2 \hat{\Delta} \sqrt{k} \bar{S}_{ij} （其中子测试尺度动能 K 与子网格尺度动能 k 的关系为 K = k + L_{ii}/2 ）（接着追踪 Germano 恒等式）。为了确定 k，我们现在使用传输方程

:
\frac{\partial k}{\partial t} + u_{j} \frac{\partial k}{\partial x_{j}} =

\tau_{ij} \bar{S}_{ij} - \frac{C_{*}}{\Delta} k^{3/2} +

\frac{\partial }{\partial x_j} \left( D \Delta \sqrt{k} \frac{\partial k }{\partial x_j} \right)

\nu \frac{\partial^{2} k }{\partial x_j \partial x_j}

其中\nu是运动粘度，C_{*},D是分别代表动能耗散和扩散的正系数。这些可以按照 DLM(+) 中的约束最小化的动态过程来确定。这种方法虽然比 DLM(+) 实施起来更昂贵，但被发现是稳定的，并且与测试的各种流的实验数据有很好的一致性。此外，DLM(k) 在数学上不可能导致计算不稳定，因为大尺度和 SGS 能量的总和不会因构造而增加。这两种结合反向散射的方法都运作良好。与 DLM(+) 相比，它们产生的模型耗散稍小，性能有所提高。 DLM(k) 模型还产生子网格动能，这可能是一个感兴趣的物理量。这些改进是在模型实现的成本有所增加的情况下实现的。

动态模型起源于 1990 年斯坦福大学湍流研究中心(CTR) 的[https://ctr.stanford.edu/ctr-summer-program 暑期项目] 。一系列“CTR-Tea”研讨会庆祝了湍流建模这一重要里程碑的[https://ctr.stanford.edu/event-type/ctr-tea 30 周年] 。

另见

流体力学
伽利略不变性——某些类型滤波器的重要属性
雷诺平均 Navier-Stokes 方程
湍流

评论 (0)