混流式水輪機特性曲線在多重邊界條件下的分區方法

馬偉超，楊桀彬，趙志高，楊威嘉，楊建東

（武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072）

0 引 言

隨著農業現代化的發展，通過修建引調水工程解決農業和生態環境用水問題是一種必然的發展趨勢[1]。以重力流輸水的長距離引調水工程首末段水頭差往往較大，常建設消能水電站，以合理運用其剩余水頭[2]，在滿足引調水工程安全的前提下，充分發揮其灌溉、發電等水資源綜合利用效益。水電站及水電機組的安全、穩定、高效運行也因此越來越受到水利部門的重視，對水電機組建模仿真精度提出了更高的要求[3-5]。水輪機綜合特性曲線是進行水電站運行仿真與調節保證分析的基本數據來源，特性曲線的完整程度、疏密程度及準確性直接影響仿真結果的可靠性。然而，水輪機廠家所提供的混流式水輪機綜合特性曲線試驗結果范圍相對較小，較難滿足水輪機各種工況過渡過程仿真需求[6]。因此，在計算前需要對水輪機特性數據進行擴展和補充，外延內插水輪機特性曲線[7-8]，目前主要處理方法有外特性法和內特性法。

外特性法充分利用水輪機綜合特性曲線和飛逸特性曲線的模型試驗數據，通過純數學方法合理預估、推測各物理量的變化規律。阮文山[9]假定力矩特性曲線將交于一點，并采用二次曲線模型擴展，此方法計算簡單有效，但在低轉速區域時誤差較大；張培等[7,10]和王珊等[11-13]利用神經網絡搭建了水輪機數學模型，有效減小了數值仿真誤差，但曲線外延精度仍存在較大的提升空間，劉冬等[14]在水輪機神經網絡模型中引入評價和修正函數，有效提高特性曲線的預測精度；張蓉生等[15]通過Delaunay三角剖分及分片線性插值的方式擴展了水輪機效率曲線；鄭源等[16-17]將水輪機特性曲線進行簡單分區處理，其區域特征描述更精確，但該方法計算量大，各分區邊界點存在不連續現象。上述外特性方法由于缺乏水輪機水力特性理論的支撐，所得結果與機組運行的真實情況存在一定差異。

內特性法通過理論推導研究了水輪機結構參數之間的內在聯系，充分反映了水輪機內部特征與水流流動規律，可以有效減小外特性法對個人經驗的依賴。常近時[18-19]推導了水輪機廣義基本方程與全特性曲線的解析表達式，方程充分反映了水輪機內部特征，但未利用外特性試驗數據，計算精度較低；趙林明等[20]基于水輪機流量調節方程，提出了水輪機線性模型；楊建東等[21-22]推導了水輪機零轉速和交點邊界及可逆式水泵水輪機的特征交點公式；門闖社等[23-24]改進了水輪機內特性模型，并結合試驗數據辨識內特性參數，但方程求解易陷入病態，參數辨識結果可能與其物理意義不符，且轉速較高時，內特性方程可能無解。

針對上述局限性，本研究綜合考慮水輪機內外特性，結合多重邊界條件、外特性試驗數據與內特性模型，提出了一種混流式水輪機特性曲線分區處理方法；對某水輪機進行實例分析，辨識其內特性參數，外延內插其特性曲線，并與典型的外特性法與內特性法對比；最后將完整的特性曲線應用于過渡過程試驗反演及對比，以期獲得較高的仿真精度，更好地應用于工程實際。

1 混流式水輪機特性曲線的邊界條件

水輪機特性曲線通常采用單位流量、單位力矩與導葉開度、單位轉速間的關系式表達?；诹髁空{節方程和能量平衡方程，可建立水輪機內特性模型[24]，如式（1）所示。

式中M11為單位力矩，N·m；Q11為單位流量，m3/s；n11為單位轉速，r/min；a1~a7為結構參數，kg/m3，具體表達式見文獻[24]，其中，a1i、a3i、a4i與開度α有關，下標i表示第i條開度線對應的參數，其他參數與開度無關。

將邊界條件作為特性曲線分區的特征點，約束各分區內特性曲線的延拓范圍。邊界條件由特性曲線外延內插的范圍確定，主要包括以下5類：開度線與n11=0軸、Q11=0軸、M11=0軸的交點、零開度線和單位力矩的交點，分別簡稱為零轉速條件、零流量條件、飛逸條件、零開度線和交點條件。

1.1 邊界條件表達式

1.1.1 零轉速條件

包括零轉速時的單位流量和單位力矩。單位流量可根據單位出力、效率和單位轉速間的關系，通過洛必達法則計算[21]，如式（2）所示。由式（1）流量調節方程推導可得單位力矩，如式（3）所示。

式中f1(n11)與f2(n11)分別為某條開度線上單位轉速關于單位出力P11(W)和水輪機效率η的函數。

1.1.2 零流量條件

零流量條件為水輪機絕對流量=0時對應的單位轉速和單位力矩。由于零流量時速度三角形的特殊性，很難直接計算對應的單位轉速，因此先推導轉輪進口處相對流速W1在進口切線方向（圖1中虛線方向）的速度分量Wb1=0時的單位轉速，并對流量系數進行修正，近似得到絕對流量=0對應的單位轉速。分析轉輪進口處的速度，可列方程組如式（4）所示。

式中下標0表示導葉出口，下標1表示轉輪進口；U為圓周速度，m/s；D為斷面直徑，m；n為水輪機轉速，r/min；V為絕對速度，m/s；W為相對速度，m/s；μ為孔口出流的流量系數，此時情況與大孔口出流近似，流量系數μ可取0.85～0.9，考慮到葉片為流線型，可取為0.9[25]；H為工況水頭，m。

將式（4）整理，零流量時的單位轉速可表達為

式中μ′為相對零流量與絕對零流量轉化的修正系數，可取 0～0.4，對于常規混流式水輪機，零流量時的單位轉速一般略大于飛逸單位轉速，以此為判斷依據，通過試算確定修正系數。

在零流量處，轉輪轉速較大，高速的水流使得葉片進口處存在低壓區，會誘導流體在低壓區形成渦流產生水力損失，故很難求出零流量處的單位力矩理論解[22]。

1.1.3 飛逸條件

由于水輪機廠家只提供某一范圍開度的飛逸工況散點數據，需要對飛逸特性曲線進行擴展。當機組飛逸且導葉完全關閉時，單位轉速和單位流量均為0，以此作為小開度區域延拓的依據。飛逸曲線表達如式（6）所示[18]。

式中nc為飛逸狀態下機組的瞬時轉速，r/min；Qc為飛逸狀態下機組的瞬時流量，m3/s；b0為導葉高度，m；r2為轉輪出口半徑，m；A2為轉輪出口斷面面積，m2；β2為水流出口角，(°)。

由于水輪機參數信息往往難以獲得，且當水輪機發生飛逸時，機組振蕩劇烈，流量與水頭等參數出現大幅振蕩[26-28]，此時遠離最優工況，方程中內特性參數與實際參數也存在偏差，故很難直接使用該式指導飛逸曲線的擴展。

補充零開度的飛逸流量與飛逸轉速后，可通過式（7）擬合飛逸條件。飛逸轉速擬合系數可利用 MATLAB的nlinfit函數求解，求解時需保證至少已知5個飛逸試驗點；飛逸流量擬合系數可通過最小二乘法求解。

式中α為導葉開度，(°)；C1Q,C2Q,B0n,B1n,B2n,B3n,B4n為擬合系數；n11c為飛逸狀態下機組的瞬時單位轉速，r/min；Q11c為飛逸狀態下機組的瞬時單位流量，m3/s。

1.1.4 零開度線

零開度時，流量為0；將Q11|α=0=0代入式（1）可知，此時力矩特性曲線為關于單位轉速的二次曲線，如式（8）所示。

1.1.5 交點條件

假定在中高轉速區域，單位力矩也為關于單位轉速的二次曲線[9]，結合零開度線的單位力矩表達式M11=a7n112，各條力矩特性曲線的交點可寫作（n11p，a7n11p2），通過最小二乘法可求出各條開度線的擬合系數和交點邊界，稱為力矩二次曲線模型，如式（9）所示。

式中BM1i，BM2i為第i條開度線的擬合系數。

1.2 簡化內特性模型與邊界條件計算

上述邊界條件中，零轉速時的流量和飛逸條件可通過外特性的試驗數據直接計算，而零轉速時的單位力矩、零流量時的單位轉速、零開度時的力矩特性曲線和交點邊界條件表達式中包含內特性模型中的參數a1i、a7、βb1、Δα0，由于原內特性模型方程求解易陷入病態，精確解難以求取[24]，不宜直接使用，故對內特性模型進行改進，選取水輪機高效率區試驗數據辨識時，葉片水流不發生脫流，可取葉片出口安放角與葉片出流角相同，即βb2≈β2，對式(1)簡化，對參數a1～a7進行合理消元、替換，統一方程量綱，如式（10）所示。

式中Y=M11/n11Q11，kg/m2；X=Q11/n11，m3·min/(r·s)；b1,b2為結構參數，kg/m3，具體表達式見文獻[24]；Δα0為導葉出流角α0與導葉開度角α的差角，(°)；ρ為水密度，kg/m3；g為重力加速度，m/s2；容積效率η0近似取為0.995[29]。

于是需要辨識的參數有：b1，b2，Δα0，a2，a3i，βb1，a7。式（10）可通過最小二乘法求解，應當注意的是，針對每一條開度線i，其參數a1i和a3i需視作獨立的變量。

參數辨識精度直接影響邊界條件計算結果的準確性。本研究給出了參數βb1的驗證方法：葉片入口安放角βb1可通過速度三角形近似計算驗證：考慮到水輪機在偏離最優工況較遠的高轉速區域的水力性能，取葉片安放角βb1比水流進口角β1略小3°～10°，而水流進口角β1表達式可根據速度三角形關系和式（4）進行推導，結果如式（11）所示，代入工況參數即可求解。

2 多重邊界條件下分區擬合特性曲線

2.1 分區處理策略

根據水輪機廠家提供的模型綜合特性曲線和飛逸特性曲線，對特性曲線進行外延和內插，其中外延是對已有試驗數據的開度線沿單位轉速軸向兩端延拓；內插是根據零開度條件和外延的開度線，插值沒有試驗數據的中間開度線。

特性曲線在小開度區域、低轉速區域和中高轉速區域的性質差異顯著，為了精確描述區域特征，以模型綜合特性曲線外圍和 5個邊界條件為界，將水輪機特性曲線分作5個區域，如表1所示。

表1 水輪機特性曲線分區Table 1 Partitioning of turbine characteristic curves

由于流量調節方程不適用于高轉速區域，力矩二次曲線模型在低轉速區域誤差明顯，能量平衡方程僅在制動工況前描述較為準確，綜合考慮，引入不同擬合方法擬合各區的特性曲線，區與區之間，用三次B樣條曲線[30]將各區的曲線段重新擬合，平滑連接，以解決分區處理在分界點處出現的插值誤差。

流量特性曲線外延內插的次序為：

1）結合零轉速條件，通過流量調節方程與能量平衡方程，將曲線外延至S2區域；

2）結合飛逸條件，通過能量平衡方程與力矩二次曲線模型，將曲線外延至S3區域；

3）結合飛逸和零開度條件，采用三次多項式分片擬合的方法，對飛逸工況前已外延的開度線進行插值；

4）結合零流量條件，通過三次B樣條曲線，同時實現S5區域的外延和全區域的平滑連接。

力矩特性曲線外延內插的次序為：

1）結合零轉速條件，通過流量調節方程與能量平衡方程，將曲線外延至S2區域；

2）結合飛逸和交點條件，通過力矩二次曲線模型，將曲線外延至S3區域和S5區域；

3）結合零轉速、飛逸和交點條件，通過三次B樣條曲線，對全區域內已外延的開度線平滑連接；

結合零開度條件，采用三次多項式分片擬合的方法對已外延的開度線進行內插。

2.2 分區處理

2.2.1 S2區域外延

在S2區域，水輪機效率較低，內特性模型與水體實際流動情況存在偏差，故在S2區域的外延中只使用內特性方程的形式，忽略部分參數的物理意義，舍棄距離S2區域較遠的飛逸試驗數據，僅使用高效率區試驗數據及零轉速條件求解。經過大量試驗與方法比選，最終僅保留原內特性方程組中容積效率項a5，將其余內特性參數替換為擬合參數，此時擬合精度較高。將式（1）簡化為式（12）：

式中A1i~A4i為第i條開度線的擬合系數；A6與A7為擬合系數；a5=30ρgη0/π，η0≈0.995。

式（12）可通過最小二乘法求解。方程求解后均可整理為M11=f(Q11,n11)的形式，聯立消去M11后，可看作是以n11為參數，Q11為變量的一元三次方程，利用盛金公式[31]進行求解。

2.2.2 基于力矩二次曲線模型的外延

基于此模型，流量特性曲線可外延至S3區域，力矩特性曲線可外延至S3與S5區域。相較于S2區域，本部分外延的不同點在于流量調節方程已不再滿足，需將式（12）替換為式（9）的力矩二次曲線模型，求解方法與S2區域類似。

2.2.3 力矩特性曲線段的平滑連接

除小開度區域外，各區域內的力矩特性曲線延拓完成后，需要連接各區域的曲線段，在全區域內使用三次B樣條曲線將力矩特性曲線段重新擬合連接。

2.2.4 特性曲線的內插

此時，需要對飛逸工況前的流量特性曲線和全區域內的力矩特性曲線進行內插?？墒褂萌味囗検椒制逯档姆绞綄崿F：將單位轉速和開度歸一化處理，以單位流量和單位力矩為變量，補充零開度條件與插值開度的零流量條件，插值計算。

2.2.5 流量特性曲線在S5區域的外延與平滑連接

S5區域偏離高效率區較遠，水力性能差，每條開度線僅包含飛逸條件和零流量條件2組數據，數據過少，很難結合理論公式外延，故結合飛逸工況前的計算結果，采用三次B樣條曲線對全區域內的流量特性曲線擬合外延。這樣將擬合范圍從 S5區域擴展到全區域，既通過補充特性曲線數據，將流量特性曲線合理外延至S5區域，又實現了飛逸工況前各區流量特性曲線的平滑連接。

3 實例計算

3.1 實例模型及方法

水輪機HLD563-F13模型轉輪進口直徑0.365 m，模型幾何比尺為19.44。其模型綜合特性曲線及飛逸特性曲線由水輪機廠家開展模型試驗得到，如圖2所示。采用Delaunay三角剖分-分片三次多項式插值[15]加密等效率線，共提取效率線與開度線的交點 215個，飛逸特性曲線數據16個，即試驗數據共231個。按本文方法、典型外特性法[9]和內特性法[24]對該特性曲線進行外延內插，并對比擬合結果以及過渡過程實測反演結果。各方法的計算流程見圖3。

3.2 內特性方程參數辨識與邊界條件

3.2.1 內特性方程參數辨識

內特性參數辨識結果及對比見表2和表3，其中典型內特性法未直接計算出葉片入口安放角βb1和參數a3i，可通過最小二乘法進一步計算得到。

參數b1和b2用于計算零轉速條件，βb1用于計算零流量條件，a7用于計算零開度線和交點條件，其余辨識的參數只包含在簡化內特性模型中，并未參與邊界條件計算。由表2和表3對比參數辨識結果可知，參數b1、b2、a3i和βb1的偏差較大，其中b1結果相差112.07 kg/m3，b2結果相差82.78 kg/m3，βb1結果相差36.45°，a3i隨著開度增加，結果偏差逐漸變小，最大偏差為5.39×105kg/m3。上述偏差較大的參數中，βb1在邊界條件計算中最為重要，對其進行驗證：以額定工況為計算工況，將相關參數代入式（11）：額定開度αr=34.74°，額定轉速nr=100 r/min，額定水頭Hr=80 m，原型轉輪進口直徑D1=7.1 m，假定水輪機導葉為標準化導葉，取D0/D1=1/1.1，可近似計算得β1=60.12°，即βb1的合理取值范圍為 50.12°～57.12°。由此可知本文方法的參數辨識結果精確性更高。

表2 不同方法下部分水輪機參數辨識結果及對比Table 2 Identification results and comparison of turbine characteristic parameters by different methods

表3 不同開度及方法的參數a3i辨識結果及對比Table 3 Identification results and comparison of a3i in different GVO by different methods(×105kg·m-3)

3.2.2 邊界條件

根據對應公式和參數辨識的結果，取流量系數為0.9，在滿足零流量時的單位轉速略大于飛逸時的單位轉速前提下，近似取修正系數為0.3，各開度的零轉速和零流量時單位參數計算結果如表4所示。

表4 零轉速條件和零流量條件下單位參數計算結果Table 4 Results of unit parameters in zero discharge and zero speed conditions

對于飛逸條件：補充零開度的飛逸數據后，代入外特性試驗結果，根據式（7）進行擬合，沿開度方向擴展的結果如圖4所示，兩曲線的擬合決定系數均為0.9997，擬合精度高，擴展結果可靠。

對于零開度條件：代入參數辨識的結果a7=?0.04 kg/m3，可得零開度時的流量特性曲線為Q11|α=0=0,力矩特性曲線為M11|α=0=?0.04n112；對于交點條件，代入模型綜合特性曲線和飛逸特性曲線試驗數據后，根據式（9），采用最小二乘法求解，可得交點條件(n11p，a7n11p2)結果為（175.31, ?1215.13）。

3.3 水輪機特性曲線外延內插結果

水輪機廠家提供了從開度10°開始，每隔2°，直到最大開度40°，共16條開度線的外特性試驗數據，圖5a與圖5b為此16條開度線外延的結果；圖5c與圖5d為從0°起，每隔1°，直到最大開度40°的開度線內插結果，其中開度線均從下至上依此增大。

特性曲線外延內插結果中開度線間均無交叉，同一開度線在分區邊界位置未出現不連續現象；且在邊界條件與試驗數據的共同約束下，擬合精度較高，各區特性曲線演化規律能較好反映水輪機過流特性[32-33]：在低轉速區域，單位流量較為平穩，單位轉速對單位流量的影響不大，單位流量主要由導葉開度決定，進入中高轉速區域后，開度線開始有下降趨勢；在制動工況區域，各條開度線上的單位流量隨著單位轉速的增加而逐漸靠近，然后急劇降低，通過零流量邊界后改變符號。

將3種方法對特性曲線的處理結果進行對比。由于經內插處理的開度線排列過于緊密，不便對比，僅節選開度為 2°、6°、10°、16°、22°、28°、34°和 40°的特性曲線對比,其中10°～40°的開度線為外延結果，有試驗數據；2°和6°的開度線為內插結果，無試驗數據。試驗數據見圖5a與5b。試驗數據處的擬合決定系數為見表5。各方法在零轉速和飛逸條件處計算的單位流量與單位力矩的誤差對比如圖6所示，其中零轉速條件參考值由式（2）與（3）確定，飛逸條件參考值由式（7）確定，對比表明：

表5 單位流量與單位力矩試驗數據處的決定系數Table 5 Determination coefficient of unit discharge and unit torque in experimental data

1）在試驗數據處，本文方法既考慮了內特性模型，又有效提升了擬合精度：與典型內特性法相比，在試驗數據處單位流量擬合決定系數提高了2.69%，單位力矩擬合決定系數提高了11.02%。

2）本文方法考慮了零轉速條件，有效降低了零轉速處的誤差：相對于典型內特性法和典型外特性法，本文方法在零轉速條件處單位流量的平均誤差分別降低0.12、0.09 m3/s，在零轉速條件處單位力矩的平均誤差分別降低192.44、215.94 N·m；

3）通過分析制動工況區域的特點和外特性法計算方法可知，典型外特性法單位力矩的飛逸位置與交點位置單位轉速明顯偏高，而典型內特性法未計算到飛逸位置與交點條件；與典型外特性法相比，本文方法在大開度飛逸條件處單位力矩的平均誤差降低89.68 N·m；

綜上分析，本文方法綜合考慮多重邊界條件、外特性試驗數據、簡化內特性模型，采取分區處理方法實現，優勢顯著。

3.4 過渡過程試驗反演

針對某電站實際情況，采用本文方法與典型外特性法計算的特性曲線進行過渡過程試驗反演與對比，典型內特性法延拓范圍不足，較難滿足過渡過程仿真計算要求，故未比較。電站為“單管單機”引水布置，引水發電系統管道總長 1 087.7 m，水輪機組飛輪力矩為127 500 t·m2，工作水頭為77 m，水輪機轉輪進口直徑為7.1 m，額定出力為367 MW，額定轉速100 r/min，甩50%負荷過渡過程，比較 2種特性曲線計算過渡過程的蝸殼壓力，其時域變化過程及 60～70 s內的精度分析與對比如圖7所示。兩種方法計算的蝸殼壓力變化趨勢與實測值基本一致；實測結果最大蝸殼壓力為1.182 MPa，本文方法計算結果為1.202 MPa，與實測結果基本吻合；典型外特性法為 1.206 MPa；最大蝸殼壓力最大誤差從0.024 MPa降低至0.020 MPa，相對誤差從2.03%降低至1.69%；隨導葉關閉，水輪機進入小開度工況區域，61.4 s時，導葉完全關死，產生較大水擊壓強，蝸殼壓力振蕩明顯。在此振蕩區域內，本文方法計算結果振蕩幅值更接近實測結果，反演效果更佳：其絕對平均誤差為0.010 MPa，典型外特性法絕對平均誤差為 0.023 MPa，降低0.013 MPa，相對誤差從3.48%降低至1.47%。

綜上分析，本研究提出的特性曲線處理方法在過渡過程試驗反演的最大蝸殼壓力及小開度區域動態過程時域響應均更加貼近實測結果。

4 結 論

本研究綜合考慮多重邊界條件、外特性試驗數據與內特性模型，提出了一種特性曲線分區處理方法。該方法有效結合了內外特性的優點，既保證了外特性試驗數據處的擬合效果，又能反映水輪機的水力特性。

特性曲線計算結果相比典型外特性法，本文結果在飛逸條件處的單位力矩平均誤差降低89.68 N·m，有效提升制動工況區域擬合精度；相比典型內特性法，在外特性試驗數據處的單位流量擬合決定系數提高了2.69%，單位力矩擬合決定系數提高了11.02%。，并可延拓特性曲線至制動工況區域。工程實例計算結果表明，本文方法與典型外特性法相比，蝸殼壓力極值最大誤差從2.03%降低至1.69%，小開度區域蝸殼壓力反演相對誤差從3.48%降低至1.47%，動態過程時域響應更接近實測結果。