不同體型調(diào)壓室的瞬態(tài)水力特性對比分析

何相慧，楊建東，楊桀彬，胡金弘，鄭貴橋，向正林

（1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072；2.南方電網(wǎng)調(diào)峰調(diào)頻發(fā)電有限公司，廣東 廣州 510630）

1 研究背景

設(shè)置調(diào)壓室是改善長輸水管道水電站運行條件的可靠措施，既可以減少管道的水錘壓力，又可以防止由于尾管真空度過高引起的水柱分離，從而提高水電站的運行安全性[1-2]。地形因素對調(diào)壓室施工有制約作用，影響了調(diào)壓室的形狀和調(diào)壓室斷面積，如何在滿足水電站安全運行的前提下，合理設(shè)計調(diào)壓室體型，是每個水電站設(shè)計中關(guān)鍵技術(shù)問題之一[3-4]。

在數(shù)值模擬方面，趙志高等[5]建立了基于電路等效理論的有壓管道系統(tǒng)數(shù)值仿真模型，有效協(xié)調(diào)了仿真精度與計算效率之間的矛盾。Harlow等[6]提出的特征線方法（MOC）被廣泛用于有壓管道系統(tǒng)的瞬態(tài)特性研究中，用于求解偏微分方程組[7-11]。陳玲等[9]對比了解析法、數(shù)值積分法和特征線法在求解調(diào)壓室涌浪過程中的區(qū)別，結(jié)果顯示特征線法考慮了水電站管路布置特性，理論精度較高，適用于大中型水電站。Yang等[10-11]將水輪機特性曲線方程改寫為空間曲面方程，并提出了基于轉(zhuǎn)速偏差函數(shù)機組邊界條件的尋根方法。盡管MOC方法可得到調(diào)壓室水位的數(shù)值解，但無法得到調(diào)壓室水位波動過程中內(nèi)部流態(tài)，該流態(tài)較為復(fù)雜多變，若產(chǎn)生旋渦、脫流、進氣等有害現(xiàn)象，會影響水電站的安全穩(wěn)定運行。于是，近年來國內(nèi)外許多學(xué)者采用CFD方法對調(diào)壓室水位波動過程及瞬變水力特性進行了研究。CAI F等[12]通過試驗和CFD數(shù)值模擬研究了立軸旋渦的產(chǎn)生原因和優(yōu)化方法，結(jié)果表明淹沒水深是影響立軸旋渦產(chǎn)生的主要因素，通過減小弗勞德數(shù)Fr和速度環(huán)量可有效抑制立軸旋渦的產(chǎn)生；鄧淞苡等[13]探討了長廊式調(diào)壓室不利流態(tài)的產(chǎn)生原因，并提出通過增加“消渦墩”的方式消除立軸旋渦；華富剛[14]討論了長上室調(diào)壓室的水面瞬態(tài)波動，發(fā)現(xiàn)CFD 計算值和測試值吻合良好，并通過特征線法與CFD方法結(jié)合進行了水電站過渡過程計算。劉飛等[15]探討了不同湍流模型和邊界條件對計算結(jié)果準(zhǔn)確性的影響，結(jié)果表明RNGk-ε湍流模型和Realizablek-ε湍流模型對強旋流的適應(yīng)性較好，預(yù)測精度能夠滿足工程要求。另外，多位學(xué)者的研究表明，CFD計算結(jié)果與試驗結(jié)果具有良好的吻合性，采用CFD 方法進行調(diào)壓室瞬變水力特性模擬具有良好的精度，可為工程實踐提供參考[16-19]。

一維特征線法計算調(diào)壓室涌浪水位已經(jīng)得到了廣泛的應(yīng)用和證實，成功地應(yīng)用于國內(nèi)外多座水電站的設(shè)計中，隨著CFD技術(shù)的發(fā)展和計算機運算能力的提升，利用CFD計算調(diào)壓室涌浪及局部流態(tài)的準(zhǔn)確性也得到了驗證，但是很少有論文對一維與三維的計算結(jié)果進行對比，本文所研究的調(diào)壓室體型也鮮有報道。因此，基于以上研究成果，本文采用一維數(shù)值仿真和三維CFD 對某水電站T型調(diào)壓室和π型調(diào)壓室兩種體型進行了數(shù)值模擬，兩種體型的調(diào)壓室斷面積相同。計算工況為雙機甩負荷工況和一臺機滿載另一臺機啟動工況（以下簡稱甩負荷工況和啟動工況），并通過對比其調(diào)壓室涌浪和內(nèi)部流態(tài)等得到較優(yōu)的調(diào)壓室斷面形狀，為工程實踐提供參考。

2 一維MOC計算模型

2.1 有壓管道的特征線法有壓管道彈性水錘的基本方程由運動過程和連續(xù)性方程組成：

式中：V為管道中的流速，由上游流向下游為正；t為時間；x為距離管道最左端的距離；g為重力加速度；f為沿程損失系數(shù)；D為管道直徑；H為水頭；a為水錘波波速。

采用MOC方法將公式轉(zhuǎn)換為兩組特征線上的常微分方程，如下所示

式中：HP為測壓水頭；SA和SB為截面周長；QP為流量。

特征線法可以解決多種邊界問題，不僅可以合理地反映水電站管路布置特點，也可以方便地考慮水流慣性、管壁彈性及摩阻的影響，便于編程實現(xiàn)計算機求解。

2.2 調(diào)壓室計算模型Topsys-TP是武漢大學(xué)開發(fā)的水電站過渡過程一維計算軟件[10-11]，已成功應(yīng)用于國內(nèi)外近百座水電站的設(shè)計。本文采用Topsys進行一維計算，其計算模型如圖1所示。T型調(diào)壓室（圖2）和π型調(diào)壓室（圖3）的算法相同，邊界條件有所差異。下面以T型調(diào)壓室為例，介紹其數(shù)學(xué)模型。

以國內(nèi)某水電站的T 型阻抗式調(diào)壓室為原型，可列出13個未知數(shù)HP1、HP2、HP3、HP4、QP1、QP2、QP3、QP4、HTP1、HTP2、QTP1、QTP2、Z，對應(yīng)的邊界條件是：

圖1 Topsys計算簡圖

圖2 T型截面調(diào)壓室模型

圖3 π型截面調(diào)壓室模型

特征線方程：

調(diào)壓室水位方程：

能量方程：

連續(xù)性方程：

式中：QP1、QP2、QP3、QP4、QTP1和QTP2為對應(yīng)位置的流量；HP1、HP2、HP3、HP4和HTP1為壓力水頭；Z為調(diào)壓室底板高程；QCP、CQP、QCM和CQM為上一時刻的已知量。

通過求解以上方程，可得到各個物理量的值。

而π型調(diào)壓室中間設(shè)置隔板，隔板上設(shè)置4個連通孔，當(dāng)水位低于94 m時，調(diào)壓室分為2個獨立的調(diào)壓室，當(dāng)水位高于94 m時，兩側(cè)通過連通孔實現(xiàn)水流交換，采用堰流公式進行計算。

3 CFD數(shù)值計算模型

3.1 有壓管道的特征線法有壓管道彈性水錘的基本方程由運動過程和連續(xù)性方程組成：

直角坐標(biāo)系下質(zhì)量守恒方程（連續(xù)性方程）的微分形式為：

若流體流速較小不考慮其可壓縮性時，流體的密度為常數(shù)，上述微分方程可改為：

動量守恒方程即N-S方程，若在流動過程中流體密度和黏性保持不變，其動量守恒方程表達式為：

Realizablek-ε模型在模擬強逆壓力梯度、射流擴散率、分離、回流、旋轉(zhuǎn)上有較高精度，在水電站過渡過程計算中與試驗數(shù)據(jù)吻合性較好[13-16]。其湍動能k及耗散率ε輸運方程為：

其中：

式中：ρ為密度；Gk和Gb為剪切產(chǎn)生項；ut為湍動黏度；YM為可壓縮湍流脈動膨脹對總的耗散率的影響；C1ε、C2、σk、σε作為默認值常數(shù)，C1ε=1.44，C2=1.92，σk=1.0，σε=1.3。

3.2 計算模型與邊界條件某水電站上游調(diào)壓室為阻抗式調(diào)壓室，水頭低，調(diào)壓室穩(wěn)定斷面積大，受地形限制無法設(shè)置雙調(diào)壓室，因此設(shè)置長上室補充調(diào)壓室穩(wěn)定斷面積，長上室底面設(shè)置一定的坡度。在滿足調(diào)壓室穩(wěn)定斷面面積的前提下，考慮地形因素和水力條件，調(diào)壓室橫截面設(shè)計為T型和π型，如圖4所示。該水電站調(diào)壓室大井由開敞段和隧洞段組成，T型調(diào)壓室大井由直方形開敞段及一條長隧洞組成，π型調(diào)壓室大井由直方形開敞段和兩條隧洞組成，開敞段中間設(shè)置隔墻分離，隔墻上設(shè)置4個連通孔。調(diào)壓室底板高程87 m，總高度35 m，連接管直徑7 m，長度31 m，計算管道長711.375 m。

采用Star-CCM+進行網(wǎng)格劃分和數(shù)值計算，其可自動生成高質(zhì)量的多面體網(wǎng)格，網(wǎng)格最大值設(shè)置為1.0 m，連接管及大井敞開段進行網(wǎng)格加密（圖3），總網(wǎng)格數(shù)67.52萬。根據(jù)上游水庫水位將入口設(shè)置為壓力入口。出口設(shè)置為質(zhì)量流量出口，質(zhì)量流量曲線由Topsys導(dǎo)出。調(diào)壓室頂部設(shè)置為壓力出口，相對壓力為0 atm。湍流模型采用Realizablek-ε，固壁邊界設(shè)置為無滑移壁面，近壁區(qū)采用標(biāo)準(zhǔn)壁面函數(shù)法處理。

調(diào)壓室涌浪水位的監(jiān)測點設(shè)置在距離阻抗孔一定距離的位置，左右對稱各設(shè)置3個監(jiān)測點，距底板3 m，監(jiān)測其靜壓值P，并取其平均值計算調(diào)壓室水位，兩種體型監(jiān)測點位置相同，監(jiān)測量通過公式Z=P/9810+Z0轉(zhuǎn)化為壓力值，其中Z0位監(jiān)測點高程。

圖4 CFD模型圖

4 計算結(jié)果與分析

4.1 兩種體型一維大波動計算結(jié)果與分析導(dǎo)葉開啟規(guī)律和導(dǎo)葉關(guān)閉規(guī)律如圖5和圖6所示，一維大波動計算結(jié)果如圖7—13所示。由圖7可以看出，甩負荷工況下，T型調(diào)壓室和π型調(diào)壓室的涌浪波動規(guī)律相同，均隨時間做周期性衰減。由圖8可以看出，啟動工況下，由于水位低于100 m，π型調(diào)壓室分為兩個單獨的調(diào)壓室，Unit2機組正常運行，因此水位不變；Unit1機組由空載啟動，受到調(diào)壓室穩(wěn)定斷面積限制，Unit1 機組側(cè)的最高涌浪比T 型調(diào)壓室高0.34 m，最低涌浪低1.78 m，最小淹沒水深為3.35 m，但其衰減速度較快。因此三維計算需重點關(guān)注π 型調(diào)壓室在啟動工況下調(diào)壓室的流態(tài)。

圖5 機組導(dǎo)葉開啟規(guī)律示意

圖6 機組導(dǎo)葉關(guān)閉規(guī)律示意

圖7 甩負荷工況調(diào)壓室涌浪波動

圖9為甩負荷工況下調(diào)壓室底板壓差，從圖中可以看出，兩種體型的調(diào)壓室底板壓差隨時間變化規(guī)律相同，向上最大為15 m 左右，向下最大為2 m 左右，且向上和向下最大壓差出現(xiàn)的時間相近。啟動工況下（圖10），T型調(diào)壓室Unit2機組側(cè)在0.02 s出現(xiàn)向下最大壓差為5.32 m。圖11—13的調(diào)保參數(shù)波動圖，兩種體型相差較小，且均在調(diào)保參數(shù)控制范圍以內(nèi)。

圖8 啟動工況調(diào)壓室涌浪波動

圖9 甩負荷工況底板壓差

圖10 啟動工況底板壓差

圖11 甩負荷工況蝸殼出口最大動水壓力

圖12 甩負荷工況尾水管最小動水壓力

圖13 甩負荷工況轉(zhuǎn)速波動

4.2 一維與三維調(diào)壓室涌浪對比分析表1為一維與三維調(diào)壓室最高/最低涌浪及周期的對比圖，圖14和圖15為一維與三維涌浪水位隨時間波動的對比。從表1可以看出：由于三維計算調(diào)壓室水面的不規(guī)則波動，曲線圖出現(xiàn)了次波峰波谷，但一維和三維的曲線圖整體波動趨勢相同。由表1初始水位的對比可知，兩種工況兩種體型下，一維計算的初始水位比三維水位高，最大差值為0.22m，是由于一維三維計算的沿程損失差值造成的，可通過調(diào)整管道糙率減小兩者差值，目前計算的差值在誤差范圍內(nèi)。通過對比表1兩種體型的一維和三維最高/最低涌浪以及周期可知：甩負荷工況下，一維和三維最高涌浪最大差值為1.41 m，最低涌浪差值最大為0.3 m；啟動工況下，一維與三維的最高涌浪和最低涌浪的差值均在1 m以內(nèi)。由于三維計算水面波動的不平穩(wěn)性，一維與三維的周期差值無明顯規(guī)律，但結(jié)合圖14和圖15可以看出，一維與三維的波動曲線變化趨勢和衰減速度吻合程度較高。因此，一維計算和三維計算的計算方法是可靠的。

表1 調(diào)壓室涌浪對比

圖14 甩負荷工況涌浪波動

圖15 啟動工況的涌浪波動

4.3 三維計算流態(tài)分析阻抗式調(diào)壓室容易產(chǎn)生立軸旋渦[20]，貫穿型立軸旋渦將會對機組及輸水管道造成嚴(yán)重影響，破壞水體正常流態(tài)，引起空化和振動等。為了研究兩種體型調(diào)壓室的液面波動形態(tài)以及是否會產(chǎn)生立軸旋渦，對其波動水面進行監(jiān)測，結(jié)果如圖16—19所示。

圖16和圖17為甩負荷工況下調(diào)壓室的液面波動情況。甩負荷開始時（即t=0 s時刻），T型調(diào)壓室和π 型調(diào)壓室水面穩(wěn)定。甩負荷開始后，機組導(dǎo)葉迅速關(guān)閉并在13 s 時全部關(guān)閉，機組引用流量減少，有壓隧洞中的水流涌向調(diào)壓室，調(diào)壓室大井阻抗孔上方出現(xiàn)鼓包形狀，由于T型調(diào)壓室上室較長，長上室內(nèi)產(chǎn)生了較大的水面波動，π型調(diào)壓室兩個長上室內(nèi)水流呈現(xiàn)梯度流動，未出現(xiàn)較大的波動。甩負荷工況下，僅在甩負荷初始階段流入調(diào)壓室流量較大時，調(diào)壓室液面擾動劇烈，甩負荷結(jié)束后，調(diào)壓室內(nèi)平面趨于平穩(wěn)。

兩種體型啟動工況下調(diào)壓室的液面波動如圖18和圖19所示，初始時刻，Unit2 機組滿負載運行，π型調(diào)壓室Unit2側(cè)水面較低，Unit1側(cè)調(diào)壓室通過下方兩個連通孔向另一側(cè)補充水流。Unit1機組啟動后，Unit1機組引用流量增加，上游調(diào)壓室補充隧洞中的水流，調(diào)壓室液面下降。T型調(diào)壓室的長上室位于兩阻抗孔中間位置，Unit1 機組啟動初始階段，長上室的水流來不及補充Unit1 側(cè)，而此時流出調(diào)壓室流量較大，在阻抗孔上方形成了吸氣旋渦，吸氣旋渦迅速蔓延到升管內(nèi)（圖18（b））。對于π型調(diào)壓室，由于長上室位于阻抗孔中心線上，流出調(diào)壓室流量較大時，長上室內(nèi)水流可迅速進行補充，水面下降較為平穩(wěn)。

圖16 甩負荷工況下T型調(diào)壓室水面波動圖

圖17 甩負荷工況下π型調(diào)壓室水面波動圖

圖18 啟動工況下T型調(diào)壓室水面波動圖

圖19 啟動工況下π型調(diào)壓室水面波動圖

圖20 啟動工況下T型調(diào)壓室水平截面速度矢量圖（t=20s）

圖21 啟動工況下π型調(diào)壓室水平截面速度矢量圖（t=20s）

為了進一步分析T型調(diào)壓室吸氣旋渦產(chǎn)生的原因，對兩種體型調(diào)壓室的水平截面速度矢量圖進行了對比，如圖20和圖21所示，水平截面距離調(diào)壓室底板3 m。由圖20可知，Unit1機組啟動后，T型調(diào)壓室長上室中的水流流向Unit1 側(cè)阻抗孔，并沿逆時針方向進入阻抗孔，由于上室較長，坡度較小，其尾部段水流基本為死水，在過渡過程后期會由于壓差緩慢流向調(diào)壓室大井。π型調(diào)壓室在過渡過程的初始階段，長上室中水流直線流向阻抗孔，調(diào)壓室大井內(nèi)水流也徑直向阻抗孔位置匯聚，未形成明顯的速度環(huán)量。

4.4 兩種體型對比與分析通過以上對比分析得出，T型調(diào)壓室在啟動工況產(chǎn)生了吸氣旋渦，而π型調(diào)壓室水流狀態(tài)較好，原因如下：（1）速度環(huán)量。T型調(diào)壓室的阻抗孔位于長上室兩側(cè)，當(dāng)流出調(diào)壓室流量較大時，長上室中的水流會順時針或逆時針流入阻抗孔中，在其上方形成速度環(huán)量。π型調(diào)壓室兩長上室位于阻抗孔中心位置，長上室中的水流可徑直進行補充，無法形成速度環(huán)量。（2）流量變化。T型和π型調(diào)壓室的長上室寬度比為0.58（如圖22），其高度相同，因此在相同條件下，π 型調(diào)壓室的過流面積遠遠大于T 型調(diào)壓室，T 型調(diào)壓室長上室的流量無法及時補充調(diào)壓室大井的流量會導(dǎo)致阻抗孔上方出現(xiàn)短時間的液面塌陷。（3）淹沒水深。啟動條件下的最小淹沒深度為3～4 m，較小的淹沒水深也是形成吸氣渦旋的重要原因。流量變化大，較強的速度環(huán)量，淹沒水深較小，都是形成吸氣漩渦的決定性因素。T 型調(diào)壓室方案在水位下降時出現(xiàn)強烈的吸氣漩渦并延伸至升管，而π型調(diào)壓室則未出現(xiàn)吸氣旋渦。因此從瞬態(tài)水力特性角度來看，π 型調(diào)壓室方案優(yōu)于T 型調(diào)壓室方案。

圖22 兩種方案過流面積比

5 結(jié)論

本文采用一維和三維數(shù)值模擬對T型調(diào)壓室和π型調(diào)壓室的過渡過程進行了對比分析，并對三維計算的調(diào)壓室局部流態(tài)進行了分析，結(jié)論如下：（1）甩負荷工況和一臺機滿載另一臺機增負荷的過渡過程中，一維和三維計算的調(diào)壓室涌浪波動曲線波動趨勢相同，均隨時間周期性衰減，涌浪極值的差值均在1.5 m以內(nèi)，驗證了CFD計算的可靠性。（2）T型截面調(diào)壓室在啟動工況水位下降時調(diào)壓室的淹沒水深較小，長上室中的水流補充阻抗孔上方的水流時呈順時針或逆時針流入，導(dǎo)致吸氣漩渦的產(chǎn)生；而π型調(diào)壓室由于長上室內(nèi)水流的即時補充，不具備產(chǎn)生吸氣旋渦的條件，因此從水力學(xué)條件上，π型調(diào)壓室方案優(yōu)于T型調(diào)壓室方案。