吴正人， 程相辉， 石祎炜， 吴 浩， 刘 梅

(1.华北电力大学 能源动力与机械工程学院，河北 保定 071003；2.华北电力大学 经济与管理学院，河北 保定 071003)

0 引 言

在流体力学研究领域中,紊流占有极其重要的地位,在日常生活中及能源利用上都十分常见。

因此,利用近壁区湍流动力学过程的减阻控制方案实现湍流摩擦阻力的降低,对工程实际中能源的充分利用具有重要意义。

自从Walsh[1]指出顺流向的对称V型沟槽减阻效果最佳后,大多数的科研人员都着眼于此。对沟槽结构探究的同时,人们对沟槽的作用机理也展开了研究。目前沟槽作用机制可以区分为横向沟槽和纵向沟槽两种。对于纵向沟槽,其主要的作用理论分别为“突出高度论”[2]和二次涡群论[3]。而横向沟槽与纵向沟槽有所不同,其横向布置的特性会阻止动量沿流向的传递,增加动量的横向交换[4]。Mariotti[5]认为横向沟槽的作用机制是回流区的无滑移边界条件的松弛。Amy[6]也发现了相同的作用机制,并且认为沟槽流域边界层的改变削弱了低速条带,延缓了流动分离。沟槽内涡结构的存在使壁面附近存在相对速度,使速度梯度减小,从而减小了粘性阻力[7]。无论是纵向沟槽还是横向沟槽,其减阻机理都与湍流相干结构密切相关。已有研究指出,湍流中高摩擦阻力的产生和内部存在的拟序结构有很大的联系[8]。

通过湍流相干结构的演化来探究沟槽减阻机理是一项重要的研究手段。但传统的涡量参数并不能很好表征真实的涡结构,于是人们试图用涡量大小来推断湍流中的拟序结构和识别涡核。然而,正如Jeong和Hussain[9]所指出的,这种方法并不总是成功的。这使人们意识到涡量并不等于整体的旋转,即旋涡。显然两者之间存在着差异,准确区分这两者是十分必要的。

从认识到不能用反对称张量或涡量来度量涡以来,人们陆续提出了以速度梯度张量的特征值为基础的涡识别方法,主要有Q、λ2、Δ、λci等方法[10-12],在涡结构的识别方面取得了很大进展。然而这些方法都存在着不同的问题[13],如不能体现漩涡方向性、物理意义并不明确、存在剪切或拉伸压缩污染、阈值敏感性、阈值选取通常依靠经验,主观性较强。采用这些识别方法进行分析时大多要辅以其他变量,增加了工作量与分析难度。最近,Liu等[14,15]将涡量进一步分解为旋转部分和非旋转部分,提出了从流体运动中提取刚性转动部分来表征涡结构的识别方法,并命名为Rortex向量(后改为Liutex参数)。刘超群等人也提出了Omega准则[16-19],该准则将涡量中的旋转量分离出来,将剪切变形量舍去,来更准确的表示涡结构,研究结果表明,当Ω为0.52时可以正确的捕获涡结构,Ω是一个近似定义涡旋的边界量。

出于传统旋涡度量方法的不足。本文使用 ANSYS FLUENT 数值计算软件,首先对Omega准则与Liutex参数在涡结构识别上的优势进行了直观展现,然后采用新型漩涡度量参数Liutex与Omega漩涡识别准则,对比分析了平板与沟槽流域的湍流相干结构的识别结果,探究了Omega准则与Liutex参数在湍流减阻研究中的适应性问题。

1 计算模型和数值计算方法

1.1 沟槽结构设定

目前对具有减阻效果的横向沟槽的具体尺寸参数并无确切定论。虽然沟槽具体结构参数并不能确定,但S. Sutardi[4]、Amy W Lang[6]、C. Y. Ching[20]等人经研究后给出了减阻沟槽结构参数的大致范围,即s+(无量纲沟槽宽度)< 50和h+(无量纲沟槽高度)< 30。

本文采用V型对称沟槽,综合考虑模拟工况与前人研究结果,将沟槽尺寸确定为s= 0.2 mm,h= 0.1 mm。表1为模拟中采用工质的物性参数。

表1 工质物性参数

1.2 数值模型

计算模型如图1所示。模拟计算采用40s×15s×50s(s= 0.2 mm)的计算模型,沟槽结构间距s=0.2 mm,高度h= 0.1 mm,顶部角度为90 °。其底部布置沟槽结构,顶部设置为光滑表面,两者均定义为 wall 条件。流向和展向部分流域具有周期性,为了确保湍流沿流向充分发展以及两个侧面对流域没有干扰,将其设置为周期性边界条件。

图1 计算模型Fig. 1 Calculation model

1.3 数值方法

本文采用大涡模拟方法对槽道流域进行计算,利用有限体积法进行空间离散,二阶迎风差分方法来处理 N-S 方程,采用 SIMPLEC 算法即压力耦合方程组的半隐式方法来求解流场。模拟初期先采用 RNG k-ε 模型进行迭代,当模拟结果符合稳定湍流特点后,将其作为大涡模拟的初场进行计算。采用大涡模拟进行非定常计算,其中大涡模拟采用 Smagorinsky 亚格子滤波模型,时间步长取10-5s,残差设定为10-5。

计算域网格采用结构化六面体网格,由于计算区域在法向上的对称性,本文只给出一半槽道的网格划分,另一半由对称性给出网格点位置。本文对上下壁面附近网格进行了加密,即在法向上采用非均匀网格,网格增长率为1.03,其中流向与展向为均匀网格。为了满足要求,初步设置网格数为316万,然后逐步提升至345万、374万、403万,在不同流速下进行数值模拟,其第一层壁面y+值如图2所示。

图2 不同流速下第一层网格y+Fig. 2 The first layer of grid y+ at different flow rates

由图2可以发现,其网格数至少到403万时,才符合大涡模拟的网格数要求。

再次对网格进行无关性验证,将网格加密到432万,在同样的条件下计算边界层内的综合湍流度,结果如图3所示。网格加密后综合湍流度并没有太大的变化,局部网格如图4所示。

图3 不同网格分辨率下综合湍流强度分布图Fig. 3 Comprehensive turbulence intensity distribution map under different grid resolutions

图4 网格局部示意图Fig. 4 Grid partial schematic diagram

为了确保数值模拟结果的准确性,论文采用了Kim[21]等人的脉动速度实验数据进行对比。对比结果如表2所示。两者可以很好的匹配,误差均在2%以内,这足以证明该模型的准确性。

表2 流向脉动速度均方根对比

2 数值结果分析

湍流作为一种流体的运动流态,是随机且非定常的,目前还没有一种具体且准确的定义。但在随机中还有着一定的规律,湍流中的拟序结构的发现使研究者对其的了解可以更进一步。

2.1 Omega准则

工程中为了探究湍流中的涡结构,一般对速度梯度等进行一定的数学处理,从而定义识别准则来进行可视化的研究。考虑到识别准则的真实性、精确性及阈值的选取问题,现针对目前应用较多的几个识别方法(Ω、涡量、Q和λ2)进行对比研究。

图5为常见识别准则的对比图。在初步确定其壁面湍流中的拟序结构时,除涡量不能准确的展示,其他3种方法均可以表现其基本特征。涡量识别方法中,壁面附近结构连成一片,6个涡核结构附近有连结情况,并且当调整阈值直到涡核结构变得不完整时,其连结情况也没有完全消失。

图5 常见识别准则对比Fig. 5 Comparison of common identification criteria

Q准则以Q>0代表涡结构,λ2准则以λ2< 0代表涡结构,但取值范围的量级均为1010,这就给涡的准确识别带来了很大麻烦,因为阈值的选取并没有统一的判定标准。换言之,不同的阈值选取会导致不同的涡结构,产生不同的物理结论。归根到底,这都是第二代涡识别方法的固有缺陷导致的。Q准则和λ2准则,均属于第二代涡识别方法,在第二代识别方法中,涡强度都是标量,有大小无方向,由于是标量,所以需要人为给定一个阈值,人为选取就具有很强的主观性,而这个阈值对涡结构的显示又具有重要影响,阈值太大,则小涡消失不见,阈值太小,强涡弱涡混合在一起,强弱不分,给研究工作带来了很大的麻烦。如上所言,阈值的选取对于Q准则和λ2准则来说,是一个具有挑战性的任务。

而对于Omega准则,虽然其也属于第二代涡识别方法,但由于其是一个正交化标量,可以同时捕捉强涡和弱涡,较其它第二代识别方法,优势明显。Ω取值范围为0到1,当其值大于0.5时,代表其旋转强度大于变形,一般情况下选择0.52即可表示涡核的基本特征,调整均在0.5～0.7的范围内,阈值选取简单方便。

Omega准则的优势不仅是阈值选取范围窄,更重要的是其对阈值敏感性低,可以更好地表征涡结构。将Omega准则与目前应用最多的Q准则进行了对比分析。通过图6可以看到,在Omega准则的阈值范围内,其6个主要的涡核结构除了变细之外,并没有发生太大的变化,其结构仍然十分完整。而右侧的Q准则的情况则不同,当阈值从3 590 000变为5 000 000时,可以看到其结构有了部分的缺失。其值增大到80 000 000时,左侧的涡头结构已经消失不见。对比两者可以发现,Omega准则的鲁棒性很好,在整个取值范围内,其主体结构并不会发生变化。但对于壁面附近的小涡捕捉效果来看,两者优劣并不容易判断。

图6 Q准则与Omega准则的对比图Fig. 6 Comparison of Q criterion and Omega criterion

对于两种识别方法来说,壁面附近的小涡都可以较为清晰的展示出来,两者仅有细微的差别。实际大涡模拟中的小涡结构较小,且其没有特定的结构,一般很难进行精确的捕捉。通过统计方法进行分析时,两种方法都可以进行参考,但综合其大涡的捕捉效果来看,Omega准则确实更具优势。

2.2 Liutex参数

涡量是经常接触到的一个物理量,在湍流的研究中曾起到了很关键的作用。与涡量相似,Liutex作为矢量参数也具有三维特征且符合右手定则。现在将两者进行对比,观察Liutex参数的优劣。

图7为涡线与Liutex线对比,图7(a)展示了涡量图中涡线的分布趋势。可以发现在主体部分的涡线大体吻合的时候,壁面附近的涡线杂乱无章,几乎无法得到有用的信息。对于图7(b)而言,无论是主体结构部分还是壁面附近的小涡结构,其Liutex线都被包含其中,两者十分的吻合。将涡线和Liutex线从等值面图中抽离出来,如图7(c)和(d)所示,这种现象表现得更加直观。从图7(c)中可以看到,涡线的分布十分的混乱,仅少数部分可以明显的看出规律,大多数的位置都被杂乱无章的涡线所覆盖,使其结构规律更难以探寻;但是对于Liutex线来说,图形显得十分的规整美观,其涡结构已经可以很明显的展现出来。就沟槽减阻研究中涡识别准则的适应性来看,涡量对于旋涡的描述效果较新的Liutex参数欠佳,Liutex参数更加简洁美观。

图7 涡量与Omega准则对比Fig. 7 Comparison of vorticity and Omega criterion

在上图中,由于Liutex为Omega准则的数学基础,所以其直观准确的线性可能是由于两者的相关性。所以进一步的采用Q准则来验证Liutex线与涡线的准确性。图8是Q= 8 000 000的局部放大图,分别采用了涡线与Liutex线对其进行描述。在画出的涡线中可以发现,除了发卡涡的涡头中较为符合之外,其他部分的涡线与识别出来的发卡涡结构有所偏差。而对于图8(b),图8(d)来说,两者几乎完美吻合。对Liutex线使用右手定则很容易发现,在发卡涡结构的内部位置(两个涡腿之间)流体从下向上流出,形成一个上抛过程,外部位置(两个涡腿以外)流体从上往下进行一个下扫过程,造成了其雷诺应力的急剧变化。

图8 Q = 80 000 000时涡线与Liutex线对比Fig. 8 Comparison of vortex line and Liutex line when Q = 80 000 000

2.3 Omega算例应用

2.3.1 平板表面的涡结构演化

图9是平板表面涡结构的时间演化过程。Ω取值为0.55,此时的涡结构已清晰的展示出来,并且没有联结情况。平板表面的涡结构多且杂乱,其发卡涡结构稳定存在6个,发卡涡结构一般处于对数区域内,其结构大,影响的空间范围广。该结构的存在会对周围区域进行扰动,发卡涡内侧“上抛”,外侧进行“下扫”运动,是湍流的一个十分重要的拟序结构。可以看到随时间的推移,其发卡涡结构并无太大的变化,但是底部小涡存在增加和破碎的情况。

图9 平板涡结构的时间演化Fig. 9 Time evolution of flat vortex structure

正如刘超群教授所言,可能小涡的形成与大涡的破碎并没有较为直接的联系,或者涡结构的等级分明,壁面附近的小涡只从壁面附近的稍大一点的涡处获取能量,并不会影响类似发卡涡这样的涡结构。总体来说,整个区域的涡结构并没有随时间有太大变化,仅仅是部分小涡发生了改变,这些变化从图中可以较清晰的展示出来,且取值范围小,在应用时便于操作,减小了很大的工作量。

2.3.2 沟槽表面的涡结构演化

图10为沟槽表面部分涡结构的时间演化过程,此时的Ω值为0.55。可以看到,沟槽表面的发卡涡数量明显减少,在此时段仅仅存在1个,表面附近的涡更多的以流向涡的形式存在,并没有向上抬升形成发卡涡,这些结构都可以清晰展示出来。据计算,流向涡的存在范围大概在y+ <60内,而发卡涡一般位于y+ = 40～100的范围内,说明沟槽结构的存在抑制了流向涡的抬升,避免其形成发卡涡结构。发卡涡结构较大,对湍流内的流动具有较大的影响,且后续涡腿附近会诱导形成新的流向涡,进而引起“上抛”和“下扫”运动,对流态产生不好的影响。发卡涡及流向涡结构没有随着时间产生大的变化,也说明了其结构的稳定性。从图中来看,Omega准则可以准确的识别沟槽流域边界层内的涡结构,对其有较好的适应性。

图10 沟槽涡结构的时间演化Fig. 10 Time evolution of groove vortex structure

沟槽表面的涡结构整体移动的速度均小于平板表面的涡结构迁移速度,其发卡涡的涡头所处的位置比起平板区域来说距离壁面的法向高度要低,这可以使边界区域的影响范围变小,是十分有利的。特别的是,与平板表面对比,沟槽底部形成了均匀的横向涡柱,这是沟槽结构的重要作用所在,这种二次涡结构并没有随着时间推移而产生太大的变化,说明其状态的稳定性,可以稳定持续的发挥“滚动轴承”的作用,这是目前横向沟槽减阻的主要理论。

2.4 Liutex算例应用

2.4.1 沟槽中心沿垂直方向的Liutex分布情况

图11为不同沟槽处的liutex的分布情况,取点位置为沟槽中心处沿y方向的一系列点位。可以看出,在沟槽内部的明显特征可以很好的显示出来,其Liutex参数在流向及法向方向数值并不大,而在展向(z)上的数值与总量Liutex-m相当,Liutex总量的大小主要是由z方向上的分量提供,而且数值为负,说明其Liutex的方向与z轴正方向相反。由右手定则可以知道,横向沟槽内的旋转方向在旋涡上部与流向的方向是相同的,这也与“滚动轴承”的现象相吻合。

图11 沟槽内部Liutex分布Fig. 11 Liutex distribution inside trench

在该沟槽内,通过Liutex分量的分布可知,发现有一个顺流向转动的旋涡存在。在底部沟槽区域的不同位置取样本点,可以发现各个图像的规律基本一致,这也表明了在底部沟槽位置,起“滚动轴承”作用的二次涡并不是某一区域独有,具有一般性。

2.4.2 沟槽垂直方向上Liutex分布

图12所表示的是沟槽对上部区域的影响。其中y= 0为沟槽底部位置,y= 0.1处为沟槽顶部也即是平板壁面处。可以看到的是,在各个方向的Liutex分量都显著的减小,波动幅度下降。Liutex-m从最大值为12 000 s-1直接下降到6 000 s-1左右,下降十分的明显。Liutex-z改变了矢量方向,光滑平板流域的数值均为正值,布置沟槽的平板上则成为了负值。除了沟槽内部的Liutex-z值为3 000 s-1左右,(根据右手定则判断沟槽内存在顺时针转动漩涡。)沟槽顶部以上区域Liutex-z波动范围均在1 000 s-1以下,较平板明显降低。至于Liutex-y,其数值最大值仅在200 s-1左右,比起其他方向的分量,数值很小,对整体流场影响不大。

图12 沟槽垂直方向上部Liutex分布Fig. 12 Liutex distribution on upper part of trench in vertical direction

在平板壁面处Liutex基本为零这是由于平板的无滑移壁面,而在沟槽顶端与平板底面相同高度处(y=0.000 1)Liutex-m接近2 000 s-1,且主要是以负的Liutex-z分量为主。主要是因为沟槽内二次涡使沟槽顶部具有了一定的漩涡强度。因此能够发现沟槽结构的作用机制:横向沟槽内顺时针转动的二次涡上部转动方向与主流方向一致,该涡结构在底部起到“滚动轴承”的作用,使沟槽壁面具有一定的滑移条件,从而改善底层流动状况及漩涡强度,也因此减小了流动的粘性阻力。

整体来看,沟槽对近壁区0.001 m高度以下的范围内的降低作用十分显著,随着高度的上升,其影响逐渐的减小。这可能与上部涡结构减小有关,光滑平板上部当高度处于0.002 m时,其Liutex已几乎为0,说明涡结构只存在于近壁区,其远离壁面的区域涡结构较少或较弱。

3 结 论

本文依托横向沟槽结构,采用大涡模拟方法,探究Liutex参数与Omega准则对算例的适应性。通过与常见准则(Ω、涡量、Q和λ2)的对比分析,以及对涡量的探究,展示了一些横向沟槽流域的实际算例,可以得出以下结论:

(1)Ω的阈值选取从0到1,对所研究的模型来说具有很好的适应性,且稳定性高,鲁棒性强,旋涡的主体结构(如发卡涡头)在调整过程中一直保持较高的完整性,但对于壁面附近的小涡结构,并没有优于Q准则的表现,通过统计方法进行分析时,两种方法都可以进行参考,但综合其大涡的捕捉效果来看,Omega准则更具优势。

(2)相较于涡线对涡结构的展现效果,新的Liutex线更加准确简洁,图形更加规整美观,涡结构展现更为清晰。

(3)Omega准则与Liutex参数对于横向沟槽结构内的流域均具有很好的适应性,可以清晰准确的显示出边界层内的涡结构及沟槽内的二次涡结构,且Omega准则的可视化与Liutex参数的准确物理意义以及方向性,充分验证了横向沟槽的宏观减阻机理(沟槽内二次涡引起的壁面无滑移边界条件的松弛)。

基于本文的研究结果,以后的工作将依据Omega准则以及Liutex涡强度参数,从壁面湍流相干结构的演化方面来探究横向沟槽结构的减阻机理。