空腔流动的拉格朗日涡动力学分析

2019-11-05韩帅斌罗勇张树海

航空工程进展 2019年5期

韩帅斌，罗勇，张树海

(1.中国空气动力研究与发展中心空气动力学国家重点实验室，绵阳 621000) (2.中国空气动力研究与发展中心计算空气动力研究所，绵阳 621000)

0 引言

空腔流动是流体力学的经典流动问题，包含丰富的涡结构相关的物理现象，例如剪切层中涡的卷起与撞击、涡与剪切层相互作用、涡声的产生与传播等。同时空腔流动具有较强的实际应用背景。在航空航天领域，剪切层内的涡与剪切层及空腔拐角相互作用产生的强烈噪声不仅影响身体健康，更会导致飞行器结构疲劳，危及飞行安全，是目前迫切需要解决的问题。自20世纪50年代起，空腔流动尤其是空腔流致噪声得到了大量研究[1-3]。空腔流致噪声与来流参数、空腔几何形态等密切相关。J.E.Rossiter[4]通过系列实验总结得到Rossiter半经验公式，能够很好地预测辐射噪声的振荡频率与来流马赫数的关系。通过H.H.Heller等[5]、C.K.W.Tam等[6]的完善，公式预测范围和精度得到进一步提高。空腔流致噪声的动力学过程为：涡扰动在剪切层中不断增长，涡撞击空腔后缘产生声波，声波扰动向上游传播，声波扰动到达前缘诱导产生新的涡扰动[3]。该动力学过程的一个关键问题是涡量扰动和压力脉动之间如何相互转化。M.V.Morkovin等[7]研究了压力信号转化为涡量扰动的物理机制；Y.P.Tang等[8]对涡-边缘相互作用的形式进行了区分，对涡波转化压力波的物理机理作出了相应的阐释；万振华[9]借助涡量及拟涡能的时空演化分析了空腔中涡结构撞击后拐角产生压力脉动的动力学过程。空腔中的涡结构作为噪声的主要来源，其动力学过程对空腔流动特性及噪声的产生和传播起着重要作用，开展空腔流动中的涡动力学分析能够为相应的流动控制和降噪提供理论依据。

识别流场中的涡结构是开展涡动力学分析的关键步骤。目前存在诸多涡判据[10-11]可用来识别流场中的涡结构，但仍没有公认的涡的严格定义。既往的涡动力学研究大多采用欧拉框架下的分析工具例如闭合流线、涡量、Q判据、准则等识别流场中的涡，然而这些方法通常仅使用某瞬时流场的速度或速度梯度提取涡结构，揭示的是流动的瞬时状态，缺失了动力学系统中与时间相关的信息，无法准确反映流动结构的历史累积效应以及相应的动力学特性。另外欧拉框架下的判据通常不具有客观性，即流场结果会随参考系改变而改变。

近年来，拉格朗日方法在复杂流动系统的动力学分析及涡识别方面取得了较大进展[12]。拉格朗日方法在时空相空间内跟踪流体粒子运动，能够客观地揭示非线性动力系统的动力学特性和流动机理。基于拉格朗日拟序结构(Lagrangian Coherent Structures，简称LCS)识别流场中的涡结构并研究其时空演化规律，能够揭示开式空腔中的拉格朗日涡的动力学特性。

本文采用拉格朗日拟序结构，并结合Q判据和涡量通量分析，对开式空腔流动中涡的生成、脱落、对流以及撞击破裂过程展开研究，分析其动力学特性，为流动控制及降噪提供理论依据。

1 数值分析方法

1.1 高精度流场计算方法

通过求解非定常可压缩Navier-Stokes方程(简记为N-S方程)对空腔流动进行直接数值模拟，守恒形式的无量纲N-S方程为

∂tρ+∂i(ρvi)=0

(1)

(2)

(3)

其中，

(4)

(5)

(6)

粘性系数的Sutherland公式为

(7)

N-S方程的对流项和粘性项分别采用五阶WENO格式和六阶中心差分格式进行离散。时间项采用三阶TVD Runge-Kutta格式求解获得高精度非定常无量纲流场数据。

1.2 拉格朗日拟序结构

拉格朗日拟序结构是基于追踪流体粒子运动而提取的流场中主导流动的吸引性、排斥性或者强剪切拟序结构，是流动的“骨架”[12]，目前已被广泛应用于流动分离[13]、涡动力学[14]、生物流体力学[15]等。诸多方法可以用来提取流场中的LCS，其中有限时间李雅普诺夫指数(Finite Time Lyapunov Exponent，简称FTLE)是被广泛应用的一种方法。对于在非自治系统

(8)

(9)

式中：λmax为矩阵的最大特征值。

在实际计算中，FTLE既可以前向时间积分也可以逆向时间积分，对应的FTLE的嵴分别捕捉流场中的前向LCS(pLCS)和逆向LCS(nLCS)，pLCS通常是流场中的排斥性结构，nLCS则通常是流场中的吸引性结构。H.Yangzi等[16]采用pLCS和nLCS的边界识别拉格朗日涡边界并追踪涡运动，本文采用该方法对空腔中的拉格朗日涡结构进行动力学分析。

1.3 Q判据

欧拉框架下常用来识别涡的一种方法是J.C.R.Hunt等[17]提出的Q判据。速度梯度张量可以分解为对称的应变率张量S和反对称的涡张量Ω，即

(10)

Q值的定义即为

(11)

在流场中Q值为正的地方，也即涡张量对流体变形的贡献大于应变率张量贡献的区域被识别为涡。Q判据具有伽利略不变性，但在旋转或加速参考系中则不能正确进行涡识别。对于本文中的开式空腔流动，参考系不存在旋转或加速，Q判据能够给出涡的正确描述，因而被用来与LCS识别的拉格朗日涡对比验证，并描述相应的涡动力学。

2 数值计算及验证

2.1 计算模型及设置

本文直接数值模拟的空腔流动计算域及空腔附近网格如图1所示，其中空腔长深比为2∶1，各项无量纲长度参数为

L=2D,Ll=7D,Lr=20D,Ly=20D

(12)

腔外网格为810×300，腔内网格为240×120，并在腔壁附近加密以捕捉边界层。空腔流动会产生强烈的声波，为避免声波在边界的反射污染流场，在计算域边界设置海绵层[18]吸收声波，从而获得准确的高精度流场数据。海绵层的长度参数为

Li=3D,Lo=10D,Lt=10D

(13)

本文针对亚声速流动开展研究，分析所用的流场物理量均为无量纲量。空腔的流动参数设置为Ma=0.8，Re=2 500，Pr=0.7，其中雷诺数基于空腔深度D。

(a) 二维空腔流动计算区域示意图

(b) 腔体附近网格分布图1 二维空腔流动计算域及网格Fig.1 Computational domain and mesh of the open cavity

2.2 计算结果及验证

本文直接数值模拟所用代码已在参考文献[13，19]中得到很好验证。针对本算例，首先进行网格收敛性验证。考虑一套加密网格，其中腔外网格为1 245×470，腔内网格为400×200，并比较两套网格下空腔底部(0.5D，0)和空腔外部(0,2D)处的压力随时间变化。两套网格下压力变化基本一致，如图2所示，可以看出：幅值和相位差别均很小，说明本计算所用网格密度已经足够。

(a) 空腔底部x=0.5D处压力随时间变化

(b) 空腔外x=0，y=2D处压力随时间变化图2 网格收敛性验证Fig.2 Validation of grid convergency

对空腔(0.995D，D)处的速度和压力脉动采样分析，结果如图3所示，流动呈现出强烈的周期性，且无量纲时间周期为T=3.75。对包含5个流动周期的468个流场片段进行快速傅里叶变换(FFT)，采样点的频谱如图4所示。

图3 速度及压力时间演化Fig.3 Time evolution of velocity and pressure

图4 空腔压力及速度脉动频谱分布Fig.4 Frequency spectrum of velocity and pressure perturbation in the open cavity

从图4可以看出：主频St=0.668。

Rossiter给出的半经验公式[20]

(14)

可见，图4得到的主频与式(14)中的第2模态即St2=0.686吻合较好，而且与K.Krishnamurty[21]的实验结果频率St=0.656吻合良好。Krishnamurty的实验纹影与直接数值模拟获得的数值纹影对比如图5所示。

(a) Krishnamurty实验纹影

(b) 数值纹影图5 实验纹影图与DNS计算所得数值纹影对比Fig.5 Comparison between experimental Schlieren and Numerical Schlieren from DNS

从图5可以看出：两者在定性上吻合良好，波系结构较为一致。

3 涡动力学分析

3.1 流场基本结构

一个周期内的流线结构如图6所示，空腔内部主要有四个回流区，它们构成了流场的主要拓扑结构。在t=0.2T时刻，空腔尾缘点处涡量等值面撕裂，在t=0.8T时刻在前缘开始形成闭合流线结构，这些特征在欧拉框架下通常被看作涡破裂和涡生成。然而流线和涡量都不具有客观性，所揭示的流场结果会随参考系的变化而变化。另外瞬时流态不能揭示涡结构的动态演化过程，其动力学特性也无法揭示。因此接下来本文采用拉格朗日方法对流场涡结构进行动力学特性分析。

(a) t=0

(b) t=0.2T

(f) t=T图6 空腔流动的流线及涡量分布Fig.6 Streamlines and vorticity distribution in the cavity

3.2 拉格朗日涡动力学分析

对空腔流动的流场进行FTLE的计算，积分时长是影响计算结果的重要参数，通常随着积分时长的增加，FTLE的嵴所提取的LCS结构更加锐利。因此在本文计算中，积分时长选取为一个时间周期T=3.75，使得FTLE的嵴能够准确捕捉到流场中的LCS。一个周期内不同时刻的pLCS(红色)和nLCS(蓝色)如图7所示，pLCS为排斥性结构，附近的流体粒子被其排斥从而流向动力学特性不同的区域；nLCS则为吸引性结构，在时空演化过程中吸引周围来自不同动力学特性区域的流体粒子向其聚集靠拢，因此pLCS和nLCS是不同动力学特性区域的边界。涡结构作为流场中的拟序结构，在涡边界内外，流动的特性存在差异，因而涡边界可用LCS进行识别。将pLCS和nLCS包络的闭合区域作为涡结构，两者交叉处的鞍点作为追踪涡结构运动的特征点来分析拉格朗日涡的动力学特性。

pLCS和nLCS包络区域形成的涡结构主要有三个：空腔前缘点剪切层卷起后形成的涡，空腔后缘点处即将发生撞击和撕裂的涡以及空腔右侧的主涡。这与流线所揭示的四个主要回流区有所不同，位于空腔右下拐角处的回流区Ⅱ和位于空腔左侧的回流区Ⅲ未被LCS揭示出来。由于LCS揭示的是拟序结构，能在一定生命周期内保持动力学及结构特征稳定，而回流区Ⅲ即使从流线观察不能保持稳定的形态且无清晰边界，因此未被LCS识别出来；回流区Ⅱ位于空腔右侧底部拐角处，速度较低，其靠近壁面的边界由于无滑移边界速度几乎为零，吸引性或排斥性较弱，因而只有与主涡相交的边界被LCS揭示出来，而完整的涡结构未被识别出来。流场中Q值大于零的区域如图7中灰黑色区域所示。Q值揭示的涡区域也主要有三个，即空腔右侧的主涡和剪切层中卷起和向下游发展的两个涡，与LCS包络所揭示的涡结构覆盖的区域基本一致。空腔右下拐角处的角涡尽管被Q值所揭示，但强度较弱，且形态和边界不够清晰稳定。流线揭示的回流区Ⅲ同样未被Q所识别如图6所示，表明空腔左侧尽管存在瞬时封闭流线，但不能保持动力学稳定性及清晰边界，该回流区不能识别为涡结构。综上，LCS与Q值都清晰地揭示了空腔流动的流场中存在的主要涡结构，而空腔右侧的角涡则由于吸引性或排斥性较弱，未被LCS识别出来。

空腔右侧的主涡在整个流动周期中形态基本维持不变，绕涡心顺时针旋转，周期与流场周期一致。在运动过程中，空腔后缘点处剪切层中的涡撕裂后，一部分流体被nLCS及pLCS包络并随主涡运动进入到主涡中；同时主涡靠近剪切层一侧流体被pLCS吸引靠近剪切层进入主流，随主流对流至下游。

剪切层中的涡结构则经历了复杂的生成、增大、脱落、对流、破裂等过程，其完整周期为流场中压力和速度脉动周期的两倍，而且剪切层内两个涡结构的完整运动动力学过程和特性是相同的。在图7中，t=0.4T时刻如图7(c)所示，剪切层中的涡结构在空腔前缘点开始生成并逐渐增大直至t=T时刻如图7(c)～图7(f)所示，从前缘点脱落进入到剪切层中，并随主流一起向下游对流。在对流过程中如图7(a)～图7(e)所示，涡结构发生形变扭曲，并在t=T+0.8T时刻如图7(e)所示，开始撞击空腔后缘点，整个撞击过程中如图7(f)、图7(a)、图7(b)、图7(c)所示，涡结构不断被空腔后缘点撕裂，最终在t=2T+0.4T如图7(c)所示，被撕裂为两部分，一部分随主流在空腔外向下游流动，另一部分则进入主涡部分随主涡在空腔内部流动。

(a) t=0

(b) t=0.2T

(d) t=0.6T

(e) t=0.8T

(f) t=T图7 LCS及Q值分布Fig.7 LCS and Q value distribution in the cavity

前缘点和后缘点附近区域的涡通量为

(15)

式中：ω为涡量。

剪切层中涡的生成、脱落、撞击和破裂时刻的量化分析通过式(15)的时间演化特性分析给出，分析结果与LCS分析给出的特征时刻相互对照验证。选取的前后缘点附近区域分别如图7中A、B所示，两个时间周期内的涡通量时间演化如图8所示，注意在两个区域内涡通量都是负值。由A区域内的涡通量时间演化可知：在t=0.35T时刻，涡通量开始逐渐积累增加，意味着涡结构开始形成；在t=0.87T时刻涡通量达到极大值并开始减小，意味着此时涡结构开始从前缘点脱落。由B区域内的涡通量时间演化可知：在t=0.76T时刻，涡通量开始迅速增加并在t=1.12T时刻达到极大值，意味着从上游而来的涡结构开始对流至后缘点处并撞击后缘点，随后涡通量开始减小，直至t=1.43T时刻涡通量减小至极小值，此时涡结构撕裂成为两部分，一部分随主流对流而下，另一部分则进入空腔中的主涡内。以上基于涡通量的涡结构动力学特性分析与LCS分析结果基本一致，印证了基于LCS分析涡动力学的准确性。