柳海濤，徐建榮，孫雙科，彭 育

（1.流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，中國(guó)水利水電科學(xué)研究院，北京 100038）（2.中國(guó)電建集團(tuán) 華東勘測(cè)設(shè)計(jì)研究院有限公司，浙江 杭州 311122）

1 研究背景

大型水電站泄洪霧化對(duì)于工程安全與周邊環(huán)境有重大影響[1-4]，如何評(píng)估與防范霧化危害，是水利水電工程設(shè)計(jì)中的重要研究課題之一。泄洪霧化的本質(zhì)是水團(tuán)、水滴以及水霧在大氣中的運(yùn)動(dòng)，其中水舌入水激濺占據(jù)主導(dǎo)地位，由此構(gòu)成了霧化降雨的主體。國(guó)內(nèi)目前主要采用隨機(jī)濺水?dāng)?shù)學(xué)模型對(duì)霧化降雨進(jìn)行分析[5-8]，該類模型通過(guò)描述每個(gè)水滴在空氣中的運(yùn)動(dòng)過(guò)程，同時(shí)考慮空氣阻力、浮力、重力等條件，求得霧化濺水區(qū)內(nèi)下墊面降雨強(qiáng)度分布。上述模型在計(jì)算中主要考慮泄洪水力學(xué)條件與地形條件的影響，對(duì)于海拔高程與氣象條件的影響尚未考慮。

目前我國(guó)水電工程建設(shè)與運(yùn)行已從低海拔地區(qū)發(fā)展到2000 m以上高海拔地區(qū)，泄洪霧化與氣象環(huán)境之間相互影響不容忽視。從定性的角度看，海拔高程對(duì)于霧化降雨的影響，主要體現(xiàn)在三個(gè)方面：（1）不同海拔地區(qū)大氣壓強(qiáng)不同，從而影響水流摻氣、擴(kuò)散過(guò)程，引起雨霧總量發(fā)生改變；（2）不同海拔地區(qū)空氣密度與黏滯系數(shù)不同，使得霧化水滴在運(yùn)動(dòng)中所受阻力、浮力等作用發(fā)生改變，引起霧化降雨分布發(fā)生改變；（3）不同地域氣壓與空氣密度的差異，使得的泄洪雨霧輸運(yùn)與沉降過(guò)程發(fā)生改變，最終影響外圍雨霧濃度分布規(guī)律。

針對(duì)上述問(wèn)題，需通過(guò)不同海拔地區(qū)對(duì)比試驗(yàn)、原型觀測(cè)、理論分析、數(shù)學(xué)模型等手段進(jìn)行系統(tǒng)研究。鑒于問(wèn)題的復(fù)雜性，本文僅針對(duì)海拔高程對(duì)于霧化降雨分布的影響，開展探索性研究。通過(guò)理論分析，在前人研究的基礎(chǔ)上，建立考慮當(dāng)?shù)睾０胃叱逃绊懙男购殪F化降雨隨機(jī)濺水?dāng)?shù)學(xué)模型，結(jié)合實(shí)際工程霧化原型觀測(cè)資料，進(jìn)行數(shù)學(xué)模型計(jì)算驗(yàn)證。然后，針對(duì)不同海拔高程條件下，泄洪霧化降雨分布的變化規(guī)律進(jìn)行敏感性分析，為今后全面考慮氣象條件對(duì)于泄洪霧化過(guò)程的影響，奠定理論研究基礎(chǔ)。

2 隨機(jī)濺水?dāng)?shù)學(xué)模型的基本理論

將泄洪水舌入水噴射過(guò)程視為一種隨機(jī)噴射現(xiàn)象，采用拉格朗日方法，對(duì)噴射水滴在空氣中的運(yùn)動(dòng)進(jìn)行跟蹤，然后通過(guò)統(tǒng)計(jì)方法求得濺水區(qū)內(nèi)下墊面上的降雨強(qiáng)度分布。

2.1 水滴運(yùn)動(dòng)的微分方程水滴在運(yùn)動(dòng)過(guò)程中，受到重力、浮力和空氣阻力的共同作用，由此建立水滴運(yùn)動(dòng)的力學(xué)微分方程：

式中：u、v、w為水滴在空間坐標(biāo)[x、y、z]處的運(yùn)動(dòng)速度，m/s；uf、vf、wf分別為水滴鄰近風(fēng)速，m/s；Cf為阻力系數(shù)；d為水滴粒徑，m；ρa(bǔ)為空氣密度，kg/m3；ρw為水的密度，kg/m3；g為重力加速度，m/s2。求解上述方程組需給定水滴初始噴射條件，然后采用4階Runge-Kutta法[9]進(jìn)行數(shù)值求解。

2.2 水滴隨機(jī)噴射條件

（1）水滴直徑d采用概率密度函數(shù)分布[8]：

（2）水滴初始噴射速度u采用概率密度函數(shù)分布：

（3）水滴出射角θ采用概率密度函數(shù)分布：

（4）水滴出射偏轉(zhuǎn)角φ采用概率密度函數(shù)分布：

式中：μ為偏轉(zhuǎn)角眾值，依據(jù)水舌平面偏轉(zhuǎn)角度取值；σ為偏轉(zhuǎn)角的均方差，σ取值在20°～30°之間，本文取22.5°。

（5）水滴噴射速度與角度的眾值采用下式表示[10]：

（6）水滴噴射顆粒流量n：

水舌入水激濺主要發(fā)生于入水前緣，其噴射厚度h可以表示為[11]：

式中：η為系數(shù)，可取25；C為含水濃度，可根據(jù)水舌入水時(shí)斷面形態(tài)判斷，一般地，對(duì)于充分發(fā)展的摻氣水舌，一般在0.03左右；v為水的運(yùn)動(dòng)黏滯系數(shù)，m2/s；R為水力半徑，m；u*為摩阻流速，m/s，其表達(dá)式為其中τ為空氣阻力

水舌入水噴射的總流量可進(jìn)一步表示為：

式中：ks為噴濺系數(shù)，ks=0.01～0.03，本文取0.01；l為水舌入水前緣總長(zhǎng)度，m。

一旦噴射總流量確定，則對(duì)應(yīng)的水滴顆粒流量可以表示為：

式中dm為整個(gè)噴射過(guò)程的水滴平均粒徑，該粒徑值未知且不同于粒徑眾值。為此，本文模型首先計(jì)算粒徑眾值對(duì)應(yīng)的顆粒流量然后，由式（2）中概率密度函數(shù)，隨機(jī)生成個(gè)噴射水滴其對(duì)應(yīng)噴射總流量為最后，式（8）中總流量對(duì)應(yīng)的顆粒流量修正為

通過(guò)上式（2）—式（9），可以隨機(jī)生成一組噴射水滴，作為初始條件，運(yùn)用式（1）計(jì)算每個(gè)水滴的運(yùn)動(dòng)過(guò)程。

2.3 海拔高程對(duì)于水滴運(yùn)動(dòng)的影響海拔高程與氣溫的變化，會(huì)引起空氣密度、氣壓、黏滯系數(shù)、阻力系數(shù)的變化，最終引起泄洪霧化降雨分布發(fā)生變化。

（1）水滴運(yùn)動(dòng)阻力系數(shù) Cf采用下式表示[12]：

（2）空氣運(yùn)動(dòng)黏滯系數(shù) va采用下式表示[13]：

式中Tr為蘭氏溫度，oR，可表示為Tr=1.8Ta+491.67，Ta為空氣溫度，℃。

（3）空氣密度ρa(bǔ)與海拔高程、氣溫的關(guān)系可用下式表示：

式中E為海拔高程，m。

（4）當(dāng)?shù)貧鈮簆與海拔高程E、氣溫Ta的關(guān)系采用下式表示：

式中：p0為海平面標(biāo)準(zhǔn)氣壓，取101 325Pa；cp為空氣定壓比熱容，取1005 J/（kg·K）；M為干燥空氣的摩爾質(zhì)量，取0.02896kg/mol；R0為普適氣體常數(shù)，取8.315 J/（mol·K）。

將上式（10）—式（12）與式（1）聯(lián)立求解，可以考慮工程當(dāng)?shù)睾０闻c氣溫對(duì)于霧雨運(yùn)動(dòng)過(guò)程的影響。由于水舌風(fēng)與霧化降雨相互伴生，不可分割，以往的研究對(duì)于水舌風(fēng)場(chǎng)多采用正態(tài)分布假定或者均一化假定處理。對(duì)此，本文通過(guò)將三維水舌風(fēng)場(chǎng)轉(zhuǎn)換為結(jié)構(gòu)化數(shù)據(jù)，然后實(shí)時(shí)求解水滴瞬時(shí)鄰近風(fēng)速，用于閉合式（1）中水滴運(yùn)動(dòng)方程，從而較全面考慮河谷中復(fù)雜風(fēng)場(chǎng)對(duì)于雨霧運(yùn)動(dòng)過(guò)程的影響。

2.4 濺水模型的數(shù)值計(jì)算方法

（1）噴射時(shí)間的離散。在濺水過(guò)程中，水滴分布于整個(gè)空間，只有當(dāng)水滴到達(dá)地面時(shí)才能形成降水。為此，在本文模型中，首先，根據(jù)水滴自由拋射后在空中的最長(zhǎng)停留時(shí)間t，確定濺水噴射總歷時(shí)T，一般應(yīng)滿足T≥2t；然后，將計(jì)算時(shí)間劃分為m個(gè)時(shí)間步長(zhǎng)dt，假設(shè)在噴射歷時(shí)T內(nèi)的水滴總量為N=Tn，n為水滴噴射顆粒流量，則每個(gè)步長(zhǎng)內(nèi)的水滴總量為Ni=N/m；第三，運(yùn)用4階龍格-庫(kù)塔法求解每個(gè)時(shí)間步長(zhǎng)內(nèi)水滴的運(yùn)動(dòng)方程，其中對(duì)于第i個(gè)步長(zhǎng)內(nèi)的Ni個(gè)水滴，其最長(zhǎng)飛行歷時(shí)為

（2）降雨強(qiáng)度的計(jì)算。降雨強(qiáng)度定義為一定時(shí)間內(nèi)穿過(guò)單位面積平面的總水量除以降雨歷時(shí)。為此，本文模型中，首先將噴射區(qū)域內(nèi)地面離散為小尺度的集雨網(wǎng)格，然后在每一個(gè)時(shí)間步長(zhǎng)dt內(nèi)，判斷Ni個(gè)水滴的垂向位置，若其在n-1時(shí)刻位于地面高程以上，而n時(shí)刻位于該平面以下，則表明水滴已降落地面，此時(shí)根據(jù)其平面坐標(biāo)，將其水量計(jì)入相應(yīng)的集雨網(wǎng)格，即有Volume=Volume0+πd36，同時(shí)該水滴的飛行終止。由于采用隨機(jī)函數(shù)模擬水滴出射條件，模型需要重復(fù)M次噴射過(guò)程，當(dāng)所有M次噴射過(guò)程計(jì)算完成后，統(tǒng)計(jì)平面上每個(gè)網(wǎng)格內(nèi)積累的水滴總體積，再除以噴濺歷時(shí)T、網(wǎng)格垂直投影面積A與重復(fù)次數(shù)M，即可得到該網(wǎng)格內(nèi)時(shí)均降雨強(qiáng)度分布P=Volume/(M AT)。

（3）水滴運(yùn)動(dòng)過(guò)程中風(fēng)速與地形的影響。模型求解水滴運(yùn)動(dòng)方程中，需要獲取水滴鄰近風(fēng)速與地面高程，前者用于求解水滴所受風(fēng)場(chǎng)加速度，后者用于判斷水滴是否降落地面，從而反映兩者對(duì)于霧化降雨分布的影響。為此，本文模型中，首先采用相對(duì)成熟的Fluent軟件求解三維河谷風(fēng)場(chǎng)，然后通過(guò)相應(yīng)的轉(zhuǎn)換程序，將非結(jié)構(gòu)化的三維風(fēng)場(chǎng)與地形數(shù)據(jù)轉(zhuǎn)換為結(jié)構(gòu)化的風(fēng)場(chǎng)數(shù)據(jù)Wind(i , j,k)及地形數(shù)據(jù)[T o po(i , j)]，其中i、j、k為三維結(jié)構(gòu)化坐標(biāo)，為節(jié)省空間，上述數(shù)據(jù)以二進(jìn)制文件存儲(chǔ)；隨機(jī)濺水模型運(yùn)行初始，先讀取上述風(fēng)場(chǎng)與地形數(shù)據(jù)，并按照其結(jié)構(gòu)化坐標(biāo)存儲(chǔ)在相應(yīng)的數(shù)組單元中。對(duì)于空間任意水滴，若要獲取其鄰近風(fēng)速，可通過(guò)對(duì)其瞬時(shí)坐標(biāo)（x,y,z）數(shù)值取整，找到其相鄰的8個(gè)結(jié)構(gòu)化坐標(biāo)[Ii,Ji,Ki]，i=1～8。根據(jù)其結(jié)構(gòu)化坐標(biāo)，直接從對(duì)應(yīng)數(shù)組單元讀取相應(yīng)的風(fēng)速數(shù)據(jù)，然后插值得到水滴的鄰近風(fēng)速Ufki為空間插值函數(shù)。同樣地，采用二維空間插值方法可以求得水滴的鄰近地面高程。

（4）數(shù)值振蕩的抑制。從水滴運(yùn)動(dòng)微分方程得知，式（1）中右端項(xiàng)為水滴所受的運(yùn)動(dòng)加速度，其阻力項(xiàng)的量值與水滴粒徑成反比，由于采用隨機(jī)函數(shù)模擬，水滴噴射粒徑為一隨機(jī)變量，當(dāng)某一水滴粒徑趨近于0時(shí)，該水滴所受阻力加速度則趨于無(wú)窮大，水滴運(yùn)動(dòng)則產(chǎn)生數(shù)值振蕩。從理論上分析可知，單位時(shí)間步長(zhǎng)內(nèi)，水滴所受的最大加速度可表示為UT為水滴終極速度，U0為水滴初始速度，dt為時(shí)間步長(zhǎng)。在水平方向上，水滴終極速度UT=U0，即為當(dāng)?shù)仫L(fēng)速；在垂線方向上，水滴終極速度為即為重力、浮力、阻力達(dá)到平衡時(shí)速度。當(dāng)水滴初始時(shí)刻所受的運(yùn)動(dòng)加速度為a0≥amax，表明水滴運(yùn)動(dòng)速度在該時(shí)間步長(zhǎng)中已經(jīng)達(dá)到穩(wěn)態(tài)，則在該時(shí)間步長(zhǎng)中，令a0=amax，下一時(shí)間步長(zhǎng)中，水滴初始速度U0=UT。

3 數(shù)學(xué)模型驗(yàn)證分析

3.1 計(jì)算條件中電建昆明院科研所于2014年8月17日，針對(duì)小灣水電站泄洪洞全開泄洪過(guò)程，水舌入水形態(tài)與泄洪霧化進(jìn)行了原型觀測(cè)。觀測(cè)過(guò)程中，壩上水位1236.32 m，實(shí)測(cè)泄洪流量3500 m3/s，當(dāng)?shù)睾０?000 m，氣溫20℃，觀測(cè)時(shí)間50 min。泄洪洞實(shí)地泄洪霧化情況見圖1，水舌入水形成的雨霧在泄洪風(fēng)場(chǎng)作用下沿兩岸爬升，高度接近300 m，通過(guò)沿岸布置霧化降雨測(cè)點(diǎn)，得到了下游降雨強(qiáng)度的實(shí)測(cè)點(diǎn)據(jù)。本次觀測(cè)工況中泄洪水舌入水條件如表1，水舌入水位置及下游河谷地形見圖2。

圖1 小灣電站泄洪洞全開霧化情況

3.2 水舌風(fēng)場(chǎng)計(jì)算結(jié)果泄洪水舌入水過(guò)程中，在河谷中形成局部風(fēng)場(chǎng)，對(duì)霧化降雨分布范圍產(chǎn)生影響。根據(jù)當(dāng)?shù)睾０闻c氣溫條件，由式（11）—式（13）求得當(dāng)?shù)卮髿怵禂?shù)1.6791×10-5m2/s，大氣密度1.0805 kg/m3，大氣壓強(qiáng)89 999.2 Pa。然后，基于Fluent軟件，求解泄洪洞水舌下游河谷風(fēng)場(chǎng)。計(jì)算中采用了RNGk-ε模型封閉連續(xù)與動(dòng)量方程[14]，水舌入水?dāng)嗝娌捎盟俣冗吔?，周圍河谷地形按照?biāo)準(zhǔn)壁面函數(shù)進(jìn)行處理，頂部為壓力出口邊界。圖3為計(jì)算得到的河谷地面附近風(fēng)速等值線圖。由圖可知，河谷泄洪風(fēng)場(chǎng)受到入水條件與地形的影響，在水舌入水附近風(fēng)速可達(dá)50 m/s量級(jí)，然后向外圍逐步衰減，縱向影響范圍可達(dá)下游1500 m，兩岸爬升高度超過(guò)300 m，其復(fù)雜形態(tài)對(duì)于霧化降雨分布將會(huì)產(chǎn)生影響。

3.3 降雨強(qiáng)度分布計(jì)算結(jié)果根據(jù)泄洪水舌入水條件、泄洪風(fēng)場(chǎng)以及當(dāng)?shù)貧庀髤?shù)，采用隨機(jī)濺水?dāng)?shù)學(xué)模型進(jìn)一步計(jì)算河谷中霧化降雨分布形態(tài)。圖4為泄洪洞下游霧化降雨強(qiáng)度分布形態(tài)，由圖可知，霧化降雨分布受到風(fēng)場(chǎng)與地形的影響，霧化雨區(qū)在兩岸爬升高度約250 m，縱向分布范圍受到下游右岸河谷地形的阻擋，邊界位于泄洪洞下游880 m。圖5為泄洪洞下游霧化降雨強(qiáng)度等值線圖，為便于比較，將原型觀測(cè)點(diǎn)據(jù)一并匯出（圖中方框數(shù)據(jù)），結(jié)果表明兩者在分布規(guī)律上甚為吻合。

表1 泄洪洞水舌入水條件

圖2 小灣泄洪洞水舌入水位置與河谷地形等高線

圖3 小灣泄洪洞下游地面附近風(fēng)速分布

圖4 泄洪洞霧化降雨分布形態(tài)

圖5 泄洪洞霧化降雨強(qiáng)度分布

表2 不同海拔地區(qū)氣象因子計(jì)算結(jié)果

圖6 海拔高程50m條件下泄洪洞下游霧化降雨分布

3.4 海拔高程對(duì)于霧化降雨影響的敏感性分析小灣電站當(dāng)?shù)睾０胃叱碳s1000 m，為分析海拔高程的影響，本文在泄洪條件、河谷地形、當(dāng)?shù)貧鉁鼐蛔兊那疤嵯?，研究不同海拔高?0 m、2000 m、3000 m情況下，泄洪霧化降雨分布的變化規(guī)律，上述海拔高程對(duì)應(yīng)的氣象參數(shù)見表2。圖6—圖8為不同海拔高程條件下，泄洪霧化降雨分布計(jì)算結(jié)果，結(jié)果分析表明：（1）隨著海拔高程的增加，在相同的泄洪條件下，泄洪洞下游霧化降雨區(qū)分布范圍會(huì)有所增加。當(dāng)海拔高程50 m時(shí)，兩岸雨區(qū)爬升高度約190～220 m，縱向邊界位于泄洪洞下游約850 m；而當(dāng)?shù)睾０胃叱?000 m時(shí)，兩岸雨區(qū)爬升高度可達(dá)230～280 m，縱向邊界位于泄洪洞下游950 m，若非受到右岸地形的阻擋，雨區(qū)縱向變化范圍會(huì)更大。（2）隨著海拔高程的增加，霧化降雨強(qiáng)度分布規(guī)律發(fā)生坦化。表3為霧化降雨強(qiáng)度沿河谷縱向的分布數(shù)據(jù)，由表可知，在大于400 mm/h的等值線區(qū)域，分布范圍有所減小，而在小于400 mm/h的區(qū)域，分布范圍有所增大。（3）海拔高程除了對(duì)霧化降雨分布產(chǎn)生影響以外，還會(huì)對(duì)泄洪雨霧總量以及雨霧輸運(yùn)過(guò)程產(chǎn)生影響，此外氣溫條件也是重要的影響因素，對(duì)此需要進(jìn)一步開展系統(tǒng)研究。

圖7 海拔高程2000m條件下泄洪洞下游霧化降雨分布

圖8 海拔高程3000m條件下泄洪洞下游霧化降雨分布

表3 泄洪霧化降雨強(qiáng)度沿河谷縱向的分布數(shù)據(jù)

4 結(jié)語(yǔ)

本文通過(guò)理論分析，建立了考慮海拔高程與氣溫條件影響的泄洪霧化隨機(jī)濺水?dāng)?shù)學(xué)模型，結(jié)合小灣電站泄洪洞霧化原型觀測(cè)數(shù)據(jù)進(jìn)行了驗(yàn)證，兩者數(shù)據(jù)吻合良好。在此基礎(chǔ)上，針對(duì)不同海拔高程條件下，泄洪洞下游霧化降雨分布變化規(guī)律進(jìn)行敏感性分析。在泄洪條件、地形、氣溫不變前提下，隨著海拔高程的增加，泄洪洞下游霧化降雨區(qū)分布范圍會(huì)有所增加，同時(shí)霧化降雨強(qiáng)度分布規(guī)律發(fā)生坦化，證明海拔高程對(duì)于霧化降雨強(qiáng)度分布的確存在一定影響。本文研究?jī)H針對(duì)泄洪霧化降雨過(guò)程，下一步擬對(duì)水舌激濺總源量以及雨霧輸運(yùn)過(guò)程，繼續(xù)開展海拔高程與氣溫條件的影響分析。