陈林烽^,2)江苏科技大学, 江苏镇江 212003

A RESIDUAL-BASED UNRESOLVED-SCALE FINITE ELEMENT MODELLING FOR IMPLICT LARGE EDDY SIMULATION ¹⁾

Chen Linfeng^,2)Jiangsu University of Science and Technology, Zhenjiang 212003, Jiangsu, China

通讯作者: 2)陈林烽, 副教授, 主要研究方向: 湍流大涡数值模拟. E-mail:chen.linfeng@hotmail.com

收稿日期:2020-02-26接受日期:2020-02-26网络出版日期:2020-09-18

基金资助:

1)江苏省自然科学青年基金.BK2018040999 江苏省高等学校自然科学基金.18KJB570001

Received:2020-02-26Accepted:2020-02-26Online:2020-09-18
作者简介 About authors


摘要
对应于湍流的大尺度与小尺度流场信息, 本文在有限元的框架下, 假设Navier-Stokes方程的解的形函数由大尺度和不可解尺度形函数叠加组成, 引入对应的权函数, 将Navier-Stokes方程的有限元变分形式分解为大尺度和不可解尺度系统. 根据不可解尺度系统, 构建基于Navier-Stokes大尺度方程残差的不可解尺度模型, 将其代入Navier-Stokes方程的大尺度系统, 进而数值求解大尺度系统得到Navier-Stokes方程的大尺度解. 该方法无需像传统的大涡模拟方法那样对方程的解进行过滤, 通过对形函数进行尺度分解实现解的尺度分解. 本文使用该方法的自编程序代码开展了槽道湍流的数值模拟. 通过与有限差分大涡模拟、DNS计算结果的比较, 发现在使用较少网格情况下该方法预测的平均流向速度在近壁区与DNS数据吻合, 在黏性外层略偏高; 该方法对雷诺应力预测偏低导致从流向向垂向方向上湍动能输运略偏低. 流向速度等值面图显示该方法有效捕捉到了大尺度旋涡结构; 同时在近壁区可以观察到明显的低速条带结构.
关键词： 大涡数值模拟;有限元方法;多尺度方法;不可解尺度模型;槽道湍流

Abstract
In consistence with large and small scales in turbulent flows, shape function space can be divided into resolved and unresolved scale spaces in a frame of finite element method. Introducing the same decomposition of the weighting function space, the variational formulations of Navier-Stokes equations can be divided into two systems of equations: resolved- and unresolved-scale equations. Generally, only the resolved-scale equation is computed, and the unresolved scales are modeled. Based on the unresolved-scale equations, an approximate residual-based unresolved-scale modeling is proposed in the present study. The large-scale equations are then computed by substituting the unresolved-scale modeling. The method is called residual-based large eddy simulation, in which unlike in the classical LES a filtering for Navier-Stokes equations is needed, multiscale decomposition is instead used. Numerical simulations of a turbulent channel flow are implemented with in-house codes of the residual-based large eddy simulation. The results show that, with a low number of elements, the mean streamwise velocity obtained using the present method is in agreement with the DNS data in the inner layer, and it is slightly overpredicted in the outer layer. Underprediction of the Reynolds stress by the present method causes a reduction of turbulence intensity transportation from the streamwise direction to the normal direction. Isosurfaces of the streamwise velocity reveals its capability of capturing the large-eddy structures. Meanwhile, low-speed streaks can be clearly observed in the sublayer near the wall.
Keywords：large eddy simulation;finite element method;variational multiscale method;unresolved-scale modelling;turbulent channel flow


PDF (6002KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
陈林烽. 基于Navier-Stokes方程残差的隐式大涡模拟有限元模型 ¹⁾. 力学学报[J], 2020, 52(5): 1314-1322 DOI:10.6052/0459-1879-20-055
Chen Linfeng. A RESIDUAL-BASED UNRESOLVED-SCALE FINITE ELEMENT MODELLING FOR IMPLICT LARGE EDDY SIMULATION ¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2020, 52(5): 1314-1322 DOI:10.6052/0459-1879-20-055

引言

作为湍流数值模拟方法的一种, 大涡模拟方法(large eddy simulation, LES)在未来的工业应用中具有极大的潜力. 大涡模拟方法的主体思想是针对湍流不同尺度涡的不同特点, 引入介于大涡尺度和小涡尺度之间的滤波器对Navier-Stokes方程进行过滤, 保留滤波器宽度以上的大涡特征, 对所保留的大尺度涡结构进行直接模拟, 对小尺度涡结构采用不可解尺度模型模拟^[1-2]. 相比于直接数值模拟(Direct numerical simulation, DNS), 它过滤了小尺度信息, 计算时网格数量需求更低, 计算耗时变得更短; 相比于雷诺平均方法(Reynolds average Navier-stokes, RANS), 它通过调节网格尺寸, 可以保留不同程度的大尺度信息, 因此除了平均尺度之外, 它还能计算得出不同程度的大尺度旋涡^[3-4]. 在过去的几十年里, 大涡模拟方法的研究一直是计算流体力学科研****们的重点关注问题^[5-6].

大涡模拟方法的核心问题是构建不可解尺度(不可解尺度)模型.~自大涡模拟方法提出以来, 已经有多种不可解尺度模型被提出. Smagorinsky^[7]在1963年提出了涡黏形式的不可解尺度模型, 随后Li-l-ly^[8]利用湍动能谱确定了Smagorinsky模型系数. 该模型最早应用于大气工程中的大涡模拟计算. 实际使用过程中发现该模型在近壁区耗散过大, 无法捕捉层流到湍流的转捩过程. Bardina等^[9]提出流动中可解尺度脉动向不可解尺度脉动输运的动量由可解尺度脉动中的最小尺度部分产生, 并且可解尺度脉动的最小尺度脉动速度和过滤掉的小尺度脉动速度相似, 从而提出尺度相似模式. 该模型在实际应用时湍动能耗散太小, 往往导致计算发散. Germano等^[10]提出二次过滤并假设二次过滤后的亚格子应力等于粗、细网格上的亚格子应力差, 从而得出动态不可解尺度模式. 动态确定模型系数时, 需要对流场信息进行统计平均. 对于湍流脉动存在统计均匀方法, 通常采用空间统计平均来代替系综平均. 对于没有统计均匀方向的流动, Maneveau等^[11-12]提出沿质点轨迹平均的动态确定模式系数方法. 实际计算结果表明, 对于规则几何边界处的湍流, 动态不可解尺度模式可以得到较好的计算结果. 对于不规则几何边界处的湍流, 二次过滤计算动态不可解尺度模型系数时往往会出现错误. 随后, Nicoud等^[13]提出了壁面自适应不可解尺度模型(WALE).

随着大涡模拟方法的发展, Hughes等^[14]在有限元方法的框架下将多尺度方法应用于大涡模拟方法. 在有限元的框架下, 假设Navier-Stokes方程的解的形函数由可解尺度和不可解尺度形函数叠加组成, 引入对应的权函数, 将Navier-Stokes方程的有限元变分形式分解为可解尺度和不可解尺度系统. Hughes等^[14]将可解尺度进一步分解为大尺度与小尺度, 采用传统的Smagorinsky不可解尺度模型模拟不可解尺度, 然后将不可解尺度模型仅作用于小尺度(或大尺度). 该方法无需使用过滤函数对Navier-Stokes方程的变量进行过滤, 通过对形函数进行尺度分解实现解的尺度分解. 不可解尺度模型仅作用于可解尺度中的小尺度(或大尺度), 更贴近湍动能输运规律. Hughes等^[15-16]通过各向同性湍流和槽道湍流的数值实验验证了基于有限元多尺度概念的大涡模拟方法较传统的大涡模拟方法更为准确, 证实了该方法的未来潜力.

Hughes等^[17]在将多尺度方法运用于对流扩散方程和Stokes方程时, 借助Green函数求出了不可解尺度系统的解析解. 通过将不可解尺度的解析解代入到可解尺度系统后, 可以使用较少的网格数量计算得到理想的数值模拟结果. 根据不可解尺度的解析解形式, 可以发现不可解尺度是可解尺度方程残余量的函数. 对于Navier-Stokes方程, 鉴于方程的复杂性, 无法通过不可解尺度系统求得解析形式. 在Hughes等^[18-22]的工作的启发下, 本文将根据Navier-Stokes方程的不可解尺度系统建立基于大尺度方程残差的不可解尺度模型, 然后将其代入Navier-Stokes方程的大尺度系统并直接求解之得到Navier-Stokes方程的大尺度解.

1 数值方法

本文将采用有限元方法对流体控制方程Navier-Stokes方程

(1) $ \begin{eqnarray} \label{eq1} &&\nabla \cdot { u}=0 \end{eqnarray} $
(2) $ \begin{eqnarray} \label{eq2} {\cal N}({ u},p)=\frac{\partial { u}}{\partial t}+({ u}\cdot \nabla ){ u}-\nu \nabla ^2{ u}+\nabla p={ 0} \end{eqnarray}$
进行离散求解^[23]. 令${ W}=\{{ w},q\}$为对应于Navier-Stokes变量${ U}=\{{ u},p\}$的权函数, 则Navier-Stokes方程的有限元变分形式可写为

(3) $ \begin{eqnarray} \label{eq3} B({ W},{ U})=({ W},{ F}) \end{eqnarray}$
其中

(4) $ \begin{eqnarray} \label{eq4} &&B({ W},{ U})= \lt({ w},\frac{\partial { u}}{\partial t} )_{\varOmega} +[{ w},({ u}\cdot \nabla ){ u}]_{\varOmega}+ \\&&\qquad ({ w},\nabla p)_{\varOmega} -({ w},\nu \nabla ^2{ u})_{\varOmega}+ (q,\nabla \cdot { u})_{\varOmega}\end{eqnarray}$
(5) $ \begin{eqnarray} \label{eq5} ({ W},{ F})=({ w},{ 0})_{\varOmega} +(q,0)_{\varOmega} \end{eqnarray}$
式中( , )$_{\varOmega }$表示前后两个变量乘积在计算域内部的有限元单元内的积分, $\varOmega $表示计算域内部空间.

2 多尺度方法

2.1 多尺度分解

鉴于湍流流动内存在多种尺度的旋涡, 可以将Navier-Stokes方程的解看成是大尺度解和小尺度解叠加而成(如图1), 即

(6) $ \begin{eqnarray} \label{eq6} { U}=\overline{ U} +{ U}' \end{eqnarray}$

图1

新窗口打开|下载原图ZIP|生成PPT
图1Navier-Stokes方程解的尺度分解示意图

Fig.1Sketch of scale decomposition of solution to Navier-Stokes equations



相应地, 可以将权函数分解为大尺度和小尺度

(7) $ \begin{eqnarray} \label{eq7} { W}=\overline{ W} +{ W}' \end{eqnarray}$
将式(6)和式(7)代入Navier-Stokes方程的变分形式(3)中可以得到分别投影到大尺度权函数和小 尺度权函数的两组方程^[24-27]

(8) $ \begin{eqnarray} \label{eq8} &&B(\overline { W} ,\overline { U} +{ U}')={ 0} \end{eqnarray}$
(9) $ \begin{eqnarray} \label{eq9} B({ W}',\overline { U} +{ U}')={ 0} \end{eqnarray}$
这里, 将式(8)称为大尺度系统, 将式(9)称为小尺度系统. 其中

(10) $ \begin{eqnarray} \label{eq10} && B(\overline { W} ,\overline { U} +{ U}')=B(\overline { W} ,\overline { U} ) + \\&&\qquad\left( {{\bar{{ w}}},\frac{\partial { u}'}{\partial t}}\right)_{\varOmega} -({\bar{{ w}}},\nu \nabla ^2{ u}')_{\varOmega} + ({\bar{{ w}}},\nabla p')_{\varOmega} + \\&&\qquad ({\bar{{ w}}},{\bar{{ u}}}\cdot \nabla { u}')_{\varOmega} + ({\bar{{ w}}},{ u}'\cdot \nabla {\bar{{ u}}})_{\varOmega} +({\bar{{ w}}},{ u}'\cdot \nabla { u}')_{\varOmega} + \\&&\qquad (q,\nabla \cdot { u}')_{\varOmega} ={ 0} \end{eqnarray}$
式中第二、三和四行为小尺度相关项. 通常情况, $\partial { u}'/\partial t$和$\nu \nabla ^2{ u}'$为较小值, 于是式(10)第二行的前两项可忽略^[28]. 为了避免在方程中出现$\nabla ^2{\bar{{ u}}}$, $\nabla { u}'$和$\nabla p'$项, 对式(10)中的对应项采用分布积分^[26], 于是

(11) $ \begin{eqnarray} \label{eq11} && B(\overline { W} ,\overline { U} +{ U}') =\left( {{\bar{{ w}}},\frac{\partial {\bar{{ u}}}}{\partial t}} \right)_{\varOmega} -(\nabla {\bar{{ w}}},{\bar{{ u}}}\otimes {\bar{{ u}}})_{\varOmega} + \\&&\qquad ({\bar{{ w}}},\nabla \bar{{p}})_{\varOmega} + (\nabla {\bar{{ w}}},\nu \nabla {\bar{{ u}}})_{\varOmega} +(q,\nabla \cdot {\bar{{ u}}})_{\varOmega} - \\&&\qquad (\nabla {\bar{{ w}}},p')_{\varOmega} -(\nabla {\bar{{ w}}},{\bar{{ u}}}\otimes { u}')_{\varOmega} -(\nabla {\bar{{ w}}},{ u}'\otimes {\bar{{ u}}})_{\varOmega} - \\&&\qquad(\nabla {\bar{{ w}}},{ u}'\otimes { u}')_{\varOmega} +(q,\nabla \cdot { u}')_{\varOmega} + \\&&\qquad ({\bar{{ w}}},({\bar{{ u}}}\cdot { n}){\bar{{ u}}})_{\varGamma} -({\bar{{ w}}},\nu \nabla {\bar{{ u}}}\cdot { n})_{\varGamma} ={ 0} \end{eqnarray}$
式中${\bar{{ u}}}\otimes {\bar{{ u}}}=\bar{{u}}_i \bar{{u}}_j$, 最后一行为分布积分得出的额外项, ( , )$_{\varGamma}$表示前后两个变量乘积在计算域边界处单元内的积分, $\varGamma$为计算域边界.

对比于传统大涡模拟方法的亚格子应力包含的Leonard, Corss和Reynolds应力项, 式(11)中包含有对应项^[30]:

(1) $-(\nabla {\bar{{ w}}},{\bar{{ u}}}\otimes {\bar{{ u}}})_{\varOmega} $为可解的Leonard应力项;

(2) $-(\nabla {\bar{{ w}}},{\bar{{ u}}}\otimes { u}')_{\varOmega} -(\nabla {\bar{{ w}}},{ u}'\otimes {\bar{{ u}}})_{\varOmega} $为Cross应力项;

(3) $-(\nabla {\bar{{ w}}},{ u}'\otimes { u}')_{\varOmega} $为Reynolds应力项.

值得注意的是, 在以上大尺度方程的推导过程中还没有使用任何模型.

对于大尺度方程式(11), 如果能找到$\{{ u}',p'\}$的精确表达式, 便可以通过求解式(11)得到精确的大尺度解. 式(11)与式(9)相互耦合, 于是可以试图通过式(9)去找到$\{{ u}',p'\}$表达式. 鉴于式(9)与式(11)同样复杂, 以下将通过分析式(9)的具体形式, 去建立$\{{ u}',p'\}$的近似模型.

2.2 不可解尺度模型

将式(9)展开后可以得到

(12) $ \begin{eqnarray} \label{eq12} &&B({ W}',\overline { U} +{ U}') =({ w}',{\cal N}({\bar{{ u}}},\bar{{p}}))_{\varOmega} +({ w}',{\cal L}_{{\bar{{ u}}}} { u}')_{\varOmega} + \\&&\qquad ({ w}',{ u}'\cdot \nabla ({\bar{{ u}}}+{ u}'))_{\varOmega} +({ w}',\nabla p')_{\varOmega} + \\&&\qquad (q',\nabla \cdot ({\bar{{ u}}}+{ u}'))_{\varOmega} ={ 0} \end{eqnarray}$
其中

(13) $ \begin{eqnarray} \label{eq13} &&{\cal N}({\bar{{ u}}},\bar{{p}})=\frac{\partial {\bar{{ u}}}}{\partial t}+({\bar{{ u}}}\cdot \nabla ){\bar{{ u}}}-\nu \nabla ^2{\bar{{ u}}}+\nabla \bar{{p}} \end{eqnarray}$
(14) $ \begin{eqnarray} \label{eq14} {\cal L}_{{\bar{{ u}}}} { u}'=\frac{\partial { u}'}{\partial t}+{\bar{{ u}}}\cdot \nabla { u}'-\nu \nabla ^2{ u}' \end{eqnarray}$
将式(12)的大尺度和小尺度项进行整理得到

(15) $ \begin{eqnarray} \label{eq15} && ({ w}',{\cal L}_{{\bar{{ u}}}} { u}')_{\varOmega} +[{ w}',{ u}'\cdot \nabla ({\bar{{ u}}}+{ u}')]_{\varOmega} + \\&&\qquad ({ w}',\nabla p')_{\varOmega} =-({ w}',{\cal N}({\bar{{ u}}},\bar{{p}}))_{\varOmega}\end{eqnarray}$
(16) $ \begin{eqnarray} \label{eq16} (q',\nabla \cdot { u}')_{\varOmega} =-(q',\nabla \cdot {\bar{{ u}}})_{\varOmega} \end{eqnarray}$
注意到式(15)和式(16)右边为大尺度相关项. 式中$\nabla \cdot {\bar{{ u}}}$和${\cal N}({\bar{{ u}}},\bar{{p}})$分别为连续性和Navier-Stokes大尺度方程残差, 即

(17) $ \begin{eqnarray} \label{eq17} &&r_{\rm c} =\nabla \cdot {\bar{{ u}}}\end{eqnarray} $
(18) $ \begin{eqnarray} \label{eq18} { r}_{\rm m} ={\cal N}({\bar{{ u}}},\bar{{p}}) \end{eqnarray} $
于是, 式(15)可简化为

(19) $ \begin{eqnarray} \label{eq19} &&({ w}',{\cal L}_{{\bar{{ u}}}} { u}'+{ u}'\cdot \nabla ({\bar{{ u}}}+{ u}')+{ r}_{\rm m} )_{\varOmega} + \\&&\qquad ({ w}',\nabla p')_{\varOmega} ={ 0} \end{eqnarray} $
假设$\nabla p'$为极小量, 式(15)可以再次简化为

(20) $ \begin{eqnarray} \label{eq20} ({ w}',{\cal L}_{{\bar{{ u}}}} { u}'+{ u}'\cdot \nabla ({\bar{{ u}}}+{ u}')+{ r}_{\rm m} )_{\varOmega} ={ 0} \end{eqnarray}$
其等效于

(21) $ \begin{eqnarray} \label{eq21} {\cal L}_{{\bar{{ u}}}} { u}'+{ u}'\cdot \nabla ({\bar{{ u}}}+{ u}')+{ r}_{\rm m} ={ 0} \end{eqnarray}$
忽略${ u}'\cdot \nabla ({\bar{{ u}}}+{ u}')$的影响, 小尺度${ u}'$的方程最终为

(22) $ \begin{eqnarray} \label{eq22} {\cal L}_{{\bar{{ u}}}} { u}'+{ r}_{\rm m} ={ 0} \end{eqnarray}$
如果式(14)中的${\bar{{ u}}}$为已知量, 则${\cal L}_{{\bar{{ u}}}} { u}'= { 0}$表示的是线性对流扩散方程. 对于对流扩散方程, Shakib等^[31-33]推导得到

(23) $ \begin{eqnarray} \label{eq23} {\cal L}_{{\bar{{ u}}}} { u}'\sim \tau _{\rm m}^{-1} { u}' \end{eqnarray}$
其中

(24) $ \begin{eqnarray} \label{eq24} \tau _{\rm m} =\left( {c_1 {\bar{{ u}}G\bar{{u}}}+c_2 \nu ^2{ G}:{ G}+\frac{c_3 }{\Delta t^2}} \right)^{-{1}/{2}} \end{eqnarray}$
式中$c_{1}$, $c_{2}$和$c_{3}$为常数, $ G$为单元局部坐标与全局坐标之间的转换系数矩阵

(25) $ \begin{eqnarray} \label{eq25} { G}=\left( {\frac{\partial { \xi }}{\partial { x}}} \right)^{\rm T}\frac{\partial { \xi }}{\partial { x}} \end{eqnarray}$
${ \xi }$为单元局部坐标, $ x$为全局坐标, ${ G}:{ G}$表示矩阵$ G$的双点积. 将式(22)代入式(21)得

(26) $ \begin{eqnarray} \label{eq26} \tau _{\rm m}^{-1} { u}'+{ r}_{\rm m} ={ 0} \end{eqnarray}$
于是解得

(27) $ \begin{eqnarray} \label{eq27} { u}'=-\tau _{\rm m} { r}_{\rm m} \end{eqnarray}$
由式(19)和式(16)并忽略式(19)中${ u}'\cdot \nabla ({\bar{{ u}}}+{ u}')$的影响可以得到

(28) $ \begin{eqnarray} \label{eq28} &&{\cal L}_{{\bar{{ u}}}} { u}'+\nabla p'+{ r}_{\rm m} ={ 0}\end{eqnarray}$
(29) $ \begin{eqnarray} \label{eq29} \nabla \cdot { u}'=-\nabla \cdot {\bar{{ u}}} \end{eqnarray} $
即

(30) $ \begin{eqnarray} \label{eq30} &&\tau _{\rm m}^{-1} { u}'+\nabla p'=-{ r}_{\rm m} \end{eqnarray}$
(31) $ \begin{eqnarray} \\&&\label{eq31} \nabla \cdot { u}'=-r_{\rm c} \end{eqnarray}$
引入算子${ g}=\sum\limits_{i=1}^e ({\partial \xi _i }/{\partial { x}}) $, 可将式(30)和式(31)写为

(32) $ \begin{eqnarray} \label{eq32} &&\tau _{\rm m}^{-1} { u}'+{ g}p'=-{ r}_{\rm m} \end{eqnarray}$
(33) $ \begin{eqnarray} \label{eq33} { g}\cdot { u}'=-r_{\rm c} \end{eqnarray}$
式(32)两边同时点乘${ g}\tau _{\rm m} $得到

(34) $ \begin{eqnarray} \label{eq34} { g}\cdot { u}'+{ g}\tau _{\rm m} \cdot { g}p'=-{ g}\cdot \tau _{\rm m} { r}_{\rm m} \end{eqnarray} $
将式(33)代入式(34), 在计算$p'$时假设${ r}_{\rm m} $已经满足${ r}_{\rm m} ={\bf 0}$, 于是得到$p'$的模型为

(35) $ \begin{eqnarray} \label{eq35} p'=-({ g}\tau _{\rm m} \cdot { g})^{-1}r_{\rm c} \end{eqnarray}$
令$\tau _{\rm c} =({ g}\tau _{\rm m} \cdot { g})^{-1}$, 则

(36) $ \begin{eqnarray} \label{eq36} p'=-\tau _{\rm c} r_{\rm c} \end{eqnarray}$
将式(27)和式(36)代入大尺度方程(11)可以得到

(37) $ \begin{eqnarray} \label{eq37} && B(\overline{ W},\overline{ U}+{ U}') =\left( \bar{ w},\frac{\partial \bar{ u}}{\partial t} \right)_{\varOmega} +(\bar{ w},(\bar{ u}-\tau _{\rm m}{ r}_{\rm m} )\cdot \nabla \bar{ u})_{\varOmega}+ \\&&\qquad({\bar{{ w}}},\nabla \bar{{p}})_{\varOmega} +(\nabla {\bar{{ w}}},\nu \nabla {\bar{{ u}}})_{\varOmega} +(q,\nabla \cdot {\bar{{ u}}})_{\varOmega} + \\&&\qquad(\nabla {\bar{{ w}}},\tau _{\rm c} r_{\rm c} )_{\varOmega} +({\bar{{ u}}}\cdot \nabla {\bar{{ w}}}+\nabla q,\tau _{\rm m} { r}_{\rm m} )_{\varOmega} - \\&&\qquad(\nabla {\bar{{ w}}},\tau _{\rm m} { r}_{\rm m} \otimes \tau _{\rm m} { r}_{\rm m})_{\varOmega}+ \\&&\qquad[{\bar{{ w}}},({\bar{{ u}}}\cdot { n}){\bar{{ u}}}]_{\varGamma} -({\bar{{ w}}},\nu \nabla {\bar{{ u}}}\cdot { n})_{\varGamma} ={ 0} \end{eqnarray}$
注意到, 此时大尺度方程仅包含大尺度变量, 即方程得到封闭.

3 时间推进方法

在数值模拟时, 将直接去求解大尺度方程式(37). 这里采用Jansen等^[34]提出的generalized-$\alpha $方法进行时间步推进. 以及将介绍generalized-$\alpha $方法应用在式(37)的详细细节.

将${ a}$定义为$\{{\bar{{ u}}},\bar{{p}}\}$的节点系数, ${\dot{{ a}}}$定义为${ a}$的时间倒数. 于是, 式(37)离散后可整理为

(38) $ \begin{eqnarray} \label{eq38} A\dot{{ a}}= B a \end{eqnarray}$
式中$ A$和$ B$分别为式(37)离散系统整理后${\dot{{ a}}}$和${ a}$的系数矩阵. 引入generalized-$\alpha $得到

(39) $ \begin{eqnarray} \label{eq39} &&{ R}({\dot{{ a}}}_{n+\alpha _{\rm m} } ,{ a}_{n+\alpha _{\rm f} } )={ A\dot{{a}}}_{n+\alpha _{\rm m} } -{ Ba}_{n+\alpha _{\rm f} } =0 \end{eqnarray}$
(40) $ \begin{eqnarray} \label{eq40} { a}_{n+1} ={ a}_{n+1} +\Delta t{\dot{{ a}}}_n +\gamma \Delta t({\dot{{ a}}}_{n+1} -{\dot{{ a}}}_n ) \end{eqnarray}$
(41) $ \begin{eqnarray} \label{eq41} {\dot{{ a}}}_{n+\alpha _{\rm m} } ={\dot{{ a}}}_n +\alpha _{\rm m} ({\dot{{ a}}}_{n+1} -{\dot{{ a}}}_{n+\alpha _{\rm m} } ) \end{eqnarray}$
(42) $ \begin{eqnarray} \label{eq42} { a}_{n+\alpha _{\rm f} } ={ a}_n +\alpha _{\rm f} ({ a}_{n+1} -{ a}_n \mbox{)} \end{eqnarray}$
式中

(43) $ \begin{eqnarray} \label{eq43} \alpha _{\rm m} =\frac{1}{2}\left( {\frac{3-\zeta }{1+\zeta }} \right),\ \ \alpha _{\rm f} =\frac{1}{1+\zeta } \end{eqnarray}$
其中$\zeta =0.5$以确保时间推进精度为二阶精度和无条件收敛

(44) $ \begin{eqnarray} \label{eq44} \gamma =0.5+\alpha _{\rm m} +\alpha _{\rm f} \end{eqnarray}$
在每一个时间步上采用迭代方法求解式(38).假设

(45) $ \begin{eqnarray} \label{eq45} &&{ a}_{_{n+1} }^{(0)} ={ a}_n \end{eqnarray}$
(46) $ \begin{eqnarray} \label{eq46} {\dot{{ a}}}_{_{n+1} }^{(0)} =\frac{\gamma -1}{\gamma }{\dot{{ a}}}_n \end{eqnarray}$
将式(45)和式(46)代回式(41)和式(42), 再将得到的${\dot{{ a}}}_{n+\alpha _{\rm m} }$和${ a}_{n+\alpha _{\rm f} } $代入式(39), 并对其使用牛顿迭代

(47) $ \begin{eqnarray} \label{eq47} { R}\lt({\dot{{ a}}}_{_{n+\alpha _{\rm m} } }^{(i)} ,{ a}_{_{n+\alpha _{\rm f} } }^{(i)} )+\frac{\partial { R}\lt({\dot{{ a}}}_{_{n+\alpha _{\rm m} } }^{(i)} ,{ a}_{_{n+\alpha _{\rm f} } }^{(i)} )}{\partial { a}_{{n+\alpha _{\rm f} } }^{(i)} }\Delta { a}_{_{n+\alpha _{\rm f} } }^{(i)} ={\bf0} \end{eqnarray}$
求得$\Delta { a}_{_{n+\alpha _{\rm f} } }^{(i)} $, 然后使用

(48) $ \begin{eqnarray} \label{eq48} &&{ a}_{_{n+\alpha _{\rm f} } }^{(i+1)} ={ a}_{_{n+\alpha _{\rm f} } }^{(i)} +\Delta { a}_{_{n+\alpha _{\rm f} } }^{(i)} \end{eqnarray}$
(49) $ \begin{eqnarray} \label{eq49} {\dot{{ a}}}_{_{n+\alpha _{\rm m} } }^{(i+1)} =\left( {1-\frac{\alpha _{\rm m} }{\gamma }} \right){\dot{{ a}}}_n +\frac{\alpha _{\rm m} }{\gamma \Delta t\alpha _{\rm f} }({ a}_{_{n+\alpha _{\rm f} } }^{(i+1)} -{ a}_n ) \end{eqnarray}$
更新${\dot{{ a}}}_{n+\alpha _{\rm m} } $和${ a}_{n+\alpha _{\rm f} } $. 当$\Delta { a}_{_{n+\alpha _{\rm f} } }^{(i)} $小于给定的小量, 迭代停止, 然后使用式(41)和式(42)计算出${\dot{{ a}}}_{n+1} $和${ a}_{n+1} $.

4 数值验证

4.1 算例设置

本文针对$Re_\tau =180$的槽道湍流问题, 对基于Navier-Stokes大尺度方程残差的大涡数值模拟方法及其自编程序^[35-37]做出验证. 其中$Re_\tau =u_\tau \delta /\nu $, $u_\tau $为壁面摩擦速度, $\delta $为槽道的半宽度, $\nu $为流体运动黏性系数. 流场的计算域为$L_x =2\delta \pi $, $L_y =2\delta $, $L_z =4\delta \pi /3$ (如图2所示).

图2

新窗口打开|下载原图ZIP|生成PPT
图2槽道流示意图

Fig.2Sketch of channel flow



计算域的网格划分采用六面体网格(如图3所示), 在流向($x$)与展向($z$)方向上采用均匀网格, 在垂向($y$)方向采用双向正切函数

(50) $ \begin{eqnarray} \label{eq50} y_j =0.5\left\{ {\frac{\tanh \left[ {\alpha _s (2j/N_y -1)} \right]}{\tanh (\alpha _s )}+1} \right\} \end{eqnarray}$

图3

新窗口打开|下载原图ZIP|生成PPT
图348 $\times $ 48 $\times $ 48个六面体单元组成的网格示意图

Fig.3Mesh sketch of 48 $\times $ 48 $\times $ 48 hexahedral elements



对壁面附近网格进行加密. 本文分别使用32 $\times $ 32 $\times $ 32和48 $\times $ 48 $\times$ 48个单元的两种网格对流场进行数值模拟. 网格的详细参数见表1, $\Delta y_{\min }^+ $表示壁面单元的无量纲尺寸.

Table 1
表1
表1计算参数
Table 1Computation parameters

新窗口打开|下载CSV

三维单元上的形函数以及有限元权函数均使用双线性(bilinear)函数, 即形函数由每一个方向上的线性函数相乘得到, 如$(1-\xi )(1-\psi )(1-\zeta )$. 图4为垂向上由双向正切函数得到的非均匀网格单元上的线性形函数示意图.

图4

新窗口打开|下载原图ZIP|生成PPT
图4垂向线性形函数示意图

Fig.4Sketch of linear shape functions in the normal direction



数值模拟中, 流向和展向采用周期性边界条件, 上、下固壁处使用无滑移边界条件(${ u}_{\rm wall} ={ 0})$. 流向进出口之间给定一个恒定的压力差来维持流体流动, 并在入口处给定一个压力参考值作为压力的边界条件. 对离散的式(37)使用gernerzlized-$\alpha $后, 可以得到式(47)的线性方程系统. 这里, 采用基于ILU的GMRES迭代算法并行求解式(47)的线性方程系统. 图5为使用48 $\times $ 48 $\times$ 48个单元计算得到的瞬时流场流向速度云图. 从入口的截面可以看到, 在壁面附近有很明显的展向旋涡.

图5

新窗口打开|下载原图ZIP|生成PPT
图548 $\times $ 48 $\times $ 48个单元瞬时流场流向速度云图

Fig.5Instantaneous contours of streamwise velocity using 48 $\times$ 48 $\times $ 48 elements

4.2 结果与分析

在流场稳定以后, 对20 000个时间步的流场数据进行了统计计算, 得到流场的平均速度和脉动速度均方根值.

图6为流向平均速度在垂向方向($y^{+})$的分布图, $\langle U\rangle $表示在$x$和$z$方向上的空间及时间综合平均统计量, $y^+=u_\tau y/\nu =yRe_\tau $. 图中分别给出了Jimenez的直接数值模拟和32 $\times$ 32 $\times$ 32, 48 $\times$ 48 $\times$ 48个单元的Smagorinsky模型大涡模拟有限差分^[3-4]和本文的计算结果. 可以看出, 在壁面附近($y^{+}<15$), 本文所用方法的计算结果与DNS数据匹配的很好, 而Smagorinsky模型过大的耗散导致该区域的计算结果低于DNS数据; 在对数率区($30<y^{+}<50$)和黏性外层($y^{+}>50$), 本文的计算结果略大于DNS计算结果. 该区域的平均速度随着网格数量增多变得更接近DNS数据. 根据这一结果, 可以判断该模型在32 $\times$ 32 $\times$ 32, 48 $\times $ 48 $\times$ 48个单元的计算中所提供的湍动能耗散略偏小.

图6

新窗口打开|下载原图ZIP|生成PPT
图6流向平均速度在垂向方向上的分布

Fig.6Mean streamwise velocity ($\langle U\rangle)$ profiles against $y^{+}$



图7给出了脉动速度均方根值($u_{\rm rms}$)在垂向方向的分布图,

$$\begin{eqnarray*} u_{\rm rms} =\sqrt{\sum\limits_{i=1}^{M\times N} {(u-\langle u\rangle )^2}}, \end{eqnarray*}$$

图7

新窗口打开|下载原图ZIP|生成PPT
图7脉动速度均方根值在垂向方向上的分布

Fig.7Root-mean-square of velocity fluctuations against $y^{+}$



其中$u$表示任一方向的流场速度, $M$为xz截面上的空间点数, $N$为时间采样数. 与DNS结果相比发现本文所用大涡数值模拟的流向脉动速度rms在壁面附近($y^{+}<15$) DNS结果匹配良好, 然而在对数率区($30<y^{+}<50$)和黏性外层($y^{+}>50$)则略大于DNS的结果; Smagorinsky模型的计算结果显示流向脉动速度均方根的峰值跟DNS结果有较大偏差. 在垂向与展向上, 两种大涡模拟的计算结果均小于DNS结果; 相比于Smagorinsky模型的计算结果, 本文所得到的结果要更接近DNS数据. 流向脉动速度均方根为湍动能的主要来源, 因此它的增加意味着流场整体的湍动能变大, 尽管其他两个方向上的脉动速度均方根略有减小, 即表明该模型所提供的湍动能耗散略偏小, 从而导致图6中的平均流向速度偏大.

以上结果说明要提高该模型在低网格数条件下的计算精确度, 需要改善该模型预测流向脉动速度rms的能力. 另外, 根据图7可以判断, 在较少网格情况下使用该模型时, 从流向向垂向及展向方向上的湍动能输运略偏低.

图8为统计平均得到的雷诺应力曲线分布图. 如图所示, 对于Smagorinsky模型, 由于它提供过大的湍动能耗散, 导致各个方向上的速度脉动均减小, 因此雷诺应力均减小. 与Smagorinsky模型相比, 本文计算得到的雷诺应力略大于Smagorinsky模型的计算结果, 尽管随着网格数的增加, 雷诺应力分布有所改进, 但结果还是偏小. 而雷诺应力值的偏小削弱了流向方向上的湍动能向垂向方向输运. 因此, 改善该模型对雷诺应力的预测可以有效提高脉动速度均方根的计算精确.

图8

新窗口打开|下载原图ZIP|生成PPT
图8平均雷诺应力项$\langle uv\rangle $在垂向方向上的分布

Fig.8Reynolds stress against $y^{+}$



图9为DNS数据和48 $\times $ 48 $\times$ 48网格计算得到的$8<U<10$ ($U$为流向速度)的等值面分布图. 如图所示, 发现DNS数据得到的等值面图上可以看到很多细小的流场结构, 而大涡数值模拟计算得到的流场仅存留了较大尺度的旋涡结构.

图9

新窗口打开|下载原图ZIP|生成PPT
图9流向速度等值面分布图: (a) DNS; (b) LES (48 $\times $ 48 $\times $ 48)

Fig.9Isosurfaces of streamwise velocity: (a) DNS; (b) present method with 48 $\times $ 48 $\times$ 48 elements



图10给出了48 $\times $ 48 $\times$ 48网格计算数据在壁面附近$y^{+}=5.4$处$xz$截面的流向速度云图. 从图中可以看出, 在近壁面的流场, 出现了明显的低速条带结构. 近壁面处出现低速条带结构为触发湍流拟序结构的第一号信息. 这说明, 尽管本文在模拟中使用的网格数较小, 该模型还是能有效触发近壁湍流的拟序结构.

图10

新窗口打开|下载原图ZIP|生成PPT
图10近壁面$y^{+}=5.4$处$xz$截面的流向速度云图

Fig.10Contours of streamwise velocity in the xz plane at $y^{+}=5.4$

5 结论

本文在有限元的框架下, 引入多尺度的概念, 提出了基于多尺度模式的大涡数值模拟方法. 论文里采用有限元形函数和权函数空间的分解来实现湍流中大尺度和不可解尺度信息的区分, 并推导得到大尺度系统的有限元变分方程. 本文通过对不可解尺度系统中各项的特点分析, 提出了不可解尺度的近似模式, 从而达到封闭大尺度方法的目的.

文中使用新提出的大涡模拟方法的自编程序实现了槽道湍流的并行数值计算. 数值计算结果验证了该方法的合理性以及自编程序的收敛性. 并得出, 该方法可以在使用较少的网格数可以得到与DNS结果接近的数值解, 并在近壁区可以观察到明显的低速条带结构. 另外, 数值结果证实该方法在网格数较少的情况下流向向另外两个方向的湍动能输运略偏低.

参考文献 View Option 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[1] 张兆顺, 崔桂香, 许春晓 等. 湍流理论与模拟. 北京: 清华大学出版社, 2017
[本文引用: 1]

( Zhang Zhaoshun, Cui Guixiang, Xu Chunxiao, et al. Theory and Modeling of Turbulence. Beijing: Tsinghua University Press, 2017 (in Chinese))
[本文引用: 1]

[2] 张兆顺, 崔桂香, 许春晓. 湍流大涡数值模拟的理论和应用. 北京: 清华大学出版社, 2008
[本文引用: 1]

( Zhang Zhaoshun, Cui Guixiang, Xu Chunxiao. Theory and Application of Large Eddy Simulation of Turbulent Flows. Beijing: Tsinghua University Press, 2008 (in Chinese))
[本文引用: 1]

[3] 董宇红, 陈林烽, 唐一敏. 存在热浮力的剪切湍流中颗粒近壁运动的大涡模拟
上海大学学报(自然科学版), 2009,15(6):600-606
[本文引用: 2]

( Dong Yuhong, Chen Linfeng, Tan Yimin. Large eddy simulation of particle motion in shear turbulent flow due to thermal buoyancy
Journal of Shanghai University $($Natural Science Edition$)$, 2009,15(6):600-606 (in Chinese))
[本文引用: 2]

[4] 陈林烽. 热剪切湍流中微细颗粒输运特性的大涡模拟研究. [硕士论文]
上海: 上海大学, 2011
[本文引用: 2]

( Chen Linfeng. An investigation on particle motion in the stratified turbulence by LES. [Master Thesis]
Shanghai: Shanghai University, 2011 (in Chinese))
[本文引用: 2]

[5] 吴霆, 时北极, 王士召 等. 大涡模拟的壁模型及其应用
力学学报, 2018,50(3):453-466
[本文引用: 1]

( Wu Ting, Shi Beiji, Wang Shizhao, et al. Wall-model for large-eddy simulation and its applications
Chinese Journal of Theoretical and Applied Mechanics, 2018,50(3):453-466 (in Chinese))
[本文引用: 1]

[6] 时北极, 何国威, 王士召. 基于滑移速度壁模型的复杂边界湍流大涡模拟
力学学报, 2019,51(3):754-766
[本文引用: 1]

( Shi Beiji, He Guowei, Wang Shizhao. Large-eddy simulation of flows with complex geometry by using the slip-wall model
Chinese Journal of Theoretical and Applied Mechanics, 2019,51(3):754-766 (in Chinese))
[本文引用: 1]

[7] Smagorinsky J. General circulation experiments with the primitive equations
Monthly Weather Review, 1963,91:99
[本文引用: 1]

[8] Lilly DK. The representation of small-scale turbulence in numerical experiment
// Proceedings of the IBM Scientific Computing Symposium on Environmental Sciences, 1967: 195-210
[本文引用: 1]

[9] Bardina J, Ferziger JH, Reynolds WC. Improved subgrid model for large-eddy simulation
AIAA meeting paper, 1980
[本文引用: 1]

[10] Germano M, Piomelli U, Moin P, et al. A dynamic subgrid-scale eddy viscosity model
Physics of Fluids, 1991,3(3):1760-1765
[本文引用: 1]

[11] Meneveau C, Lund TS, Cabot WH. A Lagrangian dynamic subgrid-scale model of turbulence
Journal of Fluid Mechanics, 1996,319(1):353-385
[本文引用: 1]

[12] Meneveau C, Katz J. Scale-invariance and turbulence models for large-eddy simulation
Annual Review of Fluid Mechanics, 2000,32(1):1-32
[本文引用: 1]

[13] Nicoud F, Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor
Flow Turbulence and Combustion, 1999,62(3):183-200
[本文引用: 1]

[14] Hughes TJR, Mazzei L, Jansen KE. Large eddy simulation and the variational multiscale method
Computing and Visualization in Science, 2000,3:47-59
DOIURL [本文引用: 2]

[15] Hughes TJR, Mazzei L, Oberai AA, et al. The multiscale formulation of large eddy simulation: Decay of homogeneous isotropic turbulence
Physics of Fluids, 2001,13(2):505-512
[本文引用: 1]

[16] Hughes TJR, Oberai AA, Mazzei L. Large eddy simulation of turbulent channel flows by the variational multiscale method
Physics of Fluids, 2001,13(6):1784-1799
[本文引用: 1]

[17] Hughes TJR. Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles, and the origins of stabilized methods
Computer Methods in Applied Mechanics and Engineering, 1995,127:387-401
[本文引用: 1]

[18] Hughes TJR, Sangalli G. Variational multiscale analysis: The fine-scale Green's function, projection, optimization, localization, and stabilized methods
SIAM on Journal Numerical Analysis, 2007,45:539-557
[本文引用: 1]

[19] Colomés O, Badia S, Codina R, et al. Assessment of variational multiscale models for the large eddy simulation of turbulent incompressible flows
Computer Methods in Applied Mechanics and Engineering, 2015,285:32-63


[20] Eikelder MFP, Akkerman I. Correct energy evolution of stabilized formulations: The relation between VMS, SUPG and GLS via dynamic orthogonal small--scales and isogeometric analysis. I: The convective--diffusive context
Computer Methods in Applied Mechanics and Engineering, 2018,331:259-280


[21] Eikelder MFP, Akkerman I. Correct energy evolution of stabilized formulations: The relation between VMS, SUPG and GLS via dynamic orthogonal small--scales and isogeometric analysis. II: The incompressible Navier--Stokes equations
Computer Methods in Applied Mechanics and Engineering, 2018,340:1135-1154


[22] Xu SZ, Gao BS, Hsu MC, et al. A residual-based variational multiscale method with weak imposition of boundary conditions for buoyancy-driven flows
Computer Methods in Applied Mechanics and Engineering, 2019,352:345-368
[本文引用: 1]

[23] 邵帅, 李明, 王年华 等. 基于非结构/混合网格模拟黏性流的高阶精度DDG/FV混合方法
力学学报, 2018,50(6):1470-1482
[本文引用: 1]

( Shao Shuai, Li Ming, Wang Nianhua, et al. High-order DDG/FV hybrid method for viscous flow simulation on unstructured/hybrid grids
Chinese Journal of Theoretical and Applied Mechanics, 2018,50(6):1470-1482 (in Chinese))
[本文引用: 1]

[24] Hughes TJR, Feijóo G, Mazzei L, et al. The variational multiscale method-- A paradigm for computational mechanics
Computer Methods in Applied Mechanics and Engineering, 1998,166:3-24
[本文引用: 1]

[25] 郑海标. 不可压缩Navier-Stokes方程变分多尺度方法研究. [博士论文]
西安: 西安交通大学, 2011


( Zheng Haibiao. Variational multiscale method for the incompressible Navier-Stokes equations. [PhD Thesis]
Xi'an: Xi'an Jiaotong University, 2011 (in Chinese))


[26] 薛菊峰, 尚月强. 非定常不可压Navier-Stokes方程基于欧拉格式的两水平变分多尺度方法
西南大学学报(自然科学版), 2018,40(9):84-90
[本文引用: 1]

( Xue Jufeng, Shang Yueqiang. A finite element variational multiscale method based on the backward Euler scheme for the time-dependent Navier-Stokes equations
Journal of Northwest University $($Natural Science Edition$)$, 2018,40(9):84-90 (in Chinese))
[本文引用: 1]

[27] 薛菊峰, 尚月强. 非定常不可压Navier-Stokes方程基于Crank-Nicolson格式的两水平变分多尺度方法
工程数学学报, 2019,36(4):419-430
[本文引用: 1]

( Xue Jufeng, Shang Yueqiang. A finite element variational multiscale method based on crank-nicolson scheme for the unsteady Navier-Stokes equations
Chinese Journal of Engineering Mathematics, 2019,36(4):419-430 (in Chinese))
[本文引用: 1]

[28] Akkerman I. Adaptive variational multiscale formulations using the discrete Germano approach. [PhD Thesis]
Delft: Delft University of Technology, 2009
[本文引用: 1]

[29] Chen LF, Wu GX. Boundary shear flow past a cylinder near a wall
Applied Ocean Research, 2019,92:101923


[30] Chen LF. Aspects of POD-based wall-layer modeling for the variational multiscale methods. [PhD Thesis]
Delft: Delft University of Technology, 2015
[本文引用: 1]

[31] Shakib F. Finite element analysis of the compressible Euler and Navier-Stokes equations. [PhD Thesis]
Stanford: Stanford University, 1989
[本文引用: 1]

[32] Shakib F, Hughes TJR. A new finite element formulation for computational fluid dynamics: IX. Fourier analysis of space-time Galerkin/least-squares algorithms
Computer Methods in Applied Mechanics and Engineering, 1991,87:35-58


[33] Shakib F, Hughes TJR, Johan Z. A new finite element formulation for computational fluid dynamics: X. The compressible Euler and Navier-Stokes equations
Computer Methods in Applied Mechanics and Engineering, 1991,89:141-219
[本文引用: 1]

[34] Jansen KE, Whiting CH, Hulbert GM. A generalized-$\alpha $ method for integrating the filtered Navier--Stokes equations with a stabilized finite element method
Computer Methods in Applied Mechanics and Engineering, 2000,190(3):305-319
[本文引用: 1]

[35] Chen LF, Edeling WN, Hulshoff SJ. Pod enriched boundary models and their optimal stabilisation
International Journal for Numerical Methods in Fluids, 2015,77(2): 92-107
[本文引用: 1]

[36] Chen LF, Hu XQ. A POD-Based variational multiscale method for large eddy simulation of turbulent channel flows
International Journal of Computational Methods, 2017,14:1750049


[37] Chen LF, Hulshoff SJ, Wang YT. 2D residual-based LES of flow around a pipeline close to a flat seabed
Computer Methods in Applied Mechanics and Engineering, 2020,363:112788
[本文引用: 1]