• 返回主页
  • CFD中的LES湍流模型

    2017.07.14:依据CFD中国用户一二的讨论增补内容


    湍流运动是目前计算流体力学中困难最多因此也最活跃的领域之一。目前可采用的数值计算方法可以大致分为三类:直接模拟(DNS)、大涡模拟(LES)和雷诺时均法(RANS)$^{1,5}$。RANS经过长期的发展,已经非常成熟。湍流中的小涡(Small Eddies)基本是各项同性(isotropic)并且表现相同(至少高雷诺数下的湍流如此)$^1$。另一方面,从主流中抽取能量的大涡(large Eddies)确是各向异性(anisotropic)并且和几何、边界、体积力高度关联。在使用RANS的时候,整个流场中必须使用同一个湍流模型对各种尺度下的湍流进行解析,但通常大涡和小涡的表现是不同的。因此研究学者对一种更完善的模型进行了探索。

    不同于RANS,LES对大涡进行解析的同时对小涡进行模化。LES认为大涡直接受边界条件的影响因此需要进行解析,但是小涡是各项同性的因此他们表现相同,可以进行模化。由于LES把小涡进行了模化,因此最小的网格单元需要大于Kolmogorov尺度$^1$,且时间步可以比DNS大的多。因此,对于给定的计算资源,相对于DNS,LES可以计算更大雷诺数的算例。

    另外,不同于RANS中的时均(Time-Averaging),LES使用的是一种空间滤波(Spatial-Filtering)技术。LES模型构成概念如下:

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

    滤波函数

    目前有不同的滤波函数可供使用。最简单的滤波函数为盒子滤波(box filter也即top-hat filter)。如下表示$^1$: \begin{equation} G(x-x')=\left\{\begin{matrix} 1/\Delta^3 & |x-x'|<\Delta/2\\ 0 & otherwise \end{matrix}\right. \end{equation} 其他比较出名的滤波函数还有斯坦福大学开发的高斯滤波函数、拉普拉斯滤波函数等。使用不同的滤波函数我们的流场变量有不同的型线,如图(1)和图(2)所示。
    盒子滤波
    图1. 盒子滤波示意图,其中纵坐标为$x$方向速度分量,横坐标为$x$方向距离。细黑线为滤波前的速度,粗黑线为滤波后的速度。

    高斯滤波
    图2. 高斯滤波示意图,其中纵坐标为$x$方向速度分量,横坐标为$x$方向距离。细黑线为滤波前的速度,粗黑线为滤波后的速度。

    在CFD代码中,有些LES模型需要植入滤波函数。滤波函数其具有如下特性:

    \begin{equation} \int_{-\infty }^{\infty} G(x)dx=1 \end{equation} \begin{equation} \overline{u+v}=\bar{u}+\bar{v} \end{equation} \begin{equation} \overline{\frac{\partial u}{\partial x}}=\frac{\partial \bar{u}}{\partial x} \label{partial} \end{equation}

    截止尺度

    截止尺度是用来表明“多大的涡才算大涡”的概念。其可以为任意大小,但是选择比网格还要小的截止尺度是没有意义的$^2$。在笛卡尔网格下,最简单的截止尺度这样计算: \begin{equation} \Delta=\sqrt[3]{\delta x \delta y \delta z} \end{equation} 其中$\delta x$等为笛卡尔网格下网格单元的边长。其他不同的截止尺度计算方法还有最大边长法、普朗特混合长法等。

    滤波NS方程

    假定我们有笛卡尔坐标系统的连续性方程为: \begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \mathbf{U})=0 \end{equation} 参考方程($\ref{partial}$),我们对$\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$分量速度进行高斯滤波之后的滤波速度分量和残余速度分量的示意图$^1$。
    残余速度
    图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} 参考方程($\ref{partial}$),同理我们有: \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} 在方程($\ref{momF}$)中,除了待求的$\bar{\mathbf{U}}$和$\bar{p}$外,我们增加了一个未知量$\overline{\mathbf{U} \mathbf{U}}$。为了将问题简化,我们把方程($\ref{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} 将方程($\ref{W}$)带入到方程($\ref{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} 对比我们原始的NS方程($\ref{mom}$),在方程($\ref{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} 把方程($\ref{T2}$)和($\ref{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} 方程($\ref{tau}$)定义的即为LES中的亚格子尺度应力(SGS),其表示模化的速度分量对解析的速度分量的影响($\tau$在OpenFOAM的代码中用$\mathbf{B}$来表示)。从数学形式来看,其起源于在对非线性对流项进行滤波的过程中。Leonard将亚格子尺度应力(SGS)进而分为3个贡献$^3$:Leonard应力,交叉应力(cross-stresses),LES雷诺应力(LES Reynolds stresses)。更进一步,LES雷诺应力部分可以分为偏应力贡献和正应力贡献。在不可压缩流体或者湍流马赫数非常小的时候,正应力贡献通常是可以忽略的。

    把方程($\ref{tau}$)带入到方程($\ref{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。

    Boussinesq假定

    可见,我们需要对SGS进行模化,最简单的模型即为通过Boussinesq假定来实现,其也为其他更高级模型的本源。首先我们将$\tau_{ij}$表达为偏应力以及法向应力部分,且其中的法向应力部分被认为是各项同性的: \begin{equation}\label{smagtau} \tau_{ij}=\left(\tau_{ij}-\frac{1}{3}\tau_{kk}\delta_{ij}\right)+\frac{1}{3}\tau_{kk}\delta_{ij} \end{equation} 其中$\tau_{ij}-\frac{1}{3}\tau_{kk}\delta_{ij}$表示偏应力,$\frac{1}{3}\tau_{kk}\delta_{ij}$表示应力中的各向同性部分。偏应力部分和解析形变率的关系可以表示为: \begin{equation} \tau_{ij}-\frac{1}{3}\tau_{kk}\delta_{ij}=-2\nu_{\mathrm{SGS}}\bar{S_{ij}}+\frac{2}{3}\nu_{\mathrm{SGS}}S_{ii}\delta_{ii} \label{smagtau2} \end{equation} 其中$\nu_{\mathrm{SGS}}$为动力亚格子尺度粘度或大涡粘度,$\bar{S_{ij}}$为解析形变率,其定义为$\bar{S_{ij}}=\frac{1}{2}\left(\frac{\partial \bar{u_i}}{\partial x_j}+\frac{\partial \bar{u_i}}{\partial x_i}\right)$。正应力部分可以和亚格子动能$k_\mathrm{SGS}$联系起来: \begin{equation} \frac{1}{3}\tau_{kk}\delta_{ij}=\frac{2}{3}\left(\frac{1}{2}\tau_{kk}\right)\delta_{ij}=\frac{2}{3}k_\mathrm{SGS}\delta_{ij} \end{equation} 即$k_\mathrm{SGS}=\frac{1}{2}\tau_{kk}$。 这样,SGS可以表示为: \begin{equation} \tau_{ij}=-2\nu_{\mathrm{SGS}}\bar{S_{ij}}+\frac{2}{3}\nu_{\mathrm{SGS}}S_{ii}\delta_{ii}+\frac{2}{3}k_\mathrm{SGS}\delta_{ij} \label{smagtau3} \end{equation} 对于不可压缩流体有: \begin{equation} \tau_{ij}=-2\nu_{\mathrm{SGS}}\bar{S_{ij}}+\frac{2}{3}k_\mathrm{SGS}\delta_{ij} \label{smagtauIncomp} \end{equation} 可以看出如果我们有$\nu_{\mathrm{SGS}}$和$k_\mathrm{SGS}$的具体值,则方程($\ref{momFFF}$)的未知变量仅为$\mathbf{\bar{U}}$,并得以封闭。随后我们看各种LES模型如何封闭$\nu_{\mathrm{SGS}}$和$k_\mathrm{SGS}$。

    · Smagorinsky亚格子模型

    在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},发现SGS的未知项变为$k_\mathrm{SGS}$,其可以通过下述模型求得: \begin{equation} \bar{S_{ij}}:\tau + C_e\frac{k^{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亚格子模型

    Smagorinsky-Lilly亚格子模型从普朗特混合长出发来这样表示$\nu_{\mathrm{SGS}}$: \begin{equation} \nu_{\mathrm{SGS}}=\rho\left(C_{\mathrm{SGS}} \Delta \right)^2 \sqrt{2\bar{S_{ij}} \bar{S_{ij}}} \label{lilly} \end{equation} 且将\eqref{smagtauIncomp}中的$\frac{2}{3}k_\mathrm{SGS}\delta_{ii}$和压力$p$结合在一起形成一个过滤压力(或者直接忽略),这样SGS也被封闭。Lilly建议$C_{\mathrm{SGS}}$的值取$0.17-0.21$。然而在某些情况下表明其会在壁面处由于过大的速度梯度而高估大涡粘度,因此有些文献建议$C_{\mathrm{SGS}}=0.1$。不同的$C_{\mathrm{SGS}}$表明了涡旋的特征是不同的因此不能用一个普适性的常数来涵盖所有算例。更高级的模型因此而生。

    · Von Driest Damping函数

    正因为Smagorinsky-Lilly亚格子模型对壁面的大涡粘度进行了高估,且由于湍流壁面附近的速度近乎为0因此大涡粘度也应该为0。则可以添加Von Driest Damping函数来实现。

    Dynamic SGS Model

    正是因为Smagorinsky模型中的$C_{\mathrm{SGS}}$不是普适性的。因此在动态亚格子尺度应力模型(Dynamic SGS Model)模型中,$C_{\mathrm{SGS}}$通过计算而来。在动态亚格子尺度应力模型中,在截断尺度$\Delta$基础之上,附加定义了一个测试截断尺度$\hat{\Delta}$。并且在俩种截断尺度上对NS方程进行滤波。Lilly通过最小二乘法来计算最终的$C_{\mathrm{SGS}}$$^7$。计算后的$C_{\mathrm{SGS}}$随时间和空间变化而变化。但是在某些情况下使用这种方法计算的大涡粘度会变为负值。更多的关于SGS模型的介绍请参考其他文献$^8$。

    最后需要提及的是,SGS模型的主要作用就是进行湍流波动进行耗散$^9$。因此SGS模型需要给定一个合适的耗散量。因此LES需要使用中心差分格式,因为这一类格式不会添加附加的耗散。所有的迎风格式都会对SGS模型添加附加的耗散特性。当然,针对这种限制,有些特殊的LES计算方法应运而生。在这些方法中(例如MILES)没有调用SGS模型。感兴趣的读者请参考相关文献$^9$看迎风格式如何在LES计算中引起附加耗散。

    LES的初始条件和边界条件

    由于LES固有的是求解非稳态NS方程,因此不同于RANS,LES的边界条件有一些其他特点。

    ·初始条件只会影响非依时类流动需要达到最终稳定的时间。最好在LES计算中添加一定的随机扰动。如果一个依时类流动和初始条件关联很大,那么推荐使用DNS或者实验数据来获取扰动;

    ·在壁面附近,可以使用壁面函数或者使用足够致密的网格以使得$y^+<1$;

    ·进口边界条件非常敏感,因为进口条件对下游流场有很大的影响。最简单的进口条件为指定一个平均速度型线并附加一个合适湍流强度的随机扰动。但是这种方法有失雷诺应力和空间的联系性。目前典型的做法为:(1).先使用RANS计算进口的雷诺应力,然后把这些雷诺应力和附加随机扰动添加在进口。(2).延长进口,这样进口可以调用一个无湍流速度。这种方法通过使流体自发发展形成湍流。但是这个上游的长度通常为水力直径的50倍以上$^6$。这种方法需要附加一个很薄的边界层。(3).从其他计算好的结果中截取剖面作为进口条件。

    ·出口边界条件非常简单,通常指定零法向梯度条件。

    ·由于湍流本身为三维的因此所有的LES和DNS均为三维的模拟(除了非常特殊的情况)。因此在计算中通常采用周期边界条件(如果平均流型是均一的)。且周期边界的距离应该大于最大的涡旋尺寸的二倍以上。

    参考文献
    [1]. Wilcox D.C. Turbulence modeling for CFD[M]. La Canada, CA: DCW industries, 1998.
    [2]. H.K.Ferziger J.H., Perić M. Computational methods for fluid dynamics[M]. Springer Science & Business Media, 2002.
    [3]. Leonard A. Energy cascade in large-eddy simulations of turbulent fluid flows[C]Turbulent diffusion in environmental pollution. 1974, 1: 237-248.
    [4]. Smagorinsky J. General circulation experiments with the primitive equations: I. the basic experiment[J]. Monthly weather review, 1963, 91(3): 99-164.
    [5]. 陶文铨. 数值传热学[M]. 西安交通大学出版社, 2001.
    [6]. Versteeg H.K, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method[M]. Pearson Education, 2007.
    [7]. Lilly D.K. A proposed modification of the Germano subgrid‐scale closure method[J]. Physics of Fluids A: Fluid Dynamics (1989-1993), 1992, 4(3): 633-635.
    [8]. Fureby C., Tabor G., Weller H.G. A comparative study of subgrid scale models in homogeneous isotropic turbulence[J]. Physics of Fluids (1994-present), 1997, 9(5): 1416-1429.
    [9]. Davidson L. Fluid mechanics, turbulent flow and turbulence modeling[J]. Chalmers University of Technology, Goteborg, Sweden (Nov 2011), 2011.


    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国