基于动态模态分解的叶道涡非定常解耦与重构

童哲铭 辛佳格 童水光

1.浙江大学流体动力与机电系统国家重点实验室，杭州，3100632.浙江大学机械工程学院，杭州，310063

0 引言

大力发展可再生能源是当前国际能源大转型的共识和行动，也是我国实现“双碳”目标的战略选择。水电是全球能源系统的重要组成部分，它的发电量在可再生资源中所占比例最高，占世界发电量的16%～17%、世界可再生电力的50%[1]。水轮机是水力发电的核心部件，其中，混流式水轮机结构紧凑、效率较高，占所有水电装机容量的60%以上，是目前应用最为广泛的机型[2]。水电的角色从提供基础负荷逐渐变成提供高度可调节负荷能源，水轮发电机组将更多地运行在偏负荷工况下以平衡电网参数[3-4]。混流式水轮机在偏负荷工况下运行时，由于转轮内有限的空间及负冲角的综合作用，水流在转轮进口处发生脱流或二次流，在转轮下环处形成叶道涡[5]。叶道涡会诱发叶片振动与机组的异常噪声，严重影响机组的运行安全。

深入研究混流式水轮机叶道涡的时空分布规律及涡结构的演化机制，对提高机组运行稳定性具有重要意义。LIU等[6]通过试验观察发现，部分负荷工况下的尾水管回流是产生叶道涡的原因。王钊宁等[7]基于数值模拟发现，水轮机转轮不能适应流量的较大变化会导致转轮内稳定连续的压力梯度被打破，这可能是叶道涡形成的因素之一。YAMAMOTO等[8]基于多相流数值计算与试验测试发现叶道涡附着在上冠表面且呈螺旋状在流道内发展，还发现旋涡中心出现的低压区诱发了旋涡空化。XIAO等[9]的研究表明，叶道涡诱发了低频压力脉动，压力脉动频率为转频的0.2～3.0倍。

通过过去的研究可以发现，叶道涡的流动结构在时空演化过程中具有显著的非定常特性，包含着尺度不一、形态各异的相干结构。为进一步理解混流式水轮机叶道涡的演变机制，需要采用低维系统分解技术提取各尺度的相干结构进行分析。本征正交分解(proper orthogonal decomposition，POD)和动态模态分解(dynamic mode decomposition，DMD)作为数据驱动的多尺度相干结构时空特征识别算法已被广泛应用。POD方法基于对整体流场的能量贡献来提取流场主导子结构，已应用于多种流动现象的研究[10-12]。SCHMID等[13]基于动力系统的Koopman分析提出了DMD方法，利用高次多项式拟合低维子空间(低维子空间的特征值和特征向量可表征流场的主导相干结构)。与POD方法相比，DMD方法能在动力学层面上描述各个物理量随时间的演化过程，并揭示不同频率下的流场时空特征，近年来已被引入流场分析领域。何昌鸿[14]应用DMD方法提取了液下式离心泵蜗壳内流场的主要模态，分析了原始流场在不同特征频率下的流动结构。LIU等[15]应用DMD方法揭示了叶轮内流场的径向不均匀性，以及扩压器内的回流是主要的相干结构，对比降阶重构、全阶重构和原始数据，发现降阶重构的流场能有效反映主流特性。HAN等[16]发现DMD方法能准确分解混流泵的叶尖泄漏涡演化的主频及其谐波频率，并能在平均流态和第一模态的基础上重构流场，对比POD方法与DMD方法后，发现 DMD 方法可应用的流动范围更加广泛。

笔者引入DMD方法对混流式水轮机在偏负荷工况下运行时的瞬态流场进行解耦重构，得到叶道涡的相干结构，研究叶道涡在时空演化过程中的非定常流场结构及其对应的物理机制，在不额外建立控制方程的前提下，降低计算成本。同时，将DMD方法与POD方法进行比较，探讨两种分解方法在分解混流式水轮机叶道涡相干结构和重构复杂流场方面的准确性和有效性。

1 物理模型和数值方法

1.1 物理模型

笔者将某混流式水轮机模型作为研究对象，如图1所示，水轮机由蜗壳、固定导叶、活动导叶、转轮和尾水管构成，叶片、固定导叶和活动导叶的数目分别为13、12、24，额定容量58MW。混流式水轮机活动导叶开度为22°(小开度)的工况为典型叶道涡工况，因此选取该工况进行计算，该工况下的水轮机转速为136.4 r/min，出力为35MW。混流式水轮机将高精度的六面体结构化网格作为计算域，并在转轮叶片和导叶叶片表面采用O型网格来匹配和适应几何曲率，如图1所示。对网格进行无关性检测，综合考虑精度和计算能力，最后确定计算域的总网格数为7 712 258。

(a)水轮机全流道计算域 (b)叶轮网格

为便于分析转轮内部叶道涡时空演变的压力脉动，并验证数值结果的有效性，分别在水轮机蜗壳进口处、固定导叶前、转轮叶片通道内叶道涡常规出现的位置设置了7个压力脉动监测点，分别为SP1、SV1、RN1～RN5，如图1所示。

1.2 数值方法

三维非定常湍流流场的数值模拟考虑了湍流剪切应力的传递，因此采用剪切应力输运k-ω湍流模型对雷诺时均N-S方程进行封闭，该湍流模型在对流体机械流场内的流动分离计算时能得到较为准确的结果。进口、出口的边界分别设为质量流量进口边界和压力出口边界，所有固体壁面均为无滑移的绝热壁面。采用二阶精度格式离散有限体积方程中的对流项，采用二阶向后欧拉格式离散瞬态项，压力-速度耦合方法采用SIMPLE算法。静止部件与旋转部件之间的界面为滑移界面的交界面，时间步长为6°。将定常结果作为初始条件，进行偏负荷工况下的非定常模拟。模拟的非定常流动总时间为10.997 07 s(约为25个旋转周期)，水轮机水力效率在第17个旋转周期收敛，本文使用最后2个旋转周期的流场。

1.3 试验验证

将基于上述数值方法的仿真结果与试验结果进行比较，以验证和分析计算流体动力学CFD求解的可靠性。

模型试验中，采用公称精度为0.15%的电磁流量计直接测量模型水轮机的流量，通过测量蜗壳进口和尾水管出口的压差，得到水轮机的静水头，进而计算出水轮机的工作水头：

(1)

通过测量施加在转轮上的力矩计算得到水轮机机械功率：

Pm= 2πnTm

(2)

式中，n为模型水轮机转速，通过脉冲法测量；Tm为力矩，通过力传感器或直接测量扭矩得到。

最终得到模型水轮机的效率

(3)

式中，Ph为水轮机水力功率。

采用公称精度为0.2级的压力脉动传感器测量固定导叶前和蜗壳进口处的压力脉动，压力脉动幅值取97%置信度的混频峰峰值，主频为快速傅里叶变换(FFT)后的第一优势频率。所有的测量都按照IEC60193标准[17]进行。

数值预测结果与实验数据一致，额定出力、水头和水力效率的误差小于2.0%。对固定导叶前(监测点SP1)、蜗壳进口(监测点SV1)的压力信号进行快速傅里叶变换，得到的压力信号主频均为0.215fn=0.49 Hz，其中，水轮机转频fn=2.2733 Hz，而数值预测的主频均为0.233fn=0.53 Hz。对比结果可知，压力脉动主频的计算值与实验值误差可接受，验证了本文应用数值方法的准确性和可靠性。

2 DMD方法

本文将DMD方法作为分析混流式水轮机叶道涡瞬态特性的重要工具，其理论基础和主要步骤如下。

首先，将数值模拟得到的流场原始数据转换为一系列连续快照矢量。相邻快照的时间步长记为Δt。这些向量可以排列成一个按时间排序的矩阵：

V1,N=[v1v2…vN]

(4)

式中，vi为流场的第i个快照，i=1，2，…，N；N为快照总数。

假设存在一个线性映射算子A使得vi+1=Avi，则根据式(2)中的线性映射关系，用快照矩阵V1,N-1和线性映射算子A表示矩阵V2,N：

V2,N=AV1,N-1=[Av1Av2…AvN-1]

(5)

对于线性系统，算子A一定存在且为常数；对于非线性系统，算子A可以看作是线性正切近似[18]。因此，明确了初始状态V1,1和线性映射算子A，即可掌握流场任意时刻的状态。算子A是反映流场内在动力学机制且具体形式未知的高维矩阵，所以用算子A的伴随矩阵S作为A的低维近似矩阵[19]，得到

AV1,N-1=V1,N-1S

(6)

对矩阵V1,N-1进行奇异值分解：

V1,N-1=UZWT

(7)

式中，U、W为正交矩阵；Z为对角矩阵。

按照奇异值从大到小的顺序从上到下排列对角元素，可得低维近似矩阵的表达式：

S=UHAU=UHV2,NWZ-1

(8)

低维近似矩阵S的特征值及其特征向量可以通过特征值分解得到：

Smi=λimi

(9)

式中，mi为低维近似矩阵特征向量；λi为特征值。

动态模态的基模态定义为

qi=Umi

(10)

动态模态的指数增长率Reωi和频率fi可通过对特征值取对数来确定：

ωi=lnλi/Δt=Reωi+Imωi

(11)

fi=Imωi/(2π)

(12)

将qi在初始时间的流场数据v1投射到原始流场即可得到模态能量：

(13)

定义模态振幅为a，其中，ai为第i个模态的振幅，表示该模态对初始流场v1的影响小大，a可表示为

a=Y-1UHv1

(14)

式中，Y为低维近似矩阵S的特征向量矩阵。

所分解的动态模态不仅可以揭示系统的相干结构，还可用来重构降价的复杂动态系统：

(15)

得到的N-1个动态模态qj按照振幅aj的顺序进行排序，而重构足够精度的原始流场所需的模态个数n将远小于N-1，因此采用上述方法将极大地减小流场重构的数据量。

3 结果与分析

3.1 叶道涡结果分析

为分析转轮内部叶道涡时空演变的压力脉动，对监测点RN1～RN5的压力脉动数据进行了快速傅里叶变换，计算出5个监测点的压力脉动主导频率和幅值，如图2所示。5个监测点处压力脉动的主频主要分布在0.5fn～1.0fn内，在靠近叶道涡起始的转轮上冠处，监测点RN1、RN2、RN3处的主频均为0.5fn，这与叶道涡呈低频特性即诱发的压力脉动频率为转频的0.2～3.0倍[9]是相符的。

(a)RN1

采用Q准则旋涡识别方法分析转轮内部叶道涡的瞬态时空特性，利用Q的大小表示旋转强度。图3所示为0.5fn对应的周期T下，Q=40 s-1时叶道涡周期性地萌生、发展、局部涡的崩塌和消失的完整过程，以及旋涡处湍动能的变化规律。t1时刻，叶片与叶轮上冠交接处出现点状涡旋，叶片背面附着片状涡(片状涡在靠近叶片出水边处的湍动能较大)。t2时刻，叶道涡在叶轮上冠处初生，叶道涡中心线垂直于叶轮上冠。t3时刻，叶道涡结构空腔呈涡管状，由叶轮上冠延伸至叶片出口附近，叶道涡的湍动能从顶端至尾端不断增大，同时，叶片与叶轮上冠交接处的点状涡旋体积增大，呈微小涡管状。t4时刻，叶道涡充分发展，在叶轮流道中间位置贯穿整个流道，整体呈弧状涡旋结构，体积达到最大值，并与叶片背面的片状涡连为一体；湍动能随叶道涡的发展不断增大，这是由于水流在离心力作用下迫使叶道涡沿出水边方向移动，导致湍流强度升高[20]。在t5至t6时间段，叶道涡尾部与片状涡逐渐脱落，叶道涡的体积逐渐减小并消失。

(a)t1=0 (b)t2=0.2T

3.2 速度场DMD分析

为进一步解耦叶道涡的相干结构，分析其模态规律，选取垂直于转轮轴向0.5跨度处的速度流场数据，使用DMD方法对该截面处的120个快照(约2个旋转周期)进行分析处理。图4a所示为DMD特征值在复平面上的分布。由于特征值是共轭复数，因此关于横轴对称。所有特征值均位于单位圆上，表明相应的模态是周期稳定的。图4b所示为DMD模态能量频谱，各阶模态对流场的重要程度可用模态能量评估。能量最高的模态出现在零频率处。零频率代表平均流场，将能量最大的动态模态记为DMD Mode 0(图4b中未显示)。根据局部峰值选择出前五阶的动力模态，一阶模态(DMD Mode 1)的频率(水轮机叶频)为29.55 Hz(13fn)，二阶模态(DMD Mode 2)～五阶模态(DMD Mode 5)的频率分别为59.11 Hz(二倍叶频)、15.92 Hz、47.73 Hz和18.18 Hz。DMD Mode 1的能量最大，反映水轮机转轮与导叶之间的动静干涉作用是影响其内部流场涡旋结构动态时空演变的主导因素。

(a)DMD特征值分布

如图5a所示，零频率的速度模态即速度场流动结构的时间平均结果是原始流场的主导流动结构。DMD Mode 0可以看作转轮轴向0.5跨度截面速度场的基本结构，原始流场可看作是在DMD Mode 0上叠加不同频率的振荡模态而形成的。转轮与导叶的动静干涉界面存在高速流体区域，这影响了转轮进口速度的周向分布。导叶区域内速度场的基本流动结构相对于旋转中心对称。图5b为DMD Mode 1的速度云图，依据图4b的能量能量谱可知，该模态为引起转轮内非定常叶道涡的主要振荡模态。在DMD Mode 1中，转轮流道主要存在两种结构即结构A和结构B，两种结构的速度为正负值且交替出现，从叶片前缘开始贯穿整个流道。结构A起源于靠近前缘的叶片压力侧；结构B由动静干涉处的脱流引发，起源于叶片吸力侧。两种结构的速度沿叶片径向向内逐渐降低，并在靠近旋转轴处重新发展，相互作用而形成叶道涡。这是由于主流方向的结构A在黏性切向力的作用下形成旋转速度分量，同时脱流形成的真空导致结构B补入，进一步增大流动负冲角，导致水流在流道内旋转、形成叶道涡。对于图5c、图5d中的DMD Mode 2和DMD Mode3，动态模态的结构A和结构B依然成对出现，且结构演变与DMD Mode 1相似，只是规模较小、数量较多，且随着频率的升高，动态模态的分布趋于碎片化。

(a)DMD Mode 0 (b)DMD Mode 1

3.3 速度场POD分析

将POD方法用于分解具有相同快照的流场，实现两种方法的对比。POD各阶模态按本征模态对应的特征值能量排序，图6所示为POD模态的特征值能量占比P。POD的平均模态是能量最大的模态，记作POD Mode 0(未在图6中显示)。POD Mode1、POD Mode 2共占50%的总能量，POD Mode3、POD Mode 4共占总能量的20%，其他模态的能量都在5%以下，因此，将POD Mode1～POD Mode4作为主导模态。图7所示为POD Mode1～POD Mode4时间序列的快速傅里叶变换分析结果。POD Mode 1的频率为29.55 Hz，与DMD的一阶模态频率相同，均为叶频。POD Mode 2以29.55 Hz为主频率，POD Mode 3、POD Mode 4以59.11 Hz为主频率。使用POD方法进行分解后，不同频率下的速度模态并未完全解耦，各阶模态有重合。

图6 POD特征值分布

(a)POD Mode 1

图8所示为POD的四阶模态结构，POD Mode 0与DMD Mode 0相似，代表了该截面速度场的平均流场特征。POD Mode 1、POD Mode 2与DMD Mode 1高度相似，也存在正负速度交替出现的结构A和结构B；POD Mode 1和POD Mode 2的相干结构相似，但存在相移，这与DMD结果不同，表明POD分解后的模态成对出现。POD Mode 3与DMD Mode 2相似，为主导模态的碎片化结构。

(a)POD Mode 0 (b)POD Mode 1

DMD和POD的相似性表明这两种方法都能有效提取水轮机流场复杂的流动结构。DMD方法与POD方法的不同之处在于，POD方法通过流动能量对空间结构进行分解，分解后的各阶模态所携带的流场信息具有空间正交性，各阶POD模态包含多个流动频率，不适合物理现象的解释；DMD方法在特定频率下对流动结构进行分解，分解后的各阶DMD模态所携带的流场信息具有时间正交性，并且单倍频率特征可使研究者更清晰地进行流动机理分析。

3.4 速度场重构

为进一步观察DMD和POD两种方法对流场特征的提取效果，利用两种分解方法得到的Mode 1～Mode 3和平均流场Mode 0对第10.630 501 s的流场截面进行重构，如图9所示。从原始流场可以看出，叶轮靠近中心轴处存在着明显的低速叶道涡旋，叶片压力侧脱流结构向转轮通道中心传播。与原始流场相比，DMD和POD都能高还原度地重构主体流场特征，但在叶片压力侧仍存在一些细微的差异。

(a)原始流场 (b)DMD重构流场

图10为通过模态叠加重构的流场与原始流场的误差云图，可以发现DMD和POD的降阶重构流场与原流场的误差均在叶道涡周围最大，但POD方法重构流场的误差在叶片压力侧、叶轮通道中心处更大。偏差分析表明，水轮机转轮内叶道涡的发生是由高阶模态分量叠加而成的。两种方法的降阶重构虽然可以捕获大多数流场细节，但对复杂涡旋的还原依然存在局限性。DMD方法重构水轮机流场的误差明显小于POD方法重构流场的误差，因此DMD方法能重构出具有更高还原度和辨识度的流场结构。

(a)DMD重构流场误差 (b)POD重构流场误差

4 结论

(1)水轮机偏负荷工况下，叶道涡的时空演化频率为0.5倍转频，叶道涡在叶轮上冠处初生，空腔结构垂直于叶轮上冠并呈涡管状，空腔结构延伸至叶片出口附近，形成片状涡并附着在叶片背面。

(2)应用DMD方法对水轮机速度流场进行分解，得到不同频率的主要模态，各阶模态中主要存在速度正负交替出现的两种结构(从叶片前缘开始，贯穿整个流道)。这两种结构分别起源于靠近前缘的叶片压力侧和叶片吸力侧，在转轮中心轴附近相互作用形成叶道涡。POD方法根据能量大小分解得到的主要模态成对出现，前两对模态即POD Mode 1～Mode 4的主频率与DMD分解的前两阶模态 DMD Mode 1、Mode 2的频率相同，因此两种方法对应的相干结构存在相似性。DMD方法在特定频率下对流动结构进行分解可得单倍频率模态，能更清晰地提供各频率下速度场的流动机理特征。

(3)DMD方法和POD方法都可对流动结构进行降阶重构，但对水轮机转轮内局部复杂叶道涡区域的重构存在误差。

(4)DMD方法可对复杂流场的内在线性动力学结构进行辨识，进而对水轮机等复杂机械结构流场进行特征分析和快速降阶重构，在不需要额外建立控制方程的前提下，大幅降低了计算成本。