李 昂

（重慶市水利電力建筑勘測設(shè)計研究院，重慶 401120）

0 引言

水利工程中，為了滿足檢修、國防（如戰(zhàn)爭）或者緊急情況（例如地震）等要求，往往需要在水庫較低高程位置修建放空建筑物，如在重力壩或者拱壩壩身設(shè)置的放空底孔，在河岸開鑿的放水隧洞等。另外，水池放空、船閘充水及泄水等均需要計算放水及充水所需時間。

對于放空時間的確定，常通過如下方法進行估算（見圖1）。

圖1 放空過程計算示意圖

為方便計算，假定上游無來流，孔口做自由出流，則dt 時段內(nèi)從孔口流出的水流體積與同一時段內(nèi)容器里面水體體積的減小量相等，即：

式中：Q 為孔口體積流率，A 為容器橫截面面積，H 為容器內(nèi)水深，μ 為孔口流量系數(shù)，b 為孔口寬度，e 為孔口高度，g 為重力加速度，取9.81 m/s2。

若孔口水頭從H1變化至H2，對（2）式進行積分得到放空所需時間為：

對于上述計算方法，有兩點不足：其一是計算過程中將孔口流量系數(shù)μ 假定為常值，這與實際情況不符，事實上，孔口流量系數(shù)不僅與孔口型式有關(guān)，還與作用水頭有關(guān)。通常情況下，隨著水頭的降低，流量系數(shù)也會逐漸減小，這就造成了計算得出的放空時間偏小。其二，式（1）中始終將孔口出流形式做閘孔出流考慮，其流量計算公式采用有壓流計算公式，然而當(dāng)水頭降至一定程度時，孔口出流形式將由閘孔出流變成堰流（即無壓流），此時上述計算過程已經(jīng)不再適用，需要采用無壓流計算公式對泄流量進行計算。根據(jù)吳持恭的研究成果[1]，當(dāng)閘底坎為平頂堰的情況，下列數(shù)值可作為判別流態(tài)的大致界線：當(dāng)為閘孔出流；當(dāng)為堰流。

當(dāng)孔口出流形式為堰流時，其泄流時間計算過程如下（式中m 為孔口為堰流時的流量系數(shù)）：

若孔口水頭從H3變化至H4，對上式進行積分得到所需時間為：

對于本文所述的放空問題，由于出流孔口長度為零，所以當(dāng)水頭變化至界限水位時，孔口出流流態(tài)將迅速由閘孔出流變?yōu)檠吡?。但是實際工程問題中的放空建筑物往往是具有一定長度的有壓孔口，此時將出現(xiàn)兩個界限水位，即上界限水位與下界限水位，此二水位之間有壓管內(nèi)將表現(xiàn)為明滿流交替狀態(tài)，常規(guī)的計算方法不再適用。

另外，隨著水頭的降低，一些不利的水力現(xiàn)象也將可能出現(xiàn)，譬如前文提到的明滿交替流將對有壓管道的結(jié)構(gòu)穩(wěn)定造成不利影響，低水位情況常常出現(xiàn)的立軸漩渦也是一個需要重點關(guān)注的現(xiàn)象。這些現(xiàn)象在常規(guī)的計算過程中都是無法預(yù)測的，隨著計算機技術(shù)以及數(shù)值計算理論的飛速發(fā)展，數(shù)值模擬的方法為該問題的解答提供了一條較為經(jīng)濟可靠的解決方案。

對于非恒定流狀態(tài)下水庫放空問題的數(shù)值模擬計算目前還未見相關(guān)文獻，其計算精度是否能夠滿足要求也不得而知。本文分別通過二維和三維的計算模型對棱柱形容器自由出流進行全時域過程的模擬，從而對數(shù)值模擬在放空問題上的應(yīng)用提供一些技術(shù)支撐。

1 數(shù)學(xué)模型

本文分別采用二維和三維RNG k-ε 雙方程模型對紊流進行模擬，采用控制體積法來對偏微分方程進行離散，采用SIMPLER 法對壓力和速度進行耦合，邊墻采用無滑移邊界條件[2]，采用VOF 法[3]對自由面進行捕捉。

本文采用的RNG k-ε 紊流模型，其連續(xù)方程，動量方程和k，ε 方程可分別表示如下：

連續(xù)方程：

動量方程：

k 方程：

ε 方程：

式中：ρ 和μ 別表示為體積分數(shù)平均的密度和分子黏性系數(shù)，p表示修正壓力，μt表示紊動黏性系數(shù)，它可以根據(jù)紊動能k 及紊動能耗散率ε 求解得出。

以上各個表達式中，i=1，2，3（對二維計算 i=1，2），即{xi=x,y,z},{ui=u,v,w}（對二維計算 {xi=x,y},{ui=u,v}）；j 表示求和下標(biāo)，方 程 中 通 用 模 型 常 數(shù)[4]η0=4.38,β=0.012,Cμ=0.0845,C2ε=1.68,σk=0.7179 和 σε=0.7179。

某水庫總庫容為1200 萬m3，死庫容110 萬m3，興利庫容1090 萬m3。放水洞尺寸為2 m×2 m（寬×高），水庫正常蓄水位時孔口水頭為30 m，其建模概化如下：

二維建模容器長 20 m（X 向），高 30 m（Y 向），孔口位于容器右下角，孔口高度2 m，網(wǎng)格間距0.1 m。三維建模容器長20 m（X 向），寬 20 m（Y 向），高 30 m（Z 向），孔口寬 2 m，高 2 m；網(wǎng)格整體間距為0.4 m，在孔口周圍（X=0～5 m;Y=-5 m～5 m；Z=0～8 m）進行局部加密，加密后間距為0.2 m。設(shè)置模型頂部為壓力進口，孔口為壓力出口。圖2 為三維模型布置圖。

圖2 三維計算模型布置圖

2 計算結(jié)果分析

圖3 為二維數(shù)值與三維數(shù)值計算的水位隨時間變化曲線及其與理論計算（按照式（3）進行計算，流量系數(shù)其中ζ 為局部損失系數(shù)，取0.5，則孔口流量系數(shù)μ=0.816）的對比?？梢钥吹剑陂_始時段，二維數(shù)值計算與三維數(shù)值計算均與理論計算具有較好的吻合度，這也驗證了數(shù)值計算結(jié)果是可靠的、可信的；隨著水位的逐漸降低，數(shù)值解逐漸偏離理論計算結(jié)果，但是二者總體變化趨勢一致，分析認為造成此種情況的原因是理論計算采用的流量系數(shù)是一個常值，而實質(zhì)上隨著水位的降低，孔口的流量系數(shù)也逐漸減小，故而同一時刻的水位數(shù)值計算結(jié)果比理論值偏大，二維數(shù)值計算與三維數(shù)值計算均具有類似的結(jié)果。

另外也可以觀察到，三維計算比二維計算具有更高的精度以及吻合度，根據(jù)前文對于孔口出流流態(tài)的判定條件可以計算得到界限水位H=e/0.65=2/0.65=3.08 m，從圖3（b）中可以看到三維數(shù)值計算結(jié)果近乎在界限水位附近才與理論計算結(jié)果逐漸出現(xiàn)相對較大的差異，而二維數(shù)值計算（圖3（a））則在比較高的水位條件下就開始出現(xiàn)偏離。使用均方根誤差（RMSE）來定性地分析數(shù)值計算結(jié)果的精度，其計算式如下：

式中：n 表示水深計算點個數(shù)，Htheoretical表示理論計算結(jié)果，Hnumerical表示數(shù)值計算結(jié)果，顯然RMSE 越小，則精度越高。計算得到二維數(shù)值計算RMSE 為1.17，而三維數(shù)值計算RMSE 為0.65，故而對于精度要求較高的放空情況，數(shù)值計算模型建議采用三維計算。

圖3 數(shù)值解與理論解對比

圖4 為孔口流量系數(shù)隨水位變化關(guān)系曲線，該圖直觀地顯示了孔口流量系數(shù)隨著庫水位的降低而減小，且水位降低越多，該趨勢越明顯，界限水位以上流量系數(shù)隨著水頭減小近似呈線性減小趨勢；當(dāng)水位降至界限水位后，孔口流量系數(shù)迅速減小，這是由于孔口流態(tài)變成堰流所致；另外，三維數(shù)值計算的流量系數(shù)比二維數(shù)值計算結(jié)果稍大，但是二者差異相對不大且變化規(guī)律一致。

圖5 為孔口單寬流量隨水位變化關(guān)系曲線，隨著庫水位的降低，單寬流量逐漸減小，當(dāng)水位降至界限水位以下后變化更為明顯，這同樣是因為孔口流態(tài)的轉(zhuǎn)換所致；同樣地，三維數(shù)值計算的單寬流量比二維數(shù)值計算結(jié)果稍大。

圖4 流量系數(shù)隨水位變化關(guān)系

圖5 單寬流量隨水位變化關(guān)系

圖6 為不同典型水位條件下（H=20 m、H=3.08 m、H=1.5 m，分別表示閘孔出流、過渡流及堰流）容器內(nèi)水流體積分數(shù)及壓強分布情況，其中左側(cè)為二維數(shù)值計算結(jié)果，右側(cè)為三維數(shù)值計算結(jié)果（取孔口中軸面，即Y=0 m 斷面）?？梢钥吹疆?dāng)水位較高（如圖7（a）所示）時，二維數(shù)值計算與三維數(shù)值計算結(jié)果幾乎沒有差異；隨著水位的降低，二者的差異逐漸趨于明顯，二維數(shù)值計算孔口上游影響區(qū)域較三維計算大，這主要是因為三維計算橫向?qū)挾龋?0 m）遠大于孔口寬度（2 m），故而孔口上游附近區(qū)域的流速迅速減小，從而造成水面降低沒有二維計算那么明顯。

為了簡化計算，建模模型完全對稱且孔口沒有長度，故而在全時域過程中不會出現(xiàn)明滿交替流，同時也沒有觀察到立軸漩渦的產(chǎn)生。但是對于實際工程問題，出流孔口往往具有一定長度，在水位逐漸下降的過程中，不可避免地會遇到明滿交替流以及立軸漩渦的產(chǎn)生，此種情況通過常規(guī)的物理模型試驗往往很難實現(xiàn)全時域過程的模擬，此時數(shù)值計算的優(yōu)勢將顯現(xiàn)無遺。

另外可以觀察到，即使是在非常低的水位條件下，容器最左側(cè)壁面（X=-20 m）附近水面也基本保持水平，且各水位條件下該區(qū)域壓強近似呈靜水壓強分布，即該區(qū)域行進流速幾乎為零。事實上，計算過程中選取的容器內(nèi)水深就處于該位置，這表明模型的建模長度滿足要求，容器水位的讀?。ㄓ捎谒畾饨唤缑婢哂幸欢ǖ母叨龋嬎阒羞x取水相體積分數(shù)為0.5 的高程位置作為容器內(nèi)水深）不需要考慮行進流速的影響。

圖6 2D 與3D 數(shù)值計算水力參數(shù)對比

3 結(jié)果與討論

本文針對實際工程中可能遇到的放空問題，將問題進行簡化后分別進行了二維和三維數(shù)值計算，結(jié)果表明：

1）運用數(shù)值計算的方法來模擬非恒定條件下的放空問題是可行的，二維與三維數(shù)值計算結(jié)果均與理論分析結(jié)果具有相對良好的吻合度，但是三維數(shù)值計算結(jié)果的計算精度要優(yōu)于二維計算，對于精度要求較高的實際問題，建議采用三維數(shù)值計算。

2）當(dāng)水位降至界限水位以下后，數(shù)值計算結(jié)果與理論計算結(jié)果出現(xiàn)較大差異，這主要是因為孔口出流類型由閘孔出流向堰流的轉(zhuǎn)換所致，這與理論分析的結(jié)論一致，數(shù)值計算較好地驗證了這一現(xiàn)象。

3）由于三維數(shù)值計算考慮了容器寬度的影響，孔口附近行進流速比二維數(shù)值計算結(jié)果偏小，其速度和壓強受影響范圍也較小，這也是三維數(shù)值計算精度更高的主要原因。