泄洪霧化水滴分檔隨機噴濺數(shù)學模型及其驗證

張 華，何貴成

(華北電力大學水利與水電工程學院，北京 102206)

水利工程在挑流泄洪時，產(chǎn)生的高速水舌與空氣摻混[1]，兩股水舌碰撞，以及水舌和下游水面的撞擊，都會形成大小不一的霧源。霧源處噴濺的水滴，在壩下游沉降，產(chǎn)生強度遠大于自然降雨的泄洪霧化降雨[2-3]。嚴重的泄洪霧化會影響電站電氣設備的正常運行、沖蝕地表，影響岸坡穩(wěn)定和廠區(qū)交通等[4-5]。研究泄洪霧化的產(chǎn)生機理、霧源強度、擴散規(guī)律、以及降雨強度，對采取合理措施減輕泄洪霧化的危害具有重要價值[6]。

泄洪霧化降雨范圍和大小的計算，常用原型觀測[7]、理論分析[8]、模型試驗[9]和數(shù)值模擬[10-11]等方法。梁在潮[12]詳細描述了泄洪霧化現(xiàn)象，推導出了水舌摻氣量、濺水影響區(qū)域和霧流擴散的理論計算公式。劉宣烈等[1]采用立體攝影和電阻式摻氣儀兩種測量方法，試驗得到了水舌斷面摻氣濃度、沿程變化與參數(shù)的關系。張華等[13]在觀測與試驗數(shù)據(jù)的基礎上，提出了水滴隨機噴濺數(shù)學模型，應用較為廣泛[14-16]。水滴隨機噴濺數(shù)學模型在初始條件中采用多個隨機變量以模擬水滴在霧源噴出時的隨機特性，隨機變量包括噴濺水滴的直徑、速度和偏移角等。

在水滴隨機噴濺數(shù)學模型基礎上，學者進行了模型的適用性驗證，以及改進工作。如在雙江口水電站數(shù)學模型噴濺試驗[17]的數(shù)值計算中，水面以上的下墊面為自然地形，驗證了水滴隨機噴濺數(shù)學模型在復雜地形情況下的適應能力。兩河口水電站的泄洪霧化計算中[18]，引入水舌風的作用，在數(shù)學模型中增加考慮了水舌風對噴濺水滴運動的影響。

水滴隨機噴濺數(shù)學模型，通常采用蒙特卡洛方法[19]進行求解，即生成大量符合假設隨機變量分布的隨機水滴，最后統(tǒng)計水滴停止運動后的分布情況，由此獲得霧化降雨大小。由于霧化降雨大小由隨機水滴統(tǒng)計產(chǎn)生，因此受到水滴樣本數(shù)量的影響較大。為得到較準確的霧化降雨數(shù)值，水滴隨機噴濺數(shù)學模型常需要生成幾百萬的水滴作為樣本，計算量大、效率低。

云物理學理論[20]為計算云滴譜隨時間和高度的演變情況[21]，采用的方法是將云滴按照大小劃分為若干檔[22]，然后隨著時間步的增加，云滴自然凝結(jié)生長，云滴大小越出原檔，各檔云滴數(shù)量也隨之增減，從而計算出云滴譜的演變情況。在燃油噴霧數(shù)值模擬[23]中，離散液滴數(shù)學模型[24]并不會計算所有油滴的運動，而是選取計算若干具有代表性的樣本油滴的運動軌跡[25]，提高了計算效率。在模擬流化床[26]的顆粒流動算法中，如OpenFOAM的DPMFoam求解器[27]，采用了顆粒微團，即多個顆粒的集合，并假定集合中的顆粒具有相同的尺寸和速度。上述方法的目的都在于減少計算霧滴、油滴和顆粒的數(shù)量，以降低對內(nèi)存和計算能力的需求。

因此，本文基于云滴分檔、代表性油滴和顆粒微團等降低計算量的思想，提出分檔水滴的概念，建立挑流泄洪霧化的水滴分檔隨機噴濺數(shù)學模型，并在假定條件下對模型進行對比驗證。

1 水滴運動微分方程及其解析解

1.1 水滴運動微分方程

如圖1所示，假設在直角坐標系Oxyz的O點處，拋出N個水滴，初始速度為v，速度與Oxy面的夾角為出射角β，速度與Oxz面的夾角為偏移角φ。從原點處拋出的水滴做斜拋運動，當水滴接觸Oxy面時，則停止下來。在某個時間間隔內(nèi)，通過統(tǒng)計水滴在水平面上的質(zhì)量分布情況，即可求得霧化降雨的數(shù)值。

圖1 水滴噴濺示意圖

考慮水滴在運動過程中受到浮力、重力和空氣阻力的作用，則水滴的運動微分方程[13]為

式中：Cf為空氣阻力系數(shù)；d為水滴直徑；ρ為水滴密度；ρa為空氣密度；vw為水滴速度；va為環(huán)境風速(vax、vay、vaz分別是3個坐標軸上的投影)；g為重力加速度。

水滴運動微分方程的初始條件為

(2)

式中：v為水滴的初始速度；β、φ分別為水滴的出射角和偏移角。

1.2 水滴運動微分方程的解析解

式(1)是一個二階非線性微分方程組，不能得到其解析解。對其做解耦合線性化的處理，轉(zhuǎn)變?yōu)槎A線性微分方程組，得到：

(3)

式中：<|vw-va|>為水滴速度與環(huán)境風速差值的系綜平均。

由初始條件(2)，并假設環(huán)境風速va為常數(shù)，可以得到式(3)的解析解：

(4)

令z(ta)=0，可以獲得水滴接觸Oxy平面的時間ta的計算公式如式(5)所示，然后可以得到水滴停留在Oxy平面上的位置為(x(ta),y(ta))。

(5)

若假定v、β和φ為隨機變量，并服從聯(lián)合概率密度分布f(v,β,φ)，則由連續(xù)型隨機變量的函數(shù)的分布理論[28]，式(4)的多元函數(shù)所對應的概率密度分布為

fp(x(t),y(t),z(t))=f(v(x(t),y(t),z(t)),

β(x(t),y(t),z(t)),φ(x(t),y(t),z(t)))|J|

(6)

2 水滴分檔隨機噴濺數(shù)學模型

在水滴隨機噴濺數(shù)學模型[13]中，水滴的初始速度v、出射角β和偏移角φ，均為隨機變量，模型通常采用蒙特卡洛方法進行求解。

為了改進蒙特卡洛方法計算量大、效率低的缺點，基于云滴分檔[22]、代表性油滴[25]和顆粒微團[27]等學術(shù)思想，提出水滴分檔隨機噴濺數(shù)學模型。主要改進內(nèi)容為，根據(jù)水滴初始參數(shù)相空間進行分檔化處理，由此形成分檔水滴，并讓一個分檔水滴攜帶概率信息，以代表多個水滴，最終減少實際的計算量，提高計算效率。

2.1 水滴初始物理參數(shù)的相空間

將水滴的初始參數(shù)v、β和φ所構(gòu)成的物理空間，稱為水滴運動初始物理參數(shù)的相空間，記為Γ(v,β,φ)，如圖2所示。相空間中的一個相點表示單個水滴的初始狀態(tài)，整體表示水滴所有可能的初始狀態(tài)。

圖2 水滴初始物理參數(shù)的相空間示意圖

2.2 相空間分檔方法

水滴初始物理參數(shù)相空間的分檔方法如下：

a.水滴初始速度的分檔。對v做均勻分布的分檔化處理[22]，這里劃分為I檔，每檔的大小為a，則每個檔的速度vi為

vi=ai(i=1,2,…,I)

(7)

b.水滴出射角的分檔。對β也同樣做均勻分布的分檔化處理：分為J檔，每檔的大小為b，則出射角βj為

βj=bj(j=1,2,…,J)

(8)

c.水滴偏移角的分檔。對φ也做均勻分布的分檔化處理：在其定義域內(nèi)均勻分為K檔，每檔的大小為c，則偏移角φk為

φk=ck(k=1,2,…,K)

(9)

2.3 分檔水滴及其概率值計算

水滴運動初始物理參數(shù)相空間通過分檔之后，形成相格(圖2中長方形線框)。一個相格中包含了多個相點，即內(nèi)含多個水滴的初始狀態(tài)信息。用一個代表性相點，來表征整個相格中所有相點，稱這個代表性相點對應的水滴為分檔水滴。分檔水滴所對應的概率密度為

Pijk=f(vi,βj,φk)

(10)

3 微分方程和數(shù)學模型的求解

為了獲得水滴運動微分方程組的解析解，并以此為基準，評價水滴隨機噴濺數(shù)學模型和水滴分檔隨機噴濺數(shù)學模型的數(shù)值解的精確度，從而驗證水滴分檔隨機噴濺數(shù)學模型的求解精度。

在不改變微分方程組本質(zhì)的情況下，為對公式(1)進行簡化，提出如下兩個假定條件：①忽略空氣阻力和浮力作用；②隨機變量v和β取為定值，φ的概率密度函數(shù)fφ滿足正態(tài)分布：

(11)

3.1 水滴運動線性微分方程的解析解

在假定條件下，得到下列參數(shù)的數(shù)值：阻力系數(shù)Cf=0，空氣密度ρa=0 kg/m3。代入式(3)，得到水滴位置隨時間的函數(shù)關系為

(12)

令z(ta)=0，則水滴經(jīng)過一段時間ta，因接觸Oxy面而停止運動，可以得到水滴停止運動時的位置與φ的函數(shù)關系為

(13)

由于fφ為正態(tài)分布，絕大部分水滴位于(-3σ，3σ)范圍內(nèi)，因此可取φ的定義域為C=(-π/3,π/3)。

y(φ)的反函數(shù)為

(14)

由假設條件②，已知φ的概率密度函數(shù)，以及式(6)，可得到水滴在y軸上投影的概率密度函數(shù)為

為便于解析解與數(shù)學模型的數(shù)值解進行對比，對解析解進行離散化處理。將值域S均勻分為M個數(shù)量的組，組距為h，即第m組的取值范圍是[hm-1,hm)，分組形式記為集合A={[hm-1,hm)|m=1,2,…,M}。則水滴處于第m組內(nèi)的概率密度的解析解為

(16)

3.2 水滴隨機噴濺數(shù)學模型的數(shù)值解

水滴隨機噴濺數(shù)學模型采用蒙特卡洛方法求解，生成的水滴為隨機產(chǎn)生，現(xiàn)生成N個水滴，其中的偏移角使用隨機數(shù)發(fā)生器產(chǎn)生，記為φn,r(n=1,2,…,N)，它滿足概率密度函數(shù)fφ。

根據(jù)式(13)，可以計算水滴偏移角為φn,r時，相應在y軸上的運動距離為yn,r。依據(jù)分組集合A，統(tǒng)計yn,r在第m個分組范圍內(nèi)的水滴的個數(shù)，記為qm,r，計算公式為

(17)

得到水滴處于第m組的概率密度為

(18)

3.3 水滴分檔隨機噴濺數(shù)學模型的數(shù)值解

水滴分檔隨機噴濺數(shù)學模型，生成的分檔水滴，在定義域內(nèi)均勻分布，并且包含概率信息，具體的計算步驟如下：

步驟1由假定條件②，只需要對φ進行分檔處理，因此在其定義域內(nèi)，產(chǎn)生N個分檔水滴，每個分檔水滴的偏移角由式(9)計算，記為φn,g(n=1,2,…,N)。

步驟2根據(jù)概率密度函數(shù)fφ，以及式(10)，計算得到偏移角為φn,g的分檔水滴所對應攜帶的概率pn,g為

pn,g=fφ(φn,g)c

(19)

步驟3由式(13)計算出分檔水滴在偏移角為φn,g時，在y軸上的運動距離yn,g。

步驟4依據(jù)分組集合A對分檔水滴在y軸的距離進行分組，得到處于第m組范圍水滴的概率qm,g為

(20)

步驟5得到水滴處于第m組范圍的概率密度為

(21)

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

式(1)是非線性微分方程組，一般情況下，只能采用數(shù)值模型進行求解。若不進行參數(shù)簡化，僅能得到兩個數(shù)學模型的數(shù)值解。沒有解析解的準確值作為基準，無法判別兩種數(shù)值模型的計算結(jié)果的誤差大小。案例是在兩個假設條件下，對微分方程組中的參數(shù)進行簡化，獲得解析解式(15)，并以此作為后續(xù)兩個數(shù)值解計算誤差的對比基準。

在兩個假定條件下，選取v=10 m/s，g=9.81 m/s2，β=π/5，M=100。根據(jù)上述3種求解方法，可以分別得到解析解和兩個模型的數(shù)值解，并對計算結(jié)果進行對比和誤差分析。

水滴隨機噴濺數(shù)學模型生成的100個水滴的偏移角φ的分布如圖3(a)所示，在概率密度較小的區(qū)域，生成的水滴個數(shù)過于稀少，甚至為0。圖3(b)為水滴分檔隨機噴濺數(shù)學模型產(chǎn)生的100個分檔水滴的偏移角分布情況，在定義域內(nèi)非常均勻。

圖3 兩個數(shù)學模型所產(chǎn)生水滴的偏移角分布

4.1 水滴隨機噴濺數(shù)學模型的數(shù)值解

當水滴總數(shù)分別為1 000、1萬、10萬和100萬時，采用水滴隨機噴濺數(shù)學模型的求解結(jié)果如圖4所示。

圖4 水滴隨機噴濺模型的數(shù)值解與解析解

當水滴總數(shù)為1 000時，水滴隨機噴濺數(shù)學模型的數(shù)值解相比解析解有較大的誤差，誤差的波動幅度也很大(圖4(a))。水滴個數(shù)達到1萬時，數(shù)值解雖然局部與解析解仍然有較大偏差，但誤差明顯減小，誤差的波動幅度也縮小(圖4(b))。隨著水滴個數(shù)繼續(xù)增加，水滴隨機噴濺數(shù)學模型的求解誤差依然比較為明顯，但繼續(xù)減小(圖4(c)(d))。

4.2 水滴分檔隨機噴濺數(shù)學模型的數(shù)值解

分檔水滴總數(shù)分別為1 000、1萬、10萬和100萬時，采用水滴分檔隨機噴濺數(shù)學模型的數(shù)值解結(jié)果如圖5所示。

圖5 水滴分檔隨機噴濺模型的數(shù)值解與解析解

分檔水滴總數(shù)為1 000時，結(jié)果如圖5(a)所示，水滴分檔隨機噴濺數(shù)學模型的數(shù)值解與解析解的差異較大，但相比圖4(a)要小許多，波動幅度相比也更小，呈現(xiàn)相對規(guī)律的鋸齒狀。當分檔水滴個數(shù)增加10倍達到1萬時的結(jié)果如圖5(b)所示，整體誤差已經(jīng)較小，已看不出數(shù)值解與解析解有明顯差異。再繼續(xù)增加分檔水滴的個數(shù)(圖5(c)(d))，對計算精度的改善作用已經(jīng)很小。

對比圖5(b)和圖4(d)可見，1萬個分檔水滴的誤差，已能達到水滴隨機噴濺數(shù)學模型使用100萬個水滴求解時的水平。

4.3 誤差分析

為量化水滴隨機噴濺數(shù)學模型和水滴分檔隨機噴濺數(shù)學模型數(shù)值解的求解精度，評估數(shù)值解與解析解之間的誤差大小。這里采用最大誤差η和平均誤差ε兩個指標進行評價，計算公式如下：

(22)

(23)

式中：pm,a為第m個分組的解析解；pm,s為第m個分組的數(shù)值解。

誤差分析結(jié)果如表1所示。隨著水滴個數(shù)的增加，水滴隨機噴濺數(shù)學模型和水滴分檔隨機噴濺數(shù)學模型的誤差均逐漸減小。在水滴總數(shù)為1萬時，平均誤差指標已比較小，分別為11.46%和3.86%；最大誤差也迅速降低，水滴隨機噴濺數(shù)學模型的誤差從143.81%降至63.58%，水滴分檔隨機噴濺數(shù)學模型的誤差從19.42%下降至12.94%。

表1 兩個數(shù)學模型的數(shù)值解誤差和計算所需時間對比

隨水滴總數(shù)的變化，在兩個模型的水滴個數(shù)相同，甚至較少的條件下，水滴分檔隨機噴濺數(shù)學模型的最大誤差會更小。如在水滴總數(shù)為1萬時，水滴分檔隨機噴濺數(shù)學模型的最大誤差為12.94%，已經(jīng)比水滴隨機噴濺數(shù)學模型在水滴總數(shù)為100萬時的15.26%更小。

進一步對兩個模型的平均誤差，也顯示出水滴分檔隨機噴濺數(shù)學模型的精度更高。當水滴分檔隨機噴濺數(shù)學模型在水滴總數(shù)為1萬時，平均誤差已經(jīng)下降到3.86%，而水滴隨機噴濺數(shù)學模型在總數(shù)達到100萬時，平均誤差才下降為3.66%，兩者誤差十分接近。

最后從收斂趨勢來看，水滴隨機噴濺數(shù)學模型，水滴總數(shù)從1 000增加到100萬時，最大誤差從143.81%，降低為15.26%，一直在逐漸下降。而水滴分檔隨機噴濺數(shù)學模型在水滴總數(shù)從1 000增加到1萬后，繼續(xù)增加水滴個數(shù)，誤差變化不大。表明了水滴分檔隨機噴濺數(shù)學模型的誤差收斂速度更快。

4.4 兩個數(shù)學模型計算所需時間分析

在i7-3 770處理器、24GB內(nèi)存的計算機配置，以及Win7系統(tǒng)、Python 3.5.2編程環(huán)境的條件下，兩個數(shù)學模型計算所需時間如表1所示。水滴總數(shù)從1 000到100萬，水滴隨機噴濺模型的計算所需時間，從0.002 8 s增加到0.405 4 s；相對應的，水滴分檔隨機噴濺模型計算所需時間從0.002 1 s增加到0.362 6 s。同一個數(shù)學模型，水滴總數(shù)越大，計算所需時間越長，呈現(xiàn)線性關系。

在相同水滴總數(shù)情況下，兩個數(shù)學模型計算所需時間的相對差異并不是很大，處于同一個數(shù)量級，但水滴分檔隨機噴濺模型所需時間更少。在水滴總數(shù)為1 000時，計算所需時間分別為0.002 8 s和0.002 1 s；而在水滴總數(shù)為100萬時，計算所需時間分別為0.405 4 s和0.362 6 s。

隨著水滴總數(shù)的增加，兩個數(shù)學模型計算所需時間相差會更大。在水滴總數(shù)為1 000時，兩個數(shù)學模型計算所需時間相差0.000 7 s；而在水滴總數(shù)為100萬時，兩個數(shù)學模型計算所需時間相差0.042 8 s。

參與計算的水滴總數(shù)越大，兩個數(shù)學模型的數(shù)值解與解析解的誤差越小。在水滴總數(shù)相同的條件下，水滴分檔隨機噴濺數(shù)學模型的誤差明顯更小，模型的計算所需時間也更少。水滴分檔隨機噴濺數(shù)學模型可以采用更小的水滴總數(shù)，減少所需計算時間，同時保證計算的準確度，從而提高了計算效率。

5 結(jié) 論

a.針對水滴隨機噴濺數(shù)學模型計算量大，效率較低的問題，提出了分檔水滴的新概念，建立了挑流泄洪霧化的水滴分檔隨機噴濺數(shù)學模型。

b.在兩個假定條件下，對水滴分檔隨機噴濺數(shù)學模型進行了驗證計算，結(jié)果顯示該模型具有更好的求解精度。其中1萬個水滴的分檔隨機噴濺數(shù)學模型的求解結(jié)果，非常接近100萬個水滴的隨機噴濺數(shù)學模型的求解精度。

c.在水滴總數(shù)相同的情況下，水滴分檔隨機噴濺數(shù)學模型計算所需時間更少，同時其數(shù)值解的誤差也更小，表明水滴分檔隨機噴濺數(shù)學模型的計算效率更高。