總體矩陣法在水電站頻率調(diào)節(jié)穩(wěn)定性分析中的應(yīng)用

馬安婷，楊建東，楊威嘉，唐韌博，馮文濤，侯亮宇

（1. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072；2. 廣東省水利電力勘測設(shè)計研究院，廣東 廣州 510635）

1 研究背景

水電站水輪機調(diào)節(jié)系統(tǒng)通常由調(diào)速器、水輪機、有壓輸水系統(tǒng)、發(fā)電機和電網(wǎng)構(gòu)成，其運行穩(wěn)定性是由管道布置方式、機組運行工況點、調(diào)速器參數(shù)等多種因素決定的，是一個有條件穩(wěn)定系統(tǒng)［1］。傳遞函數(shù)法是進(jìn)行輸水發(fā)電系統(tǒng)穩(wěn)定性分析的重要方法，在以往的研究中，常采用簡化的單管單機布置方式及簡化的一階剛性水擊模型求取系統(tǒng)閉環(huán)傳遞函數(shù)得到運行穩(wěn)定性，但該方法不適用于一洞多機等復(fù)雜管道系統(tǒng)建模及分析。

水輪機模型可分為線性水輪機模型、非線性水輪機模型。在運行穩(wěn)定性的研究中，往往采用線性水輪機模型。陳嘉謀等［2］、戴敏等［3］、馬薇［4］采用Simulink 建立了復(fù)雜管道系統(tǒng)的數(shù)學(xué)模型，分析了線性水輪機模型下調(diào)節(jié)系統(tǒng)的穩(wěn)定性和調(diào)節(jié)品質(zhì)。郭文成等［5］推導(dǎo)了3 種不同調(diào)節(jié)模式下帶調(diào)壓室水電站的傳遞函數(shù)，分析了不同因素作用下穩(wěn)定域的變化情況。Yang 等［6］研究了由調(diào)壓室涌浪引起的低頻振蕩，并建立功率調(diào)節(jié)與頻率調(diào)節(jié)相互切換的調(diào)速器模型。劉昌玉等［7］考慮到調(diào)壓室振蕩效應(yīng)和分叉管水力耦合影響，建立了水輪機及調(diào)節(jié)系統(tǒng)精細(xì)化模型。俞曉東等［8］對水電站聯(lián)合運行進(jìn)行了小波動穩(wěn)定性的研究，研究表明電站聯(lián)合運行有利于增加系統(tǒng)穩(wěn)定性并改善調(diào)節(jié)品質(zhì)。索麗生［9］、鄭向陽［10］、俞曉東等［11］采用狀態(tài)空間法分析了水力機械系統(tǒng)小波動過程，狀態(tài)空間法是通過推導(dǎo)小波動線性微分方程組，求解系數(shù)矩陣的特征值和特征向量判斷系統(tǒng)的穩(wěn)定性，但推導(dǎo)過程較為繁瑣。已有的研究中側(cè)重于對復(fù)雜管道系統(tǒng)進(jìn)行Simulink 建?；蛘哌\用狀態(tài)空間法進(jìn)行系統(tǒng)穩(wěn)定性分析，存在不能準(zhǔn)確求解系統(tǒng)穩(wěn)定域和建模推導(dǎo)過程繁瑣等問題。

為了方便管道系統(tǒng)水力振動頻率分析，Wylie 等［12］提出水力阻抗法，但管道系統(tǒng)較為復(fù)雜時，節(jié)點阻抗的推導(dǎo)也較為困難。喬杜里［13］在水力阻抗法的基礎(chǔ)上提出了傳遞矩陣法，馮文濤［14］、段煉等［15］提出了便于編程的總體矩陣法，推導(dǎo)了特定條件下水輪機的水力阻抗和調(diào)壓室阻抗，并用于復(fù)雜管道系統(tǒng)的水力振動研究。葉復(fù)萌等［16］通過阻抗法模擬閥門的自激振動，并與特征線法進(jìn)行對比。Suo 等［17］提出了將頻率響應(yīng)轉(zhuǎn)化為時域響應(yīng)的脈沖響應(yīng)方法。Kim［18-19］提出了更加通用的阻抗矩陣法（IMM），該方法可用于復(fù)雜非均勻管網(wǎng)。通常用特征線法對管網(wǎng)系統(tǒng)暫態(tài)過程進(jìn)行分析，Ranginkaman［20］使用了傳遞矩陣法對一個真實的管網(wǎng)系統(tǒng)分析其頻率響應(yīng)，并論證了傳遞矩陣法與特征線法具有同樣高的精度。Vítkovsky［21］提出了基于拓?fù)渚仃嚨姆椒?，該方法用于任意拓?fù)浣Y(jié)構(gòu)的管網(wǎng)頻率分析，并考慮了黏彈性管道材料和非彈性摩阻的影響。周建旭［22］引入機組和調(diào)速器的動態(tài)特性，得到水力阻抗表達(dá)式，但沒有將其運用于復(fù)雜管道系統(tǒng)穩(wěn)定性分析中。

鑒于以上研究現(xiàn)狀以及分析復(fù)雜管道系統(tǒng)水電站運行穩(wěn)定性的需求，本文首先推導(dǎo)了完整的機組阻抗表達(dá)式，然后借用水力系統(tǒng)振動特性分析方法，建立系統(tǒng)總體矩陣，并以最大衰減因子σmax是否小于0 作為穩(wěn)定性的判別依據(jù)。本文的方法能夠準(zhǔn)確求解復(fù)雜輸水發(fā)電系統(tǒng)的穩(wěn)定域，而且通過修改系統(tǒng)各個模塊相應(yīng)的傳遞矩陣排列順序，本方法可適應(yīng)不同的電站布置形式，極大地方便了水電站運行穩(wěn)定性的計算與分析。

2 數(shù)學(xué)模型及判別條件

2.1 機組阻抗表達(dá)式的推導(dǎo)調(diào)速器采用頻率調(diào)節(jié)［6］，PI 調(diào)節(jié)含Ty環(huán)節(jié)：

式中：Ty為接力器時間常數(shù)，s；bp為永態(tài)轉(zhuǎn)差系數(shù)；Kp為比例增益；Ki為積分增益，s-1；y 為導(dǎo)葉開度偏差相對值；x 為轉(zhuǎn)速偏差相對值；t 為時間，s。

水輪機力矩方程、流量方程［23］：

式中：mt、qt、h 分別為力矩、流量和水頭偏差相對值；eh、ex、ey分別為水輪機力矩對水頭、轉(zhuǎn)速、導(dǎo)葉開度傳遞系數(shù)；eqh、eqx、eqy分別為水輪機流量對水頭、轉(zhuǎn)速、導(dǎo)葉開度傳遞系數(shù)。

發(fā)電機加速方程［23］：

式中：Ta為機組慣性時間常數(shù)，s；mg為發(fā)電機負(fù)載力矩偏差相對值；eg為發(fā)電機負(fù)載力矩隨轉(zhuǎn)速的變化率。

引水管道二階模型［23］：

式中：G1為引水管道水擊傳遞函數(shù)；H1為引水管道末端水頭偏差相對值；H0為水輪機初始水頭，m；Q1為引水管道末端流量偏差相對值；Q0為水輪機初始流量，m3/s；hw1為引水管道特征系數(shù)；Tr1為引水管道水錘壓力波反射時間，s；s 為拉普拉斯算子， s=σ+iω；l1為引水管道長度，m；g 為重力加速度，m/s2；A1為引水管道面積，m2；a1為引水管道波速，m/s。

尾水管道二階模型［23］：

式中：G2為尾水管道水擊傳遞函數(shù)；H2為尾水管道首端水頭偏差相對值；Q2為尾水管道首端流量偏差相對值；hw2為尾水管道特征系數(shù)；Tr2為尾水管道水錘壓力波反射時間，s；l2為尾水管道長度，m；A2為尾水管道面積，m2；a2為尾水管道波速，m/s。

對基本方程（1）—（4）進(jìn)行拉普拉斯變換，得：

聯(lián)立方程（5）—（10），可得出系統(tǒng)在負(fù)荷擾動mg作用下，單管單機水輪機調(diào)節(jié)系統(tǒng)的傳遞函數(shù)為：

式中： a0=λ5λ18；a1=λ5λ19+λ6λ18；a2=λ5λ20+λ6λ19+λ7λ18+λ13λ21；a3=λ6λ20+λ7λ19+λ8λ18+λ13λ22+λ14λ21；a4=λ7λ20+λ8λ19+λ9λ18+λ14λ22+λ15λ21；a5=λ8λ20+λ9λ19+λ10λ18+λ15λ22+λ16λ21；a6=λ9λ20+λ10λ19+λ16λ22+λ17λ21；a7=λ10λ20+λ17λ22；b0=-λ1λ18；b1=-λ1λ19-λ11λ18；b2=-λ1λ20-λ2λ18-λ11λ19；b3=-λ2λ19-λ11λ20-λ12λ18；b4=-λ2λ20-λ12λ19-λ18；b5=-λ12λ20-λ19；b6=-λλ20；λ4=hw1Tr1+hw2Tr2；λ5=Taλ1；λ6=enλ1+Taeqhλ3；λ7=Taλ2+eneqhλ3+eheqxλ3；λ8=Taeqhλ4+enλ2；λ9=Ta+eneqhλ4+eheqxλ4；λ10=en；λ11=eqhλ3；λ12=eqhλ4；λ13=eyλ1；λ14=eyeqhλ3-eheqyλ3；λ15=eyλ2；λ16=eyeqhλ4-eheqyλ4；λ17=ey；λ18=Ty+bpKpTy；λ19=1+bpKp+bpKiTy；λ20=bpKi；λ21=Kp；λ22=Ki。其 中，en=eg-ex為水輪發(fā)電機組綜合自調(diào)節(jié)系數(shù)。

將式（10）代入式（8），得：

式（12）×eqy-式（9）×ey，整理得：

改寫式（7），得：

將式（13）和式（14）代入式（15），整理得：

式中：a0=λ6λ1；a1=λ6λ2+λ5λ1；a2=λ6λ3+λ5λ2-Kpλ4；a3=λ5λ3-Kiλ4；b0=Taλ1；b1=Taλ2+enλ1；b2=Taλ3+enλ2+Kpey；b3=enλ3+Kiey；c0=-λ1eqx；c1=-(λ2eqx-Kpeqy)；c2=Kieqy-eqxλ3；λ1=λ2eqx-Kpeqy；c2=Kieqy-eqxλ3；λ1=Ty+ bpKpTy；λ2=1+bpKp+bpKiTy；λ3=bpKi；λ4=eqyeh-eqhey；λ5=eqhen+eqxeh；λ6=eqhTa。

式（16）可改寫為：

式中： B=b0s3+b1s2+b2s+b3；A=a0s3+a1s2+a2s+a3；C=c0s2+c1s+c2

式中：QT為水輪機處復(fù)流量，m3/s；HT為水輪機處復(fù)水頭，m；ΔMg為發(fā)電機負(fù)載力矩變化量，N·m；Mg0為初始工況發(fā)電機負(fù)載力矩，N·m。

水輪機前后管道的復(fù)水頭HT和復(fù)流量QT之間滿足如下關(guān)系［14］：

若考慮水輪機前后管道的水體動能，式（19）可改寫為：

式中：HiD為第i段管道末端復(fù)水頭，m；HjU為第j段管道首端復(fù)水頭，m；QiD為第i段管道末端復(fù)流量，m3/s；QiU為第i 段管道首端復(fù)流量，m3/s；為第i 段末端管道平均流量，m3/s；為第j 段首端管道平均流量，m3/s；AiD為第i 段末端管道面積，m2；AjU為第j 段首端管道面積，m2；

引入恒定流條件并對上式做線性化處理，略去高階非線性項，上式化簡為：

可寫成傳遞矩陣的形式：

式（22）為完整的機組阻抗通用表達(dá)式，針對簡單或復(fù)雜的有壓輸水發(fā)電系統(tǒng)皆能適用。

2.2 管道及有關(guān)邊界條件的傳遞矩陣管道及有關(guān)邊界的傳遞矩陣詳見參考文獻(xiàn)［14］，式（23）—式（27）為本文用到的邊界條件。根據(jù)管道首末兩端復(fù)水頭和復(fù)流量的關(guān)系，可寫成傳遞矩陣，如式（23）所示：

式中：li為第i 段管道長度，m；gi為第i 段管道傳播常數(shù)；ZCi為第i 段管道特征阻抗。

對于單管單機系統(tǒng)，假定過渡過程中上、下游水庫水位不變，即H1U=0 和H2D=0，則式（24）可改寫成式（5）和式（6）的形式。

根據(jù)岔管前后復(fù)水頭和復(fù)流量的關(guān)系，可寫成傳遞矩陣的形式，如式（25）所示：

式中：xij第i 段與第j 段管道間局部損失系數(shù)；xik第i 段與第k 段管道間局部損失系數(shù)。

上游水庫邊界條件，考慮管道水體動能和局部水頭損失的管道首端面復(fù)水頭和復(fù)流量的傳遞矩陣為：

式中x 局部水頭損失系數(shù)。

下游水庫邊界條件，考慮管道水體動能和局部水頭損失的管道末端面復(fù)水頭和復(fù)流量的傳遞矩陣為：

2.3 總體矩陣的構(gòu)成及穩(wěn)定性的判據(jù)在本文研究中，忽略了系統(tǒng)水頭沿程損失和局部損失，調(diào)節(jié)系統(tǒng)為線性系統(tǒng)，水擊模型為二階水擊模型，重點分析調(diào)速器參與頻率調(diào)節(jié)時比例增益Kp與積分增益Ki對穩(wěn)定性的影響。

由管道、岔管和水輪機等傳遞矩陣可建立輸水發(fā)電系統(tǒng)總體矩陣U，參數(shù)之間關(guān)系的方程組為：矩陣B 為邊界條件列向量。由于每根管道有4 個未知數(shù)，對于n 根管道構(gòu)成的水力系統(tǒng)總體矩陣U 的大小為4n×4n，前2n 行表示管道傳遞矩陣，后2n 行表示上庫、岔管和水輪機等傳遞矩陣，矩陣B 大小為4n×1，為首末斷面復(fù)水頭和復(fù)流量構(gòu)成的未知量向量，表達(dá)式如下［14］：

假設(shè)系統(tǒng)無外部擾動，即B=0，給定一組調(diào)速器參數(shù)Kp與Ki，采用牛頓迭代法求出U=0 時對應(yīng)的頻率ω和衰減因子σ。對于總體矩陣，可求得其中最大的衰減因子σmax，當(dāng)衰減因子σmax＞0 時，系統(tǒng)處于不穩(wěn)定狀態(tài)點；當(dāng)σmax=0 時，即系統(tǒng)處于臨界穩(wěn)定狀態(tài)點；當(dāng)σmax＜0 時，即系統(tǒng)處于穩(wěn)定狀態(tài)點。進(jìn)行多次計算后可根據(jù)σmax=0 時的調(diào)速器參數(shù)值求得系統(tǒng)的穩(wěn)定域。具體計算流程圖如圖1所示。

3 對比分析

3.1 單管單機系統(tǒng)以梨園水電站為例，用總體矩陣法分析輸水發(fā)電系統(tǒng)運行穩(wěn)定性，并與傳遞函數(shù)法和Simulink 仿真結(jié)果進(jìn)行對比。該水電站基本資料：單管單機輸水發(fā)電系統(tǒng)，無調(diào)壓室，布置示意圖如圖2，管道參數(shù)如表1。機組額定出力612 MW，額定水頭106 m，額定流量621.4 m3/s，額定轉(zhuǎn)速166.7 r/min，管線水流加速時間常數(shù)TW=4.02s，機組加速時間常數(shù)Ta=9.65s。水輪機傳遞系數(shù)取eh=1.49，ex=-0.99，ey=0.67，eqh=0.5，eqx=0.0，eqy=0.77。其他參數(shù)取Ty=0.02s，bp=0.04，eg=0。

圖1 計算流程圖

圖2 梨園水電站管線示意圖

表1 管道參數(shù)

圖3所示為用總體矩陣法求解不同調(diào)速器參數(shù)對應(yīng)的最大衰減因子σmax。通過變換求取σmax=0 時對應(yīng)的Kp和Ki值，如圖4綠色線所示；用傳遞函數(shù)法進(jìn)行數(shù)值模擬，結(jié)果如圖4藍(lán)色線所示。曲線下半部分為穩(wěn)定域，上半部分為不穩(wěn)定域。可以看出，兩條線重合。

用Simulink 建立調(diào)節(jié)系統(tǒng)整體模型，取通過總體矩陣法推導(dǎo)出的穩(wěn)定域邊界上的若干組點，進(jìn)行數(shù)值仿真，可得到穩(wěn)定點與不穩(wěn)定點，如圖4所示。以兩組參數(shù)為例：①Kp=1.0，Ki=0.37；②Kp=1.0，Ki=0.39。擾動為階躍擾動，擾動值為0.1，擾動從第10 s 開始，輸出的時域響應(yīng)如圖5和圖6所示。當(dāng)Kp=1.0，Ki=0.37 時，調(diào)速器參數(shù)處于穩(wěn)定域內(nèi)，系統(tǒng)是穩(wěn)定的，響應(yīng)曲線收斂。當(dāng)Kp=1.0，Ki=0.39 時，調(diào)速器參數(shù)處于不穩(wěn)定域，系統(tǒng)是不穩(wěn)定的，響應(yīng)曲線發(fā)散。

上述3 種方法得到的計算結(jié)果，高度吻合。說明了總體矩陣法可準(zhǔn)確求解單管單機輸水發(fā)電系統(tǒng)的運行穩(wěn)定域，為進(jìn)一步驗證總體矩陣法的適用性奠定了基礎(chǔ)。

圖3 不同調(diào)速器參數(shù)對應(yīng)的最大衰減因子σmax（總體矩陣法）

圖4 總體矩陣法與傳遞函數(shù)法和Simulink模擬對比

圖5 負(fù)荷給定擾動階躍響應(yīng)（KP=1.0，Ki=0.37）

圖6 負(fù)荷給定擾動階躍響應(yīng)（KP=1.0，Ki=0.39）

3.2 一洞四機系統(tǒng)以喀臘塑克水電站為例，分別用總體矩陣法和Simulink 仿真分析電站穩(wěn)定性。該水電站基本資料：一洞四機輸水系統(tǒng)，布置示意圖如圖7，管道參數(shù)如表2。機組額定出力36 MW，額定水頭79.5 m，額定流量50.23 m3/s，額定轉(zhuǎn)速300 r/min，1#管線水流加速時間常數(shù)TW=2.98s，機組加速 時 間 常 數(shù)Ta=10.11s。 水 輪 機 傳 遞 系 數(shù) 取eh=1.4568， ex=-0.9136， ey=0.6001， eqh=0.5320，eqx=-0.0640，eqy=0.8391。其他參數(shù)取Ty=0.02s，bp=0.02，eg=0。

圖7 喀臘塑克水電站管線示意圖

表2 管道參數(shù)

圖8所示為用總體矩陣法求解不同調(diào)速器參數(shù)對應(yīng)的最大衰減因子σmax。通過變換求取σmax=0 時對應(yīng)的比例增益Kp和積分增益Ki值，如圖9綠色線所示，曲線下半部分為穩(wěn)定域，上半部分為不穩(wěn)定域。

用Simulink 建立調(diào)節(jié)系統(tǒng)整體模型，取通過總體矩陣法推導(dǎo)出的穩(wěn)定域邊界上的若干組點，進(jìn)行數(shù)值仿真，可得到穩(wěn)定點與不穩(wěn)定點，如圖9所示。以兩組參數(shù)為例：①Kp=3.0，Ki=0.53；②Kp=3.0，Ki=0.55，用Simulink 建立調(diào)節(jié)系統(tǒng)整體模型，進(jìn)行仿真模擬。擾動為階躍擾動，擾動值為0.1，擾動從第10 s 開始，輸出的時域響應(yīng)如圖10和圖11所示。當(dāng)Kp=3.0，Ki=0.53時，調(diào)速器參數(shù)處于穩(wěn)定域內(nèi)，系統(tǒng)是穩(wěn)定的，響應(yīng)曲線收斂。當(dāng)Kp=3.0，Ki=0.55 時，調(diào)速器參數(shù)處于不穩(wěn)定域，系統(tǒng)是不穩(wěn)定的，響應(yīng)曲線發(fā)散。兩種方法得到的結(jié)果完全相同，再次驗證了總體矩陣法的準(zhǔn)確性，說明了總體矩陣法可以用于一洞四機輸水發(fā)電系統(tǒng)的運行穩(wěn)定性分析中。

圖8 不同調(diào)速器參數(shù)對應(yīng)的最大衰減因子σmax（總體矩陣法）

圖9 總體矩陣法與Simulink 模擬對比

圖10 負(fù)荷給定擾動階躍響應(yīng)（KP=3.0，Ki=0.53）

圖11 負(fù)荷給定擾動階躍響（KP=3.0，Ki=0.55）

4 結(jié)論

本文將總體矩陣法應(yīng)用于水電站輸水發(fā)電系統(tǒng)運行穩(wěn)定性分析中，推導(dǎo)了機組阻抗表達(dá)式，建立了布置形式為單管單機和一洞四機的水電站輸水發(fā)電系統(tǒng)總體矩陣，以最大衰減因子σmax是否小于0 作為穩(wěn)定性的判別依據(jù)，繪制出相應(yīng)的穩(wěn)定域。本文還通過上述兩個工程實例對本方法進(jìn)行了驗證：（1）用傳遞函數(shù)法和Simulink 仿真對單管單機系統(tǒng)進(jìn)行了對比驗證，結(jié)果表明總體矩陣法所得到的穩(wěn)定域與傳遞函數(shù)法和Simulink 仿真得到的穩(wěn)定域吻合度較高， 說明了采用總體矩陣法可以計算單管單機的穩(wěn)定域。（2）建立了管道系統(tǒng)Simulink 模型對一洞四機系統(tǒng)進(jìn)行驗證，結(jié)果是總體矩陣法所得到的穩(wěn)定域與Simulink 仿真得到的穩(wěn)定域重合，說明了采用總體矩陣法可以計算復(fù)雜管道系統(tǒng)的穩(wěn)定域。當(dāng)水電站布置形式發(fā)生變化時，只需修改系統(tǒng)各個模塊相應(yīng)的傳遞矩陣排列順序即可，極大方便了運行穩(wěn)定性分析，為深入分析各種不同復(fù)雜管道系統(tǒng)水電站運行穩(wěn)定域提供了新的途徑。