考慮級配參數(shù)空間變異性的堆石壩隨機(jī)有限元分析

陸雪妮，朱 晟

(河海大學(xué)水利水電學(xué)院，江蘇 南京 210024)

0 引 言

堆石壩工程中堆石料的級配變動會在很大程度上影響材料的物理力學(xué)性質(zhì)，進(jìn)而影響壩體的變形應(yīng)力。已有研究表明，堆石料的物理力學(xué)特性與其級配分布緊密相關(guān)，李罡等[1]通過數(shù)值試驗(yàn)和室內(nèi)試驗(yàn)研究，發(fā)現(xiàn)顆粒級配對材料應(yīng)力應(yīng)變關(guān)系的影響顯著；趙婷婷等[2]基于顆粒流試驗(yàn)，發(fā)現(xiàn)表征顆粒級配分布的分形維數(shù)與堆石料的力學(xué)特性存在明顯相關(guān)關(guān)系，可將分形維數(shù)用于量化分析堆石料的力學(xué)特性；加力別克等[3]根據(jù)一系列室內(nèi)三軸壓縮試驗(yàn)結(jié)果，整理了分形維數(shù)與鄧肯E-B模型參數(shù)之間存在的二次函數(shù)關(guān)系。由此可見，堆石料的級配參數(shù)作為表征材料分布的重要參數(shù)，在堆石壩有限元計(jì)算中應(yīng)當(dāng)加以考慮，不可忽視堆石料級配效應(yīng)對壩體力學(xué)響應(yīng)的影響。朱晟等[4]結(jié)合不同級配堆石料的三軸試驗(yàn)結(jié)果，分析了堆石料級配參數(shù)與各力學(xué)特性之間的關(guān)系，基于廣義塑性模型理論，建立了本構(gòu)模型參數(shù)與分形維數(shù)之間的函數(shù)關(guān)系式，并驗(yàn)證了合理性，從而為考慮級配效應(yīng)的有限元計(jì)算提供了有效途徑。

目前，傳統(tǒng)的確定性有限元計(jì)算方法對于大壩各分區(qū)材料均采用單一級配的假定，沒有考慮到堆石料級配參數(shù)的空間變異性。即使是同一分區(qū)的堆石料，由于受到自然條件、母巖性質(zhì)、填筑過程等多種因素影響，其級配參數(shù)也會存在分布上的空間變異性，這種空間變異性是與位置和距離緊密相關(guān)的。因此，為了可以更好地反映工程實(shí)際，可在考慮級配效應(yīng)的有限元計(jì)算基礎(chǔ)上，引入隨機(jī)場理論和相關(guān)離散方法，用隨機(jī)場的統(tǒng)計(jì)特性對堆石料級配參數(shù)的空間變異性進(jìn)行表征，建立堆石料級配參數(shù)的三維隨機(jī)場模型，結(jié)合考慮級配效應(yīng)的堆石壩有限元分析方法對壩體力學(xué)響應(yīng)進(jìn)行研究。

本文以某抽水蓄能電站面板堆石壩作為工程實(shí)例，對堆石壩進(jìn)行有限元網(wǎng)格剖分，整理統(tǒng)計(jì)堆石料的級配參數(shù)和靜力計(jì)算本構(gòu)模型參數(shù)，采用Cholesky分解法建立堆石料級配參數(shù)的三維空間隨機(jī)場模型，結(jié)合考慮級配效應(yīng)的堆石壩有限元分析方法，通過多次隨機(jī)模擬，對隨機(jī)計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)分析，討論堆石料級配參數(shù)的空間變異性對堆石壩變形安全的影響。

1 對數(shù)正態(tài)隨機(jī)場模擬方法

已有眾多學(xué)者[5-11]對堆石料材料參數(shù)隨機(jī)場模型的建立展開研究，考慮到對數(shù)正態(tài)隨機(jī)場在對數(shù)值大于零的模型參數(shù)進(jìn)行模擬時的良好效果，研究中將筑壩堆石料材料參數(shù)均擬定為服從對數(shù)正態(tài)分布的隨機(jī)變量。本文采用對數(shù)正態(tài)隨機(jī)場對堆石料級配的分形維數(shù)進(jìn)行隨機(jī)模擬。

在隨機(jī)場模型中，一般采用理論自相關(guān)函數(shù)來描述堆石料同一分區(qū)中兩個不同空間位置處參數(shù)之間的自相關(guān)性[12]，由于高斯型自相關(guān)函數(shù)具有較好的平穩(wěn)性，本文選用高斯型自相關(guān)函數(shù)進(jìn)行表征，即

(1)

(2)

式中，δx、δy、δz分別為3個坐標(biāo)方向上的波動范圍；θx、θy、θz分別為3個坐標(biāo)方向上的相關(guān)距離；τx=|xi-xj|，τy=|yi-yj|和τz=|zi-zj|為任意2個空間位置點(diǎn)之間的平行于坐標(biāo)方向的間距。

對數(shù)正態(tài)隨機(jī)場模擬方法的具體步驟為：

(1)根據(jù)有限元計(jì)算網(wǎng)格將隨機(jī)場離散劃分成一系列單元，令隨機(jī)場單元與有限元單元一致，計(jì)算各單元中心點(diǎn)坐標(biāo)(xi，yi，zi)，i=1， 2，…，ne，ne為隨機(jī)場研究域內(nèi)的單元數(shù)量。

(2)采用拉丁超立方抽樣(LHS)方法隨機(jī)生成獨(dú)立標(biāo)準(zhǔn)正態(tài)空間隨機(jī)向量ξ(ξ為含ne個隨機(jī)抽樣點(diǎn)的列向量)，經(jīng)過m次隨機(jī)抽樣(m為隨機(jī)模擬次數(shù))，可以得到m個列向量。

(3)將步驟(1)得到的各單元中心點(diǎn)坐標(biāo)代入式(1)，計(jì)算得到空間任意兩點(diǎn)之間的自相關(guān)系數(shù)，構(gòu)成自相關(guān)系數(shù)矩陣

(3)

當(dāng)堆石料級配參數(shù)原始隨機(jī)場為對數(shù)正態(tài)隨機(jī)場時，需要將原始自相關(guān)矩陣C等概率轉(zhuǎn)換成相應(yīng)的標(biāo)準(zhǔn)正態(tài)空間的等效自相關(guān)系數(shù)C0，計(jì)算公式為

(4)

(4)采用Cholesky分解方法對上一步得到的等效自相關(guān)矩陣C0進(jìn)行分解，得到下三角矩陣L(ne×ne)，乘上步驟(2)得到的獨(dú)立標(biāo)準(zhǔn)正態(tài)隨機(jī)向量樣本ξ，即能夠求得標(biāo)準(zhǔn)參數(shù)高斯隨機(jī)場

(5)

然后結(jié)合均值μ、標(biāo)準(zhǔn)差σ等概率分布統(tǒng)計(jì)特征值，進(jìn)行等概率變換，將上述得到的標(biāo)準(zhǔn)參數(shù)高斯隨機(jī)場轉(zhuǎn)換成相應(yīng)的對數(shù)正態(tài)隨機(jī)場，具體變換公式為

圖1 大壩典型剖面示意

(6)

(5)以上過程實(shí)現(xiàn)了一次對數(shù)正態(tài)隨機(jī)場的建立，對于步驟(2)中的m次隨機(jī)抽樣結(jié)果，將步驟(3)～步驟(4)進(jìn)行m次重復(fù)，即實(shí)現(xiàn)了對數(shù)正態(tài)隨機(jī)場的m次模擬與建立。

2 工程實(shí)例

2.1 單元網(wǎng)格剖分與荷載分級

根據(jù)所研究堆石壩工程的壩址地形及壩體分區(qū)資料，圖1為大壩分區(qū)示意，模擬大壩填筑施工過程，逐級施加荷載，共分32級模擬大壩填筑過程，建立大壩三維有限元計(jì)算網(wǎng)格模型，如圖2所示，共剖分網(wǎng)格節(jié)點(diǎn)59 298個，總單元57 359個，壩體底部施加固定約束。

圖2 堆石壩有限元網(wǎng)格模型

2.2 堆石料級配參數(shù)空間隨機(jī)場模擬

采用分形分布模型[13]擬合主堆石區(qū)、下游堆石區(qū)等主要分區(qū)壩料的顆粒級配曲線。對于堆石料，顆粒質(zhì)量的分形分布模型為

(9)

式中，P(di)為小于di的顆粒質(zhì)量分?jǐn)?shù)；di為顆粒直徑；dmax為最大粒徑；D為分形維數(shù)。

結(jié)合施工填筑檢測級配資料，其中主堆石料的實(shí)測級配數(shù)據(jù)共計(jì)147組，下游堆石料級配數(shù)據(jù)共計(jì)99組，擬合各分區(qū)壩料的顆粒級配曲線，如圖3所示，相關(guān)系數(shù)均大于0.95，認(rèn)為堆石料級配基本滿足分形分布。

圖3 堆石料填筑級配曲線

繪制分形維數(shù)的頻率分布直方圖和對數(shù)正態(tài)概率分布圖，如圖4、5所示，正態(tài)概率圖用于檢查一組數(shù)據(jù)是否服從正態(tài)分布，對主堆石區(qū)和下游堆石區(qū)的分形維數(shù)D取對數(shù)后繪制正態(tài)概率分布圖，發(fā)現(xiàn)散點(diǎn)近似呈一條直線，說明主堆石區(qū)和下游堆石區(qū)的分形維數(shù)D滿足對數(shù)正態(tài)分布。對各區(qū)分形維數(shù)D進(jìn)行統(tǒng)計(jì)分析，根據(jù)3σ準(zhǔn)則對各區(qū)分形維數(shù)D中含粗大誤差的數(shù)據(jù)予以剔除，得到的各區(qū)參數(shù)分布范圍為主堆石區(qū)2.34～2.71，下游堆石區(qū)2.44～2.71，見表1。

圖4 主堆石區(qū)分形維數(shù)D概率統(tǒng)計(jì)規(guī)律

圖5 下游堆石區(qū)分形維數(shù)D概率統(tǒng)計(jì)規(guī)律

表1 各區(qū)級配分形維數(shù)分布統(tǒng)計(jì)分析

表2 堆石料分形維數(shù)隨機(jī)場的數(shù)字特征

在確定隨機(jī)場特征值的基礎(chǔ)上，根據(jù)前文所述對數(shù)正態(tài)隨機(jī)場模擬方法，分別對主堆石區(qū)和下游堆石區(qū)的分形維數(shù)D進(jìn)行隨機(jī)模擬，采用中點(diǎn)映射法將各分區(qū)隨機(jī)場中的分形維數(shù)D賦值給對應(yīng)的壩體有限元網(wǎng)格單元，即可得出壩體堆石料分形維數(shù)的三維隨機(jī)場。本文對堆石料分形維數(shù)隨機(jī)場進(jìn)行了200次隨機(jī)模擬，根據(jù)其中一次隨機(jī)模擬得到的堆石料分形維數(shù)D的三維隨機(jī)場，繪制堆石料三維整體和最大斷面的分形維數(shù)D的空間分布云圖，如圖6所示，發(fā)現(xiàn)主堆石區(qū)和下游堆石區(qū)的分形維數(shù)值分布沒有明顯的分界，主要原因是兩個分區(qū)分形維數(shù)的統(tǒng)計(jì)特征值之間差異不大，且變異系數(shù)均較小。在建立三維隨機(jī)場模型后，即可利用隨機(jī)場模型中的分形維數(shù)值進(jìn)行壩體有限元計(jì)算。

圖6 主堆石區(qū)和下游堆石區(qū)的分形維數(shù)D的空間分布云圖

圖7 竣工期典型斷面應(yīng)力變形分布云圖(單位：位移cm；應(yīng)力MPa)

2.3 靜力計(jì)算本構(gòu)模型參數(shù)

計(jì)算采用可以反映堆石料剪脹性的統(tǒng)一廣義塑性模型[16]。考慮級配效應(yīng)對大壩有限元計(jì)算的影響，需要研究筑壩料級配參數(shù)與其力學(xué)特性之間的關(guān)系。已有研究[4]以大量堆石壩工程的試驗(yàn)結(jié)果為依據(jù)，分析了表征抗剪強(qiáng)度、剪脹性、壓縮性等特性的參數(shù)與級配參數(shù)之間的關(guān)系，其建立的廣義塑性模型參數(shù)與分形維數(shù)之間的函數(shù)關(guān)系式可將級配參數(shù)直接應(yīng)用于堆石壩有限元計(jì)算中。本文主要研究考慮堆石料級配參數(shù)空間變異性時壩體的結(jié)構(gòu)響應(yīng)，因此對于庫盆回填區(qū)、過渡區(qū)、墊層區(qū)等體積較小、對大壩主體影響較小的分區(qū)筑壩料，堆石料的分形維數(shù)按照單一級配的假定作均一化處理，僅探究對壩體變形應(yīng)力影響較大的主堆石區(qū)和下游堆石區(qū)的分形維數(shù)隨機(jī)分布。

對于主堆石區(qū)、下游堆石區(qū)的堆石料，需要考慮級配參數(shù)的空間變異性，對主堆石區(qū)、下游堆石區(qū)的堆石料級配參數(shù)隨機(jī)場進(jìn)行離散，實(shí)現(xiàn)足夠多次數(shù)的對數(shù)正態(tài)隨機(jī)場的模擬與建立，然后根據(jù)單元中心點(diǎn)的坐標(biāo)，將所建隨機(jī)場模型中的分形維數(shù)值一一映射到相應(yīng)位置的網(wǎng)格單元中，根據(jù)分形維數(shù)與廣義塑性模型參數(shù)之間對應(yīng)的函數(shù)關(guān)系式[4]，完成堆石料單元級配參數(shù)與本構(gòu)模型參數(shù)之間的映射，從而直接將級配參數(shù)直接應(yīng)用于壩體有限元計(jì)算中。

2.4 考慮級配參數(shù)空間變異性的堆石壩三維有限元應(yīng)力變形分析

以與圖6隨機(jī)場對應(yīng)的隨機(jī)有限元計(jì)算結(jié)果作為典型案例，選取壩體最大斷面作為典型斷面，整理竣工期堆石壩應(yīng)力變形隨機(jī)響應(yīng)結(jié)果，如圖7所示。圖7中，位移豎直向上為正，豎直向下為負(fù)；水平位移向下游為正，向上游為負(fù)；壓應(yīng)力為正，拉應(yīng)力為負(fù)。

由圖7a可知，壩體沉降分布規(guī)律符合一般規(guī)律，壩體沉降較顯著的位置主要集中于主堆石區(qū)、下游堆石區(qū)以及二者交界處。隨機(jī)性有限元方法算得的竣工期壩體最大豎直沉降量為98.5 cm，位于下游次堆石區(qū)183.7 m高程附近，約1/2壩高處。

由圖7b可知，壩體的水平向位移基本上主要指向下游，這主要是受到了向下游傾斜的壩基地形和堆石體的泊松效應(yīng)兩方面的影響。隨機(jī)性有限元方法計(jì)算得到竣工期的水平向位移最大值為16.3 cm，指向下游方向，位于下游次堆石區(qū)的170.0 m高程附近。

由圖7c、7d可知，竣工期的堆石壩大、小主應(yīng)力的等值線分布均呈現(xiàn)出和壩坡近似平行的趨勢，并且應(yīng)力值隨高程自上而下逐漸增大，在壩基處達(dá)到最大值。隨機(jī)性有限元方法算得壩體大主應(yīng)力最大值為2.92 MPa，位于主堆石區(qū)底部地基處；小主應(yīng)力最大值為1.40 MPa，也位于主堆石區(qū)底部位置。

圖8 隨機(jī)有限元計(jì)算結(jié)果的統(tǒng)計(jì)穩(wěn)定性(竣工期)

3 隨機(jī)有限元計(jì)算結(jié)果分析

3.1 計(jì)算結(jié)果的統(tǒng)計(jì)穩(wěn)定性

對壩體主堆石區(qū)和下游堆石區(qū)堆石料的分形維數(shù)分別進(jìn)行200次隨機(jī)模擬后，完成相應(yīng)的隨機(jī)有限元計(jì)算，在對200次隨機(jī)有限元計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)分析前，需要先判斷經(jīng)過多次隨機(jī)模擬后的輸出結(jié)果是否達(dá)到統(tǒng)計(jì)穩(wěn)定狀態(tài)。通過繪制隨機(jī)計(jì)算結(jié)果的滑動平均值以及滑動標(biāo)準(zhǔn)差隨模擬次數(shù)的波動曲線，作為判斷隨機(jī)模擬是否穩(wěn)定的依據(jù)，若曲線波動性減弱、波動趨于平穩(wěn)，則判定隨機(jī)模擬計(jì)算結(jié)果達(dá)到穩(wěn)定狀態(tài)，認(rèn)為隨機(jī)模擬次數(shù)已足夠。

滑動平均值Runava(n)和滑動標(biāo)準(zhǔn)差Runsd(n)分別定義為

(10)

(11)

根據(jù)以上計(jì)算公式，將竣工期壩體豎直沉降、水平向位移、大主應(yīng)力及小主應(yīng)力等變形應(yīng)力指標(biāo)的最大值作為壩體響應(yīng)值，繪制各響應(yīng)值的滑動平均值Runava和滑動標(biāo)準(zhǔn)差Runsd隨模擬次數(shù)變化的波動示意圖，如圖8所示。可以發(fā)現(xiàn)當(dāng)隨機(jī)模擬次數(shù)達(dá)到60次時，各響應(yīng)值的波動曲線波動性均已減弱，曲線趨于平穩(wěn)，當(dāng)模擬次數(shù)達(dá)到140次以上曲線波動性已經(jīng)非常穩(wěn)定，認(rèn)為曲線達(dá)到收斂狀態(tài)。因此對壩體實(shí)現(xiàn)200次隨機(jī)場模擬，認(rèn)為隨機(jī)模擬次數(shù)滿足穩(wěn)定要求，由此得到的統(tǒng)計(jì)分析結(jié)果足以很好地代表總體的統(tǒng)計(jì)特性。

3.2 統(tǒng)計(jì)分析

對200次隨機(jī)有限元計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)分析，統(tǒng)計(jì)特性如表3所示，由均值和標(biāo)準(zhǔn)差計(jì)算得到各響應(yīng)值的變異系數(shù)分別為1.60%、1.56%、3.27%、1.85%、2.52%，各響應(yīng)值的概率分布情況如圖9、10所示。對壩體各隨機(jī)響應(yīng)值分別繪制正態(tài)概率分布圖，用于檢驗(yàn)是否服從正態(tài)分布，如圖9、10所示，各響應(yīng)值的正態(tài)概率分布圖上的點(diǎn)可以較好地用直線擬合，認(rèn)為滿足正態(tài)分布，由此表明，堆石料級配的分形維數(shù)滿足對數(shù)正態(tài)分布時，壩體應(yīng)力變形基本滿足正態(tài)分布。

3.3 與確定性有限元比較

如果不考慮堆石料級配參數(shù)的空間變異性，大壩各分區(qū)按照單一級配計(jì)算，即主堆石區(qū)、下游堆石區(qū)的分形維數(shù)取其平均值計(jì)算，作為確定性有限元計(jì)算結(jié)果，列于表4中。根據(jù)各壩體響應(yīng)值的均值表征隨機(jī)性有限元的計(jì)算結(jié)果，計(jì)算得到竣工期最大豎直沉降為97.8 cm，比確定性有限元方法算得的最大豎直沉降95.1 cm多出2.84%；最大水平向位移為16.3 cm，比確定性有限元方法算得的最大水平向位移16.0 cm多出1.88%；大主應(yīng)力最大值為2.98 MPa，比確定性有限元方法算得的大主應(yīng)力最大值2.85 MPa多出4.56%；小主應(yīng)力最大值為1.43 MPa，比確定性有限元方法算得的小主應(yīng)力最大值1.34 MPa多出6.72%。根據(jù)表4的計(jì)算結(jié)果比較情況可以看到，與忽略級配參數(shù)空間變異性的確定性有限元相比，隨機(jī)性有限元的變形和應(yīng)力均偏大，計(jì)算結(jié)果偏于安全，兩者相對誤差在2%～7%左右，由此可見，如果不考慮堆石料級配參數(shù)的空間變異性，很可能會在一定程度上低估壩體變形和應(yīng)力。

表3 壩體200次隨機(jī)模擬響應(yīng)值的統(tǒng)計(jì)特性(竣工期)

圖9 竣工期壩體變形概率統(tǒng)計(jì)規(guī)律

圖10 竣工期壩體應(yīng)力概率統(tǒng)計(jì)規(guī)律

表4 考慮級配參數(shù)空間變異性的隨機(jī)有限元計(jì)算結(jié)果與確定性有限元對比(竣工期)

4 結(jié) 論

本文以某面板堆石壩工程為例，發(fā)展了基于Cholesky分解的堆石壩隨機(jī)有限元計(jì)算方法，建立堆石料級配參數(shù)的三維隨機(jī)場模型，結(jié)合考慮級配效應(yīng)的堆石壩有限元分析方法對壩體力學(xué)響應(yīng)進(jìn)行研究，主要結(jié)論如下：

(1)堆石料級配曲線基本滿足分形分布模型，分形維數(shù)D可較好地表征級配曲線，主堆石區(qū)、下游堆石區(qū)等分區(qū)的分形維數(shù)頻率分布近似服從對數(shù)正態(tài)分布。

(2)基于Cholesky分解法實(shí)現(xiàn)堆石料分形維數(shù)的對數(shù)正態(tài)隨機(jī)場模擬，并將隨機(jī)參數(shù)映射至有限元分析模塊，計(jì)算考慮級配參數(shù)空間變異性時的壩體結(jié)構(gòu)響應(yīng)分析。計(jì)算結(jié)果表明：當(dāng)隨機(jī)模擬至200次時，壩體各響應(yīng)值達(dá)到統(tǒng)計(jì)穩(wěn)定，且堆石料級配的分形維數(shù)滿足對數(shù)正態(tài)分布時，算出的壩體應(yīng)力變形基本滿足正態(tài)分布。

(3)與壩料各分區(qū)參數(shù)采用單一級配假定的確定性有限元相比，隨機(jī)性有限元的變形和應(yīng)力均偏大，計(jì)算結(jié)果偏于安全，兩者相對誤差在2%～7%左右。如果不考慮堆石料級配參數(shù)的空間變異性，會在一定程度上低估壩體變形和應(yīng)力。