水电是重要的可再生能源之一[1]，在电网调峰调频中起着至关重要的作用[2]．国内对电力资源的需求不断增加，确保电网的安全稳定运行变得越来越重要[3]．随着国内水电工程的不断发展，出现了越来越多具有长引水隧洞的水电站[4]．若在长引水隧洞电站上游只设置一个调压室，由于上游引水管道水流惯性过大，因此调压室面积会很大[5]，这可能会导致地形地质及施工条件无法满足安全稳定的要求．为了解决上述问题，现在采用设置上游串联双调压室的布置方式[6]．李振轩等[6]分析了双调压室面积变化对系统稳定性的影响；滕毅等[7]分析了调压室面积和管道特性对系统稳定性的影响；基于托马假定，陈岚等[8]推导了调压室临界断面面积的稳定域；郭文成等[9]运用开环模型分析了上游串联双调压室系统甩负荷时调压室的涌浪特性．现有研究大多根据现代控制理论[10]，采用传递函数法[11]或者状态空间法[12]，根据特征根实部的正负来判断上游串联双调压室系统的稳定性，对于系统的振荡特性及失稳机理尚缺乏理论研究．本研究根据相关因子法得到系统特征值和状态变量的相关关系，结合实际算例分析了系统发生扰动时的振荡特性及管道参数和调压室面积对系统稳定性的影响，为带上游串联双调压室水电站的设计提供理论依据．1 数学模型带上游串联双调压室的水力发电系统布置图如图1所示，图中：Fu2和Fu1为主副调压室面积；Zu2和Zu1分别为主副调压室水位．10.13245/j.hust.238687.F001图1带上游串联双调压室的水力发电系统布置长引水隧洞动量方程Tw1dq1dt=-2α1Q02H0q1-Z1，(1)式中：Tw1为长引水隧洞水流惯性常数；α1为长引水隧洞水头损失系数；q1=(Q1-Q0)/Q0为长引水隧洞流量相对变化率，Q1为长引水隧洞流量，Q0为稳定工况时水轮机工作流量；t为时间；H0为稳定工况时水轮机工作水头；Z1=(Zu1-Z10)/H0为副调压室水位相对变化，Z10为副调压室初始水位．副调压室连续方程Fu1dZ1dt=Q0H0q1-Q0H0q2，(2)式中q2=(Q2-Q0)/Q0为短引水隧洞流量相对变化率，Q2为短引水隧洞流量．短引水隧洞动量方程Tw2dq2dt=Z1-2α2Q02H0q2-Z2，(3)式中：Tw2为短引水隧洞水流惯性常数；α2为短引水隧洞水头损失系数；Z2=(Zu2-Z20)/H0为主调压室水位相对变化，Z20为主调压室初始水位．主调压室连续方程Fu2dZ2dt=Q0H0q2-Q0H0q3，(4)式中q3=(Q3-Q0)/Q0为压力管道流量相对变化率，Q3为压力管道流量．压力管道动量方程Tw3dq3dt=-2α3Q02H0q3-h+Z2，(5)式中：Tw3为压力管道水流惯性常数；α3为压力管道系数；h=(H-H0)/H0为水轮机水头相对变化率，H为水轮机水头．水轮机力矩和流量方程mt=ehh+exφ+eyμ；(6)q3=eqhh+eqxφ+eqyμ，(7)式中：mt=(Mt-Mt0)/Mt0为水轮机动力矩相对变化率，Mt为水轮机动力矩，Mt0为水轮机动力矩初始值；φ=(n-n0)/n0为机组转速相对变化率，n为机组转速，n0为机组转速初始值；μ=(τ-τ0)/τ0为机组开度相对变化率，τ为机组开度，τ0为机组开度初始值；eh，ex，ey，eqh，eqx，eqy分别为水轮机特性系数．发电机方程Tadφdt=mt-mg-spφ，(8)式中：Ta为机组惯性时间常数；mg=(Mg-Mg0)/Mg0为水轮机阻力矩相对变化率，Mg为水轮机阻力矩，Mg0为水轮机阻力矩初始值；sp为电网自调节系数．PI调速器方程(bt+bp)Tddμdt+bpμ=-Tddφdt-φ，(9)式中：bt，bp和Td分别为调速器暂态调差系数、永态调差率和缓冲时间常数．将式(6)和(7)代入式(5)、(8)和(9)，消去mt和h，然后利用式(1)~(5)、式(8)和(9)建立上游双调压室水力发电系统小波动模型，根据状态空间法[13]可以得到dX/dt=AX+R，(10)式中：A为系统矩阵；X=[q1,Z1,q2,Z2,q3，φ,μ]为状态矩阵；R为系统扰动量．通过识别矩阵A的特征值实部的正负，可以判断系统的稳定性．当A的特征值实部全部为负时，系统处于稳定状态．此外，采用四阶龙格库塔法求解式(10)，可以得到系统受到扰动后的动态过程，如机组转速波动和调压室水位波动[14]．值得说明的是，进行小波动过渡过程计算时，一般假定波动是微小的，因而可以忽略基本方程的非线性项，将方程线性化，应用现代控制理论对系统的稳定性进行理论分析．2 相关因子法对于矩阵A，其特征值λi对应的右特征向量可表示为AUi=λiUi，左特征向量Vi可表示为ViTA=ViTλi，即ATVi=λiVi，其中Ui和Vi均为列向量．量度第k个状态变量和第i个特征值的相关性的相关因子表达式为[15]Pki=vkiuki/(ViTUi)，(11)式中：uki和vki分别为右特征向量Ui和左特征向量Vi的第k行元素．Pki反映了状态变量和特征值之间的相关关系．Pki=0，表示第k个状态变量与第i个特征值无关；Pki≠0，表示第k个状态变量与第i个特征值相关．Pki越小，相关性越小；Pki越大，相关性越大．3 算例分析3.1　工程概况以某上游串联双调压室电站为例，电站布置如图1所示，基本参数如下：L1=16 206.85 m，α1=4.41×10-5 s2/m5，Tw1=25.27 s；L2 = 621.58 m，α2=1.53×10-6 s2/m5，Tw2=0.97 s；L3=649.69 m，α3=6.23×10-5 s2/m5，Tw3=3.14 s；Q0=241.9 m3/s；H0=129 m；eh=1.653，ex=-1.15，ey=0.811；eqh=0.347，eqx=-0.21，eqy=0.85；Sp=0．根据上述参数，分别取4种方案对系统进行分析，即：方案1，Fu1=600 m2，Fu2=800 m2，bt=0.5，Td=8；方案2，Fu1=300 m2，Fu2=800 m2，bt=0.5，Td=8；方案3，Fu1=600 m2，Fu2=600 m2，bt=0.5，Td=8；方案4，Fu1=600 m2，Fu2=800 m2，bt=0.4，Td=5．将上述方案的参数代入矩阵A，采用遍历法求解全部调速器参数组合(bt和Td)下系统状态矩阵的特征值实部，通过判断特征值实部是否全部小于0可以得到系统的稳定域，如图2所示．图中每个方案均只有一条稳定域边界线，在边界线上方系统是稳定的．10.13245/j.hust.238687.F002图2不同方案下的稳定域方案1和4的调压室面积组合完全一样，两个方案的稳定域是相同的．和方案1相比，方案2的副调压室的面积较小，方案3的主调压室的面积较小，因此这两个方案的稳定域都比方案1小．此外，方案1的调速器参数组合在稳定域内，因此系统发生扰动时是稳定的；而方案2~4的调速器参数组合均处于失稳区域，因此系统发生扰动时会发散．用相关因子法对系统的特征值和状态变量的相关因子进行求解，分析系统的振荡特性．3.2　基于相关因子法的振荡特性分析根据方案1~4的参数及电站参数，利用式(11)计算4种方案下系统各个特征值与状态变量的相关因子，结果如表1所示．10.13245/j.hust.238687.T001表14种方案下系统各个特征值与状态变量相关因子方案特征值Pkiq1Z1q2Z2q3φμ1λ1，λ2-0.000 1±0.007 1i0.5040.2120.0070.2960.0020.0000.006λ3，λ4-0.000 2±0.075 4i0.0060.2930.4940.2050.0020.0000.004λ5，λ6-0.082 2±0.321 9i0.0000.0000.0000.0000.0050.8470.148λ7-0.112 00.0000.0000.0000.0060.8520.0960.7582λ1，λ20.000 1±0.008i0.5020.1330.0100.3750.0030.0000.007λ3，λ4-0.000 4±0.095i0.0100.3700.4900.1280.0010.0000.003λ5，λ6-0.082 2±0.321 9i0.0000.0000.0000.0000.0050.8470.148λ7-0.112 00.0000.0000.0000.0060.8520.0960.7583λ1，λ20.000 4± 0.00 76i0.5060.2490.0060.2600.0030.0000.006λ3，λ40.000 1±0.080 4i0.0050.2570.4950.2400.0030.0000.006λ5，λ6-0.082 6±0.322 4i0.0000.0000.0000.0000.0050.8470.148λ7-0.112 00.0000.0000.0000.0080.8510.0960.7554λ1，λ2-0.000 1±0.007 1i0.5040.2130.0070.3000.0020.0000.003λ3，λ4-0.000 2±0.075 5i0.0060.2920.4940.2080.0020.0000.002λ5，λ60.007 7±0.340 5i0.0000.0000.0000.0030.2080.8620.067λ7-0.211 30.0000.0000.0000.0060.8220.1780.770矩阵A的特征值存在三对共轭复根，即整个系统有3个振荡环节，分别为低频振荡环节S1、中频振荡环节S2和高频振荡环节S3，共轭复根的虚部为振荡环节的振荡频率，实部反映了振荡环节的稳定性．以方案1为例，低频振荡环节S1的振荡频率为0.007 1 Hz，对应的特征值为λ1和λ2，且q1，Z1和Z2对应λ1和λ2的相关因子分别为0.504，0.212和0.296，其他状态变量对应λ1和λ2的相关因子均约为0，这说明S1对应的状态变量为q1，Z1和Z2，即低频振荡环节S1主要受长引水隧洞、副调压室和主调压室的影响．中频振荡环节S2的振荡频率为0.075 4 Hz，对应的特征值为λ3和λ4，且Z1，q2和Z2对应λ3和λ4的相关因子分别为0.293，0.494和0.205，其他状态变量的相关因子均约为0，说明S2对应的状态变量为Z1，q2和Z2，即中频振荡环节S2主要受副调压室、短引水隧洞和主调压室的影响．高频振荡环节S3的振荡频率为0.321 9 Hz，对应的特征值为λ5和λ6，φ和μ对应λ5和λ6的相关因子分别为0.847和0.148，其他状态变量的相关因子均约为0，说明S3对应的状态变量为φ和μ，即高频振荡环节S3主要受机组和调速器的影响．根据其他方案也可以得出相同的结论．由于方案1的特征值实部均小于0，因此方案1的系统稳定．与方案1相比，方案2减小了Fu1，低频振荡环节S1的特征值实部大于0，说明了方案2中系统失稳主要是由于低频振荡环节S1失稳造成的，减小Fu1会降低低频振荡环节S1的稳定性．与方案1相比，方案3减小了Fu2，低频振荡环节S1和中频振荡环节S2的特征值实部均大于0，说明了方案3中系统失稳主要是由于低频振荡环节S1和中频振荡环节S2失稳造成的，减小Fu2会降低低频振荡环节S1和中频振荡环节S2的稳定性．与方案1相比，方案4减小了调速器参数，高频振荡环节S3的特征值实部大于0，说明了方案4中系统失稳主要是由于高频振荡环节S3失稳造成的，减小调速器参数会降低高频振荡环节S3的稳定性．3.3　数值验证为了验证上述结论，采用四阶龙格-库塔法分别求解方案1~4的机组转速的时域过程，如图3所示．10.13245/j.hust.238687.F003图3方案1~4的机组转速时域过程分析图3可得：在机组转速的时域中，存在3种不同频率的波，分别对应低频振荡环节S1，中频振荡环节S2和高频振荡环节S3．此外，机组转速时域过程还存在明显的主波与尾波之分，主波主要发生在机组转速振荡的刚开始，具有振荡频率高和收敛速度快的特点，因此主波对应的是高频波(高频振荡环节S3)．而尾波发生在主波振荡之后，表现为低频波(低频振荡环节S1)和中频波(中频振荡环节S2)的叠加．方案1中的机组转速的主波快速收敛，而后尾波也逐渐收敛，因而系统是稳定的．方案2中机组转速的主波快速收敛，当计算时间大于5 000 s后转速的时域过程只存在低频波，且低频波逐渐发散，因此方案2的系统失稳是由低频振荡环节S1引起的．方案3中机组转速的主波快速收敛，而后低频波和高频波都逐渐发散，因此方案3中的系统失稳是由低频振荡环节S1和中频振荡环节S2共同引起的．方案4中的高频波快速发散，此时系统转速表现为只存在高频波，因此方案4中的系统失稳是由高频振荡环节S3引起的．综上所述，以上运用状态空间法模拟机组转速的时域过程得到的结论与上述相关因子法分析得到的结论完全一致，验证了用相关因子法进行系统的振荡特性分析的正确性．3.4　敏感性分析由于尾波发生在主波之后，且调速器参数越大主波越稳定这一结论与单调压室系统相一致，因此主要研究系统尾波(即低频振荡环节和中频振荡环节)稳定性的影响因素．根据3.2节的分析可知：低频振荡环节S1主要受长引水隧洞、副调压室和主调压室的影响，中频振荡环节S2主要受副调压室、短引水隧洞和主调压室的影响，因而选取不同的L1，L2，Fu1和主调压室面积Fu2对低频振荡环节S1和中频振荡环节S2的稳定性进行分析．根据现代控制理论可知：特征值实部反映了系统的稳定性，当特征值实部大于0时系统失稳，且特征值实部越大系统稳定性越差．采用方案1的参数和电站参数，分别取不同的L1，L2，Fu1和Fu2，对矩阵A的特征值实部进行求解，从中选取系统低频振荡环节和中频振荡环节的特征值实部，得出系统参数对低频环节和中频环节特征值实部(E)的影响如图4所示．10.13245/j.hust.238687.F004图4系统参数对低频环节和中频环节特征值实部的影响随着L1增大，低频振荡环节的特征值实部先快速增大，而后增大速度明显变慢，随后当L1增大到大于1.86×104 m时，低频振荡环节的特征值实部大于0，此时低频振荡环节失稳．当L1增大时，中频振荡环节的特征值实部缓慢增大．随着L2的增大，低频振荡环节的特征值实部缓慢增大，当L2增大到大于2 744 m时，低频振荡环节的特征值实部大于0，此时低频振荡环节失稳．当L2增大时，中频振荡环节的特征值实部先迅速增大，而后增大速度明显变慢，当L2增大到大于897 m时，中频振荡环节的特征值实部大于0，此时中频振荡环节失稳．总体而言，随着L1和L2的增大，管道中的水流惯性增大，不利于系统的稳定．随着Fu1的增大，低频振荡环节的特征值实部迅速减小，说明增大副调压室的面积可以有效提升低频振荡环节的稳定性，当Fu1458 m2时，低频振荡环节特征值的实部大于0，低频振荡环节失稳．当Fu1增大时，中频振荡环节的特征值实部迅速增大，说明增大副调压室的面积会降低中频振荡环节的稳定性，当Fu11 150 m2时，中频振荡环节的特征值实部大于0，表明中频振荡环节失稳．随着Fu2的增大，低频振荡环节和中频振荡环节的特征值实部均快速减小，说明增大Fu2可以有效提升低频振荡环节和中频振荡环节的稳定性，当Fu2647 m2时，低频振荡环节和中频振荡环节的特征值实部均大于0，表明低频振荡环节和中频振荡环节均失稳．Fu2比Fu1对低频振荡环节和中频振荡环节的特征值实部影响更大，说明系统的稳定性主要受主调压室的影响．综上所述，减小L1或者增大Fu1和Fu2可以提升系统低频振荡环节的稳定性，减小L2和Fu1或者增大Fu2可以提升系统中频振荡环节的稳定性．4 结论a．通过相关因子分析法，可以对系统的特征值和状态变量的相关因子进行求解，得到系统低频、中频和高频振荡环节的主要影响因素．b．机组转速的时域过程存在多频振荡特性，且有主波与尾波之分，主波对应的是高频波，尾波表现为低频波和中频波的叠加．c．减小长引水隧洞的长度或者增大副调压室和主调压室的面积可以提升系统低频振荡环节的稳定性，减小短引水隧洞的长度和副调压室的面积或者增大主调压室的面积可以提升系统中频振荡环节的稳定性，增大调速器参数可以提升系统高频环节的稳定性．