CFD中的LES大涡模拟


如果打不开图像,请右键在新标签页打开图像后刷新几次

1. 引言

湍流运动是目前计算流体力学中困难最多因此也最活跃的领域之一。当湍流存在,则住在其他相关的流动现象,并引致能量耗散、混合以及传热。没有三维的涡,则没有真正的湍流,因为只有在三维的流动中,涡旋才能进行伸展并产生新的涡旋。目前可采用的数值计算方法分为三类:直接模拟(Direct Numerical Simulation,DNS)、大涡模拟(Large Eddy Simulation,LES)和雷诺时均法(Reynolds-averaged Navier–Stokes,RANS)。RANS经过长期的发展,已经非常成熟。但RANS通过将速度进行平均后,并不能捕获湍流中的小涡结构。同时,这些小涡基本是各项同性的。另一方面,从主流中抽取能量的大涡却是各向异性,并且其和计算域的几何、边界、体积力高度关联。在使用RANS的时候,整个流场中必须使用同一个湍流模型对各种尺度下的湍流进行解析,但通常大涡和小涡的表现是不同的。因此研究学者对一种更完善的模型进行了探索。

不同于RANS,LES对大涡进行解析的同时对小涡进行模化。LES认为大涡直接受边界条件的影响因此对其解析,但小涡是各项同性的因此他们表现相同,可以进行模化。由于LES把小涡进行了模化,因此最小的网格单元需要大于Kolmogorov尺度(最小的涡旋尺度)。同时LES的时间步可以比DNS大的多。因此,对于给定的计算资源,相对于DNS,LES可以计算更大雷诺数的算例。另外,不同于RANS中平均的概念,LES使用的是一种空间滤波技术。LES模型的概念如下:

  1. 首先要确定一种滤波函数和截止尺度$\Delta$。这样,就可以对所有大于截止尺度的涡进行非稳态计算;
  2. 使用滤波函数对依时变量进行空间滤波操作,在这一步,小于截止尺度的涡被过滤掉;
  3. 在解析大涡和模化小涡的数学操作中,会产生一个亚格子尺度应力项(Sub-grid-scale Stress,SGS),亚格子尺度应力需要通过SGS模型来模化;

在LES中,截止尺度$\Delta$是用来表明“多大的涡才算大涡”的概念。其可以为任意大小,但是选择比网格还要小的截止尺度是没有意义的。在笛卡尔网格下,最简单的截止尺度这样计算:

\begin{equation} \Delta=\sqrt[3]{\delta x \delta y \delta z} \end{equation}

其中$\delta x$等为笛卡尔网格下网格单元的边长。其他不同的截止尺度计算方法还有最大边长法、普朗特混合长法等。

2. 方程与模型
2.1. 滤波N-S方程

笛卡尔坐标下的连续性方程为:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \mathbf{U})=0 \end{equation}

对$\nabla \cdot(\rho \mathbf{U})$进行滤波后有滤波连续性方程:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \mathbf{\bar{U}})=0 \end{equation}

其中$\mathbf{\bar{U}}$为滤波后的速度。且有:

\begin{equation} \mathbf{U}=\mathbf{\bar{U}}+\mathbf{U'} \end{equation}

其中$\mathbf{U'}$为残余速度,且$\mathbf{U'} \neq 0$。图(3)表示的是一个对$x$分量速度进行高斯滤波之后的滤波速度分量和残余速度分量的示意图。


残余速度
图3. 滤波速度和残余速度示意图,其中上端的细线为$x$分量速度,上端粗线为高斯滤波后的$x$分量速度。下端细线为$x$分量残余速度,下端粗线为高斯滤波后的$x$分量残余速度。其中$\Delta=0.35$

另外,有不可压缩动量方程为:

\begin{equation} \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \mathbf{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\nu \nabla \mathbf{U}) \label{mom} \end{equation}

同理有:

\begin{equation} \frac{\partial \mathbf{\bar{U}}}{\partial t}+\nabla \cdot (\overline{\mathbf{U} \mathbf{U}})=-\nabla \frac{\bar{p}}{\rho}+\nabla \cdot(\nu \nabla \mathbf{\bar{U}}) \label{momF} \end{equation}

在方程\eqref{momF}中,除了待求的$\bar{\mathbf{U}}$和$\bar{p}$外增加了一个未知量$\overline{\mathbf{U} \mathbf{U}}$。为了将问题简化,把方程\eqref{momF}的第二项进行变化:

\begin{equation} \nabla \cdot (\overline{\mathbf{U} \mathbf{U}})=\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}})+\left(\nabla \cdot (\overline{\mathbf{U} \mathbf{U}})-\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}}) \right) \label{W} \end{equation}

将方程\eqref{W}带入到方程\eqref{momF}中有:

\begin{equation} \frac{\partial \mathbf{\bar{U}}}{\partial t}+\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}})=-\nabla \frac{\bar{p}}{\rho}+\nabla \cdot(\nu \nabla \mathbf{\bar{U}}) -\left(\nabla \cdot (\overline{\mathbf{U} \mathbf{U}})-\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}}) \right) \label{momFF} \end{equation}

对比最初的的N-S方程\eqref{mom},方程\eqref{momFF}中的最后一项$-\left(\nabla \cdot (\overline{\mathbf{U} \mathbf{U}})-\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}}) \right)$为滤波操作产生的特殊项。对其展开有:

\begin{equation} \nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}})=\nabla \cdot \left[\begin{matrix} \bar{u}_1\\ \bar{u}_2\\ \bar{u}_3 \end{matrix}\right][\bar{u}_1, \bar{u}_2, \bar{u}_3]=\nabla \cdot \left[ \begin{matrix} \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{matrix} \right] \label{T1} \end{equation} \begin{equation} \nabla \cdot (\overline{\mathbf{U} \mathbf{U}})=\nabla \cdot \overline{\left[\begin{matrix} u_1\\ u_2\\ u_3 \end{matrix}\right][u_1, u_2, u_3]}=\nabla \cdot \left[ \begin{matrix} \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{matrix} \right] \label{T2} \end{equation}

把方程\eqref{T2}和\eqref{T1}代入到$-\left(\nabla \cdot (\overline{\mathbf{U} \mathbf{U}})-\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}}) \right)$有:

\begin{equation} -\left(\nabla \cdot (\overline{\mathbf{U} \mathbf{U}})-\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}}) \right) = -\nabla \cdot \left[ \begin{matrix} \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{matrix} \right] \label{T3} \end{equation}

定义一个应力$\tau$为

\begin{equation} \tau= \left[ \begin{matrix} \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{matrix} \right] \label{tau} \end{equation}

方程\eqref{tau}即为LES中的亚格子应力,其表示模化的速度分量对解析的速度分量的影响($\tau$在OpenFOAM的代码中用$\mathbf{B}$来表示)。从数学形式来看,其起源于在对非线性对流项进行滤波的过程中。Leonard将亚格子尺度应力分为3个贡献$^3$:Leonard应力、交叉应力(cross-stress)和LES雷诺应力(LES Reynolds stresses)。更进一步,LES雷诺应力部分可以分为偏应力贡献和正应力贡献。在不可压缩流体或者湍流马赫数非常小的时候,正应力贡献通常是可以忽略的。把方程\eqref{tau}带入到方程\eqref{momFF}有:

\begin{equation} \frac{\partial \mathbf{\bar{U}}}{\partial t}+\nabla \cdot (\mathbf{\bar{U}} \mathbf{\bar{U}})=-\nabla \frac{\bar{p}}{\rho}+\nabla \cdot(\nu \nabla \mathbf{\bar{U}}) -\nabla \cdot \tau \label{momFFF} \end{equation}

下一步需要对$\tau$进行模化。需要提及的是$\mathbf{\bar{U}}$是依时的,因此LES是一种非稳态的计算。并且,由于湍流本身的固有特性,LES通常为三维的(除了非常特殊的情况)。同时,还可以注意到,在截止尺度$\Delta \rightarrow 0$的时候,$\mathbf{\bar{U}} \rightarrow \mathbf{U}$,在这种情况下,所有尺度的涡旋均被解析,因此LES也就转变为了DNS。

2.2. 亚格子模型

现在需要对SGS进行模化,最简单的模型即为通过Boussinesq假定来实现,其也为其他更高级模型的本源。首先将$\tau_{ij}$表达为偏应力以及法向应力部分,且其中的法向应力部分被认为是各项同性的:

\begin{equation}\label{smagtau} \tau=\tau-\frac{1}{3}\tau\bfI+\frac{1}{3}\tau\bfI \end{equation}

其中$\tau-\frac{1}{3}\tau\bfI$表示偏应力,$\frac{1}{3}\tau\bfI$表示应力中的各向同性部分。偏应力部分和解析形变率的关系可以表示为:

\begin{equation} \tau-\frac{1}{3}\tau\bfI=-\nu_{\mathrm{SGS}}(\nabla\bar{\bfU}+\nabla\bar{\bfU}^{\mathrm{T}})+\frac{2}{3}\nu_{\mathrm{SGS}}(\nabla\cdot\bar{\bfU})\bfI \label{smagtau2} \end{equation}

其中$\nu_{\mathrm{SGS}}$为动力亚格子粘度或大涡粘度。正应力部分可以和亚格子动能$k_\mathrm{SGS}$联系起来:

\begin{equation} \frac{1}{3}\tau\bfI=\frac{2}{3}k_{\mathrm{SGS}}\bfI \end{equation}

这样,SGS可以表示为:

\begin{equation} \tau=-\nu_{\mathrm{SGS}}(\nabla\bar{\bfU}+\nabla\bar{\bfU}^{\mathrm{T}})+\frac{2}{3}\nu_{\mathrm{SGS}}(\nabla\cdot\bar{\bfU})\bfI+\frac{2}{3}k_{\mathrm{SGS}}\bfI \label{smagtau3} \end{equation}

对于不可压缩流体有:

\begin{equation} \tau_{ij}=-\nu_{\mathrm{SGS}}(\nabla\bar{\bfU}+\nabla\bar{\bfU}^{\mathrm{T}})+\frac{2}{3}k_{\mathrm{SGS}}\bfI \label{smagtauIncomp} \end{equation}

如果有$\nu_{\mathrm{SGS}}$和$k_\mathrm{SGS}$的值,则方程\eqref{momFFF}的未知变量仅为$\mathbf{\bar{U}}$,并得以封闭。在Smagorinsky亚格子模型中,$\nu_{\mathrm{SGS}}$可以表示为:

\begin{equation} \nu_{\mathrm{SGS}}=\rho C_{\mathrm{SGS}} \Delta \sqrt{k_\mathrm{SGS}} \label{Smag} \end{equation}

将方程\eqref{Smag}带入到方程\eqref{smagtauIncomp},发现未知项变为$k_\mathrm{SGS}$,其可以通过下述模型求得:

\begin{equation} \frac{1}{2}(\nabla\bar{\bfU}+\nabla\bar{\bfU}^{\mathrm{T}}):\tau + C_e\frac{k_\mathrm{SGS}^{3/2}}{\delta} = 0 \end{equation}

其可以进一步的演变为一个关于$k_\mathrm{SGS}$的二元一次方程,其最终的解为$k_{SGS}=2\frac{C_{k}}{C_{e}}\Delta^{2}\overline{S_{ij}^{2}}$(参考链接)。在开源CFD软件OpenFOAM中,$C_\mathrm{SGS}=0.094$,$C_\epsilon=1.048$。需要注意的是,$k_\mathrm{SGS}$同样可以采用传输方程的方式来求解,这种模型对应的即为一方程亚格子动能传输方程。在Smagorinsky-Lilly亚格子模型中,从普朗特混合长出发来这样表示$\nu_{\mathrm{SGS}}$:

\begin{equation} \nu_{\mathrm{SGS}}=\rho\left(C_{\mathrm{SGS}} \Delta \right)^2 \sqrt{\frac{1}{2}(\nabla\bar{\bfU}+\nabla\bar{\bfU}^{\mathrm{T}}):(\nabla\bar{\bfU}+\nabla\bar{\bfU}^{\mathrm{T}})} \label{lilly} \end{equation}

且将\eqref{smagtauIncomp}中的$\frac{2}{3}k_\mathrm{SGS}\bfI$和压力$p$结合在一起形成一个过滤压力(或者直接忽略),这样SGS也被封闭。Lilly建议$C_{\mathrm{SGS}}$的值取$0.17-0.21$。然而在某些情况下表明其会在壁面处由于过大的速度梯度而高估大涡粘度,因此有些文献建议$C_{\mathrm{SGS}}=0.1$。不同的$C_{\mathrm{SGS}}$表明了涡旋的特征是不同的因此不能用一个普适性的常数来涵盖所有算例。更高级的模型因此而生。

正是因为$C_{\mathrm{SGS}}$不是普适性的。因此在动态亚格子模型中,$C_{\mathrm{SGS}}$通过计算而来。在动态亚格子模型中,在截断尺度$\Delta$基础之上,附加定义了一个测试截断尺度$\hat{\Delta}$。并且在俩种截断尺度上对N-S方程进行滤波。Lilly通过最小二乘法来计算最终的$C_{\mathrm{SGS}}$。计算后的$C_{\mathrm{SGS}}$随时间和空间变化而变化。但是在某些情况下使用这种方法计算的大涡粘度会变为负值。更多的关于SGS模型的介绍请参考其他文献。最后需要提及的是,SGS模型的主要作用就是进行湍流波动进行耗散$^9$。因此SGS模型需要给定一个合适的耗散量。因此LES需要使用中心差分格式,因为这一类格式不会添加附加的数值耗散。

3. 初始条件和边界条件

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

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

2018.01.11大修

东岳流体 2014 - 2019
勘误、讨论、补充内容请前往CFD中文网