紅石巖堰塞壩力學參數隨機場模擬研究

賈 宇 峰，葛 培 杰，相 彪，陳 永 康

(1.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024； 2.大連理工大學 工程抗震研究所，遼寧 大連 116024； 3.云南農業大學 水利學院，云南 昆明 650201)

0 引 言

2014年8月3日，在云南省魯甸縣和巧家縣交界處發生了6.5級地震。地震導致魯甸縣火德紅鄉紅石巖村的牛欄江干流兩岸發生大面積的山體崩塌和滑坡。地震引發的滑坡形成了順河向長約1 807 m，高103 m，寬307 m，體積約1×107m3的堰塞體，如圖1所示。從圖中可以看到，堰塞體左岸側為古滑坡體，河谷中心處為新堆積體，底部為河床沖積層。由于左岸古滑坡體穩定性較差，采用爆破方式拆除堰塞體極易引發次生災害。并且，堰塞體體積巨大，采用常規方式進行拆除不僅需要對上部古滑坡體進行加固處理、設計合適的排渣場地，還需要承擔高額的運輸費用和較長的處理時間。因此，通過在堰塞體中設立混凝土防滲墻和帷幕灌漿區，可將其改建為永久性擋水建筑物，從而建成紅石巖堰塞湖多功能大型綜合水利樞紐工程。其中，堰塞壩最大壩高103 m，庫容2.6億m3，達到大(2)型水庫規模。由圖1可以看到，地震引起滑坡形成的新堆積體位于河道中央，是堰塞壩防滲墻所在區域。為了準確分析新堆積體的力學參數特性，預測壩體變形，分析混凝土防滲墻的變形和應力分布特性，對堰塞壩新堆積體進行了詳細的現場勘察。

圖1 堰塞體現場Fig.1 Photo of the barrier dam

現場勘察結果表明，紅石巖堰塞壩新堆積體主要由弱風化、微新風化和新鮮的白云巖和灰巖組成。新堆積體上部18～40 m為孤石和塊石，夾雜少量碎石和砂；下部以塊石、碎石混合砂土為主。新堆積體上部細顆粒較少，存在架空現象，密度較低，下部顆粒級配基本連續，孔隙率低，密實度較高，變形模量較高。表1中列出了7次現場密度試驗的測量結果。圖2給出了現場鉆孔試驗測定的不同位置新堆積體的顆粒級配和相應的采樣位置分布和采樣點高程。

表1 堆積體現場密度采樣結果Tab.1 Field results of accumulation body density g/cm3

圖2表明新堆積體上游側表面區域(ZK133，埋深9.5～12.6 m)粗顆粒較多，1 mm以上顆粒百分比超過74%。在新堆積體中部(ZK109-2，埋深46.3～47.9 m)的中間區域，小粒徑顆粒含量明顯提高，1 mm以下顆粒含量超過40%，在下游側底層區域(ZK101-2，埋深 54.7～56.3 m)小顆粒含量更高，1 mm以下顆粒含量超過 50%。為了確定新堆積體顆粒的力學特性，進行了15組隨機取樣的巖石單軸抗壓強度試驗和4組直剪試驗。試驗結果表明，由于巖石風化程度不同，室內試驗的單軸抗壓強度差異性較大。干抗壓強度在27.0～192.8 MPa，濕抗壓強度在18.0～169.1 MPa，巖石抗剪角度在53.2°～61.7°，黏聚力在5.0～5.61 MPa。

圖2 現場鉆孔試驗采樣位置與試樣顆粒分布Fig.2 Sampling location and particle distribution of field borehole test

新堆積體的現場試驗結果表明：新堆積體的密度分布規律具有明顯的離散性；新堆積體的顆粒組成呈現顯著的空間變異性，并且其小顆粒百分比遠高于土石壩的堆石料，其大顆粒粒徑又明顯超出心墻摻粒料范圍；并且土顆粒由白云巖和灰巖混合構成，顆粒風化程度差異導致其強度存在較大的離散性。傳統的均質材料數值模型無法描述新堆積體力學參數的空間變異性。因此，本文根據空間隨機場理論，結合現場試驗和室內三軸試驗構建了紅石巖新堆積體力學參數的隨機場有限元模型，為紅石巖堰塞壩安全分析奠定了基礎。

1 隨機場理論的基本概念

空間變異性是天然土體固有的一種屬性。地震引起山體滑坡形成的紅石巖新堆積，其土體參數分布具有明顯的空間變異性。已有研究結果表明，自然土體的密度和強度參數的分布規律可以采用標準正態分布或自然對數正態分布描述[1]。在此基礎上，根據現場觀測數據，通過相關函數來描述土體內部不同空間坐標下的土體參數之間的相關性，建立描述土體參數空間變異性的隨機場有限元模型。

1.1 隨機場的構建

目前，建立隨機場有限元模型的成熟方法有Karhunen-Loeve級數展開法[2]、協方差矩陣分解法[3]等。上述離散方法在邊坡穩定性、盾構隧道變形分析方面應用較廣。如李典慶[4]、蔣水華[5]等學者采用Karhunen-Loeve級數展開法構建了土坡強度參數的隨機場，進行了土坡可靠度分析。袁葳等[6]學者則采用Karhunen-Loeve級數展開法描述土體內摩擦角和黏聚力的空間變異性。程紅戰等[7]學者則采用協方差矩陣分解法模擬土體彈性模量變異性，進行盾構隧道數值分析。李健斌等[8]學者則對采用Karhunen-Loeve級數法建立的盾構隧道隨機有限元模型的力學響應進行了敏感性分析。然而上述Karhunen-Loeve級數離散方法在針對不規則區域時，第二類Fredholm方程求解相對較難實現。而基于Cholesky分解的協方差矩陣分解法可適用于不規則區域且編程較為簡單，便于工程應用。另一方面，隨機場理論在土石壩安全評價分析方面應用較少。楊鴿等[9]學者提出考慮力學參數的空間變異性的土石壩動力分析結果明顯有別于傳統的確定性方法。

綜上所述，本文根據紅石巖堰塞壩現場試驗確定新堆積體力學參數的相關距離，采用高斯函數擬合壩體內各個單元之間的相關系數，通過室內三軸試驗確定力學參數的統計規律。采用協方差分解法建立堰塞壩新堆積體隨機有限元模型。

1.2 協方差分解法

堰塞壩新堆積體有限元模型中，任意兩個單元i和單元j之間的土體參數相關性可以采用三維高斯型相關函數描述，其表達式為

(1)

式中：τx、τy、τz表示i單元和j單元形心在x、y、z3個方向的距離；δx、δy、δz為3個方向的相關距離。由式(1) 可以求出壩體內任意兩個單元參數之間的相關系數(aij=ρij)，組成對稱正定的n階相關系數矩陣A，對A矩陣進行三角分解，可得：

A=LLT

(2)

其中，L為下三角矩陣。新堆積體隨機場有限元模型中，n個單元土體參數為

H(X)=μ+σLZ

(3)

式中：H(X)是n階參數向量，且X={(xi，yi，zi)，i=1，…，n}代表n個單元的形心坐標，Z是代表不同位置單元土體參數分布概率的n階獨立標準正交隨機向量，μ和σ分別為新堆積體參數的平均值和標準差。因此，由(3)式確定的新堆積體參數向量H符合各個單元之間的相關性，并且其統計參數也收斂于現場試驗和室內三軸試驗確定的統計值。

2 新堆積體參數相關距離的確定

2.1 空間遞推法

土體參數的相關距離δ是描述其空間分布規律影響范圍的參數，反映了不同方向上土體結構的變化程度[10]。本文通過空間遞推法計算相關距離[11]。空間遞推法通過方差折減函數Γ2(h)與土體參數樣本波動范圍h的關系計算相關距離。其基本步驟是將h取為取樣間距Δz的倍數，即h=iΔz(i=1，2，3…)，然后將相鄰i個數據點構成一組新的樣本，計算新樣本的方差為γ2(h)，S2為初始樣本數據的方差，那么方差折減函數為

(4)

2.2 根據現場試驗確定相關距離

紅石巖新堆積體由地震引發山體滑坡形成，無法通過壩體內各個高程不同位置測量的密度或孔隙比等參數計算相關距離。因此，本文通過現場勘察試驗測量的剪切波速、表面波速和視電阻率計算相關距離。圖3給出了部分測點的豎直向剪切波速沿深度的分布圖。圖4～5給出了兩個順河向剖面的表面波速和視電阻率分布情況。以豎向相關距離計算為例，主要計算步驟如下：

(1) 先對圖3中豎向剪切波速的樣本數據進行趨勢項分解，建立趨勢項參數沿豎向分布的經驗公式。并在此基礎上，從實測剪切波速中扣除趨勢項分量，確定由隨機項分量產生的方差S2。

(2) 根據現場實際采樣數據點的間距值Δz估算數據的波動范圍h，即h=iΔz。以臨近的i個數據點的剪切波速作為一組新樣本，求出該組數據的方差γ2(h)。

(4) 繪制方差折減函數與波動范圍的分布曲線，即Γ2(h)～h圖。

(5) 根據所建立的Γ2(h)～h分布曲線，確定方差折減函數值趨于平穩的波動范圍值h=i*Δz，根據δu=i*ΔzΓ2(i*)確定剪切波速的相關距離。

圖4 順河向剖面表面波速分布(單位：m/s)Fig.4 Surface wave velocity distribution of a section along the river

圖5 順河向剖面視電阻率分布(單位：Ω·m)Fig.5 Apparent resistivity distribution of a section along the river

已有研究結果表明，土體參數的相關距離可用來表征采樣范圍內不同位置土體參數之間的互相關性。因此，相同土層中剪切波速確定的相關距離同樣也可以用于描述其它土體參數的互相關性[12]。新堆積體的水平相關距離也可以采用上述方法由圖4～5中的實測表面波速和視電阻率確定。各個方向上的相關距離如表2所列，水平向相關距離取兩種現場數據的計算結果的平均值(33.94 m)。

根據豎向相關距離和水平相關距離，結合式(1)可以求出新堆積體內各個位置點的力學參數在豎直向和水平向的相關系數。進而構建符合現場地勘實驗數據分布規律的力學參數隨機場。

表2 紅石巖新堆積體物理參數的相關距離Tab.2 Correlation distance of Hongshiyan new accamulation material m

3 新堆積體力學參數的統計規律

3.1 試驗方案

新堆積體力學參數的統計規律根據室內常規三軸試驗確定。將紅石巖新堆積體的現場采樣級配通過等量替代法進行縮尺，擬定具有代表性的3條級配曲線，即上包線級配、平均線級配和下包線級配作為室內試驗級配，如圖6所示。

圖6 三軸試驗級配曲線Fig.6 Triaxial test grading curve

通過現場7次密度試驗，擬合得到新堆積體的密度分布曲線(均值為1.98 g/cm3，標準差為0.21)。結合密度的正態分布特征和室內試驗實際情況，確定室內試驗密度為2.28，2.17，2.06，1.95，1.84 g/cm3，如圖7所示。在試驗級配和密度確定之后，分別進行圍壓為400，800，1 200 kPa的多組室內固結排水三軸壓縮試驗，試驗方案見表3。

圖7 室內試驗密度確定Fig.7 Determination of laboratory test density

表3 室內三軸壓縮試驗方案

3.2 力學參數統計分析

按照表3的方案進行三軸試驗，典型的三軸試驗結果如圖8所示。鄧肯E-B模型能夠較好地描述土體的非線性力學特性，廣泛應用于土石壩數值分析，是目前SL 274-2020《碾壓式土石壩設計規范》推薦的土石料本構模型。根據三軸試驗結果整理了各個試驗方案下的鄧肯張E-B模型參數列入表4中。

圖8 試驗方案7的三軸試驗曲線Fig.8 Triaxial test curve of Test 7

由于E-B模型參數中，K、n、Kb和m是影響土石壩靜力分析的主要參數。因此，本文僅構建了上述4個參數的隨機場有限元模型。為了確定參數分布規律，對上述4種參數的試驗數值進行了K-S檢驗[13]。K-S檢驗法將數據樣本的累積頻數函數與理論分布函數差異的最大值D作為衡量樣本分布規律的統計量。表5給出了4個參數的D值，而D7，0.05為0.483 4。可以看到，根據對數正態分布計算的D值均小于分位值0.483 4，參數服從對數正態分布。這主要是因為鄧肯E-B模型參數在數值上都是正值，并且模型數值的分布離散性較大，其原因是粗粒含量對模型的力學特性存在較大影響，相同密度條件下隨著粗顆粒含量的增大模型參數變化幅度較大。

表4 新堆積體E-B模型參數Tab.4 E-B model parameters of new accumulation body

表5 模型參數統計規律Tab.5 Parameter statistical characters

4 新堆積體力學參數隨機場的構建

紅石巖堰塞壩三維有限元模型如圖9所示。有限元模型中單元數為10 083，節點數為8 579。

圖9 紅石巖堰塞體三維有限元模型Fig.9 Three-dimensional finite element model of Hongshiyan barries dam body

參數K、Kb、n、m均服從對數正態分布，而協方差矩陣分解法是在標準正態空間下進行的。因此，需要將標準正態空間變換到對數正態空間。對H取對數[8]，便可得出對數正態分布下的空間隨機場He：

He=exp(μln+σlnLZ)

(5)

(6)

式中：μln和σln為對數正態分布下參數的平均值和標準差。

根據式(5)，結合所建立的紅石巖堰塞壩三維有限元模型，生成新堆積參數的隨機場有限元模型：

(1) 采用拉丁超立方抽樣方法隨機產生一組獨立標準正態隨機向量Z。

(2) 根據有限元模型中各單元形心點的空間坐標Xi(i=1，2，…，10 083)，采用表2中的相關距離，根據式(1)求出新堆積體有限元模型的相關系數矩陣C，并進行三角分解，得到系數矩陣L。

(3) 根據式(6)求出對數正態分布下的平均值μln和標準差σln。

(4) 根據式(5)求出新堆積體單元各個參數的隨機場He。

圖10和圖11為參數K的某次隨機抽樣結果在典型斷面的分布情況。參數Kb、n、m的隨機場實現與參數K類似，由于篇幅有限暫未列出。

圖10 順河向剖面K參數隨機場分布Fig.10 Distribution of K parameter random field along the river section

圖11 橫河向剖面K參數隨機場分布Fig.11 Distribution of K parameter random field along the dam axial section

5 結 論

本文根據紅石巖堰塞壩新堆積體現場試驗資料，結合7組室內三軸試驗結果構建了K、n、Kb、m參數的隨機場模型，從中得出一些有益的結論：

(1) 通過空間遞推法依據現場剪切波速求出紅石巖堰新堆積體豎向相關距離為2.84 m，根據現場表面波速和視電阻率求出新堆積體水平相關距離平均值為33.94 m。

(2) 室內試驗參數結果經K-S檢驗后，發現新堆積體鄧肯張E-B模型參數K、Kb、n和m符合對數正態分布規律。

(3) 采用平穩性較好的高斯相關函數構建的協方差矩陣能夠有效地描述紅石巖新堆積體力學參數的空間變異性。