我们知道，现阶段数值模拟湍流的方法主要有三个层次：

1. 雷诺平均（RANS）应用最广泛，但雷诺平均改变NS方程的瞬时性，湍流模型的封闭一直是其难点；
2. 大涡模拟（LES）方法建立在人们对湍流的基本认识：大尺度涡由流场边界、流体粘性等决定，具有明显的各向异性和非均匀性，小尺度涡基本是均匀各向同性的。通过滤波建立大涡方程，通过亚格子应力来模化小尺度湍流。
3. 直接数值模拟（DNS），精度最高，但计算量巨大，工程应用方面具有算法和计算资源的困难。

下面记录自己学习LES的一些收获。

滤波函数

常见的滤波函数

Name	Filter function	Transfer function
General	\(G(r)\)	\(\hat{G}(k)=\int_{-\infty}^{\infty}e^{ikr}G(r)dr\)
Box	$\frac{1}{\Delta}H(\frac{1}{2}\Delta-	r
Gaussian	\((\frac{6}{\pi \Delta^2})^{1/2}exp(-\frac{6r^2}{\Delta^2})\)	\(exp(-\frac{k^2\Delta^2}{24})\)
Sharp spectral	\(\frac{sin(\pi r / \Delta)}{\pi r}\)	$H(k_c-
Cauchy	\(\frac{a}{\pi\Delta[(r/\Delta)^2+a^2]}, a=\frac{\pi}{24}\)	$exp(-a\Delta
Pao		$exp(-\frac{\pi^{2/3}}{24}(\Delta

假设有一个一维的场变量\(u=u(x)\)，则其滤波后表达式如下：

\[\bar{u}(x)=\int_{-\infty}^{\infty}G(r)u(x-r)dr \]

可以看到，滤波后函数其实是原函数与滤波器的卷积。

下面我们已最为常见的帽形函数（top-hat）为例：

\[G(x-x^{'})=\left\{ \begin{array}{ll} 1/\Delta & |x-x^{'}|\le\Delta/2 \\ 0& |x-x^{'}|>\Delta/2 \end{array} \right. \]

假设有一个一维的变量：

\[u(x)=10sin(x)+sin(wx),x\in[-\pi,\pi],w=\frac{\pi}{\Delta} \]

这是一个高频振荡的函数，图像如下：

对该函数进行滤波：

\[\Delta=\frac{2\pi}{N} \]

\[\bar{u}(x_i)=\int_{-\infty}^{\infty}(10sin(x)+sin(wx))G(x_i-x)dx \]

\[\bar{u}(x_i)=\int_{x_i-\Delta/2}^{x_i+\Delta/2}(10sin(x)+sin(wx))\frac{1}{\Delta}dx \]

\[\bar{u}(x_i)=\frac{1}{\Delta}[-10cos(x)-\frac{1}{w}cos(wx)]|_{x_i-\Delta/2}^{x_i+\Delta/2} \]

绘制滤波前后图像如下：

可见滤波就是会将小尺度（高频）脉动量过滤掉。滤波尺度不一样，滤波后的形态也不一致。

对NS方程做滤波运算

NS方程如下：

\[\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho \vec{u})=S_m \]

\[\frac{\partial (\rho \vec{u})}{\partial t}+\nabla\cdot(\rho \vec{u}\vec{u})=\nabla \cdot \vec{\vec{\sigma}}+\vec{f}+S_u \]

\[\frac{\partial (\rho E)}{\partial t}+\nabla\cdot(\rho \vec{u}E)=\nabla \cdot (\vec{u}\cdot\vec{\vec{\sigma}})+\vec{u}\cdot\vec{f}+\nabla\cdot\vec{q}+S_e \]

针对无源不可压缩动量方程，如下：

\[\frac{\partial \mathbf u}{\partial t}+\nabla\cdot(\mathbf u \mathbf u)=-\nabla\frac{p}{\rho} + \nabla \cdot (\nu \nabla \mathbf u) \]

滤波后有：

\[\frac{\partial \bar{\mathbf u}}{\partial t}+\nabla\cdot(\overline{\mathbf u \mathbf u})=-\nabla\frac{\bar{p}}{\rho} + \nabla \cdot (\nu \nabla \bar{\mathbf u}) \]

注意到：

\[\nabla\cdot(\overline{\mathbf u \mathbf u})=\nabla\cdot(\bar{\mathbf u}\bar{\mathbf u})+(\nabla\cdot(\overline{\mathbf u \mathbf u})-\nabla\cdot(\bar{\mathbf u}\bar{\mathbf u})) \]

所以有：

\[\frac{\partial \bar{\mathbf u}}{\partial t}+\nabla\cdot(\bar{\mathbf u}\bar{\mathbf u})=-\nabla\frac{\bar{p}}{\rho} + \nabla \cdot (\nu \nabla \bar{\mathbf u})-(\nabla\cdot(\overline{\mathbf u \mathbf u})-\nabla\cdot(\bar{\mathbf u}\bar{\mathbf u})) \]

其中：

\[\nabla\cdot(\overline{\mathbf u \mathbf u})=\nabla\cdot\overline{ \left[ \begin{array}{l} u_1 \\ u_2 \\ u_3 \end{array} \right] \left[u_2\ u_2 \ u_3 \right] }=\nabla\cdot\left[ \begin{array}{lll} \overline{u_1 u_1} & \overline{u_1 u_2} & \overline{u_1 u_3} \\ \overline{u_2 u_1} & \overline{u_2 u_2} & \overline{u_2 u_3} \\ \overline{u_3 u_1} & \overline{u_3 u_2} & \overline{u_3 u_3} \\ \end{array} \right ] \]

\[\nabla\cdot(\bar{\mathbf u}\bar{\mathbf u})=\nabla\cdot \left[\begin{array}{l}\bar{u}_1 \\ \bar{u}_2 \\ \bar{u}_3\end{array}\right ] \left[\bar{u}_1\ \bar{u}_2\ \bar{u}_3\right]=\nabla\cdot\left[ \begin{array}{lll} \bar{u}_{1} \bar{u}_{1} & \bar{u}_{1} \bar{u}_{2} & \bar{u}_{1} \bar{u}_{3} \\ \bar{u}_{2} \bar{u}_{1} & \bar{u}_{2} \bar{u}_{2} & \bar{u}_{2} \bar{u}_{3} \\ \bar{u}_{3} \bar{u}_{1} & \bar{u}_{3} \bar{u}_{2} & \bar{u}_{3} \bar{u}_{3} \end{array} \right] \]

所以有：

\[-(\nabla\cdot(\overline{\mathbf u \mathbf u})-\nabla\cdot(\bar{\mathbf u}\bar{\mathbf u}))=-\nabla\cdot \left[ \begin{array}{lll} \overline{u_{1} u_{1}}-\bar{u}_{1} \bar{u}_{1} & \overline{u_{1} u_{2}}-\bar{u}_{1} \bar{u}_{2} & \overline{u_{1} u_{3}}-\bar{u}_{1} \bar{u}_{3} \\ \overline{u_{2} u_{1}}-\bar{u}_{2} \bar{u}_{1} & \overline{u_{2} u_{2}}-\bar{u}_{2} \bar{u}_{2} & \overline{u_{2} u_{3}}-\bar{u}_{2} \bar{u}_{3} \\ \overline{u_{3} u_{1}}-\bar{u}_{3} \bar{u}_{1} & \overline{u_{3} u_{2}}-\bar{u}_{3} \bar{u}_{2} & \overline{u_{3} u_{3}}-\bar{u}_{3} \bar{u}_{3} \end{array} \right] \]

定义亚格子应力如下：

\[\tau=\left[\begin{array}{lll} \overline{u_{1} u_{1}}-\bar{u}_{1} \bar{u}_{1} & \overline{u_{1} u_{2}}-\bar{u}_{1} \bar{u}_{2} & \overline{u_{1} u_{3}}-\bar{u}_{1} \bar{u}_{3} \\ \overline{u_{2} u_{1}}-\bar{u}_{2} \bar{u}_{1} & \overline{u_{2} u_{2}}-\bar{u}_{2} \bar{u}_{2} & \overline{u_{2} u_{3}}-\bar{u}_{2} \bar{u}_{3} \\ \overline{u_{3} u_{1}}-\bar{u}_{3} \bar{u}_{1} & \overline{u_{3} u_{2}}-\bar{u}_{3} \bar{u}_{2} & \overline{u_{3} u_{3}}-\bar{u}_{3} \bar{u}_{3} \end{array}\right] \]

亚格子模型

可以看到，对NS方程做滤波之后出现了亚格子应力项，其描述小尺度涡对于大尺度涡的影响。需要利用模型将其封闭。

计算注意事项

由于LES求解非稳态N-S方程，因此不同于RANS，LES的边界条件有一些其他特点：

1. 初始条件只会影响非依时类流动需要达到最终稳定的时间。最好在LES计算中添加一定的随机扰动。如果一个依时类流动和初始条件关联很大，那么推荐使用DNS或者实验数据来获取扰动；
2. 在壁面附近，可以使用壁面函数或者使用足够致密的网格以使得y+<1；
3. 进口边界条件非常敏感，因为进口条件对下游流场有很大的影响。最简单的进口条件为指定一个平均速度型线并附加一个合适湍流强度的随机扰动。但是这种方法有失雷诺应力和空间的联系性。目前典型的做法为：（1）先使用RANS计算进口的雷诺应力，然后把这些雷诺应力和附加随机扰动添加在进口。（2）延长进口，这样进口可以调用一个无湍流速度。这种方法通过使流体自发发展形成湍流。但是这个上游的长度通常为水力直径的50倍以上。这种方法需要附加一个很薄的边界层。（3）从其他计算好的结果中截取剖面作为进口条件；
4. 由于湍流本身为三维的因此所有的LES和DNS均为三维的模拟（除了非常特殊的情况）。因此在计算中通常采用周期边界条件（如果平均流型是均一的）。且周期边界的距离应该大于最大的涡旋尺寸的二倍以上；