任万龙，赵志范，娄珊珊，高原，郝宗睿

(齐鲁工业大学(山东省科学院)山东省科学院海洋仪器仪表研究所，山东 青岛 266001)

空腔流场易产生振动和噪声，是航空航天领域一类受到广泛关注的流动问题[1]。在高速水下航行器及水力空化领域，空腔流动由于其速度较高，压力降低，极易诱发空化，影响水下航行器的航行性能。

孙亚琴等[2]通过实验和理论分析，研究了水平空腔旋流压强的产生机理及构成，分析了其影响因素，探究了压强的变化规律。韩帅斌等[3]运用涡动力学理论，结合拉格朗日拟序结构识别流场中旋涡结构，研究了空腔流动中的涡旋结构演化规律，包括碰撞、合并等过程。万振华等[4]分析了空腔中涡结构碰撞后的压力脉动分布，运用拟涡能的概念分析了空腔流动。曲俊超等[5]在空腔前部添加了微射流，改变了空腔流动结构，一方面减小了水下航行阻力，另一方面降低了航行器的水动力噪声。刘哲等[6]采用动态模态分解的方法，探究了不同空腔结构的流态分析，重点解析了后壁面倾斜角对流动的影响。尽管有诸多学者开展了空腔流动的研究，然而针对高宽比对空腔流动的研究较少。

在空腔流动的数值计算中，常用的数值模拟方法有雷诺时均方程(Reynolds-averaged Navier-Stokes,RANS)、大涡模拟(large edy simulation,LES)和直接数值模拟(direct numerical simulation,DNS)。其中DNS能最有效地获得流场内的湍流细节，但是消耗资源，计算成本高，计算时间长；而RANS方法通常会过滤掉很多流场细节，导致计算流场不够精确[7]。因此，通常采用LES方法开展空腔流动流场数值模拟。

本文以空腔流场为研究对象，采用LES开展空腔流动的数值计算，重点分析高宽比对空腔流动的影响，探究监测点压力脉动的时频特性，获得空腔内固定位置的速度分布和压力分布。结合涡量分布，获得高宽比对空腔流动的具体影响，为后续开展空腔流动的应用奠定基础。

1 模型及方法

大涡模拟方法的控制方程如下[8]：

(1)

(2)

式中，ρ为液体密度，uj为不同方向的速度，μ为黏性系数，p为流体压力。

大涡模拟方程对质量和动量方程进行Favre滤波得到：

(3)

(4)

亚格子应力为：

(5)

涡黏度模型是一种广泛使用的亚格子网格尺度模型，其假设是子网格尺度应力与过滤后的大流量应变率张量Sij的模量成正比[9]。

(6)

(7)

湍黏性系数μt是应用壁面局部涡黏度模型(WALE模型)封闭的亚格子尺度湍流黏度，WALE模型可以从层流剪切返回值[10]。

(8)

式中，LS和Sijd计算公式如下

LS=min(kd,CSV1/3)，

(9)

(10)

式中，k为冯卡门数，d是距离最近壁的距离，V是计算单元体积，Cs为默认WALE常数，通常取0.5。

2 空腔模型及边界条件

空腔计算流域如图1所示。空腔的长度为L，高度D，宽度W，空腔前缘长度LI，空腔后部LA。为研究后续高宽比对空腔性能的影响，选择3个典型工况：D/L分别为0.5、1.0、1.5。来流速度选择30 m/s。基于腔长的来流参数的雷诺数ReL=ρ∞U∞L/μ∞=1.2×105，其中U∞为空腔来流速度，网格分布在各个壁面及空腔开口处的剪切层附近加密，壁面第一层网格保证Y+≤1。进口设置为速度进口，出口设置成为压力出口。

图1 空腔流动计算流域Fig.1 Calculation domain of cavity flow

3 计算结果及验证

3.1 网格验证

由于在数值计算过程中，网格的划分直接影响到计算结果，本文开展了网格无关性验证。以D/L=0.5的空腔为模型，将计算网格划分了粗、中、细三套网格，Case1空腔网格数为100×200，Case2为150×240，Case3为170×270。网格划分形式如图2所示，比较了不同网格数量下，分析空腔内(0，0)的压力随时间的变化，如图3所示。通过计算发现，Case2和Case3的网格压力脉动变化基本一致，网格的数量对计算结果影响较少。考虑计算成本的问题，本文选择Case2的网格作为计算网格，开展后续研究。后续针对不同D/L的空腔，按比例对空腔网格进行划分，使其满足相关精度要求。

图2 网格划分形式Fig.2 Grid generation form

图3 时域压力脉动Fig.3 The pressure fluctuation of time domain

3.2 流场内压力脉动

在空腔中心位置P1(0，0)设置监测点，通过计算，获得不同高宽比空腔中非定常流动，记录压力脉动，获得如图4所示的压力脉动分析。通过对比分析，发现高宽比D/L=0.5时，压力脉动呈现明显的周期性分布，当D/L=1.0和D/L=1.5时，其压力脉动周期性不明显，在大周期范围内，还存在一定的小周期脉动。这说明高宽比的不同，对空腔内P1的压力脉动影响较大，当高宽比增大时，空腔流动在压力脉动范围内增大，且周期性不明显。通过对压力脉动的傅里叶变换，分析得到频域曲线，从图4中发现，当高宽比增大时，峰值频率增加，这说明在空腔口位置处的压力脉动更加剧烈。

图4 P1点的压力脉动及其频域分析Fig.4 The pressure fluctuation at P1 and the frequency domain analysis

3.3 时均速度分布

在计算过程中，截取y=0.001位置的曲线，线上的时均速度分布、时均压力分布，如图5所示。在y=0.001位置距离空腔口较近，截取此曲线位置的时均压力脉动和时均速度，能够有效反映空腔范围内的流场特征。图5显示，低高宽比的空腔在y=0.001位置处的压力脉动振荡剧烈，在靠近来流方向压力脉动低，考虑去流方向的位置，压力脉动大，说明此处涡旋结构复杂。高宽比等于1.0和1.5的空腔压力脉动较为平缓，波动较小。图5(b)展示了y=0.001位置的速度分布，从图中可以看出，D/L=0.5的速度波动较大，空腔速度分布呈现一定的周期分布，而其他高宽比的空腔速度波动较小，与前面压力脉动的分布正好对应。

图5 y=0.001流场特性分析Fig.5 Analysis of flow field characteristics at y=0.001

3.4 速度场分布

计算过程中获得不同高宽比空腔的速度场分布情况，如图6所示。从图中可以看出，随着高宽比的增大，流场内速度分布呈现不同趋势。高宽比为0.5时，由于在流向方向上距离较大，速度分布沿流线分布，在腔内分布存在明显的低速区，表明了在高宽比较小时，空腔内速度变化较大，有明显的低速区，内部流动非定常特性比较明显；高宽比为1.0时，垂直于流向方向上，呈现明显的两个大小相等，方向相反的速度分布，表明空腔内部存在两个明显的涡旋结构，内部流动规律明显。随着高宽比的继续增加，速度分布区域呈现3个明显交替的区域，这表明空腔区域内部具有速度梯度分布，呈现涡心数量增加。

图6 x方向速度场对比Fig.6 Comparison of x-direction velocity field

3.5 涡量分布

通过仿真计算获得不同高宽比下空腔的时均涡量分布情况，具体如图7所示。从图中可以看出，D/L=0.5的空腔整个腔内，涡量分布较大，与前面的速度分布和压力分布对应。当高宽比较小时，空腔内流动活跃，内部出现很多涡旋结构，流动较为湍急，非线性特征明显。D/L=1.0时，涡量值集中在空腔开口域和空腔内部，形成固定区域，空腔开口域涡量值较大，随着空腔高度的增加，其涡量值减小，集中于空腔中心区域。流体流过空腔后，在空腔的后缘形成系列涡旋，且P1位置的压力脉动形成波动较小的周期性流动。D/L=1.5的空腔其涡量分布与D/L=1.0相似，随着空腔高度的增加，空腔内的涡量值和区域都减小了，空腔内的涡量高的值集中在空腔口附近及空腔内部的上边缘位置，由于空腔口的作用在后缘同样形成一系列涡量，相比其他两种空腔流动，D/L=1.5的空腔后缘涡量值明显小于前面几种，这说明空腔高宽比对空腔内部流动影响较大。空腔开口位置是涡量的集中区，腔内形成明显的圆形涡量分布，意味着在内部形成较为规则的环形涡流，其涡量分布与前述x方向速度分布相对应。

图7 涡量场对比分析Fig.7 Comparison of vorticity field

4 结论

通过数值计算，对比了ReL=1.2×105情况下，高宽比对空腔流动的影响，通过分析获得如下结论：

(1)高宽比的增加，空腔开口位置的压力脉动周期性降低，周期性不明显。空腔内压力脉动峰值降低，减弱了压力脉动。

(2)高宽比增加，速度和压力的波动减小，内部涡旋结构较为稳定，呈现独立较大的涡旋结构，而高宽比较小的空腔流动，内部流动复杂，具有多个小的涡旋结构，脉动性较强。

(3)时均涡量分布表明，随着高宽比的增加，涡量分布较为单一，集中在空腔开口域附近及空腔内部，且整体的脉动值较小，而高宽比较小的空腔流动，涡量值较大，且分布区域广，说明其内部流动脉动性较强。

经过本文的计算分析，发现不同高宽比的空腔流动特性变化较大，高宽比较大的流动较为平稳。因此在水下航行器及船舶设计过程中尽量避免出现小的空腔，以免增加航行阻力，该研究为后续船舶水下开孔等系列空腔流动奠定了一定的理论基础。