基于蒙特卡羅法的挑流消能風(fēng)險(xiǎn)分析

尤嘉惠，王長(zhǎng)新

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

基于蒙特卡羅法的挑流消能風(fēng)險(xiǎn)分析

尤嘉惠1，王長(zhǎng)新2

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

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

關(guān)鍵詞：挑流消能；風(fēng)險(xiǎn)率；Matlab；蒙特卡羅；隨機(jī)數(shù)

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

1影響挑流消能風(fēng)險(xiǎn)的不確定性因素

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

正因?yàn)榇嬖谥@些不確定性因素，風(fēng)險(xiǎn)才會(huì)產(chǎn)生。我們將挑流消能的風(fēng)險(xiǎn)不確定性大致分為以下兩種：水文不確定性、水力不確定性[3]。

1.1　水文不確定性

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

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

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

(1)

1.1.2洪水過程不確定性

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

(2)

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

1.2　水力不確定性

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

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

σ=1.25θ

(3)

(4)

1.2.2水位不確定性

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

(5)

1.2.3水力計(jì)算模型的不確定性

在對(duì)泄水建筑物進(jìn)行設(shè)計(jì)時(shí)所運(yùn)用到的水力學(xué)方面計(jì)算公式也存在不確定性，可通過一個(gè)不確定性系數(shù)λ來表示，規(guī)定其均值為1，可適當(dāng)?shù)淖们檫x擇其標(biāo)準(zhǔn)差。

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

2挑流消能風(fēng)險(xiǎn)計(jì)算模式

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

將挑流消能的風(fēng)險(xiǎn)定義為：在正常運(yùn)行情況下，挑坎末端至沖刷坑最深點(diǎn)的平均坡度大于最大臨界坡度的概率。計(jì)算模式為

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　水下射程(水面以下水舌長(zhǎng)度水平投影)L1

(10)

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

(11)

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

(12)

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

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

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

3.1　計(jì)算步驟[12]

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

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

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

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

3.2　隨機(jī)數(shù)抽樣

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

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

(1) P-Ⅲ型分布

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

(2) 極值Ⅰ型分布

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

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

(13)

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

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

4算例

4.1　工程概況

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

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

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

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

4.2　挑流風(fēng)險(xiǎn)計(jì)算

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

表3　各隨機(jī)變量取值表

(13)

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

表4　風(fēng)險(xiǎn)計(jì)算成果表

(4) 對(duì)比考慮洪水過程不確定性因素與否對(duì)風(fēng)險(xiǎn)的影響，見表5。

表5　洪水過程不確定性對(duì)風(fēng)險(xiǎn)的影響對(duì)比表

4.3　靈敏度分析

表6　靈敏度分析表

圖1　各不確定性因素對(duì)風(fēng)險(xiǎn)的影響(100000次)

圖2各不確定性因素對(duì)風(fēng)險(xiǎn)的影響(1000000次)

4.4　計(jì)算結(jié)果

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

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

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

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

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

5結(jié)論

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

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

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

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

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

參考文獻(xiàn)：

[1]李亮.基于混凝土大壩的泄洪建筑物風(fēng)險(xiǎn)指標(biāo)體系及泄洪風(fēng)險(xiǎn)研究[D].西安:西安理工大學(xué)，2007.

[2]徐祖信,郭子中.二灘大壩表孔跌流中孔挑流的沖刷風(fēng)險(xiǎn)分析[J].水動(dòng)力學(xué)研究與進(jìn)展，1990,5(3):102-109.

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

[4]王長(zhǎng)新.施工導(dǎo)流風(fēng)險(xiǎn)分析與計(jì)算[J].八一農(nóng)學(xué)院學(xué)報(bào),1994，12(17):11-15.

[5]陳鳳蘭,王長(zhǎng)新.施工導(dǎo)流風(fēng)險(xiǎn)分析與計(jì)算[J].水科學(xué)進(jìn)展，1996，7(4):361-366.

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

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

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

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

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

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

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

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

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

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

DOI:10.3969/j.issn.1672-1144.2015.04.029

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

作者簡(jiǎn)介：尤嘉惠(1990—)，女，河南項(xiàng)城人，碩士，研究方向?yàn)樗煽慷壤碚?。E-mail:873560879@qq.com 通訊作者：王長(zhǎng)新(1957—)，男，遼寧撫順人，教授，博士生導(dǎo)師，主要從事水工水力學(xué)、河流泥沙和水力可靠度理論方面的研究和教學(xué)工作。E-mail: 420273169@qq.com

中圖分類號(hào)：TV135.2+3 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 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