劉德偉 ，李連俠 ，廖華勝 ，羅樹

（1.水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室，四川 成都 610065；2.Michigan State University, East Lansing,USA；3.西南電力設(shè)計院有限公司，四川 成都 610021）

堰塞體浸潤線的Monte Carlo模擬

劉德偉1，李連俠1，廖華勝2，羅樹火昆3

（1.水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室，四川 成都 610065；2.Michigan State University, East Lansing,USA；3.西南電力設(shè)計院有限公司，四川 成都 610021）

浸潤線的高低對堰塞體的滲透穩(wěn)定至關(guān)重要，而堰塞體材料滲透系數(shù)往往具有很強(qiáng)的隨機(jī)性，如果采用確定性模型則不能準(zhǔn)確地模擬實際的滲流場及浸潤線位置。本文利用Monte Carlo方法在堰塞體材料滲透系數(shù)滿足對數(shù)正態(tài)分布的條件下，對唐家山堰塞體滲流場進(jìn)行了300次二維隨機(jī)模擬。結(jié)果表明，受材料隨機(jī)性的影響，水頭場、浸潤線分布及滲透坡降場等均呈現(xiàn)隨機(jī)特性，但和滲透系數(shù)場的隨機(jī)性并非簡單的線性關(guān)系，水頭方差空間變化特性較為明顯，越往下游其值越大；對比確定性模型，隨機(jī)性模型中出逸點的計算結(jié)果更具統(tǒng)計意義，對工程防護(hù)中有一定指導(dǎo)意義；由于受邊界條件影響較小，堰體中間部位滲流特性基本滿足正態(tài)分布；浸潤線逸出點處受邊界條件影響較大，逸出位置分布集中，不服從正態(tài)分布，其余部位浸潤線位置高程隨機(jī)過程基本服從正態(tài)分布。

滲流；堰塞體；Monte Carlo方法；浸潤線；滲透系數(shù)

1 研究背景

2008年5月12 日，汶川大地震引發(fā)了34座大型堰塞湖和數(shù)以百計的小型堰塞湖，對下游地區(qū)人民的生命財產(chǎn)將構(gòu)成巨大的威脅。堰塞體有不同于人工堆筑壩體的地方：堰塞體是由新鮮的、不均勻的、欠固結(jié)的滑坡土體堆積而成，級配較廣，土體組成粒徑從黏性顆粒到塊石不等[1]。透水性強(qiáng)和弱的材料相互摻，造成滲透系數(shù)分布的非均質(zhì)性。這種非均質(zhì)性使得應(yīng)用確定性方法處理非確定的堰塞湖滲流問題具有很大的局限性，很難客觀地反映實際的滲流特性。

劉海峰等[2]以汶川地震中綿遠(yuǎn)河堰塞湖群排險工程為例，定性分析堰塞湖群潰決的風(fēng)險影響因素；嚴(yán)祖文等[3]就根據(jù)既往的研究成果，主要闡述了堰塞湖天然壩壩體安全應(yīng)急評估的基本方法和原則，基于唐家山堰塞體的巖體結(jié)構(gòu)穩(wěn)定性分析、宏觀現(xiàn)象監(jiān)控與地表位移監(jiān)測，研究堰塞體的整體穩(wěn)定性；崔銀祥等[4]通過三維滲流計算評價堰塞體滲透穩(wěn)定性，在進(jìn)行穩(wěn)定性評價時先評價壩體內(nèi)水力坡度最大方向上格點的滲流穩(wěn)定性，進(jìn)而評價整個壩體的滲透穩(wěn)定性。

從前人的研究成果來看，大多采用確定性模型進(jìn)行模擬，對堰塞湖的材料特性對滲流影響的研究很少；很少考慮到堰塞體材料滲透系數(shù)具有隨機(jī)性這一特點；對堰塞湖浸潤線模擬較少，而浸潤線的高低直接影響到堰塞湖的滲流穩(wěn)定。而隨機(jī)地下水模擬方法則考慮了介質(zhì)的空間變異性和由此引起的因變量的不確定性，能夠較為客觀地反映現(xiàn)實規(guī)律[5]。Monte-Carlo方法就是一種使用隨機(jī)樣本技術(shù)，通過計算機(jī)隨機(jī)模擬而獲得問題的近似解的方法，它特別適用于解具有一定概率分布的問題。本文依據(jù)現(xiàn)有堰塞湖資料，建立適當(dāng)?shù)母呕Ｐ?，采用自編的飽?非飽和滲流計算程序[6]，并引入Monte Carlo模擬進(jìn)行二維立面滲流數(shù)值模擬。

2 計算模型及參數(shù)的確定

目前，已經(jīng)有許多隨機(jī)數(shù)值方法用于分析地下水輸運問題方法，包括Monte-Carlo方法、擾動方法及基于Taylor展開的一次二階矩法等[7]。其中，Monte-Carlo方法使用最為廣泛。

2.1 隨機(jī)地下水模型的建立

二維飽和/非飽和滲流的基本方程[8]如下：

式中：kmh m為滲透系數(shù)，是壓力水頭h的函數(shù)，h= H-y；H為水頭；x為水平方向坐標(biāo)；y為垂向坐標(biāo)；n為土體孔隙度；c′mh m為相對飽和度；λ是飽和度，λ=茲/n，θ為體積含水率；SS為貯水率。

本文模擬堰塞體內(nèi)穩(wěn)定滲流的浸潤線，此時方程（1）的右邊項為零，求解時只需要給定k-h的關(guān)系即土水特征曲線，堰塞體兩種材料土水特征曲線采用Van Genuchten模型。

國內(nèi)外學(xué)者對于含水層的滲透系數(shù)做出了大量研究[5，6，9]，發(fā)現(xiàn)無論是在微觀顆粒還是宏觀區(qū)域來講k值都具有隨機(jī)性。Freeze[10]統(tǒng)計了大量飽和滲透系數(shù)資料，結(jié)果表明，飽和滲透系數(shù)可以用對數(shù)正態(tài)分布來描述，其概率密度分布如下：

式中：滋為飽和滲透系數(shù)均值；滓為飽和滲透系數(shù)標(biāo)準(zhǔn)差。

2.2 模型尺寸邊界條件的確定

如圖1所示，壩體底寬為堰塞壩順河向長度為803.4 m，壩頂高程為759.55 m，壩底高程為669.55 m。壩體主體為厚度90 m的碎裂巖。壩基為苦竹壩庫區(qū)沉積的8 m厚度含泥粉細(xì)砂，下伏風(fēng)化基巖。上游壩坡被滑坡擠壓隆起的含泥粉細(xì)砂覆蓋，其壩頂高程為735 m。上游壩坡坡度為20°。下游壩坡分為三段：上部坡比為 1∶0.7，坡高為50 m；下部坡比為 1∶0.5，坡高為 20 m；中部為緩坡，坡比為 1∶2.5。

由于模型底部粉細(xì)砂下臥風(fēng)化基巖，邊界作為不透水設(shè)置。順河向兩側(cè)為自由臨空面，上游為定水頭邊界（第一類邊界），下游自由面邊界，堰塞體腳處為定水頭邊界。據(jù)堰塞體工程地質(zhì)評價以及壩體組成物質(zhì)和經(jīng)驗判斷，模型中粉細(xì)砂和碎裂巖飽和滲透系數(shù)均值分別取8.64×10-2m/d和 8.64 m/d，Van Genuchten模型（式（2））中參數(shù)m均取1.58，參數(shù)α取0.001 7 m-1，其相應(yīng)的允許坡降分別為0.1和0.6。

圖1 唐家山堰塞體二維模型

在對比不同尺度的計算結(jié)果的情況下，取均值K1區(qū)域的相關(guān)尺度為2 000，取均值K2區(qū)域的相關(guān)尺度為1 000。Freeze曾系統(tǒng)地分析過一些含水量，認(rèn)為砂礫石含水層方差滓=0.92；淤泥層方差滓=2.14。本文在比較分析后，取方差為滓=3。為了對Monte-Carlo方法隨機(jī)計算過程進(jìn)行檢驗，在堰塞體內(nèi)部加入監(jiān)測井（見圖1）。監(jiān)測井不為源匯項，不影響計算結(jié)果，僅為結(jié)果輸出設(shè)立。

3 Monte-Carlo模擬結(jié)果分析

3.1 隨機(jī)過程檢驗

為檢驗Monte-Carlo方法樣本個數(shù)的統(tǒng)計代表性，共進(jìn)行了300次樣本模擬計算。圖2和圖3為不同監(jiān)測井不同樣本容量統(tǒng)計計算得到的對數(shù)滲透系數(shù)和水頭的樣本均值。由圖2，3可見，隨著樣本容量的增加，樣本均值快速趨近于定值。說明樣本容量300可到統(tǒng)計所需樣本容量的要求。

圖2 不同井監(jiān)測的對數(shù)滲透系數(shù)的樣本均值

圖3 各監(jiān)測井水頭均值

本文視飽和滲透系數(shù)隨機(jī)場為標(biāo)準(zhǔn)對數(shù)正態(tài)分布場，計算是否能保持應(yīng)有的統(tǒng)計特性直接影響到隨機(jī)場模擬結(jié)果的質(zhì)量。下面選取編號為W5，W6，W7，橫坐標(biāo)都為510 m，縱坐標(biāo)差值為20 m的3個監(jiān)測井考察統(tǒng)計性質(zhì)。

圖4 對數(shù)飽和滲透系數(shù)概率密度分布直方圖

從圖4對數(shù)飽和滲透系數(shù)概率密度分布直方圖看出，編號為W5，W6，W7等3個井的對數(shù)飽和滲透系數(shù)概率密度基本上關(guān)于lnk1=2.15，lnk2= 2.4，lnk3=2.16對稱分布，符合正態(tài)分布特征。

從圖5水頭概率分布直方圖看出，水頭概率密度基本上關(guān)于H水頭=730.6 m對稱分布，符合正態(tài)分布特征。

3.2 典型樣本計算結(jié)果分析

由飽和滲透系數(shù)KS的隨機(jī)性，可以推斷計算結(jié)果的隨機(jī)性，選取樣本號為50和250的兩個不同隨機(jī)樣本計算所得的水頭等值線分布、滲流速度分布、滲透坡降進(jìn)行結(jié)果分析。

圖5 水頭概率密度分布直方圖

圖6對比不同隨機(jī)場的計算結(jié)果，發(fā)現(xiàn)水頭分布沿程降低，降低幅度基本相同。受材料隨機(jī)性的影響，各個樣本之間等水頭線分布不盡相同，等水頭線在垂直方向有扭曲現(xiàn)象。

圖6 不同隨機(jī)場的水頭等值線及滲透速度分布圖

圖7 不同隨機(jī)場的滲透坡降等值線計算結(jié)果

圖8 確定性模型滲透坡降等值線圖

由圖7堰塞體內(nèi)滲透坡降等值線分布可以更明顯地看出材料隨機(jī)性對滲流的影響。各個隨機(jī)場計算的滲透坡降分布范圍大致相同，但是等值線幾乎不同，滲透坡降等值線破碎，不連續(xù)，隨機(jī)特性明顯，堰塞體內(nèi)最大滲透坡降均未超過允許坡降。對比圖8確定性模型滲透坡降等值線圖，可以看出隨機(jī)性模型中滲透坡降大于0.2的區(qū)域相對較大，不利于堰體的滲流穩(wěn)定。

3.3 浸潤線計算結(jié)果分析

為了研究分布規(guī)律，在本文的Monte Carlo數(shù)值模型中選取 X=410，510，560，660，760 m 等 5 個斷面及逸出點的浸潤線位置作為分析對象。各個斷面的統(tǒng)計結(jié)果如表1所示。

表1 不同斷面浸潤線高程統(tǒng)計規(guī)律計算表

圖9為300個樣本的浸潤線分布，圖中粗實線為采用確定性模型（滲透系數(shù)取均值）得到的浸潤線。

由圖9可以看出，受材料隨機(jī)性影響，浸潤線分布隨樣本不同而不同，隨位置向下游移動浸潤線分布逐漸分散，大約在X=650 m斷面處分散程度達(dá)到最大，由于受下游壩坡邊界影響，再往下游浸潤線出現(xiàn)聚攏現(xiàn)象。

圖9 浸潤線分布圖

現(xiàn)對300個樣本的逸出點高程做出統(tǒng)計結(jié)果見圖10。

圖10 逸出點高程分布直方圖

逸出點高程分布在674.28~684.25 m之間。其中有85.7%的逸出點高程分布在676.04~680.15 m之間。圖9黑色線為采用確定性模型計算的浸潤線，其出逸點高程為680.8 m，高于隨機(jī)性模型的浸潤線的平均高程678.603 m。說明采用均值的K進(jìn)行計算得到的結(jié)果與采用隨機(jī)K場進(jìn)行多樣本模擬后得到的結(jié)果的均值并非一致，隨機(jī)性模型中出逸點的計算結(jié)果更具統(tǒng)計意義，對在工程防護(hù)中擬訂防護(hù)措施和范圍具有明確的指導(dǎo)意義，減少盲目性。

4 結(jié) 論

1）Monte-Carlo模擬結(jié)果表明，受材料隨機(jī)性的影響，水頭場、浸潤線分布及滲透坡降場等均呈現(xiàn)隨機(jī)特性，不同樣本之間水頭分布不盡相同，水頭均值場與確定性模型計算結(jié)果非常接近；水頭方差場與滲透系數(shù)K的方差場的關(guān)系并非簡單的線性關(guān)系：水頭方差空間變化特性較為明顯，越往下游其值越大，說明最大的不確定性發(fā)生在堰塞體下游部位，其變化程度更為激烈。

2）采用均值的滲透系數(shù)K進(jìn)行計算得到的結(jié)果與采用隨機(jī)滲透系數(shù)場進(jìn)行多樣本模擬后得到的結(jié)果的均值并非一致，隨機(jī)性模型中出逸點的計算結(jié)果更具統(tǒng)計意義，對在工程防護(hù)中擬訂防護(hù)措施和范圍具有更明確的指導(dǎo)意義，減少盲目性。

3）典型采樣位置的模擬結(jié)果表明，水頭標(biāo)準(zhǔn)差在上游側(cè)均較小，隨位置向下游移動，由于下游逸出點的非線性，其值進(jìn)一步增大，使得逸出點附近水頭的不確定性程度增加，位于堰體中部的監(jiān)測點的概率密度分布更接近正態(tài)分布，其余各處由于受邊界條件影響其分布均不同程度偏離正態(tài)分布。說明輸入?yún)?shù)的正態(tài)特性并不能保證輸出結(jié)果的正態(tài)分布特性，輸入和輸出之間的非線性關(guān)系和邊界條件的影響會破壞這種概率分布的線性傳遞。

［1］劉海峰，周宏偉.地震堰塞湖風(fēng)險影響評價定性分析［J］.東北水利水電，2013（2）55—57.

［2］李守定，李曉，張軍，赫建明，李世海，汪陽春.唐家山滑坡成因機(jī)制與堰塞壩整體穩(wěn)定性研究［J］.巖石力學(xué)與工程學(xué)報，2010，29（S1）：2908—2915.

［3］嚴(yán)祖文，魏迎奇，蔡紅.堰塞壩滲透穩(wěn)定性評估［J］.長江科學(xué)院院報，2009，6（10）：122—125.

［4］崔銀祥,聶德新，等.通過三維滲流計算評價某滑坡壩滲透穩(wěn)定性［J］.水土保持研究，2005，12（2）：97—100.

［5］王亞軍，張我華，陳合龍.長江堤防三維隨機(jī)滲流場研究［J］.巖石力學(xué)與工程學(xué)報，2007，26（09）：1824—1824.

［6］廖華勝，李連俠，LI Shu-guang.地下水非平穩(wěn)隨機(jī)模型及空間變異性與非均勻性相互關(guān)系研究的展望［J］.水利學(xué)報，2004，35（10）：13—21.

［7］孫瑩潔，等.滑坡體滲流場與庫水位波動響應(yīng)關(guān)系研究［J］.人民長江，2014，45（7）：66—69.

［8］ SUDICKY E A.A natural gradient experiment on solute transport in a sand aquifer:spatial variability of hydraulic conductivity and its role in the dispersion process［J］.Water Research，1986，22（13）：2017—2029.

［9］ K.Binder,D.W.Heerman.MonteCarloSimulationinStatistical Physics An Introduction［M］.New York:Springer verlag, Berlin Heideberg，2010，5th ed.

［10］ R.A.Freeze.A stochastic-conceptual analysis of onedimensional groundwater flow in nonuniform homogeneous media.Water Resource Research，1975，11（5）：725—741.

TV131.61；P315.9

A

1002－0624（2016）04－0043－05

國家自然科學(xué)青年基金（51209154）。

李連俠，博士，副教授，主要從事工程水力學(xué)和地下水研究。

2015-10-10