尤嘉惠，王長新

(1.烏魯木齊水業(yè)建設投資有限公司， 新疆 烏魯木齊 830049；2.新疆農(nóng)業(yè)大學 水利與土木工程學院， 新疆 烏魯木齊 830052)

基于蒙特卡羅法的挑流消能風險分析

尤嘉惠1，王長新2

(1.烏魯木齊水業(yè)建設投資有限公司， 新疆 烏魯木齊 830049；2.新疆農(nóng)業(yè)大學 水利與土木工程學院， 新疆 烏魯木齊 830052)

摘要：在現(xiàn)今實際的挑流工程中，既要保證挑流消能的局部沖刷不會影響壩身安全和岸坡穩(wěn)定，同時又要求更經(jīng)濟的造價。為了在安全與經(jīng)濟之間取得平衡，綜合考慮洪峰流量和洪水過程的不確定性，用蒙特卡羅方法對挑流消能進行了風險率計算，并分析了下泄流量、沖坑系數(shù)、洪峰消減系數(shù)等不確定性因素對風險的影響。結(jié)果表明在考慮洪水不確定性時，挑流消能的風險有所減小，下泄流量、沖坑系數(shù)和洪水過程不確定性對挑流消能風險的影響比較顯著。

關鍵詞：挑流消能；風險率；Matlab；蒙特卡羅；隨機數(shù)

挑流消能是工程泄水消能方式中常見的銜接消能方式，泄水建筑物在下泄洪水時，勢能轉(zhuǎn)化成動能，水舌在鼻坎挑出過程中消耗的能量一般不超過10%～20%[1]，下游河床承受絕大多數(shù)的能量。因此，挑流消能要求下游河床基巖要有較高的抗沖能力，使下游的局部沖刷不會危及壩身安全和岸坡穩(wěn)定。一般要求挑流挑距應不小于最大沖坑深度的3～4倍，這樣才能保證工程的安全[2]。但為了節(jié)省工程投資，還有工程布置的要求，建筑物需要盡可能減小泄流寬度，以至于增加了單寬流量。為了工程在安全與經(jīng)濟之間取得平衡，本文運用蒙特卡羅法對挑流消能進行風險分析。

1影響挑流消能風險的不確定性因素

挑流消能系統(tǒng)中影響挑流局部沖刷的因素較多，諸如：工程布置、下泄流量、挑坎尺寸、沖坑后下游水深、地質(zhì)條件和運行情況等，這些因素都具有不確定性。

正因為存在著這些不確定性因素，風險才會產(chǎn)生。我們將挑流消能的風險不確定性大致分為以下兩種：水文不確定性、水力不確定性[3]。

1.1　水文不確定性

1.1.1洪峰流量不確定性[4]

由于自然因素(地貌、氣溫、暴雨時間和區(qū)域等)的隨機性引起洪峰流量的不確定性，使得天然河道的實際來流量及水位超過設計值[5]。

洪峰流量的不確定性在挑流消能風險的分析計算中是相當重要的。根據(jù)資料，我國最大洪峰流量服從的是P-Ⅲ型分布[3]。其概率密度函數(shù)為：

(1)

1.1.2洪水過程不確定性

在進行挑流消能風險分析與計算時，洪水過程不確定性對它是有影響的，所以必須要考慮到。這是因為：實際發(fā)生的洪水過程與進行設計時調(diào)洪演算后得到的洪水過程不可能完全一致，這就使洪水在調(diào)洪后所得的最大下泄流量以及最高水位與設計值不一致。為了描述水庫調(diào)洪的作用，我們引入洪峰消減系數(shù)η這個概念：

(2)

式中：Q0為經(jīng)過水庫調(diào)蓄作用后能夠下泄的最大出庫流量；Qf為河道中洪水的洪峰流量。

1.2　水力不確定性

1.2.1過流結(jié)構(gòu)不確定性

在工程施工過程中，由于各個方面的原因，如施工技術、工藝、施工人員以及工程管理等，使得建筑物的實際過流結(jié)構(gòu)尺寸與設計時的過流結(jié)構(gòu)尺寸存在一定的誤差，從而工程實際的下泄流量也與設計下泄量不一致。這在工程上是允許的。根據(jù)誤差理論，對于過流結(jié)構(gòu)尺寸的不確定性進行研究，認為結(jié)構(gòu)尺寸是服從正態(tài)分布，將其設計值作為均值，而標準差按下式進行估算[7]：

σ=1.25θ

(3)

(4)

1.2.2水位不確定性

由于洪峰流量具有不確定性，洪水過程及工程運行情況也具有不確定性，這就導致了水位不確定性的產(chǎn)生。

(5)

1.2.3水力計算模型的不確定性

在對泄水建筑物進行設計時所運用到的水力學方面計算公式也存在不確定性，可通過一個不確定性系數(shù)λ來表示，規(guī)定其均值為1，可適當?shù)淖们檫x擇其標準差。

在用到的水力學計算公式中，選取相關變量及各種參數(shù)時也存在著不確定性，如：流速系數(shù)及糙率等。由中心極限定理可知，可以將這些參數(shù)當作是服從正態(tài)分布的，其設計值定為均值，在可變范圍內(nèi)根據(jù)式(3)和式(4)來估算其標準差[8]。

2挑流消能風險計算模式

綜合以上各不確定性因素，本文主要從挑流消能沖刷坑對挑流(泄水)建筑物影響的角度來分析計算挑流消能的風險。

將挑流消能的風險定義為：在正常運行情況下，挑坎末端至沖刷坑最深點的平均坡度大于最大臨界坡度的概率。計算模式為

R=P(ik

(6)

用極限狀態(tài)方程表示為

Z*=ik-i

(7)

2.1　沖刷坑深度ts[10]

ts=Kq0.5Z0.25-ht

(8)

式中：K為沖坑系數(shù)；q為下泄單寬流量；Z為上下游水位差；ht為下游水深。

2.2　挑流鼻坎至下游水面間的挑流距離，即空中挑距L0

(9)

2.3　水下射程(水面以下水舌長度水平投影)L1

(10)

式中：T為沖刷坑深度(從下游水位起算的)，T=Kq0.5Z0.25；β為挑流射出水舌的外緣與下游水面間的夾角；

(11)

將以上的各因素帶入功能函數(shù)，即：

(12)

3蒙特卡羅法(Monte Carlo Method法)

蒙特卡羅法是從目標函數(shù)對應的分布中，人為的生成一個含有特殊隨機變量的系列，然后對其進行模擬與計算，最終的風險是對得到的計算結(jié)果通過檢驗估算出來的[11]。具體步驟如下：

利用偽隨機數(shù)生成的方法，抽樣生成一組隨機數(shù)樣本值，然后將其轉(zhuǎn)化為服從各變量概率分布的隨機數(shù)組，再將其代入到極限狀態(tài)方程(Z*=R-S)中，得到大量的Z*值，統(tǒng)計Z*值小于0的次數(shù)，最后小于0的次數(shù)與總抽樣次數(shù)的比被當作該系統(tǒng)的風險值。

3.1　計算步驟[12]

(2) 根據(jù)Qf所屬的分布，代入相應數(shù)據(jù)，生成Qf的隨機數(shù)(用舍選法)；

(3) 模擬生成挑流消能計算中不確定因素(Qf、η、K、a、B、R、θ)的隨機數(shù)；

(4) 根據(jù)水力學公式，計算出平均坡度i，即沖刷坑深度與挑流射程的比值；

(5) 將(4)中求得的最大臨界坡度ik與平均坡度i相減，得到Z值。若Z大于0，則失效次數(shù)N=N+1，重復步驟(3)至(5)，當N等于M時結(jié)束重復；

3.2　隨機數(shù)抽樣

Matlab軟件其自身配備一些隨機數(shù)的發(fā)生器，它能夠直接產(chǎn)生服從某些特定分布的隨機數(shù)，如均勻分布、正態(tài)分布等[8]。但是，與本文計算相關的P-Ⅲ型分布和極值Ⅰ型分布，在軟件中沒有能直接產(chǎn)生其隨機數(shù)的發(fā)生器。

3.3　P-Ⅲ型分布與極值Ⅰ型分布的隨機數(shù)抽樣

(1) P-Ⅲ型分布

前文已經(jīng)說明，洪峰流量Qf服從的是P-Ⅲ型分布，其對應的概率密度函數(shù)已知。根據(jù)舍選抽樣法[13-14]，并利用Matlab軟件編制相關程序來產(chǎn)生該不確定性因素的隨機數(shù)。

(2) 極值Ⅰ型分布

求解極值Ⅰ型分布時，必須首先對x=F-1(r)進行求解，則服從該分布的隨機數(shù)x即可得到。文獻[15]已經(jīng)提供了產(chǎn)生該分布隨機數(shù)的相關公式。如下：

xi=μx-0.45σx-0.7797σxln(-lnui)

(13)

式中：ui為服從均勻分布的隨機數(shù)；xi為服從該極值Ⅰ型分布的隨機數(shù)。

本文計算中所涉及到的洪峰消減系數(shù)η經(jīng)過假設檢驗后，認為服從極值Ⅰ型分布。

4算例

4.1　工程概況

一水利工程主要由擋水建筑物大壩、泄水建筑物溢流壩及底孔組成。溢流壩段位于原河槽部位，兩側(cè)為非溢流壩段。設計洪水重現(xiàn)期為500 a一遇，相應洪峰流量為4 318.75 m3/s，校核洪水重現(xiàn)期為5 000 a一遇，相應洪峰流量為5 640.64 m3/s。從水文站給出的實測資料中可以看出，最大的洪峰流量為Q=3 040 m3/s。通過對歷史洪水資料進行頻率分析，得到的洪水特征參數(shù)及頻率結(jié)果見表1。經(jīng)過調(diào)洪驗算后的成果表見表2。

表1　洪水洪峰流量統(tǒng)計特征參數(shù)及頻率

表2　調(diào)洪計算成果表

重力壩高140.33 m，泄流段總長約為116 m，堰頂高程為594 m。溢流段采用鼻坎挑流消能，鼻坎高程494 m，反弧半徑17.3 m，挑角33°。

4.2　挑流風險計算

(1) 建立極限狀態(tài)方程：Z=ik-i(Qf、η、K、a、B、R、θ)，其中7個不確定性因素的取值見表3。

表3　各隨機變量取值表

(13)

(2) 分別產(chǎn)生服從P-Ⅲ型分布的下泄流量Qf的隨機數(shù)，服從極值Ⅰ型分布的洪峰消減系數(shù)η的隨機數(shù)，服從正態(tài)分布的挑坎尺寸(坎高a、坎寬B、反弧半徑R、挑角θ)和沖坑系數(shù)K的隨機數(shù)，將以上隨機數(shù)代入到極限狀態(tài)方程中，最后統(tǒng)計ik

表4　風險計算成果表

(4) 對比考慮洪水過程不確定性因素與否對風險的影響，見表5。

表5　洪水過程不確定性對風險的影響對比表

4.3　靈敏度分析

表6　靈敏度分析表

圖1　各不確定性因素對風險的影響(100000次)

圖2各不確定性因素對風險的影響(1000000次)

4.4　計算結(jié)果

(1) 由表4可知，考慮了洪峰流量和洪水過程不確定性后，挑流建筑物風險率平均為0.00214。

(2) 由表4可知，挑流消能的風險隨著模擬次數(shù)的改變而變化。模擬次數(shù)由1萬次逐漸增加到1000萬次時，挑流風險值是逐漸向?qū)嶋H風險值趨近的，但是實際值現(xiàn)階段還是未知的，這就需要今后對此進行大量的實驗或分析大量的樣本來得到。

(3) 由表5可知，考慮洪水過程不確定性時的風險小于不考慮它時的風險，風險減小了33.75%。由此可知，在分析計算挑流消能風險時，有必要考慮洪水過程(洪峰消減系數(shù)η)的不確定性。

(4) 由表6、圖1及圖2可知，下泄流量Qf(風險率平均為0.00204)、洪峰消減系數(shù)η(風險率平均為0.001915)和沖坑系數(shù)K(風險率平均為0.00145)相對來說，對挑流風險的影響較大；挑坎尺寸(坎高、坎寬、反弧半徑、挑角)對挑流風險影響較小。

(5) 通過Matlab軟件計算風險，模擬次數(shù)為1000萬時，耗時為1 179 s，提高了蒙特卡羅方法的效率。

5結(jié)論

(1) 用蒙特卡羅法計算挑流風險率的關鍵是產(chǎn)生服從已知分布的隨機數(shù)，這項工作的計算量是龐大的，而通過Matlab軟件能快速得到，解決了這一問題。

(2) 通過靈敏度分析可知，下泄流量、洪峰消減系數(shù)和沖坑系數(shù)對挑流風險的影響較大。

(3) 本文利用Matlab軟件進行編程，成功解決了軟件不能直接生成P-Ⅲ型分布及極值Ⅰ型分布的隨機數(shù)的這一問題。

(4) Matlab軟件直接運算矩陣的功能，能夠使編程不再用非常繁瑣的If語句，從而使編程工作得到大大的簡化，提高了效率。

(5) 本文只是利用Matlab軟件計算了挑流風險率，對于挑流如何優(yōu)化和系數(shù)的選取等問題還需要研究。

參考文獻：

[1]李亮.基于混凝土大壩的泄洪建筑物風險指標體系及泄洪風險研究[D].西安:西安理工大學，2007.

[2]徐祖信,郭子中.二灘大壩表孔跌流中孔挑流的沖刷風險分析[J].水動力學研究與進展，1990,5(3):102-109.

[3]徐祖信,郭子中.混凝土高壩泄洪可靠度研究[J].河海大學學報，1990,18(2):76-82.

[4]王長新.施工導流風險分析與計算[J].八一農(nóng)學院學報,1994，12(17):11-15.

[5]陳鳳蘭,王長新.施工導流風險分析與計算[J].水科學進展，1996，7(4):361-366.

[6]趙經(jīng)華.泄洪風險計算軟件包的研究及其在工程風險計算中的應用[D].烏魯木齊:新疆農(nóng)業(yè)大學,2005．

[7]楊惠蓮,張濤.誤差理論與數(shù)據(jù)處理[M].天津:天津大學出版社,1992.

[8]張虎.新疆某水利樞紐工程施工導流風險分析和導流方案多目標決策[D].烏魯木齊:新疆農(nóng)業(yè)大學,2013.

[9]邱秀云.水力學[M].烏魯木齊:新疆電子出版社,2004.

[10]李煒.水力計算手冊[M].北京:中國水利水電出版社,2006.

[11]吳世偉.結(jié)構(gòu)可靠度[M].北京:人民交通出版社,1990.

[12]王麗學,林鳳偉,汪可欣，等.基于蒙特卡洛模擬的泄洪風險率計算[J].人民長江，2008，39(19):20-22.

[13]叢樹錚.水文學的概率統(tǒng)計基礎[M].北京:水利出版社,1981．

[14]王長新,王惠民,徐祖信，等.泄洪風險計算方法的比較[J].水力發(fā)電，1996，12(1):13-16.

[15]呂滿英.考慮洪水過程不確定性的泄洪風險分析[D].烏魯木齊:新疆農(nóng)業(yè)大學,2002.

DOI:10.3969/j.issn.1672-1144.2015.04.029

收稿日期：2015-03-01修稿日期：2015-04-18

作者簡介：尤嘉惠(1990—)，女，河南項城人，碩士，研究方向為水力可靠度理論。E-mail:873560879@qq.com 通訊作者：王長新(1957—)，男，遼寧撫順人，教授，博士生導師，主要從事水工水力學、河流泥沙和水力可靠度理論方面的研究和教學工作。E-mail: 420273169@qq.com

中圖分類號：TV135.2+3 文獻標識碼： A 文章編號： 1672—1144(2015)04—0146—05

The Risk Analysis of the Trajectory Energy Dissipation Based on Monte Carlo Method

YOU Jiahui1, WANG Changxin2

(1.UrumqiWaterIndustryConstructionInvestmentCo.,Ltd.,Urumqi,Xinjiang830049,China;2.CollegeofHydraulicandCivilEngineering,XinjiangAgriculturalUniversity,Urumqi,Xinjiang830052,China)

Abstract：In actual hydraulic projects, it is required that the local scour of the trajectory energy dissipation should not affect the safety of the dam body and the stability of the slope, meanwhile a satisfactory low cost should be upheld. In order to achieve the balance between safety and economy, the risk of the trajectory energy dissipation was calculated using Monte Carlo method, with the consideration of the uncertainties of flood peak discharge and flood process. And then the effects of uncertainty factors on the risk of trajectory energy dissipation were analyzed, such as the discharge, the scour coefficient, the peak attenuation coefficient and etc. The results showed that the risk decreased when the flood process uncertainty was in consideration and the discharge, the scour coefficient and the flood process uncertainty significantly influenced the risk of the trajectory energy dissipation.

Keywords:trajectory energy dissipation; risk rate; Matlab; Monte Carlo method; random number