抽水蓄能电站能够灵活地开停机，跟踪并适应剧烈的负荷变化，并且通过调制调频保持电压稳定性，具有很好的峰值负荷削减功能，已经成为提供大规模电网储能的重要商业手段[1]．由于抽水蓄能电站运行工况复杂多样，因此当运行偏离最优工况并发生大波动过渡过程，特别是甩负荷工况时，机组的运行参数变化剧烈，容易激发水泵水轮机内部流场的快速瞬变，破坏稳定运行，甚至造成机组与人员的伤害[2]．研究甩负荷过程中水泵水轮机内部流场的瞬变性能对诊断易损结构、优化运行过程具有重要价值．目前有多种方法分析甩负荷过渡过程的机组瞬变性能，一种是通过一维动力学数值模型反映流道中的压力、流量等关键参数[3-4]，该方法计算时间短，但可靠性较低，无法全面地反映整个设备的空间参数变化．还有一种研究设备过渡过程的方法是运用三维计算流场动力学(CFD)技术对水轮机流道的压力流速分布与流态变化进行研究，该方法能更好地反映水力机械过渡过程中流道空间的三维空间参数分布，其试验结果与实际效果更加符合[5-7]，计算得到的结果接近于真实的物理模型，但模型求解时间往往长达数日甚至数月．模型降阶可以有效地解决效率与精度间的矛盾，通过将大规模甚至超大规模的系统或模型投影至低阶空间，用较少的自由度拟合原系统或模型，从而降低运算难度，提高时间效率．对于简单的线性系统，降阶模型中最常用的方法为Krylov子空间法[8-11]、平衡截断法[12]和函数正交分解．函数正交分解主要根据系统运算获取的数据样本通过构造一组基函数对系统进行降阶，被称为本征正交分解(POD)，是主流降阶方式．文献[13]将POD与LTI系统降阶模型结合，分析了流体冷却热电池的动态单输入单输出温度瞬态变化．文献[14]提出了基于序贯POD代理模型自适应优化，对二维翼型气动模型进行设计优化．文献[15-16]基于POD实现了两相流场的模型降阶，采用POD代理实现了叶型对水力性能的优化．根据以上研究可知基于POD的降阶模型对流场模型的预测可靠．综合上述研究，本研究以中国江西省的抽水蓄能电站实测数据为依据，基于POD的瞬态降阶模型(reduced order model，ROM)方法，建立以高水头水泵水轮机计算流体动力学(CFD)流场的降阶模型，融合水泵水轮机部位瞬态ROM改进的CFD运算模型，比较传统MOC模型与MOC-POD模型的精度与计算效率．1 相关理论1.1　快照POD降阶理论快照POD降阶模型用于CFD流场仿真以解决因体积元过多而导致的空间矩阵求解过慢或无法求解的问题．由CFD结果提取网格节点参数在各个时刻的数据合成快照输出矩阵A，并转化为平均值偏差矩阵A˜，即A(xi,tj)=u(x1,t1)u(x1,t2)⋯u(x1,tn)u(x2,t1)u(x2,t2)u(x2,tn)⋮⋱⋮u(xm,t1)u(xm,t2)⋯u(xm,tn)(i=1,2,⋯,m;j=1,2,⋯,n); (1)a˜j=aj-A¯, (2)式中：u(xi,tj)为通过数值计算获得的参数值；xi (i=1,2,⋯,m)为流量、压力、转矩等m个流场参数；tj (j=1,2,⋯,n)为n个时刻，通常情况下n≫m；aj为A的第j个行向量，a˜j为A˜的第j个行向量，A¯中的各元素为所有时刻的平均值，A¯=1n∑j=1naj．快照矩阵是一个半正定矩阵，满足AAT∈Rn×n．为找到一组使1n∑j=1n|〈a˜j,ϕj〉L2|2达到最大值且∑j=1n|〈a˜j,ϕj〉L2|2=1的POD正交基函数ϕi(x)，其中1〈∙,∙〉L2为L2范数中的典型内积，对A˜进行奇异值分解(singular value decomposition，SVD)转换，即A˜=UΣVT，(3)式中：U∈Rn×n和V∈Rm×m分别为由A˜A˜T和A˜TA˜的正交向量组成的正交矩阵；Σ为大小为n×m的对角矩阵．由于n≫m，对角元素Σii由m个非负数σi组成并按降序排列，即σ1≥σ2≥…≥σm≥0，被称为A˜(或A˜T)的奇异值．对于式(3)，若令UΣ=Q，则Q∈Rn×m，且A˜=QVT．设Q={q1,q2,…,qm}，V={v1，v2,⋯,vm}，则A˜=QVT=∑i=1mqiviT．(4)输出变量函数z(x,t)由A˜表示，模态系数函数λi(t)由qi表示，正交基函数ϕi(x)由viT表示．根据文献[17]中的证明，POD正交基函数及模态系数为ϕi(x)=A˜viσi1/2；(5)λi(t)=ϕi(x)TA˜ϕi(x)Tϕi(x)．(6)对A˜进一步实行低秩近似，对于任意rm，在Σ中设σr+1=σr+2=⋯=σm=0并得到矩阵Σr，用于计算A˜的最优r秩的近似值A˜≈A˜r=UΣrVT．(7)通常将矩阵U和V的前r列替代原矩阵，并将Σr替换为r×r的主子式．矩阵V的前r列给出了逼近数据的POD最优正交基，满足最优性条件，同时误差最小，任意时刻的流场变量可以表示为z(x,t)=z¯(x)+∑k=1rλk(t)ϕk(x)，(8)式中：z¯(x)=A¯为流场空间变量的平均值；∑k=1rλk(t)ϕk(x)为由一组r个基函数线性组合而成的降维空间．1.2　抽水蓄能电站MOC-POD模型本研究提出了一种将特征线法(method of characteristics，MOC)有压管道和快照POD水泵水轮机模型耦连的全流道降阶模型构建方法．MOC模拟有压管道和水轮机控制系统的方法已比较成熟．通过等效电路模型提取(equivalent circuit extraction，ECE)的方法求解高精度的发电电动机响应模型，ECE通过预先计算电机内部的参数构建响应面，利用电机性能补全并拓展其状态空间[18]．对于全流道数值模型，受到水泵水轮机和发电电动机作用影响的动力学模型在第τ个时间步上的迭代过程如下．步骤1 确定初始化时刻的导叶开度yτ、转速aτ,α和流量qτ,β，其中α和β分别为关于转速和流量的两个迭代数．步骤2 由水轮机POD模型，通过进出口节点平均压力得到相对水头hτ,β，同时运用MOC求解有压过水管道各节点的压力Hsβ，其中s表示第s个节点．步骤3 根据全特性曲线和相似理论求解水泵水轮机有功功率Pτ．步骤4 相对水头hτ,β作为过水系统状态方程的输入，求解水泵水轮机流量qτ+1,β=f(hτ,β)．步骤5 判断|qτ+1,β-qτ,β|δ1是否满足，其中δ1为最大流量绝对误差．若满足，则进行步骤6；若不满足，则判断β是否达到最大迭代数，若是则进行步骤6，否则设置权重参数λ计算qτ+1,β+1=(1-λ)qτ+1,β+λqτ,β，此时迭代数α加1，并返回步骤2．步骤6 由水轮机POD降阶模型得到相对转矩μτ,α=fPOD(τ,  yτ,  aτ,α,  qτ,β)．步骤7 运用ECE降阶模型计算发电电动机的相对转速aτ+1,α=fECE(μτ,α)．步骤8 判断|aτ+1,α-aτ,α|δ2是否满足，其中δ2为最大转速绝对误差．若满足，则进入下一步计算；若不满足，则判断α是否达到最大迭代数，是则进入下一步计算，否则设置权重参数η计算aτ+1,α+1=(1-η)aτ+1,α+ηaτ,α，迭代数α加1，并返回步骤2．2 全流道数值建模与应用2.1　研究对象本研究的标的抽水蓄能电站位于江西省，配备有四台完全相同的机组，为一管双机型结构．每台机组的额定功率为306.1 MW，比转速为35.15，额定转速为500 r/min，额定水头为540.0 m，额定流量为62.09 m3/s，安装高度为93 m．由此可知该电站的水轮机为高水头、高转速的可逆式水轮机．基本设计参数有：蜗壳进口直径Din=2 100.0 mm；尾水管出口直径Dout=4 400.0 mm；转轮高压侧直径D1=3 850.1 mm；转轮低压侧直径D2=1 934.9 mm；水轮机工况下100%导叶开度ymax=22.4°；转轮叶片数为9；双列叶栅组数为20．通过布置在水泵水轮机上的传感器获得甩负荷工况下设备关键位置的指标量数据．主要传感器有水力测量盘和压力脉动变送器(瑞士KELLER公司的21Y系列)和水力状态控制器(水泵水轮机系统内置，可在线测量主轴转速与导叶开度)等．在实测参数的基础上，构建高水头水泵水轮机的CFD数值模型，湍流模型选用SST k-ω模型．建立水泵水轮机简化等比例三维流场域模型，采用混合网格的方式进行网格划分，对位于双列叶栅流域与转轮叶片流域的液固交界面处的网格进行细化．水泵水轮机流体计算域模型见图1．10.13245/j.hust.220812.F001图1水泵水轮机流体计算域模型研究运用商业软件Fluent 19.0对水泵水轮机进行流场仿真，采用流量进口和压力出口的边界条件，通过编写用户定义函数(UDF)文件嵌入CFD模型．运用动网格方法，实现了20个活动导叶的刚体旋转运动[19]．转轮网格设置为滑移网格，流体域为显式旋转，转速通过UDF脚本实现设置．通过网格无关性验证，最终选取网格数为1.130 540 8×107个，网格质量均在0.4以上，符合流场瞬态分析的运算稳定性．根据实测数据，验证水泵水轮机的瞬态模型精度，边界条件参照现场实测数据，对甩负荷工况下的瞬态流场进行数值计算，总迭代时间为42 s．考虑转轮处滑移网格的刷新及动静干涉的频率分析，选取时间步长为0.001 s．通过均方根误差(root mean square error，RMSE)和平均绝对误差(mean absolute error，MAE)衡量无量纲化数值结果的精确度，表2为CFD模拟与实测无量纲数据误差表，表中：Hin为蜗壳进口压力；M为转轮转矩；Cp为进口压力脉动强度；HDTC为尾水锥管侧壁压力．10.13245/j.hust.220812.T001表2CFD模拟与实测无量纲数据误差表误差种类Hin/mM/(N•m)CpHDTC/m均方根误差0.017 00.023 70.027 30.029 0平均绝对误差0.015 80.021 10.026 20.034 0同时，可以通过在CFD计算事先设置测点或事后提取节点信息的方式实现参数对象的拓展．如在水泵水轮机的蜗壳、双列叶栅、水轮机流道或尾水管内部设置数个测点，获得甩负荷过程中的压力曲线，以此构造快照矩阵，然后进行模型降阶，以起到拓展参数范围的效果．由表2可知：通过CFD数值模拟得到的参数值与实测值的变化趋势相同，极值点对应的大小与时间也基本符合，同时RMSE与MAE处于工程允许的范围内．2.2　POD降阶模型构建完整地运行了两组甩负荷工况的瞬态计算，第一组为一段式导叶关闭规律(wicket gate closing law，WGCL)，总关闭时间为42 s(即关闭方式为42 s-0 s-97%)，第二组为两段式WGCL，关闭方式为10 s-35 s-20%．每隔0.5 s提取一组完整的节点参数矩阵，共有174组快照矩阵．为简化ROM的构建，去除了数据中的噪声，通过r阶模态的近似拟合获得ROM响应矢量Xred．进一步计算出不同模态阶数下的均方根误差及对应参数的时域变化特征．具有r阶模态预测的RMS误差为εRMS=∑i=1m||Xi-Xr,i||2/∑i=1m||Xi||2，(9)式中：Xi为第i个输出参数；Xr,i为r阶降阶压缩下第i个输出参数．不同模态数下均方根误差见图2，图中εEMS为均方根误差．可以看到：随着选择模态数的增加，曲线的拟合程度变优，最大绝对误差值逐渐减小．当选择模态数r时，若模态数较少，则整体误差较大；当考虑地模态数过多时，尽管预测误差减小，但会导致ROM构建过程过慢甚至无法实现，因此须要合理选择适当模态数来平衡误差与构建效率，并对学习点进行激励误差评估．10.13245/j.hust.220812.F002图2不同模态数下均方根误差保证ROM拥有低于1%的RMS误差，模态数不应小于5．当模态数为5时，RMS误差为0.87%．本研究构建了模态数为5的动态ROM，构建过程耗时15h 34min，为一次性过程．2.3　MOC-POD模型应用为验证降阶模型对水泵水轮机压力外特性的预测精度，选取两组一段式WGCL (42.03 s-0 s-97%，34.89 s-0 s-97%)和一组两段式WGCL (4.23 s-41.88 s-17.6%)作为结果对照，均为甩100%负荷工况，时间步长为0.001 s．从电站现场获取了实测数据(采样频率为120 Hz)，并将MOC-POD与传统一维MOC模型进行对比，实测数据与仿真数据对比见图3，图中：t为时间；Hout为蜗壳出口压力．10.13245/j.hust.220812.F003图3实测数据与仿真数据对比从运算时间角度对比，在56核心112线程的服务器上，对于三组瞬态甩负荷数值实验，运用传统MOC模型的计算耗时44min 45s、41min 35s和51min 7s，而运用MOC-POD模型计算仅耗时56 s，52 s和63 s，降低运算时间约98%，运算效率得到极大地提升．三组甩负荷工况实验各指标参数对比如表3所示．由图3和表3可知：采用三种不同WGCL的水轮机甩负荷工况，一维MOC模型的均方根误差与最大相对误差均较大，无法满足性能曲线的高精度预测．而通过MOC-POD预测关键性能参数，RMS误差与最大相对误差有了很大程度上的降低，同时曲线波形与实测参数曲线拟合性更强，证明MOC-POD模型具有较高的曲线预测可靠性．10.13245/j.hust.220812.T002表3三组甩负荷工况实验各指标参数对比WGCL参数MOC仿真无量纲压力MOC-POD仿真无量纲压力均方根误差最大相对误差均方根误差最大相对误差一段式I(42.03s-0s-97%)Hin7.2315.880.835.53Hout5.109.431.963.56一段式II(34.89s-0s-97%)Hin4.6715.471.414.48Hout8.9932.385.099.16二段式(4.23s-41.88s-17.6%)Hin3.079.110.982.17Hout8.6128.392.467.29%3 结论a. 在对高水头水泵水轮机CFD数值模型的基础上，运用快照POD理论对水泵水轮机在甩负荷工况下的瞬态模型进行了降阶压缩，实现了空间监测点性能参数的拓展与提取．b. 将MOC有压管道与水轮机控制系统建立了MOC-POD耦连模型，实现了抽水蓄能电站的甩负荷过渡过程运行的精确模拟，模态预测的均方根误差不高于1%．c. 通过将实测数据与MOC模型和MOC-ROM模型所求的外特性曲线与流体激励力进行对比，发现新方法较传统一维全流道模型模拟甩负荷工况外特性预测，均方根误差与最大相对误差普遍降低70%以上，运算时间减少98%．