基于动力学模态分解法的绕水翼非定常空化流场演化分析1)

谢庆墨 陈 亮 张桂勇 孙铁志

(大连理工大学船舶工程学院,大连 116024)

(哈尔滨工业大学机电工程学院,哈尔滨 150001)

(工业装备结构分析国家重点实验室,大连 116024)

引言

当局部压力降低至饱和蒸气压时，液体发生由液相转变为气相的现象称为空化现象[1].空化现象常常出现在高速船舶的舵、螺旋桨、水翼等高速运动的流体机械表面以及高速出入水的航行体表面[2-4].空泡的周期性脱落和溃灭会产生压力脉冲和辐射噪声[5].一方面，这会使流体机械表面出现明显的振动和材料剥蚀[6]，另一方面，当压力脉动频率与桨自身的固有频率接近时，还会引发共振，对流体机械的水动力性能带来很大的影响[7].

空化的产生及其流体动力学特性一直是研究的重点[8-10].Knapp 等[11]首先对片云空化进行了广泛的研究，观察了片云空化从片云空化腔中剥离/脱落的过程.Leroux 等[12]通过压力传感器试验研究了NACA66 水翼周围空化的不稳定性，并进行了数值模拟，结果表明回射流是导致空泡脱落的主要原因.王雅赟等[13]采用不同的湍流模型模拟二维水翼空化流动，发现空泡尾部大尺度旋涡的形成导致了空泡的断裂，为后人进行数值模拟时应该选用哪种湍流模型提供了有效的参考.除了理想流体，流体的可压缩性与热力学效应也会对流场产生影响[14].王畅畅等[15]模拟了可压缩流体的空化现象，对空泡溃灭时的激波进行了重点研究，详细说明了大尺度云空泡的溃灭过程.孙铁志等[16]研究了不同空化模型在空化流场数值模拟中的影响，并考虑了热力学效应对空化流动带来的影响，为低温流体的空化现象提供参考.随着计算机运算能力的发展，更加复杂的泄漏涡空化也开始成为研究的热点[17].彭凯等[18]研究了轴流泵顶隙附近的涡流结构，采用滑移网格技术模拟了浆叶的旋转，并对旋涡与旋涡、旋涡与壁面之间的相互作用进行了分析.Cheng 等[19]使用大涡模拟(LES)捕获不稳定的尖端泄漏空化流，成功模拟出了尖端泄漏涡、分离涡、诱导涡这些复杂的涡结构，并将尖端泄漏空化流的空间演化分为3 个阶段,揭示了空化对涡度和湍流的影响机理.

由于空化具有强烈的非定常性，并伴随有复杂的旋涡结构[20].空泡与壁面、空泡与旋涡之间具有强烈的相互作用，是空泡水动力学研究的热点[21-23].空化流场这一复杂的流动结构具有多尺度特性.目前的研究大多都在对整体流场进行分析，若可以将空化流场分解成若干个不同尺度的流动结构，将对进一步分析空化流动机理提供参考.近年来，动力学模态分解(dynamic mode decomposition，DMD)已逐渐应用于计算流体动力学，被认为是评估非线性波动数据的重要分析工具[24].DMD 方法可以有效地将流场分解成若干个模态，并识别出主导模态进行分析.且具有模态频率和增长率单一的特点，在对周期性流动的分析中有较大优势.目前形式的DMD 算法最初由Schmid[25]提出，Rowley 等[26]进一步对其进行扩展，他们将DMD 解释为Koopman 分析的一种形式.寇家庆和张伟伟[27]首次对现有的DMD 方法与类似的本征模态分解(POD)进行了比较，阐述了DMD 方法在流体力学分析领域的优势.并总结了现有的DMD 方法以及应用领域，对DMD方法的现状与前景作了深入讨论.对于空化现象的DMD 分解，韩亚东和韩磊[28]使用DMD 方法研究了文丘里管的空化流动，讨论了各阶模态的选取原则以及模态的特征.Liu 等[29]使用DMD 方法研究了高雷诺数下的水翼空化流动特征.

目前的研究工作来看，DMD 方法用于空化流场研究的取得了令人鼓舞的进展，但基于DMD 方法探究非定常空化流场特性仍需要进一步研究.本文采用RANS 方法，结合Schnerr-Sauer 空化模型求解非定常的水翼空化流场，获取流场数据，并利用动力学模态分解方法，研究了水翼空化流场各阶模态的主要特征.同时对前四阶模态进行重构，得到前四阶模态的对应流场，研究了各阶模态的流动特征.最后对不同空化数下模态2 的对应流场进行了比较，研究了空化数对空泡脱落产生旋涡对流场的影响机理.

1 计算模型

1.1 基本控制方程

均相流模型是空化流场数值模拟使用的常规模型，其质量守恒和动量守恒控制方程为

其中，p是压力，ui是沿第i方向的速度.混合密度ρ 为

其中，αl是水体积分数，αv是水蒸气体积分数.

1.2 空化模型

空化模型通常由输运方程表征，本文采用Schnerr-Sauer 模型[30]，其输运方程为

其中，αv是水蒸气体积分数，源项m+和m−表示为

其中，m+为蒸发源项，m−为凝结源项，pv为当地饱和蒸气压.Rb为空泡半径，表示为

其中Nb=1.0×1013.

1.3 湍流模型

本文拟采用SSTk-ω 湍流模型，SSTk-ω 模型由Menter 提出[31]，具体方程如下

其中，α，β为经验系数，α=0.44，β=0.082 8；k代表湍流的动能，ω表示比耗散率，Pk表示湍动能的产生项，Dω为扩散项，µt为涡黏系数.它们分别定义为

其中，F1，F2为混合函数；S为剪应力张量的常数项；σω,2为经验系数.

1.4 动力学模态分解

采用DMD 方法进行模态分解时，首先需要将非定常流场的时间序列以快照序列矩阵的方式呈现，快照数据来源为物理试验或数值仿真，且任意两个快照之间的时间间隔均为∆t，即

线性动力学矩阵A被定义为非线性映射的近似，该矩阵控制系统从一个时刻到下一个时刻的演变

由于A的维度很高，难以直接计算，因此需要通过降维的方式从数据序列中计算出A.因此，通过将最终快照表示为前几张快照的线性组合，建立一个伴随矩阵S作为所需的近似

通过友矩阵变换，S的形式为

由于S中的未知量仅为a向量，即可以计算使残差r最小的a来求得S

当残差r很小时，S的特征值就近似于A的特征值.S矩阵可以看作是A矩阵的低维形式，因此，S的特征值能够代表A矩阵中的主要特征值.对S矩阵进行特征值分解

然后给出DMD 模态ϕi

则任意时刻的流场快照可以通过前m个流场快照表示

其中，m为DMD 模态数，µj为模态D的特征值，ϕj为模态ϕi的列向量.

在流体力学的分析过程中，往往只能够获得数据的实数部分，如速度，压力等，对于实数数据进行DMD 分解，所得模态的特征值是共轭的，他们所对应的特征向量也是共轭的.流场数据的虚部，即流场变化的相位信息，在标准DMD 中并没有体现，因此会造成一定误差[25].因此，在进行标准DMD 方法前，使用了希尔伯特--黄变换[32]对原始数据进行预处理.得到原始数据波动的相位信息，与原始数据叠加得到复数形式的数据.通过希尔伯特--黄变换可以得到更精确的原始数据，使DMD 分解结果更加准确.其中数据的实部代表流场的物理信息，例如速度、压力，虚部代表这些物理信息的瞬时相位

其中，xi(t)为第i个探测点的时间进程,IMFi(t)为Huang 经验模态分解法中的本征模态函数.

1.5 模型设置

本文模拟的模型为NACA66 水翼，该翼型的最大厚度位于距前缘45%处，厚度与弦长的比值为12%.最大拱度比位于距前缘50%处，弧度与弦长的比值为2%.来流攻角α=6◦，弦长C=0.150 m.计算域与网格如图1 所示，设水翼弦长为C，进流段长度为2C，出流段长度为5C，上下壁面相距2C.采用C 型网格，为了提高计算精度，在水翼周围及尾流段进行加密，网格数为38 750.边界条件设置为：入口处为速度进口，出口处为压力出口，在水翼表面与上下壁面均施加无滑移条件.流动速度为V=5.33 m/s，其对应的雷诺数为Re=8.93×105.出口压力由空化数σ 来确定，初始温度设置为298 K，时间步长为1.0×10−4s.

图1 计算域及网格划分Fig.1 Computing domains and mesh generation around the hydrofoil

2 数值计算方法验证

为验证数值模拟的有效性，将空化数σ=1.00的结果与试验[12]进行对比，所有条件设置与试验一致.图2 为一个脱落周期内水翼空泡形态的试验与数值模拟对比.可以看出，数值模拟可以有效模拟出空泡周期性演变中发展、脱落、溃灭的特性，且数值模拟的形态与试验结果吻合较好.

图3 为翼型吸力侧中部x/C=0.6，x/C=0.7 两个压力监测点压力值的数值预报与试验数据对比.从图中可以看出，各观测点的压力存在周期性波动，这与空泡周期性生成和脱落有关，其中压力值较低的区域处在空化区域，压力值高主要是空化溃灭引起的.但每个周期压力波动的峰值都不相同，体现了空化现象的强非定常性.对比试验和数值结果可以看出，试验和数值结果在振荡频率和压力波动的大小上都有很好的一致性.综上所述，本文的空化流场数据与试验结果吻合较好，从而验证数值方法的有效性.

图2 数值预报的空泡形态周期演化与试验结果对比Fig.2 Comparison of the cavity evolution between numerical and experimental results

图3 数值预测的水翼吸力面中部压力值与试验结果[6]对比Fig.3 Comparison of the pressure fluctuation on the suction surface of hydrofoil between numerical and experimental results

3 结果与分析

3.1 水翼空化流场动力学模态分解

本文对空化数σ=1.00 的空化流场中的速度场进行DMD 分解，图4 表示特征值在复域上的分布情况.与标准DMD 有所不同，改进后的DMD 方法的不以共轭的形式存在，而是单独存在.黑色圆圈为单位圆，特征值点与单位圆的关系揭示了对应模态的非稳态特征.其中特征值位于单位圆外代表该模态是发散的，位于单位圆内代表该模态是收敛的，位于单位圆上代表该模态稳定.由于原始数据经过了希尔伯特--黄变换，还原出的原始数据虚部为负值，因此流场数据的有效信息大部分集中在虚部为负值的区域.虚部为正值的区域为对应频率大，振幅小，这是由于希尔伯特--黄变换存在一定的误差，这个误差存在于数据的开始部分，所占的比例很小，在此可以忽略[32].结果表明，在虚部为负值的区域，大多数特征值点都在单位圆的范围内，这说明该流动条件下的空化流动是周期性的.图4 给出了选取的前四阶模态，均在单位圆上，表1 给出了前四阶模态对应的频率和放大率，其中模态的频率和放大率可以由下式求出

其中，gj为模态的增长率，ωj为模态的对应频率.由于前四阶模态的放大率接近于0，因此选取的前四阶模态都是稳定的.前四阶模态的选取方法将在下文进行详细论述.

表1 前四阶模态对应频率和放大率Table 1 Frequency and magnification of the first fourth modes

图4 模态特征值的分布Fig.4 Distribution of mode eigenvalues

动力学模态中包含的能量高低是评价模态在流场中主导作用的重要依据.模态中的能量高低可由模态ϕi的范数表示[26]

图5 动力学模态的能量分布Fig.5 Energy distribution of dynamic mode

图6 表示DMD 分解得到的前四阶动力学模态.模态值只取实部，虚部与实部主要体现在相位上的差别，在分析结果时一般可忽略.模态1 对应于频率为零的平均流场.模态在水翼表面以及水翼尾部的上方，体现了片状空泡与云状脱落空泡产生的区域.模态2 的特征主要集中在水翼的前端，结构与数值模拟中的空泡演化形态相似，体现了空泡脱落的特征结构.同时也表明片状空泡的周期演化是空化现象呈现周期性的主要原因.因此模态2 可以描述出空泡的发展趋势以及周期性演化特征.模态3 的结构具有相反的正负值，且呈现交替分布的特征.这表明流场中存在明显的相干结构，图中正值区域为空泡的脱落点，负值区域描述了空泡的溃灭阶段的流动特性，这意味着当一个结构由于空化而生长，另一个结构由于传输和凝结而减弱.模态4 的相干结构尺度更小，这说明模态4 捕捉了能量更少，占比较小的结构，对应的频率为空泡脱落频率的3 倍，体现了空泡云脱落现象中的高频特征.随着频率增加，模态的振幅减小，相干结构的尺寸也变小.

图6 DMD 前四阶动力学模态Fig.6 The first 4 dynamic modes of DMD

3.2 基于模态分解的空化流场分析

为了进一步揭示DMD 分解中各阶模态表达的物理意义，对空化数σ=1.00 的工况进行DMD 分解，取前四阶模态特征值使用重构公式，还原出单阶模态对应的流场.其中模态1 为平均流场，对应流场不随时间变化而变化.图7 表示模态2,3,4 对应的流场，由于不同模态对应的特征频率不同，因此对应的周期表示为T1,T2,T3.对于模态2、模态3、模态4，相对应的，流场描述了空化在生长的过程中水翼后缘的大尺度旋涡与回射流的形成，旋涡的周期脱落，以及水翼尾缘处形成的旋涡.观察流线图，在模态2 中，旋涡Ⅰ先从水翼吸力面头部开始形成，在向尾部移动的过程中不断生长，移动到接近尾缘时，尾缘处开始产生新的旋涡Ⅱ.旋涡Ⅰ在脱离水翼后停止移动，并开始向下游水平分裂出新的旋涡Ⅲ，在旋涡Ⅲ完全形成后，旋涡Ⅰ退化成节点，旋涡Ⅲ继续向下游移动.旋涡Ⅱ形成后不断向下游运动.在模态3 中，水翼吸力面形成的旋涡Ⅰ在脱离水翼后，与水翼尾缘形成的旋涡Ⅱ之间存在融合行为.融合的过程中，旋涡Ⅰ逐渐缩小，旋涡Ⅱ在下游移动的同时向上移动，形成新的旋涡Ⅱ继续向下游移动.模态4 中同时存在多个旋涡，其旋涡行为与模态2 和模态3 基本相似，这揭示了两阶模态的旋涡在运动的过程中存在高频行为.

图7 模态2，3，4 的对应流场流线图Fig.7 Flow field streamline of mode 2,3 and 4

3.3 不同空化数下的重构流场分析

为了对不同空化数下的空化流场进行比较，分别将DMD 方法应用于σ=1.25，σ=0.75，σ=0.50的流场数据.图8 表示3 种不同空化数条件下模态2 的对应流场.对比不同空化数流场可以看出，在模态2 中，不同空化数下的旋涡运动模式基本一致.都包含主旋涡向下游运动，主旋涡分裂形成新的小旋涡.不同的是，随着空化数的减小，脱落的主旋涡的尺寸增大，这说明空化数的减小使流场变化更加剧烈，对水翼的影响更强.

图8 σ=1.25，σ=0.75，σ=0.50 工况下模态2 的对应流场流线图Fig.8 Flow field streamline of mode 2 for σ=1.25，σ=0.75，σ=0.50

表2 表示各阶模态对应的频率值，对于不同空化数下的流场，模态1 的频率均为0 Hz，这表明平均流场均为流场中的主要结构，占据着流场中的大部分能量，从模态2 的频率变化可以看出，频率随着空化数的减小而减小.这主要是由于随着空化数的降低，空化现象逐渐剧烈，空穴长度逐渐增长，导致空化的生长、脱落时间变长.模态3、模态4分别为模态2 的2,3 倍频，这是对流场中的主要信息进行补充描述.

表2 各阶模态频率Table 2 Frequency of each mode

4 结论

本文针对水翼开展非定常空化流动数值预报，并利用动力学模态分解方法对空化流场进行分析，得到如下结论：

(1)利用DMD 方法对速度场进行分解，得到各阶流场模态.第一阶模态表示平均流场，其频率为0 Hz，是流场的主要结构.第二阶模态与空泡周期性脱落频率相同，与空泡脱落的动力特性有关.第3,4 阶模态频率分别为第二阶模态的2,3 倍频，这表示流场中的高频行为.

(2)通过观察各阶模态对应流场一个周期内的流场演化，结果表明，第二阶模态描述了旋涡1在水翼吸力面生长，脱离水翼后形成新的旋涡并向下游移动.旋涡2 在水翼尾缘形成并向下游移动.第三阶模态揭示了旋涡1 与旋涡2 在向下游移动的过程中存在融合行为.第四阶模态中存在更多的旋涡，旋涡运动较为复杂，包含旋涡分裂、融合、退化.

(3)各模态对应流场中的旋涡结构与DMD 方法得到的模态特征值呈现的相干结构相对应.与空泡脱落频率相同的第二阶模态反应了流场的主要旋涡结构.旋涡结构的尺寸随着模态阶数的增大而减小.

(4)随着空化数的降低，空泡脱落频率逐渐减小.模态2 的主旋涡尺寸随着空化数的降低而增大，说明空化数的降低加剧了流场的演化，形成的旋涡更大.

本文通过采用动力学模态分解方法对绕水翼空化流动进行了分解、分析以及重构，重点通过DMD 方法的运用，进一步获得了空化流场流动结构特性.虽然取得了一定的研究进步，但是在结合空泡形态演化和动力学模态方法之间的内在关联机制仍略显不足.此外，在低空化数条件下，RANS 方法存在过度预报湍流黏性问题，导致对空化瞬变过程的预报存在一定的不足.在今后的研究中将提升对非定常空化流动的预报，从而进一步深入开展相关研究.