用三維曲面插值技術(shù)求解水輪機實時運行參數(shù)

結(jié)少鵬，虞 鴻

(浙江省水利水電勘測設(shè)計院，浙江 杭州 310002)

在水利工程規(guī)劃設(shè)計和實際運行調(diào)度中一些關(guān)鍵參數(shù)或結(jié)果的獲得常常需要查找特定的已知數(shù)據(jù)[1- 2]、參數(shù)表或曲線圖，如實用堰淹沒系數(shù)曲線[3]、水輪機模型綜合特性曲線[4]等；另一方面，一些復雜方程的求解也需要求助于圖解法，如矩形斷面明渠底流消能水力計算、梯形斷面渠道正常水深、臨界水深等。雖然一些可以用復雜公式(如隱式公式)表達的問題目前已經(jīng)能通過程序進行查表和插值[5- 6]，甚至可以通過試算、迭代等數(shù)值方法求解，簡化了大量設(shè)計人員的工作，但還有一些問題因為無法找到對應(yīng)的數(shù)學表達式，必須查找實驗或經(jīng)驗圖形曲線才能獲得對應(yīng)參數(shù)。具體來說，模型試驗繪制的曲線，如水輪機轉(zhuǎn)輪綜合特性曲線、實用堰淹沒系數(shù)曲線等，受到眾多參數(shù)的影響，往往難以用簡單函數(shù)或其組合進行良好的擬合，每次參數(shù)查找基本都需要由設(shè)計人員根據(jù)曲形進行取值，有較大的隨意性，而且往往需要根據(jù)多條曲線上的鄰近數(shù)據(jù)進行單向或雙向線性插值計算，在種種限制條件下查詢的過程不僅效率不高，獲得的結(jié)果也有較大的隨機誤差，也難以用程序?qū)崿F(xiàn)計算過程。

在當前水利數(shù)字化和信息化蓬勃發(fā)展的背景下[7- 12]，傳統(tǒng)的參數(shù)曲線處理方法已經(jīng)不能很好地適應(yīng)程序化的處理方式和實時精細化控制對參數(shù)的要求。目前用于計算的圖表通常表現(xiàn)為笛卡爾坐標系中的等值曲線簇，每條線上的各個點都可以認為是分別對應(yīng)縱、橫軸自變量的函數(shù)值，因此上述曲線簇可以認為是來自于一個復雜三維曲面的成果；如水電站工程中關(guān)鍵設(shè)備水輪機的運轉(zhuǎn)綜合特性曲線圖，包含橫坐標“功率”、縱坐標“水頭”及等效率值曲線簇。對水輪機實時工況點效率值的計算，在數(shù)學上就構(gòu)成三維曲面問題，可以通過克里格方法將類似曲面數(shù)據(jù)網(wǎng)格化，使之能用程序處理。

1 克里格(Kriging)方法簡介

克里格法是比較靈活和準確的網(wǎng)格化插值方法之一，它可以通過在整體預測中給予較小的權(quán)重來補償集群數(shù)據(jù)。

Z(p)為區(qū)域Ω上的一個隨機過程，p∈Ω。Z(p)在p，p+h兩點處的值之差的方差之半為Z(p)在p方向上的變差函數(shù)γ(h)。

(1)

Z(p)在p1，p2處的兩個隨機變量Z(p1)和Z(p2)的二階混合中心距，即：

Cov{Z(p1)，Z(p2)}=E[Z(p1)·Z(p2)]-E[Z(p1]·E[Z(p2)]

(2)

Ω上有n個樣本點，Zi=Z(pi)在pi處的值，則p0處的最優(yōu)線性估計值為

(3)

(4)

求解方程組，可得

(5)

2 算例分析

2.1 工程概況

某水庫工程以供水、防洪為主，兼顧灌溉、發(fā)電及改善河道水環(huán)境等綜合利用。水庫配套水電站型式為地面引水式，發(fā)電廠房位于大壩左岸下游側(cè)。電站裝設(shè)3臺臥式混流式水輪發(fā)電機組，單機容量1200kW。水輪機運轉(zhuǎn)綜合特性曲線如圖1所示。

工程正在建設(shè)水庫智慧運行管理平臺，基于設(shè)備在線監(jiān)測系統(tǒng)及平臺計算分析功能，精確掌握各類機電設(shè)備的運行狀態(tài)數(shù)據(jù)，以便進行智能運行管理。為了滿足電站精細化管理的需求，需要實時計算水輪機運行工況點效率，同時計算水輪機的過流量，并將工況點實時展示在運轉(zhuǎn)綜合特性曲線圖中。由于機組容量較小，廠房空間有限，機組未裝設(shè)流量和效率監(jiān)測設(shè)備，因此水輪機流量和效率數(shù)據(jù)的采集是智慧管理平臺開發(fā)的一個難點。

在電站的自動化監(jiān)測系統(tǒng)中，可以實時讀取的參數(shù)有：發(fā)電機功率Pg，蝸殼進水口壓力PW，下游尾水位ZT等。

水輪機的過流量QT可通過式(6)計算：

(6)

(7)

式中，Pg—發(fā)電機功率，kW;ηg—發(fā)電機效率，可近似為常數(shù)。

H=PW+ZW-ZT

(8)

式中，PW—蝸殼進水口壓力，m;ZW—蝸殼進水口高程，為常數(shù)，m；ZT—電站尾水位，m。

經(jīng)以上初步分析可知，通過電站現(xiàn)有的監(jiān)測系統(tǒng)參數(shù)，可以間接計算得到水輪機運轉(zhuǎn)綜合特性曲線圖中工況點的橫、縱坐標值(即：水輪機出力PT和凈水頭H)，只要能計算出水輪機的實時效率值，平臺開發(fā)的數(shù)據(jù)問題就可迎刃而解。而水輪機效率值的計算必須采用程序化的處理方法，以滿足流量數(shù)據(jù)的實時計算和實時顯示。

2.2 等值線的處理及選點

水輪機運轉(zhuǎn)綜合特性曲線是由模型試驗數(shù)值制成的。對于模型試驗難以模擬的過渡過程工況，曲線圖上有些區(qū)域是空白的。傳統(tǒng)的人工計算方法下，水輪機效率值是根據(jù)出力PT和凈水頭H從運轉(zhuǎn)綜合特性曲線近似查得，遇到的也是穩(wěn)態(tài)常規(guī)運行工況，故不需要對等效率線進行特殊處理。

而智慧管理平臺在計算水輪機效率時，會涉及水輪機的開機、停機等過渡過程工況，如果在這個區(qū)域內(nèi)沒有曲線數(shù)據(jù)，則會顯示錯誤。因此，需要對運轉(zhuǎn)綜合特性曲線圖進行處理，用專業(yè)經(jīng)驗按照曲線簇的變化趨勢，將外圍的等效率曲線補全，保證水輪機運行范圍內(nèi)的效率數(shù)據(jù)是完整的，以便計算機能夠正常計算處理。處理后的運轉(zhuǎn)綜合特性曲線圖如圖2所示。

圖1 水輪機運轉(zhuǎn)綜合特性曲線

圖2 處理后的水輪機運轉(zhuǎn)綜合特性曲線

為了將等效率曲線簇生成三維曲面，需要將等效率曲線的坐標點數(shù)據(jù)導出至Civil 3D軟件中，因此需要先進行選點。每條等效率線上的點的選擇應(yīng)該能反映曲線的走勢，并保證曲線的精度。在曲線的平直段，點的設(shè)置可以適當稀疏一些，而在曲率半徑小的位置，點的設(shè)置應(yīng)當加密。每條曲線上的點應(yīng)將其Z坐標設(shè)置為效率值。水輪機運轉(zhuǎn)綜合特性曲線的點的設(shè)置如圖3所示。

圖3 等效率線的點設(shè)置

通過list命令將CAD中綜合特性曲線上點的x(對應(yīng)PT)、y(對應(yīng)H)、z(對應(yīng)ηT)坐標值導出提取至Excel表中，進行適當?shù)淖鴺俗兓纬?PT，H，ηT)數(shù)據(jù)表，見表1。

表1 從CAD中提取的綜合特性曲線數(shù)據(jù)

2.3 三維曲面生成及柵格數(shù)據(jù)的輸出

將表1中的特性曲線等值線數(shù)據(jù)導入至Civil 3D軟件中，以便生成三維曲面。曲面生成過程中，選擇自然鄰近內(nèi)插法或克里格方法對曲面進行光滑處理以提高曲面精度。由此，綜合特性曲線圖中離散的等值線簇轉(zhuǎn)變?yōu)檫B續(xù)的三維曲面。如圖4所示。

圖4 克里格方法處理的三維曲面圖(運轉(zhuǎn)綜合特性曲線)

三維曲面生成后利用指定的坐標生成點，可以直接取得曲面Z坐標的數(shù)值，就完成了一次數(shù)據(jù)查詢，完全免去了插值計算的過程。

利用Civil 3D軟件的按指定柵格處生成點命令，將離散成滿足使用精度要求的網(wǎng)格化數(shù)據(jù)(如圖5所示)，并按指定的格式輸出“X-Y-Z”的數(shù)據(jù)文件供軟件開發(fā)使用(實際文件是按地形數(shù)據(jù)“東距-北距-高程”等格式導出的)。原有的三維曲線數(shù)據(jù)點是粗線條的、不均勻的、不規(guī)律的，經(jīng)過克里格方法轉(zhuǎn)換以后，數(shù)據(jù)點變?yōu)榫毜那揖鶆蚍植嫉木W(wǎng)格化數(shù)據(jù)，便于計算機進行程序化的查詢計算。Civil 3D軟件導出的數(shù)據(jù)見表2。

圖5 水輪機運行范圍內(nèi)的柵格設(shè)置

表2 Civil 3D導出的網(wǎng)格化數(shù)據(jù)表

將Civil 3D軟件導出的特性曲線網(wǎng)格化數(shù)據(jù)上傳至智慧管理平臺，作為水輪機實時效率計算的基礎(chǔ)數(shù)據(jù)。

2.4 水輪機實時參數(shù)計算

以表2中的一系列柵格數(shù)據(jù)為基礎(chǔ)，根據(jù)水輪機的實時工況參數(shù)(出力PT(坐標值x)、實時凈水頭H(坐標值y))，計算水輪機實時工況點p與每個格柵點i的距離Li：

(9)

式中，Li—工況點p與格柵點i的距離；xp—實時工況點p的橫坐標；yp—實時工況點p的縱坐標；xi—格柵點i的橫坐標；yi—格柵點i的縱坐標。

通過式(9)的排序計算，找到離實時工況點p最近的4個格柵點A、B、C和D，然后選擇最近4個格柵點的效率平均值，作為實時工況點p的效率值ηp。適當提高格柵矩陣的密度會提高計算精度。

ηp=(ηA+ηB+ηC+ηD)/4

(10)

得到水輪機實時效率值ηp后，即可按式(1)計算水輪機實時過流量QT。如此可得到水輪機的所有實時運行參數(shù)，滿足工程運行管理的需求。

3 結(jié)語

(1)水利工程的建設(shè)中涉及較多的復雜曲線圖，這些參數(shù)曲線難以找到對應(yīng)的數(shù)學表達式或擬合函數(shù)。本文提出的程序化的插值計算方法，解決了傳統(tǒng)人工讀數(shù)和插值計算精度不高且效率低下的問題。

(2)在綜合特性曲線處理時，補全等效率曲線簇需要專業(yè)經(jīng)驗，依賴執(zhí)行人的判斷力，精度難以控制。但考慮到水輪機在此區(qū)域內(nèi)運行的時間很短，處理誤差對整體的監(jiān)控效果影響十分有限。

(3)現(xiàn)有的水輪機綜合特性曲線都是根據(jù)試驗數(shù)據(jù)繪制而成的。而實際上，模型試驗時工況點的位置很難剛好處于等效率線上，等效率線也是由試驗人員借助經(jīng)驗描繪，本身就存在一定的誤差。若采用本文所闡述的方法，可以直接將所有試驗工況點的原始數(shù)據(jù)導入Civil 3D軟件，不需繪制等效率線，可以避免人為因素導致的隨機誤差。