亚临界区圆柱绕流相干结构壁面模化混合RANS/LES 模型

季梦 尤云祥3)† 韩盼盼 邱小平 马乔 吴凯健

1) (上海交通大学,海洋工程国家重点实验室,上海 200240)

2) (上海君昱信息科技有限公司,上海 201800)

3) (上海交通大学三亚崖州湾深海科技研究院,三亚 572000)

本文发展了一种具有壁面模化大涡模拟能力的雷诺平均纳维-斯托克斯 (RANS)和大涡模拟(LES)方法的混合模型(简称WM-HRL 模型),致力于对亚临界区雷诺数钝体绕流相干结构这类复杂流动现象进行高置信度的CFD 解析模拟研究.该方法通过一个仅与当地网格空间分布尺寸有关的湍动能解析度指标参数rk 即可实现从RANS 到LES 的无缝快速转换,并且RANS/LES 混合转换区的边界位置及其各个分区(包括RANS区、LES 区及RANS/LES 混合转换区)对湍动能的解析能力均可通过两个指标参数 nrk1-Q和nrk2-Q 准则进行预先设定.通过对雷诺数Re=3900 下圆柱绕流场的系列数值模拟研究,获得了能够高置信度解析并捕捉其绕流场中三维时空瞬态发展相干结构特性的湍动能解析度指标参数 nrk1-Q和nrk2-Q 准则的组合条件.研究表明,该WM-HRL 模型不仅能够准确获取圆柱绕流场中剪切层小尺度K-H 不稳定性结构的精细谱结构,而且在同一套网格系统下通过变化湍动能解析度指标参数 nrk2-Q和nrk1-Q 准则的组合条件,还可以精细解析圆柱绕流场中两类不同回流区的长度结构特征,及其对应的圆柱尾部近壁面处V 和U 形两个平均流向速度剖面的分支结构特性.

1 引言

钝体绕流问题无论在理论上还是在工程实践中,都有着重要的研究价值.所谓亚临界区雷诺数钝体绕流,即边界层为层流状态,而尾流则为湍流状态的一类流动现象.在理论上亚临界区雷诺数钝体绕流场研究中,最具挑战性的难题当属其复杂三维时空瞬态发展的相干结构特性问题,包括上游低强度湍流、自由剪切层中相干结构的产生与发展、高强度湍流溃灭,以及随后产生的涡泄现象等[1].

对圆柱绕流问题,其亚临界区雷诺数的范围一般地定义为 400 ≤Re≤2.0×105.在这个雷诺数区圆柱绕流场中除了会出现大尺度涡泄这类相干结构外,还会出现其他一些更为复杂的三维时空瞬态发展的相干结构现象,包括剪切层小尺度Kelvin-Helmholtz 不稳定性(K-H 不稳定性)结构,在圆柱尾部会出现两类不同长度的回流区结构现象,以及在圆柱尾部近壁面区的平均流向速度剖面会出现V 形和U 形两个流动分支等[2].

研究表明,在圆柱绕流雷诺数位于亚临界区的情况,当雷诺数Re从400 增大到约1200 时,从圆柱表面分离的剪切层开始出现不稳定性现象[3],称为K-H 不稳定性.研究进一步表明,圆柱绕流场中周期性涡泄(Karman 涡街)现象发生于雷诺数Re～190时[4],但当雷诺数Re达到5000 附近时,圆柱绕流涡泄出现突然转变现象[3],其主要特征表现为泄涡结构开始变得不再具有周期性,这种现象可以维持到雷诺数Re=2.0×105.

所谓圆柱绕流场中的K-H 不稳定性,是指一条速度不连续的切变线上产生涡度集中而导致的流动不稳定性现象.Karman 涡街属于一类频率相对较低(频率记为fvs)的大尺度相干结构,而K-H不稳定性则属于一类频率相对较高(频率记为fkh)的小尺度相干结构,其主要特征表现为宽频的信号特性,且其峰值频率受雷诺数Re的影响而显著变化.在亚临界区雷诺数不大于5000 的情况,这两种不稳定性模式通常可以共存,且两者的频率近似满足[3]fkh/fvs=0.023Re0.67.

在圆柱绕流泄涡为周期性大尺度Karman 涡街的情况,利用RANS 模型通常能够获取其较为准确的涡泄频率fvs的值,以及Karman 涡街的流态结构.由于RANS 模型只能提供圆柱绕流场的时均量信息,而不能获得其三维空间瞬态发展的信息,因此对圆柱绕流大规模流动分离及剪切层小尺度K-H 不稳定性等这类复杂流动问题,RANS 模型并不适用[5].

直接数值模拟(DNS)[6-12]和大涡模拟 (LES)[6-8,13-23]以及部分平均N-S(PANS)[2,24-26]等尺度解析模拟(SRS)方法则可以弥补RANS 的这种缺陷,并已成为研究这类复杂常钝体绕流问题的主要手段.总体上,通过基于DNS,LES 及PANS 等大量CFD 数值模拟研究工作,并结合相关的模型实验[18,27]研究工作,对亚临界区雷诺数下圆柱绕流场涉及的层流分离、层流-湍流转捩、周期性涡脱落及剪切层不稳定性等复杂流动现象的形成机理及其特征等问题,有了较为深入的认识.

在雷诺数Re=3900 下的圆柱绕流是典型的亚临界区雷诺数流动,在其流动中除了存在大尺度涡泄及剪切层小尺度K-H 不稳定性这两类相干结构外,无论在模型实验中还是在CFD 数值模拟中,都发现在其流动中还会出现两类特殊的流动结构现象.其中,第一个特殊流动现象为在其绕流场中会出现两种不同长度的回流区结构,而第二个特殊流动现象是在圆柱尾部近壁面区的平均流向速度剖面会出现V和U 形两个流动分支结构.

Lourenco 和Shih[27]的实验发现,该雷诺数下圆柱绕流的回流区长度为Lr/D=1.18,并且在x1/D=1.06处的平均流向速度剖面呈V 形.然而,在同样的雷诺下,Parnaudeau等[18]的实验发现,圆柱绕流的回流区长度为Lr/D=1.51,且在x1/D=1.06处的平均流向速度剖面呈U 形.同时,两个实验所测得的不同站位处平均流向和横向速度及各向同性和异性雷诺应力剖面特性等均不相同.为此,众多学者采用CFD 数值模拟方法对此问题进行了长时间的研究.

Tremblay[8]采用DNS 和LES,Breuer[15]及Wong 和Png[19]采用LES,Pereira等[2]及Luo等[24]采用PANS,均复现了Lourenco 和Shih[27]的实验结果.Parnaudeau等[18]及Franke等[17]采用LES,Song等[11]及Ooi等[12]采用DNS,则均复现了Parnaudeau等[18]的实验结果.Tremblay[8]认为Lourenco 和Shih[27]的实验结果有误,Afgan等[20]认为Lourenco 和Shih[27]的统计周期不够进而导致结果未收敛,Kravchenko等[16]认为Lourenco 和Shih[27]的实验由于受到了外部干扰而导致剪切层过早转捩.

从目前的文献资料看,许多学者认为Parnaudeau等[18]的实验结果更具有可信度,因此大部分学者都采用Parnaudeau等[18]的实验结果作为CFD 数值模拟的基准参考数据,然而这种说法很难令人信服.考虑到从实验、DNS、LES 及PANS等诸多方面均可获得V 形和U 形结果,出现这两种结果的背后必然有其更深层次的理论和物理方面的机制,如圆柱体长径比、展向网格分辨率及湍流模型对绕流场湍动能的解析能力等.

Ma等[7]采用DNS 及LES,对圆柱体展向长度的影响进行了研究,发现当展向长度为2πD时产生V 形结果,而当展向长度为πD时产生U 形结果,并将产生两种结果的原因归于展向长度设置问题.Dong等[9]采用DNS 得出了类似结论.然而,Kravchenko等[16]基于LES 的研究发现,当展向网格分辨率为πD/8 时,产生V 形结果,而当展向网格分辨率为πD/48 时,产生U 形结果,但在保持足够密的展现网格分辨率时,如果将展向长度从πD增大至2πD时并不会产生明显区别.

Xia等[6]采用展向很疏的网格分辨率开展DNS 研究(在展向仅布置四层网格),却惊奇地复现了V 形结果.Kim[13]基于LES 研究获得了与Kravchenko等[16]一致的结论,即展向网格分辨率是导致U 形和V 形两个分支结构的主要原因.然而,Breuer[15]采用LES 的研究发现,在设置展向πD/32 及πD/64 两种网格分辨率的情况下均出现V 形结果.

对PANS 方法,其对钝体绕流场湍动能的解析能力通过一个滤波参数fk(=ku/k) 进行控制.其中,k为总的湍动能,ku为不可解湍动能.Pereira等[2]及Kořínek等[26]采用PANS 对V 形和U 形问题进行了研究.研究发现,当fk大于某个值fk0时,会产生V 形结果,而当fk小于某个值fk0时,会产生U 形结果,他们将此现象归因于湍流模型对湍动能的解析能,即具有高湍动能解析能力的湍流模型产生U 形结果,而具有低湍动能解析能力的湍流模型产生V 形结果.

为准确获取亚临界区雷诺数圆柱绕流中剪切层小尺度K-H 不稳定性结构高频成分的频率fkh,PANS 模型中的滤波参数一般地需要满足fk≤0.25,即PANS 对湍动能的解析能力至少需要达到75%[2].这意味着,对PANS 而言,为准确获取亚临界区雷诺数下圆柱绕流场中剪切层小尺度K-H 不稳定性结构高频成分的频率fkh,其所需网格量应当与LES 的网格量相当.

DNS 需要解析边界层内所有尺度的湍流,对网格分辨率的要求特别高,需要的网格量特别巨大,不适合高雷诺数钝体绕流CFD 计算.LES 直接解析大尺度的湍流,而小尺度湍流则用亚格子应力模化,虽然与DNS 相比网格量要少很多,但对高雷诺数钝体绕流的工程化应用仍很遥远[28].由此可见,对PANS 而言,其工程化应用前景亦是如此.

在诸多湍流模型中,兼顾计算精度与资源消耗的混合RANS/LES(HRL)模型已成为当今CFD领域研究与应用的热点,包括脱体涡模拟(DES)、延迟脱体涡模拟(DDES)及增强版脱体涡模拟(IDDES)等[29].下文将这三类模型统称为DES 类模型.

D'Alessandro等[30]基于DES,对不同网格分辨率及湍流模型的能力进行了研究与评估,认为V 形及U 形两种结果与网格分辨率及相应湍流模型性能关系密切.研究表明,标准SA-DES 模型在不同加密程度的网格下仅能预测V 形结果,而SA-IDDES 模型在很密网格下可预测U 形结果,在较疏网格下则可预测V 形结果.

综上所述,目前的研究仍存在诸多未解问题,主要为对V 形及U 形两个平均流向速度剖面分支结构产生的机理尚不能形成统一的认识.特别地,为何Ma等[7]的结果与其他相关文献的结果矛盾? 第二,为何Breuer[15]所用展向网格分辨率与Kravchenko等[16]相同,但获得的结果并不一致?第三,SA-DES 与SA-IDDES 模型所采用的RANS和LES 模型一样,其区别主要在于RANS 与LES之间的转换方式不同,为何网格分辨率对其结果却有如此大的影响[30]?

另一方面,根据目前所能查阅到的文献资料,利用混合RANS/LES 模型来研究亚临界区雷诺数圆柱绕流的模型主要为DES 及DDES[24,30].虽然这两类模型能够较为准确地获得与实验结果一致的一阶统计量(压力系数、流向及横向速度剖面分布等)的计算结果,但对二阶统计量(各向同性及异性雷诺应力剖面分布等)的计算结果与实验结果仍有较大差异.同时,由于这类模型对边界层中小尺度流动结构的解析能力有限,因此难以准确获取圆柱绕流剪切层小尺度K-H 不稳定性结构的精细信息,特别是其频谱结构及频率fkh的准确值.

有鉴于此,本文发展了一种壁面模化混合RANS/LES 模型(WM-HRL),该方法也属于一类HRL 方法.WM-HRL 模型与传统DES 类模型的主要不同之处在于,可实现自剪切层小尺度K-H不稳定性结构发生区域即做具有至少80%湍动能解析能力的完全LES 计算,不仅可有效地减少计算网格的数量,而且还可以有效解析剪切层小尺度K-H 不稳定性结构特征,并准确地捕捉其频谱结构及特征频率等信息.

在此基础上,本文以亚临界区雷诺数Re=3900 下的圆柱绕流问题为对象,对该WM-HRL模型的能力进行系列数值模拟和评估研究,包括对圆柱绕流场中剪切层小尺度K-H 不稳定性和两类不同长度回流区精细结构解析与捕捉能力的评估,以及对圆柱尾部x1/D=1.06 处平均流向速度剖面出现V 形和U 形两个流动分支结构形成机理的进一步认识等.

2 理论模型

考虑不可压缩流体的钝体绕流问题.设流体密度为ρ0,运动黏度系数为ν.建立直角坐标系为o-x1x2x3,其中ox3轴垂直 向上为 正,u=(u1,u2,u3) 为流体运动的速度矢量.对各类混合RANS/LES 模型,虽然其在钝体近壁面区取采用RANS进行计算,而在远离钝体近壁面的区域采用LES 进行计算,但两者均采用如下RANS 统一框架下的控制方程对流场进行计算:

其中,顶标“—”表示雷诺时均.τki为Reynolds 应力,采用如下Boussinesq 近似进行计算:

其中,k为湍动能,ω为比耗率,ε=β*kω为耗散率,S¯ki为形变率张量,vt为涡黏系数.

为封闭RANS 方程(1),需要引入相应的湍流模型,本文采用SSTk-ω 模型.基于该湍流模型的混合RANS/LES 模型,可通过修改其k方程中的色散项而建立,具体如下[31]:

其中,α1为模型系数,取值为0.31;F2为混合函数为形变率张量的模.

为保证迭代求解的稳定性,Menter等[32]加入了湍流黏度的限制性条件,对湍流动能生成项Pk进行了如下的修正:

其中,β*为SST 的模型常数,取值为0.09.

对DES 模型[29],定义为

其中,lRANS为RANS尺度,lLES为LES尺度,可分别表示为

其中,CDES为模型常数,Δ为LES 滤波宽度.

在(8)式中,lRANS和lLES与当地网格中湍流特征尺度相关.当lRANS＜lLES时,˜lhyb=lRANS,此时DES 为RANS 模式;当lRANS＞lLES时,˜lhyb=lLES,此时DES 为LES 模式.由此可见,DES 模型没有明确的RANS/LES 混合转换界面.该模型从RANS到LES 的转换完全由RANS 和LES 尺度的相对大小决定,或者说主要由网格尺度的空间分布决定.因此,在使用该模型时,通常要求流向和展向网格尺度不能同时小于边界层的厚度.

然而,当边界层内流向网格和展向网格加密至某个中间情况时,即当网格尺度小于RANS 尺度而又远大于LES 计算壁湍流所需的网格尺度时,此时DES 中的LES 尺度会提前进入边界层,导致边界层内模化的雷诺应力不足(MSD),缺失的湍流脉动又不足以被解析出来,当流向遇到一定的逆压梯度时,则产生非物理的流动分离现象,即网格诱导非物理分离(GIS),甚至发生对数律层不匹配(LLM)现象[33].

有鉴于此,Spalart等[33]提出可在混合长度尺度中引入边界层的识别函数来延迟RANS 到LES 的转换,从而防止LES 提前进入边界层,进而得到DDES 方法,具体如下:

其中,fd为一个经验性混合函数,具体形式如下:

其中,dw为网格计算点到壁面的距离,κ为冯卡门常数,取值为0.41.

函数fd的分布特点为,在近壁层的某个范围内(与网格及当地流场有关)等于0,即在此区域内混合长度=lRANS,此时DDES 模型还原为RANS模式.在此区域外的计算区域,混合长度决定于两个尺度lRANS和lLES的相对大小,此时DDES 与DES的表现是一样的.

DDES 模型虽然解决了原始DES 模型存在的MSD 及GIS 等缺陷,但其仍继承了DES 存在的其他问题,如LLM 缺陷.为此,Shur等[34]提出了增强版的DDES 模型(IDDES),但在原始IDDES定义的混合长度中,含有一个所谓的提升函数fe,这是一个纯人工构造的函数.Gritskevich等[31]指出,该函数的作用仅为增加涡黏系数,但这种人为增加涡黏系数的作用是否合理有待商榷,因此建议采用如下的简化版本:

hmax为计算单元的最大网格步长.

对DES,DDES 和IDDES 模型,当其进入LES后,在局部平衡流条件下,SST 模型k方程中的生成项与色散项相等[35],即

此外,由文献[35]可得

由(14)式和(15)式可得

由(16)式可得

由(16)式和(17)式,可得:

(19)式正好与经典Smagorinsky 亚格子涡黏模型一致.因此,DES,DDES 和IDDES 模型都是利用SST 模型k方程中产生项与色散项平衡这种极限情况下来间接等效LES 模式而实现的.

在(8)式中,滤波尺寸Δ控制LES 能否在Kolmogorov 能量谱尺度下解析尽可能多含能尺度的湍流场.在经典DES,DDES 及IDDES 模型(DES类模型)中,Δ一般取为最大网格尺寸Δmax.根据这3 类DES 类模型中“局部平衡流”的假定,即在边界层外的区域,要实现DES 类模型从RANS 模式转换为LES 模式,其中SST 湍流模型k方程中生成项与色散项需要达到平衡,此时DES 类模型相当于经典的Smagorinsky 型亚格子涡黏模型.

Breuer等[36]认为,在DES 类模型中采用最大网格尺寸并不合适,进而建议采用如下的体积立方根尺度:

其中,Δ1,Δ2和Δ3分别为3 个坐标方向的网格尺寸.

Reddy等[37]针对DDES 模型建议了一个新的网格混合形式如下:

Shur等[34]则建议采用如下的定义:

其中,Cw=0.15,hwn为壁面法向网格步长.

对两方程的SST 湍流模型,(8)式中CDES的取值可采用如下加权形式[31]:

在(23)式中,CDES,in为IDDES内层RANS分支的系数,CDES,out为IDDES 类模型外层LES 分支的系数,一般地可按下式取值:

对IDDES 模型,在边界层内一般为完全RANS模式,而其完全LES 模式一般发生在边界层外[38].由此可见,IDDES 除了可以避免MSD 及GIS 的问题外,也可以避免LLM 的问题.本文研究的亚临界雷诺数圆柱绕流问题,其剪切层小尺度K-H不稳定性发生在对数律层区内,由于RANS 模式难以准确地捕捉到这类小尺度K-H 不稳定性结构的三维时空瞬态发展流动的精细结构,因此IDDES 同样难以高置信度地解析这类非定常、非平衡流动现象的精细结构特性.

克服IDDES 模型缺陷的一种有效途径是使其具有壁面模化大涡模拟的能力.有鉴于此,本文构造一种新的混合函数fnd,使新所构造的混合RANS/LES 模型,除了具有延迟脱体涡模拟(DDES)的能力外,同时具有壁面模化大涡模拟(WMLES)的能力,将其称为WM-HRL 模型.

为此,首先引入湍动能解析度指标概念,即rk=ku/k,其中ku为未解湍动能.当rk=1 时,可得ku=k,即湍动能被完全模化,此时WM-HRL 为完全RANS 模式.当rk=0 时,即ku=0,即湍动能被完全解析,此时WN-HRL 为完全DNS 模式.对LES 来说,为准确地捕捉剪切层小尺度K-H 不稳定性结构这类特殊流动现象,其对湍动能的解析能力至少需要达到75%[2],即rk≤0.25 .

定义kc为截断波数,它表示LES 的滤波宽度,可由当地网格尺寸Δ*确定如下:

根据Kolmogorov的-5/3 谱定律,当kc位于惯性亚区时,ku可由下式进行计算:

其中,Ck(≈1.5) 为Kolmogorov 常数.由此可得:

其中,Lsgs为亚格 子滤波尺度,Ltur为湍流 含能尺度,它们可分别表示如下:

由(27)式和(28)式可得:

对当地网格尺度Δ*,一种合理的选择[39]为Δ*=hmax.此时,rk可以改写为

根据Nyquist-Shannon 样本定理,为了能够准确解析相关尺度的湍流结构,湍流含能尺度Ltur需要满足如下条件[39]:

其中,Nr为某个长度尺度结构能够被解析的网格单元数量,它与单位波长的波能够被重构所需的样本点数量一致.

对RANS 模式,由于其不能对湍流结构进行直接解析,因此Ltur≤hmax.由(30)式可得,此时rk≥1.对LES 模式,Larsson等[39]指出,湍流大尺度脉动场信息能够被准确捕捉的条件是Nr至少为2,即Ltur≥2hmax,将其代入(30)式,可得

有鉴于此,设rk1≤.本文将当地网格满足条件rk＞rk1的计算区域定义为WM-HRL 模型的完全RANS 区域,即

其中,P为计算区域Dcomp中的网 格单元,rk|P为网格单元P上的湍动能解析度指标值.

对WM-HRL 模型来讲,其第2 个关键是如何合理地确定RANS/LES 混合转化区域Dhyb.为此,设rk2为WM-HRL 模型进入完全LES 模式时人们期望的大尺度湍流之解析度指标值.据此,可以将WM-HRL 的RANS/LES 混合转换区域Dhyb定义为

对PANS 模型,rk也是衡量其对钝体绕流场大尺度湍流解析能力的关键性控制参数.Pereira等[2]对Re=3900 下圆柱绕流问题进行了系列数值模拟,研究发现,当rk＞0.5 时,PANS 模型尚不足以充分解析绕流场中大尺度涡结构的信息.同时,研究表明:只有当rk≤0.5 时,LES 才具有较好地解析大尺度涡结构信息的能力.因此,把rk2取为小于0.5 的某个值将是一种合理的选择.

当rk2=0.5 时,由(30)式可知,Ltur≥2.83hmax.再结合(31)式,可取Nr=3 .此时,Ltur≥3hmax,将其代入(30)式可得rk≤0.48 .由此可见,可将r¯k2的取值修正为r¯k2=0.48 .这意味着,当WM-HRL模型进入完全LES 模式后,其在区域Dhyb中对大尺度涡结构的解析能力至少可达52%.

至此,对本文将构造的一种新的WM-HRL 模型中如何合理地确定两个关键区域DRANS和Dhyb进行阐述,分别给出确定WM-HRL 模型进行完全RANS 模式的区域DRANS之rk1准则,以及确定WM-HRL 模型进行RANS/LES 混合转换模式的区域Dhyb之rk2准则.其中,rk1≤r¯k1且rk2≤r¯k2.后文将其分别称为rk1-Q准则和rk2-Q准则.

对如上 所述的rk1-Q和rk2-Q准则,需要利 用(30)式进行计算,其中RANS 含能尺度Ltur在进行CFD 计算之前属于未知量.因此,这两个准则无法用于在进行CFD 计算网格设置时来具体确定两个关键区域DRANS和Dhyb的边界位置.

为此,下面继续构造rk的一种仅依赖于当地网格尺寸的计算公式.Pope[40]指出,在对数律层区,Ltur与dw成正比关系,即Ltur=Cwdw.在高雷诺数的情况,Han等[41]指出,Cw可近似取为Cw≈2.5 .由此可见,(30)式可以改写为:

对(35)式定义的rk,其值有可能会出现大于1 的情况.为避免这种情况,将 (35) 式修改为如下形式:

在(36)式的形式下,所定义的湍动能解析度指标参数rk将始终不会超过1.0,即rk≤1.0 .根据(36)式可得如下结论.

首先,当rk≥rk1时,dw/hmax≤0.4(rk1)-3/2.在rk≥时,可得dw/hmax≤0.8,此时LES 将不能被真正激活.一般地,可将WM-HRL 模型为完全RANS 的区域定义为

其次,当rk2＜rk＜rk1时,0.4(rk1)-3/2＜dw/hmax＜0.4(rk2)-3/2.在rk2≤时,可得dw/hmax≥1.2,此时LES 正好被激活.一般地,可将WM-HRL 模型为RANS/LES 混合转换的区域定义为

将由(37)式确定的WM-HRL 模型为完全RANS 的区域称为nrk1-Q准则,而将由(38)式确定之WM-HRL 模型为RANS/LES 混合转换的区域称为nrk2-Q准则.由此可知,与前面所述之原rk1-Q和rk2-Q准则相比,新的nrk1-Q和nrk2-Q准则仅与当地网格尺寸相关,因此可以很容易地根据这两个准则来进行WM-HRL 模型中两个关键区域DRANS和Dhyb的确定及其相应的网格设置.

根据如上所述两个nrk1-Q和nrk2-Q准则,可构造一种新的混合函数fnd如下:

其中,fs为待定函数.同时,将混合长度修改为

由(39)式可知,在当地网格满足nrk1-Q准则的情况,即rk≥rk1此时mk=rk1,可得fnd=1 .再结合(40)式可知,此时WM-HRL 模型为完全RANS模式.另一方面,当rk＜rk1时,由(39)式可知,mk=rk＜rrk1,此时fns=1-fs.

根据nrk2-Q准则,函数fs需要满 足如下条 件:第一,当rk2＜mk＜rk1时,fs需要满足 0＜fs＜1 ;第二,当mk≤rk2时,fs需要满足fs=1 ;第三,当mk=rk1时,fs需要满足fs=0 .

能够同时满足上述3 个条件的函数可用一个关于mk的指数函数表示,具体如下:

其中,rk1＞rk2为自定义参数,且取值在0 和1之间.

下面证明,在(42)式的条件下,由(41)式构造的函数fs,满足其所需的3 条要求:首先,由(42)式可知,当mk=rk1时,α=-0.75,由(41)式可知,此时fs=0 ;其次,当rk2＜mk＜rk1时,-0.75＜α＜-0.25,由(41)式可知,此时 0＜fs＜1 ;第三,当mk≤rk2时,α=-0.25,由(41)式可知,此时fs=1.

目前,常用的一类WMLES 模型为所谓的代数壁面模化大涡模拟(Alg WMLES)[34],其基本思想是对Smagorinsky 型LES 的亚格子涡黏系数乘以一个阻尼系数Fr,即:

其中,vsgs为亚格子涡黏系数;CSMAG为模型常数,取值在0.1—0.18 之间;A为常数,一般取为25;y+为无量纲壁面距离.

由(43)式可知,当y+约大于60 时,阻尼函数Fr趋于1,此时Alg WMLES 即为完全Smagorinsky 型亚格子涡黏模型.对y+≤60 的这个近壁区域,正好为黏性底层和过渡层,由于黏性底层很薄,其范围约为y+≤10.0,此时Alg WLMES 的亚格子涡黏系数vsgs≈0.0,这相当于DNS.在过渡层内,Fr从0 开始增大到1,此时Alg WLMES 相当于准DNS.由此可见,对高雷诺数钝体绕流问题,为了捕捉黏性底层及过渡层这两个区域内足够精细的湍流信息,此时Alg WMLES 所需的网格量几乎与DNS 的网格量相当.

另一类常用的WMLES 模型为所谓的壁面应力模化大涡模拟(WRMLES)[39],该模型根据湍流边界层速度剖面的对数律来计算壁面剪切力,并输入到LES 边界网格作为边界条件.在WRMLES中,网格间距(Δ)与局部边界层厚度(δ)通常需要满足条件δ/Δ≈20—30 .这种较粗的近壁网格有可能会导致缺乏湍流应力的现象,同时基于最近邻LES 速度对壁面应力进行建模的壁面应力边界条件会增大边界层内部区域的总应力,因此可能会致壁面应力被低估或高估,进而发生LLM 问题[40].

WM-HRL 模型与Alg WMLES 和WRMLES模型均不相同,该模型在黏性底层内一般为完全RANS 解析模式,在过渡层内的某个区域内为RANS/LES 混合解析模式,在对数律层区域则为完全LES 模式.由于RANS 模式和RAN/LES 混合解析模式所需网格数量远小于DNS 和准DNS模式的网格量,而且在混合函数(39)式—(42)式的双重保护下,不仅可以避免MSD 及GIS 的问题,同时也可避免LLM 的问题.因此,该WM-HRL模型有望成为高雷诺数壁湍流三维时空瞬态发展湍流场高置信度解析的一种实用化的CFD 工具.

本文利用作者团队自研CFD 软件NUWA:FLOWUV系统,对如上所述WM-HRL 模型开发了相应的植入式程序.该自研CFD 软件系统,采用有限体积法(FVM)求解RANS 方程.其中,压力速度求解采用PISO 算法,对流项及扩散项的空间离散采用二阶中心差分格式,时间离散格式为二阶隐式格式.

3 计算模型与设置

本文对雷诺数Re=3900 下圆柱绕流场问题,采用WM-HRL 进行数值模拟与分析.其中,Re=U∞D/v,U∞为来流速度,D为圆柱直径.计算区域为矩形区域,如图1 所示.具体设置如下:圆柱底面中心位于坐标原点 (0,0,0),入口边界位于x1/D=-10 ;出口边界位于x1/D=15 ;前后两个边界分别位于x2/D=±4 ;展向高度为L3/D=π .

图1 计算区域设置Fig.1.Computational domain schematic.

边界条件设置如下:在入口边界,速度入口设置为自由来流条件,即=(U∞,0,0) ;湍动能按湍流强度I=5% 设定,其中k=3(U∞I)2/2 ;比耗率按ω=k/vt设定,其中vt/v=10.0 .在出口边界,设置为零压梯度出口,即∇p¯=0 .在两个垂直侧面及上下边界,均设置为对称边界条件.在圆球表面边界,设置为无滑移条件,即==0,湍动能设置为k=0,比耗率设置为ω=

采用ANSYS ICEM 软件进行分块网格划分,如图2 所示.单元网格为六面体,圆柱体壁面处第一层网格位于y+=0.66 .沿圆柱体展向网格的尺寸为Δ3/D=π/64,水平方向的网格平均尺寸为Δ1(Δ2)/D=0.02.自第一层起,各层网格在径向的增长率为1.07,在边界层内总计布置了45 层网格,整个计算域的网格数量为143 万.

图2 计算网格剖面Fig.2.Computational grid configuration.

对SSTk-ω 模型,在钝体近壁区通常采用所谓的壁面函数模型[42].该壁面函数模型对湍动能k、比耗率ω 及壁面剪切应力等,在黏性底层(y+＜5)均给出了明确的边界条件.在本文所构建的WM-HRL 模型中,在钝体近壁区采用的是SSTk-ω 模型,为使所构建的WM-HRL 模型在钝体近壁区发挥出实际的壁面模化作用,第一层网格一般需要设置在黏性底层的y+＜1 内.

对雷诺数Re=3900 下的圆柱绕流问题,本文通过系列数值模拟研究表明,将第一层网格设置在y+=0.66 处,可同时兼顾计算精度及网格总量控制等要求.另外,Lehmkuhl等[10]采用DNS 的研究表明,在湍流核心区Kolmogorov 尺度的平均值为D=0.02,为保证网格密度能够捕捉到最小尺度之圆柱绕流的相干结构,网格尺寸需要满足:=0.9.有鉴于此,本文将径向的增长率设置为1.07.

经统计,对雷诺数Re=3900 的情况,以圆柱体展向高度L3/D=π 为对象进行CFD 计算的主要文献如表1 所示.其中,Pereira等[2]采用的是L3/D=3.0 .表中还列出了展向网格分辨率Δ3/D以及所用总网格量等信息,并与本文所采用的相关网格信息进行比较.

表1 在雷诺数Re=3900 下圆柱绕流文献中所用计算模型与网格参数设置情况比较Table 1.Comparisons of computational models and grid parameters in references for flow around a circular cylinder at Reynolds number Re=3900.

由表1 可知,Lehmkuhl等[10]采用的展向网格分辨率最高,所用网格量也最大.Tremblay[8]和Breuer[15]采用的展向网格分辨率一致,但水平方向的网格分辨率不同.Luo等[24]采用的展向网格分辨率略低于前三篇文献的情况,但水平方向的网格分辨率高于Breuer[15].Pereira等[2]和D'Alessandro等[30]采用的展向网格分辨率最低,两者的总网格量接近.本文所采用的展向网格分辨率与Tremblay[8]和Breuer[15]的一致,但水平方向的网格分辨率均要低于这两篇文献的情况.总之,本文所用计算网格数量均少于DNS[10],LES[8,15],PANS[2,24]及DES类模型[24,30]的计算网格数量.

在数值计算中,无量纲化时间步长Δt*(=ΔtU∞/D) 取值为 6.8×10-3,库朗数CFL＜1,计算时长为70 个泄涡周期,并对后45 个泄涡周期的数据做统计平均,以获取圆柱绕流场统计量等信息,并与文献中相关实验结果及CFD 数值模拟结果进行比较分析,如表2 所示.

表2 文献中雷诺数Re=3900 下圆柱绕流场相关统计量的实验和数值结果Table 2.Experimental and numerical results for flow statistical characteristics from references for flow around a circular cylinder at Reynolds numbers Re=3900.

在表2 中,对PANS 模型,滤波控制函数fk定义为fk=ku/k.Lehmkuhl等[10]采用了两种DNS模式.其中,“Mode H”表示“高能模态”,利用该模态的DNS 得到V 形结果.而“Mode L”表示“低能模态”,利用该模态的DNS 得到U 形结果.

由表2 可知,对Pereira等[2](PANS) (fk=0.25,0.5)的情况,与表中所列实验结果[18,27]相比,分离角的计算结果均明显偏小.对Luo 等(PANS)[24](fk=0.5)的情况,分离角ϕs、阻力系数Cd及基线压力系数 -Cpb的计算结果与表中所列实验结果结果[27]相比均明显偏大,相对误差分别达到9.176%,37.76%和63.333%.对D'Alessandro 等(SA-DES)[30]的情况,阻力系数Cd的计算结果与表中所列实验结果[27]相比均明显偏大,相对误差达22.7%.

由表2 可进一步发现,在Parnaudeau等[18]和Lourenco 和Shih[27]的实验文章中未给出剪切层小尺度K-H 不稳定性的频谱特征及其无因次频率的实验结果,而在Tremblay[8],Breuer[15],Luo等[24]及D'Alessandro等[30]的CFD 计算中,也未给出剪切层小尺度K-H 不稳定性的频谱特征及其无因次频率的计算结果,但Lehmkuhl等[10]和Pereira等[2]给出了剪切层小尺度K-H 不稳定性的频谱特征及其无因次频率的计算结果.

Parnaudeau等[18]的实验测得回流区长度为Lr/D=1.51,并得到U 形流向速度剖面.Lourenco和Shih[27]的实验测得回流区长度为Lr/D=1.18,并得到V 形流向速度剖面.Lehmkuhl等[10]利用DNS 之Mode H 计算得到的回流区长度Lr/D=1.26,与Lourenco 和Shih[27]实验结果的相对误差为6.78%.而Lehmkuhl等[10]利用DNS 之Mode L计算得到回流区长度Lr/D=1.55,与Parnaudeau等[18]实验结果的相对误差为2.65%.

Tremblay[8]及Breuer[15]采用LES 均得到V形剖面,但是他们的CFD 计算所得回流区长度与Lourenco 和Shih[27]实验测得的回流区长度有较大差异,相对误差分别达到-11.86%和16.27%.Pereira等[2]采用PANS在fk=0.25,0.5 下分别得到U 形和V 形剖面.当fk=0.5 时,CFD 计算所得回流区长度与Lourenco 和Shih[27]实验测得的回流区长度接近,相对误差为-5.08%.当fk=0.25时,CFD 计算所得回流区长度与Parnaudeau等[18]实验测得的回流区长度相比偏大,相对误差达14.56%.此外,对Pereira等[2](PANS) (fk=0.5)的情况,与表中所列DNS[10]结果相比,K-H 不稳定性频率的计算结果明显偏大,相对误差达15.67%.

Luo等[24]采用PANS在fk=0.1 和0.5 下均得到V 形剖面.当fk=0.1 时,CFD 计算所得回流区长度与Lourenco 和Shih[27]实验测得的回流区长度的相对误差为7.63%.当fk=0.5 时,CFD计算所得回流区长度与Lourenco等[27]实验测得的回流区长度相比明显不符.Luo等[24]采用SSTDES 得到V 形剖面,但计算所得回流区长度与Lourenco 和Shih[27]实验测得的回流区长度相比也明显不符.

D'Alessandro等[30]采用SA-DES 也得到V 形剖面,但计算所得回流区长度与Lourenco 和Shih[27]实验测得的回流区长度相比明显不符.D'Alessandro等[30]采用SA-IDDES 及v2-f DES 均得到U 形剖面,采用SA-IDDES 计算得到的回流区长度与Parnaudeau等[18]实验测得的回流区长度的相对误差为-5.50%.但采用v2-f DES 计算得到的回流区长度与Parnaudeau等[18]实验测得的回流区长度相比偏大,相对误差高达11.13%.

对钝体绕流问题,其边界层包括黏性底层、过渡层及对数律层.对本文WM-HRL 模型来讲,这三类边界层范围的信息将是至关重要的.为此,利用图2 所示计算网格系统,首先进行RANS 数值模拟,结果表明:在雷诺数Re=3900 下圆柱绕流场的边界层厚度为y+≤180.0,黏性底层厚度为y+≤10.0,过渡层厚度为 10.0＜y+＜30.0 .

对雷诺数Re=3900 下的圆柱绕流场问题,其相干结构主要包含两类结构.其中,第一类为与尾流中涡泄相关的大尺度不稳定性结构,其无因次特征频率为,而第二类为与分离剪切层脉动相关的小尺度不稳定性结构,称为K-H 不稳定性,其无因次特征频率为.对本文WM-HRL 模型来讲,圆柱绕流剪切层小尺度K-H 不稳定性结构发生的位置信息将是至关重要的.

剪切层的速度梯度较大,属于一个位置狭小的窄带,且对监测点的分布位置敏感[23].为获取圆柱绕流剪切层小尺度K-H 不稳定性结构发生的准确边界位置,结合相关文献中DNS[10],PANS[2]及LES[23]等CFD 数值模拟结果,本文设置了如图3所示的18 个监测点,其坐标位置如表3 所示.其中,监测点P1—P14 位于剪切层区域,监测点P15 位于剪切层脱落区,而监测点P16—P18 位于尾流中心线上.表3 中所列的这些监测点均位于圆柱体底部位置.同时,对表3 中所列的每个监测点,沿圆柱体展向同时均匀布置了其他4 个相应的监测点.因此,实际上总计布置了90 个监测点.

表3 监测点坐标信息Table 3.Coordinate information of the probes.

图3 剪切层小尺度K-H 不稳定性结构监测点分布Fig.3.Location configuration of the probes for the small scale K-H instability structure in the shear layer.

在此基础上,利用本文WM-HRL 进行系列数值模拟,对所得各监测点处的横向脉动速度进行Lomb 频谱分析,详细分析见4.3 节.在做Lomb谱分析时,用的是表3 中所列各监测点相对应的5 个监测点处横向脉动速度做展向平均所得,其做法与Lehmkuhl等[10](DNS)的做法相同.

结果表明,在各个监测点的Lomb 频谱中,均出现一个频率相对较低的谱峰,该谱峰所对应频率与涡泄频率一致.同时,对各种工况的系列CFD数值模拟结果中,在测点P4 的Lomb 频谱中,均观察到一个频率相对较高的谱峰,该谱峰所对应频率与剪切层小尺度K-H 不稳定性结构的频率一致.由此可见,雷诺数Re=3900 下圆柱绕流场中剪切层小尺度K-H 不稳定性结构发生在对数律层中y+≥89.4的区域.

在图2 所述网格系统设置下,利用(36)式进行计算,结果表明:在对数律层的y+≥67.4 的区域中,rk均不大于0.2,即rk≤0.2 .由于圆柱绕流剪切层小尺度K-H 不稳定性结构发生在y+≥89.4的区域,因此本文WM-HRL 模型对该剪切层小尺度K-H 不稳定结构为完全LES 解析模式,且其对湍动能的解析度至少为80%.

4 数值结果与分析

当利用WM-HRL 进行数值模拟时,影响圆柱绕流场计算置信度的主要因素包括两个边界位置和3 个区域的湍动能解析度.其中,两个边界分别为RANS 结束边界(记为ΓRANS)和LES 启动 边界(记为ΓLES),3 个区域分别为RANS 区、RANS/LES 混合转换区和LES 区.为下文陈述简便计,记为RANS 结束边界ΓRANS的位置,而记为LES 启动边界ΓLES的位置.

4.1 WM-HRL 模型参数影响规律

首先讨论LES 启动边界ΓLES位于剪切层小尺度K-H 不稳定性结构发生区域的情况,并考虑两种情况:第一种情况为=105.8,rk2=0.1556 ;第 二种情况为=113.9,rk2=0.1484 .在每种情况下,均设置11 种RANS 结束边界ΓRANS的位置及其相应的rk1,如表4 所示.其中,=7.9位于黏性底层,=13.3,14.9,18.4,20.4,27.1 位于过渡层,=38.4,41.7,49.2,72.7 位于对数律层但在剪切层K-H 不稳定区外,而=91.2则位于剪切层K-H 不稳定性结构发生区域内.

表4 当 ΓLES 位于剪切层小尺度K-H 不稳定性结构发生区域内时,相关流场统计量的数值结果Table 4.Numerical results for flow statistic characteristics when ΓLES is located in the K-H instability region of the shear layer.

下面讨论LES 启动边界ΓLES位于剪切层小尺度K-H 不稳定性结构发生区域外,但仍位于对数律层区的情况,并考虑3 种情况:第一种情况为=49.2,rk2=0.2546 ;第二种 情况为=72.4,rk2=0.1983 ;第三种情况为=84.6,rk2=0.1713.在每种情况下,均设置若干种RANS 结束边界ΓRANS的位置及其相应的rk1,其设置原则与表4 一致,如表5 所示.在各种及rk1和及rk2的组合下,利用WM-HRL 进行数值模拟所得相关流场统计量的结果如表5 所示.

表5 当 ΓLES 位于剪切层小尺度K-H 不稳定性结构发生区域外且在对数律层内时,相关流场统计量的数值结果Table 5.Numerical results for flow statistic characteristics when ΓLES is located in the log-law layer and outside the K-H instability region of the shear layer.

由表5 可知,只有当RANS 结束边界ΓRANS位于过渡层,且rk1≤=0.63 时,才能同时准确获得和Lr/D的计算结果.在此条件下,当计算得到V 形速度剖面时,Lr/D的计算结果与Lourenco和Shih[27]的实验结果相比,相对误差最大为-9.32%,而的计算结果与Lehmkuhl等[10]的DNS 之Mode H 计算结果相比,相对误 差最大为12.68%.当计算得到U 形速度剖面时,Lr/D的计算结果与Parnaudeau等[18]的实验结果相比,相对误差最大为-12.58%,而的计算结果与Lehmkuhl等[10]的DNS 之Mode H 计算结果相比,相对误差最大为-20.15%.

最后讨论LES 启动边界ΓLES位于过渡层的情况,并考虑5 种情况,分别为=29.6 .在每种情况下,均 设置若干种RANS 结束边界ΓRANS的位置及其相应的rk1.在各种的组合下,利用WM-HRL 进行数值模拟所得相关流场统计量的结果如表6 所示.

表6 当 ΓLES 位于过渡层时,相关流场统计量的数值结果Table 6.Numerical results for flow statistic characteristics when ΓLES is located in the buffer layer.

需要指出的是,在上文统计中,对ϕs之CFD计算值,已去掉了Pereira等[2](PANS) (fk=0.25和0.5)的异常计算结果,以及Luo等[24](PANS)(fk=0.5)的异常计算结果.对Cd之CFD 计算值,已去掉Luo等[24](PANS)(fk=0.5)的异常计算结果,以及D'Alessandro等[30](SA-DES)的异常计算结果.对-Cpb之CFD 计算值,已去掉Luo等[24](PANS)(fk=0.5)的异常计算结果.

由于剪切层小尺度K-H 不稳定性发生的位置是未知的,但位于对数律层区中.系列数值模拟结果表明,对本文所构造的WM-HRL 模型,为了能够准确解析并捕捉剪切层小尺度K-H 不稳定性的结构特征及其频谱特性,可将其LES 模式的启动边界均设置在过渡层内,而RANS 模式的结束边界则既可设置在过渡层内也可设置在黏性底层内,并且使rk2≤0.2,即在对数律层区的网格具有至少80%的湍动能解析度.

在圆柱体展向长度及其网格系统已经确定的条件下,对DNS[6-12]及LES[6-8,13-23]来讲,通过CFD计算只能得到回流区长度及流向速度剖面形状的其中一种分支结构.对DES 类模型[24,30]来讲,不同RANS 湍流模型的类型及RANS/LES 之间不同的混合转换方式等,均可能会影响其对回流区长度及流向速度剖面形状分支结构的计算结果.

对PANS 模型[2,24-26]而言,其核心思想是通过引入一个fk参数来调控流场可解/不可解湍流量的比例而实现.对基于SSTk-ω 的PANS 模型,fk参数可调控其ω 方程中耗散项之值的变化.一般地,fk值越小,耗散项的值也越小,从而使求解得到的比耗率ω 增大,进而使湍流涡黏系数vt减小,可解尺度就释放的越多[25].

对PANS 模型,Pereira等[2]指出:为准确获取圆柱绕流回流区结构及其长度信息,fk的值不能大于0.5,即fk≤0.5 .同时,当fk从0.5 减小到某个值(文中为0.25)后,圆柱尾部近壁面处流向速度剖面的形状从V 形转变为U 形.但Luo等[24]在fk=0.5 时通过CFD 计算并没有得到准确的回流区长度,而且在fk=0.1 时,得到的是V 形剖面.导致这种相矛盾结果的原因之一,可能与两位作者所使用网格结构及其空间分辨率分布不同有关(具体见表2).

本文采用WM-HRL 模型的系列数值模拟结果表明,在同一套网格系统下,通过改变WM-HRL模型中两个边界位置及3 个区域的湍动能解析度的组合,既可以通过数值模拟获得与Parnaudeau等[18]实验结果一致的回流区长度及圆柱近壁面处平均流向速度的U 形剖面,也可以通过数值模拟获得与Lourenco 和Shih[27]实验结果一致的回流区长度及圆柱近壁面处平均流向速度的V 形剖面.

这表明对本文所建立的WM-HRL 模型,可以在同一套网格系统下通过变化湍动能解析度指标参数nrk1-Q和nrk2-Q准则的组合条件,精细地解析圆柱绕流场中两类不同回流区长度结构特征,及其对应的圆柱尾部近壁面处U 形和V 形两个流向速度剖面的分支结构.

4.2 一阶和二阶统计量特性

选取表6 中相关数值模拟结果案列,讨论雷诺数Re=3900 下圆柱绕流场的一阶和二阶统计量特性的计算结果,并与Parnaudeau等[18]及Lourenco和Shih[27]的相关实验结果进行比较分析.对U形流向速度剖面的情况,选取4 个组合工况如下.Case AU:=18.4;Case BU:=7.9,=29.6.其中,对Case AU 工况,其RANS 结束边界位于黏性底层和过渡层交界处,而LES 启动边界位于过渡层内;对Case BU 工况,其RANS 结束边界和LES 启动边界均位于过度层内;对Case CU 工况,其RANS 结束边界位于黏性底层内,而LES启动边界位于过渡层内;对Case DU 工况,其RANS 结束边界位于过渡层内,而LES 启动边界位于过渡层与对数律层交界处.

对V 形流向速度剖面的情况,选取4 个组合工况如下.Case AV:=13.3;Case BV:=14.9,=29.6.其中,对Case AV 工况,其RANS 结束边界位于黏性底层内,而LES 启动边界位于黏性底层与过渡层交界处;对Case BV 工况,其RANS 结束边界位于黏性底层和过渡层交界处,而LES 启动边界位于过渡层内;对Case CV 工况,其RANS结束边界和LES 启动边界均位于过渡层内;对Case DV 工况,其RANS 结束边界位于过渡层内,而LES 启动边界位于过渡层与对数律层交界处.

在针对一阶和二阶统计量的分析中,对Case AU—Case DU 这4 种工况,除了将其与Parnaudeau等[18]的实验结果进行比较外,同时还与Lehmkuhl等[10]利用Mode H 之DNS 计算结果进行了比较.对Case AV—Case DV 这4 种工况,除了将其与Lourenco 和Shih[27]的实验进行比较外,同时还与Lehmkuhl等[10]利用Mode L 之DNS 计算结果进行了比较.

4.2.1 一阶统计量特性

在图4 中,分别给出了上述8 种工况下圆柱表面周向系数Cp分布特性的数值结果.由于Parnaudeau等[18]和Lourenco 和Shih[27]的实验没有给出圆柱表面压力系数的数据,此处采用Norberg[43]的实验结果进行比较分析.

图4 圆柱表面周向压力系数 Cp 分布特性Fig.4.Azimuthal distribution characteristics for pressure coefficient along the circular cylinder surface.

由图4 可知,无论是对Case AU—Case DU这4 种工况,还是对Case AV—Case DV 这4 种工况,利用本文WM-HRL 模型所得沿圆柱表面周向压力系数Cp分布的数值结果之间的差异均很小.对图4(a)的情况,本文计算结果与Norberg[43]实验结果相比的相对误差最大为10.3%,而Lehmkuhl等[10]利用DNS(Mode L)计算结果与Norberg[43]实验结果相比的相对误差最大为3.1%.对图4(b)的情况,本文计算结果与Norberg[43]实验结果相比的相对误差最大为14.0%,而Lehmkuhl等[10]利用DNS(Mode H)计算结果Norberg[43]实验结果相比的相对误差最大为0.6%.

在图5 中,分别给出了前述8 种工况下沿尾流中心线平均流向速度剖面特性的数值结果.由图5可知,在圆柱表面正后方中心线上的点P(0.5D,0)处,平均流向速度的值为0,在达到一个负的最小值后开始增大,并在中心线上的点Q(Lr+0.5D,0)处又变为0.其中,Lr/D为回流区长度.Parnaudeau等[18]的实验测量获得Lr/D=1.51,Lehmkuhl等[10]利用Mode L 之DNS 计算得到Lr/D=1.55,如图5(a)所示.而Lourenco 和Shih[27]的实验测量获得Lr/D=1.18,Lehmkuhl等[10]利用Mode H之DNS 计算得到Lr/D=1.26,如图5(b)所示.此即在雷诺数Re=3900 下圆柱绕流两个实验及DNS 计算中出现两类不同回流区分支结构的情况.

图5 沿尾流中心线平均流向速度剖面特性Fig.5.Distribution characteristics of mean stream-wise velocities along the wake centerline.

由图5 进一步可发现,无论是对Case AU—Case DU 这4 种工况,还是对Case AV—Case DV这4 种工况,利用本文WM-HRL 模型所得沿尾流中心线平均流向速度分布的数值结果之间的差异均很小.对图5(a)的情况,在回流区外,本文数值结果与Parnaudeau等[18]的实验结果之间的相对误差最大为5.5%,而与Lehmkuhl等[10]利用Mode L 之DNS 的计算结果之间的相对误差最大为3.2%.

对图5(b)的情况,本文数值结果与Lehmkuhl等[10]利用Mode H 之DNS 的计算均吻合良好,两者之间的相对误差最大为4.2%.但无论是本文计算结果还是Lehmkuhl等[10]利用Mode H 之DNS的计算结果,在x1/D∈[2.36,3.26]的范围内,它们与Lourenco 和Shih[27]的实验结果均有一定的差异.

图6 中,分别给出了前述8 种工况下在3 个不同站位(x1/D=1.06,1.54,2.02)处平均流向速度剖面特性的数值结果.站位x1/D=1.06 位于剪切层K-H 不稳定性结构发生区域,同时也位于回流区内.Parnaudeau等[18]的实验测量获得U形速度剖面,而Lehmkuhl等[10]利用Mode L 之DNS 计算也得到U 形剖面,如图6(a)所示.而Lourenco 和Shih[27]的实验测量获得V 形速度剖面,而Lehmkuhl等[10]利用Mode H 之DNS 计算也得到V 形剖面,如图6(b)所示.此即在雷诺数Re=3900 下,圆柱绕流两个试验中所测量得到的平均流向速度剖面出现U 形和V 形两个流动分支结构的情况.

图6 圆柱后方不同站位处平均流向速度剖面特性Fig.6.Distribution characteristics of mean stream-wise velocities at different locations in the backside of the circular cylinder.

由图6 可知,无论是对Case AU—Case DU 这4 种工况,还是对Case AV—Case DV 这4 种工况,利用本文WM-HRL 模型所得该站位处平均流向速度分布的数值结果之间的差异均很小.在位于剪切层K-H 不稳定性结构发生区域的站位x1/D=1.06 处,对Case AU—Case DU 这4 种工况,本文计算所得平均流向速度剖面均为U 形,与Parnaudeau等[18]的实验及Lehmkuhl等[10]利用Mode L 之DNS 计算所得速度剖面形状均一致.对Case AV—Case DV 这4 种工况,本文计算所得平均流向速度剖面均为V 形,与Lourenco 和Shih[27]的实验及Lehmkuhl等[10]利用Mode H 之DNS 计算所得速度剖面形状均一致.

站位x1/D=1.54 和2.02 位于回流区中.由图6可见,Parnaudeau等[18]和Lourenco 和Shih[27]的实验结果均获得V 形流向速度剖面,而Lehmkuhl等[10]利用Mode L 和H 之DNS 计算也均获得V形流向速度剖面.同时,由图6 还可观察到,对Case AU-Case DU 这4 种工况,利用本文WMHRL 模型计算所得平均流向速度剖面均为V 形,与Parnaudeau等[18]的实验及Lehmkuhl等[10]利用Mode L 之DNS 计算所得速度剖面形状一致.对Case AV—Case DV 这4 种工况,利用本文WMHRL 模型计算所得平均流向速度剖面也均为V 形,与Lourenco 和Shih[27]的实验及Lehmkuhl等[10]利用Mode H 之DNS 计算所得速度剖面形状均一致.

图7 分别给出了前述8 种工况下在3 个不同站位(x1/D=1.06,1.54,2.02)处平均横向速度剖面特性的数值结果.

图7 圆柱后方不同站位处平均横向速度剖面特性Fig.7.Distribution characteristics of mean cross-flow velocities at different locations in the backside of the circular cylinder.

由图7 可见,无论是对Case AU—Case DU 这4 种工况,还是对Case AV—Case DV 这4 种工况,利用本文WM-HRL 模型所得圆柱后方不同站位处平均横向速度分布的数值结果差异均很小.

对图7(a)的情况,在x1/D=1.06 站位处,本文计算结果与Parnaudeau等[18]实验结果之间的相对误差最大为3.5%,而与Lehmkuhl等[10]利用Mode L 的DNS 计算结果之间的相对误差最大为3.8%.在x1/D=1.54 站位处,本文计算结果与Parnaudeau等[18]实验结果之间的相对误差最大为9.9%,而与Lehmkuhl等[10]利用Mode L 的DNS计算结果之间的相对误差最大为10.3%.在x1/D=2.02 站位处,本文计算结果与Parnaudeau等[18]实验结果之间的相对误差最大为15%,而与Lehmkuhl等[10]利用Mode L 之DNS 计算结果之间的相对误差最大为16.1%.

对图7(b)的情况,本文计算结果与Lehmkuhl等[10]利用Mode H 之DNS 计算结果一致,但两者与Lourenco 和Shih[27]的实验结果均有一定的差异.在x1/D=1.06 站位处,本文计算结果与Lehmkuhl等[10]利用Mode H 的DNS 计算结果之间的相对误差最大为7.8%.在x1/D=1.54 站位处,本文计算结果与Lehmkuhl等[10]利用Mode H 之DNS 计算结果之间的相对误差最大为5.4%.在x1/D=2.02 站位处,本文计算结果与Lehmkuhl等[10]利用Mode H 之DNS 计算结果之间的相对误差最大为10.0%.

Tremblay[8]采用DNS 和LES 均复现了Lourenco 和Shih[27]对雷诺数Re=3900 下,圆柱绕流在其后方不同站位处平均流向及横向速度剖面特性的实验结果.其中,LES 亚格子模型分别采用了经典Smagorinsky 模型(Smag)及动态模型(Dyn).结果表明,DNS的结果要优 于LES的结果,而Dyn 亚格子模型的结果则要优于Smag 亚格子模型的结果.

对PANS 模型,其对圆柱后方各站位处平均流向及横向速度剖面的计算精度依赖于滤波控制参数fk的取值方式.对DES 类模型,其对圆柱后方各站位处平均流向及横向速度剖面的计算精度则依赖于所用RANS 模型的类型及RANS/LES的混合转换方法等.总体上,在网格空间分辨率分布及相关的模型参数设置合理的条件下,这两类模型对圆柱后方各站位处平均流向及横向速度剖面的计算精度,一般地可以达到与LES 相当的计算精度[2,24-26,30].

本文采用WM-HRL 模型的系列数值模拟结果表明,在同一套网格系统下,通过改变WM-HRL模型中两个边界位置及3 个区域的湍动能解析度的组合,既可以通过数值模拟获得与Parnaudeau等[18]实验所得一致的圆柱后方各站位处的平均流向及横向速度剖面,也可以通过数值模拟获得与Lourenco 和Shih[27]实验所得一致的圆柱后方各站位处的平均流向及横向速度剖面.

4.2.2 二阶统计量特性

图8 分别给出了前述8 种工况下在3 个不同站位(x1/D=1.06,1.54,2.02)处各向同性流向雷诺应力剖面特性的数值结果.

图8 圆柱后方不同站位处各向同性流向雷诺应力剖面特性Fig.8.Distribution characteristics of isotropic stream-wise Reynolds stresses at different locations in the backside of the circular cylinder.

由图8 可见,对Case AU—Case DU 的情况,在位于剪切层小尺度K-H 不稳定性结构发生区域的站位x1/D=1.06 处,本文计算结果与Parnaudeau等[18]的实验结果及Lehmkuhl等[10]利用Mode L 之DNS 计算结果相比,主要差异在两个“峰”处.在位于回流区的站位x1/D=1.54,2.02处,本文计算结果与Parnaudeau等[18]的实验结果及Lehmkuhl等[10]利用Mode L 的DNS 计算结果相比,在两个“峰”处的值小于相应的实验值及DNS 计算值.

对Case AV—Case DV 的情况,本文计算结果均与Lourenco 和Shih[27]实验结果及Lehmkuhl等[10]利用Mode H 之DNS 计算结果具有良好的一致性.在位于剪切层小尺度K-H 不稳定区的站位x1/D=1.06 处,本文对BV 和DV 工况的计算结果与Lourenco 和Shih[27]实验结果之间的最大相对误差分别为1.7%和5%,而与Lehmkuhl等[10]利用Mode H 之DNS 计算结果之间的最大误差分别为0.3%和3.7%.在x1/D=1.54 站位处,本文计算结果与Lourenco 和Shih[27]实验结果之间的最大相对误差分别为16.1%,3.4%,3.8%,6.4%.在x1/D=2.02 站位处,本文计算结果与Lourenco和Shih[27]实验结果之间的最大相对误差分别为4.6%,11%,3.5%,6.7%.

图9 分别给出了前述8 种工况下在3 个不同站位(x1/D=1.06,1.54,2.02)处各向同性横向雷诺应力剖面特性的数值结果.

图9 圆柱后方不同站位处各向同性横向雷诺应力剖面特性Fig.9.Distribution characteristics of isotropic cross-flow Reynolds stresses at different locations in the backside of the circular cylinder.

由图9 可见,对Case AU—Case DU 的情况,在位于剪切层小尺度K-H 不稳定区的站位x1/D=1.06 处,本文计算结果与Parnaudeau等[18]实验结果相比的最大相对误差分别为23.6%,7.7%,15.7%,9.1%,而Lehmkuhl等[10]利用Mode L 之DNS 计算结果与Parnaudeau等[18]实验结果相比的最大相对误差为14.5%.

在位于回流区的站位x1/D=1.54 处,本文计算结果与Parnaudeau等[18]实验结果相比的最大相对误差分别为13.7%,25.3%,6.4%,13.6%,而Lehmkuhl等[10]利用Mode L 之DNS 计算结果与Parnaudeau等[18]实验结果相比的最大相对误差为15%.

在位于回流区的站位x1/D=2.02 处,本文计算结果与Parnaudeau等[18]实验结果相比的最大相对误差分别为3.3%,7.2%,12.9%,0.5%,而Lehmkuhl等[10]利用Mode L 的DNS 计算结果与Parnaudeau等[18]实验结果相比的最大相对误差为7%.

对Case AV—Case DV 的情况,无论是本文计算结果还是Lehmkuhl等[10]利用Mode H 的DNS 计算结果一致,两者均与Lourenco 和Shih[27]的实验结果有一定的差异.在位于剪切层小尺度K-H 不稳定区的站位x1/D=1.06 处,本文对Case BV 和Case DV 工况的计算结果与Lehmkuhl等[10]利用Mode H 的DNS 计算结果之间的相对误差较小,分别为6.8%和3.7%,而对Case AV 和CaseCV工况,本文计算结果与Lehmkuhl等[10]利用Mode H之DNS 计算结果的相对误差较大,分别为29.5%和29.3%.

在位于回流区的站位x1/D=1.54 处,对Case AV—Case DV 工况,本文计算结果与Lehmkuhl等[10]利用Mode H 的DNS 计算结果的相对误差分别为7.75%,15.2%,14%和1%.在位于回流区的站位x1/D=2.02 处,对Case AV—Case DV工况,本文计算结果与Lehmkuhl等[10]利用Mode H 之DNS 计算结果的相对误差分别为1.7%,1.2%,0.5%,2.2%.

图10 分别给出了上述8 种工况下在3 个不同站位(x1/D=1.06,1.54,2.02)处各向异性雷诺应力剖面特性的数值结果.对此情况,Lehmkuhl等[10]没有给出相关的DNS 之计算结果.

图10 圆柱后方不同站位处各向异性雷诺应力剖面特性Fig.10.Distribution characteristics of anisotropy cross-flow Reynolds stresses at different locations in the backside of the circular cylinder.

由图10 可见,对图10(a)的情况,在位于剪切层小尺度K-H 不稳定性发生区域的站位x1/D=1.06 处,本文在Case AU—Case DU 这4 种工况下的计算结果均与Parnaudeau等[18]的实验结果具有良好的一致性.在位于回流区的站位x1/D=1.54 处,本文对Case AU 和Case BU 的计算结果与Parnaudeau等[18]的实验结果也具有良好的一致性,但在 |x2/D|＜0.48 的范围内本文对Case CU和Case DU 的计算结果,与Parnaudeau等[18]的实验结果有一定的差异,最大相对误差分别为38.7%和14.8%.在位于回流区的站位x1/D=2.02处,本文对Case BU—Case DU 的计算 结果 与Parnaudeau等[18]的实验结果具有良好的一致性,但在 |x2/D|＜0.48 的范围内本文对Case AU 的计算结果,与Parnaudeau等[18]的实验结果有一定的差异,最大相对误差为30.8%.

对图10(b)的情况,在位于剪切层小尺度KH 不稳定性发生区域的站位x1/D=1.06 处,以及在位于回流区的站位x1/D=1.54 处,本文对Case AV—Case DV 这4 种工况下的计算结果均与Lourenco 和Shih[27]的实验结果具有良好的一致性.在位于回流区的站位x1/D=2.02 处,本文对Case AV—Case DV 的计算结果与Lourenco和Shih[27]的实验结果相比,在“峰”和“谷”处有一定的差异.其中,在“峰”处的相对误差分别为30.2%,29.5%,38%,43.8%,而在“谷”处的相对误差分别为32.7%,31.8%,41.1%,40.8%.

Tremblay[8]采用DNS 和LES 同时也均复现了Lourenco 和Shih[27]对雷诺数Re=3900 下,圆柱绕流后方不同站位处雷诺应力剖面特性的实验结果.其中,LES 亚格子模型分别包括经典Smagorinsky 模型(Smago)及动态模型(Dyn).结果表明,DNS 的结果要优于LES 的结果,而Dyn 亚格子模型的结果要优于Smago 亚格子模型的结果.

对PANS 模型,其对圆柱后方各站位处雷诺应力剖面的计算精度同样地依赖于滤波控制参数fk的取值方式.对DES 类模型,其对圆柱后方各站位处雷诺应力剖面的计算精度则依赖于所选用RANS 模型的类型及RANS/LES 之间的混合转换方法等.

由于这两类模型均采用了相关的RANS 模型,而RANS 模型对湍流尺度的解析能力是基于对涡黏系数的调控而实现,但这种调控与计算网格的空间分布特性对湍流尺度的系统性分辨率直接相关.通常情况下,虽然可以通过加密网格来调控这两类模型对湍流尺度的解析能力,但在近壁面网格较密的条件下,会出现对网格空间分布特性过于敏感的缺陷,导致近壁面RANS 区域被破坏,而网格空间分布尺度还没有小到足够进行大涡模拟的程度[25].

本文采用WM-HRL 模型的系列数值模拟结果表明,在同一套网格系统下,通过改变WM-HRL模型中“两个边界位置”及“三个区域的湍动能解析度”的组合,既可通过数值模拟获得与Parnaudeau等[18]实验结果一致的圆柱后方各站位处的雷诺应力剖面,也可以通过数值模拟获得与Lourenco 和Shih[27]实验结果一致的圆柱后方各站位处的雷诺应力剖面.

研究发现,对本文所构建的WM-HRL 模型,在RANS/LES 混合转换区边界位置及其两个指标参数nrk1-Q和nrk2-Q准则取值的某些组合下,计算所得二阶统计量与相关文献的实验结果相比,存在较大的误差,其主要原因如下.

第一,在本文所构建的WM-HRL 模型中,其RANS 模型采用的是SSTk-ω模型.该湍流模型在边界层内为标准的Wilcoxk-ω两方程模型(简称SKW 模型).对WM-HRL 模型,其RANS 模型必须直接作用于圆柱近壁面处,但对SKW 模型,会存在如下缺陷:包括可能会过大地预测回流区的长度,还可能会过大地预测壁面剪切应力的值及过小地预测壁面湍动能的值[44].

第二,在本文所构建的WM-HRL 模型中,LES模型采用的是经典的Smargorinsky 型亚格子模型.该亚格子模型虽然概念简单且数值求解稳定性高,但不仅会在近壁面产生过大的数值耗散,而且定常的Smagorinsky 常数难以描述复杂时空变化下的湍流场.

为克服上述两个缺陷,对RANS 模型可改用Peng等[44]提出的低雷诺数Wilcoxk-ω两方程模型的修正版本.Peng等[44]的研究表明,该修正版本可有效地提高其对近壁面流的解析能力,包括可减小近壁面处ω的值及增大近壁面处k的值等,并具有更好地模拟分离、回流及再附等复杂流动现象的能力.

另一方面,对LES 模型可改用动态(Dyn)亚格子模型.Tremblay[8]的研究表明,对雷诺数Re=3900 下的圆柱绕流场问题,与基于经典Smagorinsky 型亚格子涡黏系数的LES 模型相比,采用基于动态(Dny)亚格子涡黏系数的LES 模型,可以更好地获得圆柱绕流场相关统计量的计算结果.这些正是本文作者团队将在下一阶段重点开展的相关研究工作.

4.3 相干结构频谱及流场特征

首先研究剪切层小尺度K-H 不稳定性结构发生区域的范围问题.为此,选取Case CU 和Case CV 两种工况进行分析.研究对象为这两种工况下,在图3 中P1—P12 这12 个监测点处横向脉动速度的Lomb 谱特征,结果如图11 和图12 所示.

图11 在Case CU 工况下,在P1—P12 监测点处横向脉动速度的Lomb谱Fig.11.Lomb spectrums of the cross-stream fluctuation velocities at different probes P1-P12 for the Case CU.

图12 在Case CV 工况下,在P1—P12 监测点处横向脉动速度的Lomb谱Fig.12.Lomb spectrums of the cross-stream fluctuation velocities at different probes P1-P12 for the Case CV.

由图11 和图12 可见,在12 个监测点P1—P12处横向脉动速度的Lomb 谱中,均出现一个频率相对较低的谱峰,该峰值所对应频率与涡泄频率一致,其无因次频率≈0.22.同时,在12个监测点P1—P12处横向脉动速度的Lomb 谱中,还会出现一个与该频率峰值相对应的倍频()成分.前者为圆柱体升力脉动的主频成分,而后者为圆柱体阻力脉动的主频成分.

对Case CU 工况,在测点P4—P10 处横向脉动速度的Lomb 谱中(见图11),还可以观察到一个频率更高的谱峰,该峰值所对应频率与剪切层小尺度K-H不稳定性的频率一致,其无因次频率在1.33—1.47之间.进一步的研究表明,在此工况下圆柱绕流场中剪切层小尺度K-H 不稳定性结构发生在对数律层中 89.4 ≤y+≤253.5 的区域.

对Case CV 工况,在测点P4—P7 处横向脉动速度的Lomb 谱中(见图12),也可以观察到一个频率更高的谱峰,该峰值所对应频率也与剪切层小尺度K-H 不稳定性的频率一致,其无因次频率也在1.43—1.44 之间.进一步的研究表明,在此工况下圆柱绕流场中剪切层小尺度K-H 不稳定性结构发生在对数律层中 89.4 ≤y+≤167.4 的区域.

讨论其他6 个监测点P13—P18 处流向和横向脉动速度的Lomb 谱特征,计算工况选取Case CU 和Case CV 两种情况,结果如图13 所示.其中,监测点P13 和P14 分别位于剪切层偏上及偏下的位置,监测点P15 位于剪切层脱落区,监测点P16 位于回流区内的尾流中线上,而监测点P17 和P18 位于回流区外的尾流中线上.

图13 在Case CU (第1 和第2 列)和Case CV (第3 和第4 列)工况下,在P13—P18 监测点处流向(第1 和第3 列)及横向(第2 和第4 列)脉动速度的Lomb谱Fig.13.Lomb spectrums of the stream-wise (from the first to third rows) and cross-stream (from the second to fourth rows) fluctuation velocities at different probes P13-P18 for the Case CU (from the first to second rows) and the Case CV (from the third to fourth rows).

由图13 可见,在Case CU 和Case CV 两种工况下,对P13—P18 这6 个监测点,无论是在其流向脉动速度的Lomb 谱中,还是在其横向脉动速度的Lomb 谱中,均清晰地出现两个频率较低的谱峰,一个谱峰对应涡泄频率(),另一个谱峰对应涡泄频率的倍频().

对位于剪切层脱落区的监测点P15 的情况,在Case CU 和Case CV 两种工况下,无论是在其流向脉动速度的Lomb 谱中,还是在其横向脉动速度的Lomb 谱中,除了均清晰地出现单频频率及其倍频的谱峰外,还可观察到一个三倍频率的谱峰.

对位于其他区域的监测点P16—P18 的情况,在Case CU 和Case CV 两种工况下,从其流向脉动速度的Lomb 谱中只能清晰地看到单频频率及其倍频的谱峰,而从其横向瞬态速度的Lomb 频谱图中则只能清晰地看到单频频率及其三倍频率的谱峰,但没有出现明显的二倍频率的谱峰.

由图13 进一步可见,在Case CU 和Case CV两种工况下,对位于剪切层偏上方的监测点P13,无论是在其流向脉动速度的Lomb 谱中,还是在其横向脉动速度的Lomb 谱中,除了均清晰地出现与涡泄相关的单频频率及倍频频率的谱峰外,还可清晰地观察到一个频率更高的谱峰,其对应频率与剪切层小尺度K-H 不稳定性结构的频率一致,其无因次频率在1.43—1.44 之间.

在Case CU 和Case CV 两种工况下,对位于剪切层偏下方的监测点P14,在其流向脉动速度的Lomb 谱图中,只清晰地出现了与涡泄相关的单频频率及倍频频率的谱峰.然而,在其横向脉动速度的Lomb 谱图中,除了均清晰地出现与涡泄相关的单频频率及倍频频率的谱峰外,还可清晰地观察到一个与剪切层小尺度K-H 不稳定性结构频率一致的谱峰,其无因次频率也在1.43—1.44 之间.

最后,从图13 可见,在Case CU 和Case CV 两种工况下,对监测点P13 和P14,无论是在其流向脉动速度的Lomb 频谱图中,还是在其横向脉动速度的Lomb 谱中,均没有出现明显的三倍频率的谱峰.

对雷诺数Re=3900 下圆柱绕流场中出现三倍频率的现象,在Parnaudeau等[18]的DNS 及Pereira[2]等的PANS 计算中都有发现,而在Lehmkuhl等[10]采用DNS 及Tremblay[8]和Breuer[15]采用LES 的计算中则没有发现.Ong 和Wallace[45]指出,流体行为中的峰值现象可以通过线性稳定性理论进行预测,用于研究各种不稳定性问题.然而,这个理论是否能够适用于本文所观察到的这种复杂现象,或者本文算例所出现的这类复杂频谱现象是否还存在其他影响因素,值得在后续工作中做进一步的研究.

最后讨论圆柱绕流场中两类相干结构的流场特征,计算工况选取Case AU—DU 及Case AV—DV 两类情况,结果如图14 和图15 所示.其中,图14 的左列为x3/D=π/2 断面上展向涡量ω3=∂u2/∂x1-∂u1/∂x2的云图,右列为流向速度云图,且图中白色垂直虚线为回流区末端位置(即u1/U∞=0位置).图15 为沿圆柱体长度的三维展向涡量云图.

图14 在Case AU—DU(前4 行)和Case AV—DV(后4 行)工况下,圆柱绕流涡量(左)及流向速度(右)云图Fig.14.Contours of the span-wise vorticity (left) and stream-wise velocity (right) for both Case AU-DU (the first four lines) and Case AV-DV (the last four lines).

图15 在Case AU—DU (左)和Case AV—DV (右)工况下,圆柱绕流展向三维涡量云图Fig.15.Contours of the three-dimensional span-wise vorticities both Case AU-DU (left) and Case AV-DV (right).

在Case AU—DU 及Case AV—DV 两类工况下,从图14 (左)均清晰地可见附着于圆柱上下侧壁的两个层流分离剪切层结构,表现为两条位置狭小的窄带结构.同时,从图14 (右)则可清晰地观察到圆柱尾部因流动分离而形成的两个回流区分支结构.具体来说,第1 个分支结构发生在x1/D=1.06 处平均流向速度剖面为U 形的情况,回流区长度较长;而第2 个分支结构发生在x1/D=1.06 处平均流向速度剖面为V 形的情况,回流区长度较短.

结合图11 和图12 的结果进一步可知,在Case AU—DU 及Case AV—DV 两类工况下,两个剪切层中小尺度K-H 不稳定性结构的触发位置都基本相同,均约位于y+=89.4 处,但两个剪切层中小尺度K-H 不稳定性结构的消逝位置不同.在Case AU—DU 工况下,其消逝位置约为y+=253.5.而在Case AV—DV 工况下,其消逝位置约为y+=167.4 .

这意味着,当x1/D=1.06 处平均流向速度剖面为U 形时,剪切层中小尺度K-H 不稳定性结构的长度较大,而且相应的回流区长度也较长.而当x1/D=1.06 处平均流向速度剖面为V 形时,剪切层中小尺度K-H 不稳定性结构的长度较小,而且相应的回流区长度也较短.因此,对雷诺数Re=3900 下圆柱绕流中出现两类回流区分支结构,以及在圆柱尾部近壁面处平均流向速度出现U 形和V 形两类剖面结构形状的形成机理,可能与WMHRL 模型中两个边界位置及3 个区域的湍动能解析度发生变化时,该湍流模型对圆柱绕流剪切层中小尺度层流分离流动的解析能力不同有关.

对DNS 模型,在同一套网格系统下,其对圆柱绕流剪切层中小尺度层流分离流动的解析能力是确定的,因此只能获得某一种回流区分支结构及与之对应的圆柱尾部近壁面平均流向速度剖面的某一种形状.但对不同的网格系统,由于其对圆柱绕流剪切层中小尺度层流分离流动的解析能力相应地会不同,进而可能会出现在不同网格系统下获得不同的回流区分支结构及与之对应的圆柱尾部近壁面流向速度剖面的U 形或V 形结构.

对PANS 及DES 类模型,它们可选用不同的RANS 湍流模型,也可以改变湍流模型中的相关模型参数.对DES 类模型,其执行RANS/LES 之间转换过渡的方式还可以不同.另一方面,对LES 模型,其可以选用经典Smagorinsky 及Dyn等[9]不同的亚格子涡黏模型.这些复杂的因素都可能会导致即使在同一套网格系统下,使得这些湍流模型对圆柱绕流剪切层中小尺度层流分离流动的解析能力不同,进而可能会获得不同的回流区分支结构及与之对应的圆柱尾部近壁面平均流向速度剖面的U 形或V 形结构.

亚临界雷诺数区圆柱绕流中主要有两类相干结构.其中,一类为剪切层小尺度K-H 不稳定性结构,而另一类为大尺度的Karman 涡街结构.为进一步阐述它们的形成机理及其特征,在图15中,给出了 在Case AU—DU 和Case AV—DV两类工况下,圆柱绕流展向三维涡量云图的数值结果.

由图15 并结合图14 可见,由于剪切层两侧的速度差诱发K-H 不稳定性,使剪切层失稳并出现滚卷现象.一方面,这种现象会在剪切层中导致大量小尺度的旋涡产生于此区域,在横向瞬态速度的Lomb 谱中表现为宽频信号的凸起(见图11 和图12).另一方面,这种现象还会在圆柱尾迹区形成另一类大尺度的旋涡结构,即Karman 涡街结构,在横向瞬态速度的Lomb 谱中表现为频率相对较低的窄频信号的凸起(见图11—图13).

由图15 还可以清晰地观察到,附着于圆柱壁面的剪切层伴随着其失稳、滚卷等复杂现象在圆柱尾流区形成流向与展向相嵌套的三维涡体结构.由于本文所构建的WM-HRL 模型具有能够精细解析湍流谱结构的能力(见图11—图13),进而不仅能够精细地解析各种小尺度涡结构(如K-H 不稳定性结构等),又能够精细地解析各种大尺度涡结构(如Karman 涡街等).因此,该WM-HRL 模型表现出预期的相当于LES 甚至DNS 对各种小尺度及大尺度运动解析的能力.

5 结论

本文通过修改SST-k-ω湍流模型中k方程的色散项,构造了一种新的具有壁面模化能力的钝体绕流混合RANS/LES 模型,即WM-HRL 模型.该模型在钝体近壁面的某个计算区域内采用RANS 模式,在经历一个RANS/LES 混合转换过渡区后,在其余计算区域则采用LES 模式,且该LES 模式等效为经典Smagorinsky 亚格子模型.

通过构造一个仅与当地网格空间分布尺寸相关的湍动能解析度指标参数rk,提出了确定该WM-HRL 模型两个边界位置和3 个区域湍动能解析度的nrk1-Q和nrk2-Q准则.其中,nrk1-Q准则用于确定WM-HRL 模型中RANS 模式结束边界位置,及其在RANS 区湍动能的最大解析度(即rk1),而nrk2-Q准则用于确定WM-HRL 模型中LES 模式启动边界位置及其在LES 区对湍动能的最小解析度(即rk2).

在此基础上,本文构造了一个新的具有双重保护功能的混合函数fnd.第一重保护通过nrk1-Q准则实现,保证在rk≥rk1时WM-HRL 为完全RANS模式.通过第一重保护层的设置,可使WM-HRL既具有延迟脱体涡模拟(DDES)的能力,又具有壁面模化(WM)的能力.在RANS 区最大湍动能解析度指标rk1的合理设置下,可避免MSD 和GIS的问题.

第二重 保护通过nrk2-Q准则实 现,可保证在rk≤rk2时WM-HRL 为完全LES 模式.在此第二重保护层设置的条件下,可通过rk2值的合理设置,一方面在进行计算网格设置时即可实现对LES 启动边界位置的预先设定,另一方面还可自定义LES 区中WM-HRL 对钝体绕流各种复杂湍流场的解析能力.

在合理设置nrk1-Q和nrk2-Q准则的条件下,既可保证WM-LES 模型在RANS/LES 混合转换区即可激活LES 模式,同时还可保证在完全LES 区即可完全激活LES 模式,进而避免LLM 的问题.

基于该WM-HRL 模型,本文以雷诺数Re=3900 下圆柱绕流为对象,开展了系列数值模拟分析与评估,得到如下主要结论.

1)即使当LES 的启动边界位于圆柱绕流剪切层小尺度K-H 不稳定性区内时,在合理设置的nrk2-Q准则(一般需要rk2≤0.5)的条件下,该WMHRL 模型也能获取准确的涡泄频率、分离角、阻力系数及基准线压力系数等统计量信息.

2)为使WM-HRL 模型能够准确获取圆柱绕流剪切层小尺度K-H 不稳定性结构特征及其频谱特性,以及回流区两个分支结构的长度信息,可将LES 模式的启动边界位置均设置在过渡层内,而RANS 模式的结束边界位置则既可在黏性底层也可在过渡层内,并且使rk2≤0.2,即在对数律层区LES 对湍动能具有至少具有80%的解析能力.

3)在上述2)的条件下,WM-HRL 模型能够获取与相关实验精度相当的圆柱绕流一阶和二阶统计量的信息,以及剪切层小尺度K-H 不稳定性结构的无因次频率f¯kh及回流区两个分支结构的长度信息,并且当回流区长度的数值结果与Parnaudeau等[18]的实验结果一致,在x1/D=1.06 处的平均流向速度剖面为U 形,而当回流区长度的数值结果与Lourenco 和Shih[27]的实验结果一致,则为V 形.

特别地,在本文的系列数值模拟中,通过两个边界位置及3 个区域的湍动能解析度的各种不同组合,在仅用一套网格系统的条件下,利用新构造的WM-HRL 模型,针对雷诺数Re=3900 下圆柱绕流场问题,获得了两类长度不同的剪切层小尺度K-H 不稳定性结构特性,具体如下.

第一类剪切层小尺度K-H 不稳定性结构的长度较长,此时回流区长度也较长,与之相应的圆柱尾部近壁面处的平均流向速度剖面为U 形.第二类剪切层小尺度K-H 不稳定性结构的长度较短,此时回流区长度也较短,与之相应的圆柱尾部近壁面处的平均流向速度剖面为V 形.

关于雷诺数Re=3900 下圆柱绕流场中回流区的两个分支结构,以及与之相应的圆柱尾部近壁面处U 形和V 形两类平均流向速度剖面形状等问题,已经成为困扰学界三十余年的难题.本项研究为后续更深层次认识这一难题提供了一种新的手段.