多場耦合作用下天然土質(zhì)邊坡安全性數(shù)值仿真研究*

謝承煜， 劉承波， 朱佳妮， 鐘紅軍， 謝洪高， 盛美玉

(湘潭大學(xué) 環(huán)境與資源學(xué)院，湖南 湘潭 411105)

0 引言

隨著經(jīng)濟高速發(fā)展，各項基礎(chǔ)建設(shè)蓬勃興起，如水利水電、房屋建設(shè)、能源與交通以及礦山開采等，都不可避免地遇到各類邊坡安全問題.在工程實際中，邊坡往往承受多物理場耦合作用，邊坡體的存在狀態(tài)更是一個復(fù)雜的體系，包括壓力(如重力、邊坡開挖后的應(yīng)力重新分布以及地震活動引起的動態(tài)載荷)、水力(如地下水流滲入和雨水入滲)和邊坡體本身力學(xué)結(jié)構(gòu)等，在多種影響因素的耦合下，導(dǎo)致了邊坡體的失穩(wěn)破壞[1-4]，其所造成的人員和財產(chǎn)損失備受社會關(guān)注，其中，以天然土質(zhì)結(jié)構(gòu)為主的邊坡相較于巖質(zhì)邊坡而言更易于發(fā)生失穩(wěn)破壞，主要的災(zāi)害類型為滑坡災(zāi)害，是工程建設(shè)中的重點防治對象，且我國又是一個滑坡等自然災(zāi)害頻發(fā)的國家[5-8]，根據(jù)國家統(tǒng)計局地質(zhì)災(zāi)害數(shù)據(jù)顯示(如圖1)，以2019年為例，全年地質(zhì)災(zāi)害共發(fā)生6 181次，其中滑坡災(zāi)害占4 220次，累計占全年災(zāi)害次數(shù)的68.3%，是崩塌災(zāi)害的3倍.

圖1 我國近五年地質(zhì)災(zāi)害數(shù)據(jù)柱形圖Fig.1 Column chart of geological hazard data in recent five years in China

物理場耦合作用于邊坡體的穩(wěn)定性分析經(jīng)過了較長的發(fā)展歷程，由兩場耦合到三場耦合再到四場耦合的研究，隨著數(shù)學(xué)理論方法的不斷發(fā)展，科學(xué)技術(shù)的不斷進(jìn)步，對邊坡多場耦合作用的數(shù)學(xué)模型建立及數(shù)值模擬計算方式也得到了不同程度的改善.蔡亞飛等[9]運用ABAQUS對滲流-應(yīng)力耦合及降雨入滲作用下的邊坡進(jìn)行了穩(wěn)定性分析；張良以等[10]結(jié)合數(shù)值模擬運用對滲流場-應(yīng)力場-膨脹應(yīng)變場耦合下非飽和膨脹土邊坡逐漸破壞演化規(guī)律進(jìn)行了系統(tǒng)的研究；朱和玲等[2]運用數(shù)值模擬方法研究了基于降雨-地震耦合作用的土質(zhì)邊坡穩(wěn)定性；Wu等[11]對溫度場-滲流場-機械荷載耦合下的邊坡變化特性進(jìn)行了數(shù)值模擬，并與傳統(tǒng)的極限平衡計算結(jié)果進(jìn)行了對比；Qu等[12]運用COMSOL Multiphysics 對應(yīng)力場-滲流場-溫度場-化學(xué)場耦合下的巖質(zhì)邊坡進(jìn)行了穩(wěn)定性研究.盡管對多場耦合作用于邊坡體的穩(wěn)定狀況已取得了一些研究成果，但仍不夠完善，作用于邊坡體的多種物理場之間的相互作用、相互影響是一種非常復(fù)雜的關(guān)系，不同的研究手段、內(nèi)外在因素的不同等都會造成研究結(jié)果的差異，故仍需加以改進(jìn)或完善，使其更貼近于實際應(yīng)用.基于此，本文對在多場耦合作用下邊坡的安全性展開研究，利用數(shù)值模擬分析方法，揭示天然土質(zhì)邊坡在降雨-滲流-應(yīng)力場耦合作用下的位移變形規(guī)律，以期達(dá)到邊坡安全防護(hù)的目的.

1 多場耦合作用下土質(zhì)邊坡安全性數(shù)值模擬

運用COMSOL Multiphysics軟件進(jìn)行數(shù)值模擬計算，通過建立二維有限元模型，設(shè)置降雨-滲流-應(yīng)力耦合工況參數(shù)，并分階段進(jìn)行求解，從應(yīng)力、塑性、位移等云圖中提取數(shù)據(jù)，以此分析土質(zhì)邊坡的變化特性.

1.1 建立幾何有限元模型

(1) 幾何模型

為方便操作，導(dǎo)入在CAD軟件中構(gòu)建的二維邊坡模型，如圖2所示，其邊界范圍：左右邊界距離取54 m，上下邊界距離取30 m，坡面頂點距左邊界的距離取20 m，坡角距右邊界的距離取10 m、距下邊界的距離取10 m.同時，為了不影響邊界設(shè)定和求解域的完整性，對CAD導(dǎo)入模型作修復(fù)操作，相對修復(fù)容差為1.0×10-5[13].

圖2 二維土質(zhì)邊坡模型尺寸圖(單位：m)Fig.2 2D dimension drawing of soil slope model(unit:m)

(2) 網(wǎng)格剖分

為了更為精確高效地對幾何模型進(jìn)行網(wǎng)格剖分，將網(wǎng)格的序列類型切換為用戶控制網(wǎng)格，網(wǎng)格剖分方式為映射，共剖分1 699個單元(如圖3所示).其中，在映射設(shè)置中的“域選擇”依次點擊幾何模型的域，直至選中整個模型，且在尺寸設(shè)置中，將“預(yù)定義”設(shè)置為超細(xì)化，單元大小參數(shù)：最大單元大小為1.08 m，最小單元大小為0.004 03 m，最大單元增長率為1.2，曲率因子為0.25，狹窄區(qū)域分辨率為1[14].

圖3 二維土質(zhì)邊坡網(wǎng)格剖分圖Fig.3 2D mesh subdivision of soil slope

(3) 材料參數(shù)

材料參數(shù)的設(shè)定決定了材料的變化特性，邊坡土體的材料為彈塑性材料，服從Drucker-Prager準(zhǔn)則和莫爾-庫倫準(zhǔn)則[15-16]，其參數(shù)如表1所示.

表1 二維土質(zhì)邊坡土體的材料參數(shù)表

(4)邊界條件

邊坡模型的下邊界設(shè)置為固定約束，左右邊界為輥支承，則其余為自由邊界；同樣，在“體積力”的擴展選項中選擇重力，且在重力設(shè)置中選擇所有域，為邊坡模型添加自重力.

1.2 工況分析設(shè)計方案

土質(zhì)邊坡在受物理場作用的情況下，其應(yīng)力及孔壓隨時間在不斷變化，為確保土質(zhì)邊坡的變形在物理場作用下產(chǎn)生，而非自重力，故分四步層層遞進(jìn)(如圖4所示)，將物理場接口分為四組(如表2所示)；先計算土質(zhì)邊坡在彈性狀態(tài)下的初始應(yīng)力和孔壓，進(jìn)而將所得結(jié)果導(dǎo)入后計算土質(zhì)邊坡在塑性狀態(tài)下的應(yīng)力和孔壓，此時土質(zhì)邊坡在降雨滲流條件下為動態(tài)變化，即為瞬態(tài)研究，最后求算出土質(zhì)邊坡在強度折減法下的位移變形和塑性應(yīng)變.

圖4 工況分析步驟圖Fig.4 Working condition analysis step diagram

表2 工況分析設(shè)計方案表

表2(續(xù))

2 計算結(jié)果分析

為每一個研究步驟選擇好物理場接口后，通過求解得到各個研究步驟的計算結(jié)果.

2.1 初始應(yīng)力與孔壓

為了確保土質(zhì)邊坡的變形是受物理場作用下所產(chǎn)生的，而非自重力，在進(jìn)行多場耦合計算之前，需要先對土質(zhì)邊坡進(jìn)行初始應(yīng)力與孔壓的計算，如圖5和圖6所示.

從圖5中可以得出初始主應(yīng)力分布整體呈波浪形，由邊坡內(nèi)部向外輻射，層狀邊界清晰，可明顯得出邊坡主應(yīng)力分布走向，層次感較強，并在坡頂和坡腳處出現(xiàn)了剪應(yīng)力集中現(xiàn)象；從圖6中可以得出邊坡受自重影響，其壓力主要集中在邊坡底部，分布均勻，層狀清晰，并在邊坡底部出現(xiàn)壓力最大值，壓力值隨高度的增大而減小.

圖5 初始主應(yīng)力等值云圖 圖6 壓力云圖Fig.5 Contour map of initial principal stress Fig.6 Pressure nephogram

2.2 添加塑性變化

將第一組物理接口的解導(dǎo)入第二組物理接口中，并為模型添加塑性變化，則得到如圖7～10結(jié)果.

圖7 彈塑性狀態(tài)下主應(yīng)力等值云圖 圖8 彈塑性狀態(tài)下壓力云圖Fig.7 Isogram of stress in elastoplastic state Fig.8 Pressure nephogram in elastoplastic state

從圖7～圖10可以得出，將研究1的結(jié)果導(dǎo)入第二組物理接口，并為邊坡土體材料添加土壤塑性，邊坡土體材料由原來的彈性添加塑性變化為彈塑性本構(gòu)模型后，對比圖5和圖7，其應(yīng)力集中表現(xiàn)在彈塑性狀態(tài)下相比較有緩和；對比圖6和圖8，其壓力分布范圍無明顯變化，壓力分布水平稍有點下移；土質(zhì)邊坡有效塑性應(yīng)變?nèi)鐖D9所示，主要集中在坡角處，呈橢球狀；如圖10所示，土質(zhì)邊坡在彈塑性狀態(tài)下的最大位移主要表現(xiàn)在坡頂至下坡面，呈弧形狀，并在距坡面不遠(yuǎn)處有集中體現(xiàn).

圖9 彈塑性狀態(tài)下有效塑性應(yīng)變等值云圖 圖10 彈塑性狀態(tài)下位移云圖Fig.9 Equivalent nephogram of effective plastic strain in elastoplastic state Fig.10 Displacement diagram in elastoplastic state

2.3 降雨-滲流-應(yīng)力場耦合

根據(jù)上一組物理接口得到的結(jié)果導(dǎo)入第三組物理接口中，并從應(yīng)力、飽和度、塑性、位移四個方面選取具有顯著特點的三個不同時間段的特征進(jìn)行分析.

(1) 應(yīng)力特征

圖11中(a)為降雨滲流前主應(yīng)力等值云圖，(b)為降雨滲流12 h后主應(yīng)力等值云圖，(c)為降雨滲流120 h后主應(yīng)力等值云圖.土質(zhì)邊坡在降雨滲流后的主應(yīng)力變化不明顯，層次清晰，尤其在降雨滲流12 h后與之相比無明顯變化，但從圖中可以得出，在坡頂至坡角處有應(yīng)力集中體現(xiàn)，尤其在坡角處應(yīng)力集中明顯，且有緩慢向邊坡內(nèi)部延伸現(xiàn)象.

圖11 主應(yīng)力等值云圖Fig.11 Principal stress contour map

(2) 飽和度特征

采用Richards方程模型，水力傳導(dǎo)率為0.2 m/h，并在右邊界和左下邊界設(shè)置10 m的壓力水頭，壓力水頭與壓力水頭時間變化關(guān)系如圖12所示：其隨時間變化呈圓弧形下降.圖13中(a)為降雨滲流前飽和度云圖，(b)為降雨滲流72 h后飽和度云圖，(c)為降雨滲流120 h后飽和度云圖.隨著時間的推移，由于降雨滲流作用，水流向邊坡體內(nèi)滲流，土壤由非飽和狀態(tài)轉(zhuǎn)變?yōu)榫植匡柡蜖顟B(tài)，或由非飽和到飽和漸變向內(nèi)擴張.

圖12 壓力水頭與壓力水頭時間變化關(guān)系曲線圖Fig.12 Pressure head and pressure head time variation curve

圖13 飽和度云圖Fig.13 Saturation nephogram

(3)塑性特征

圖14中(a)為降雨滲流2.4 h后有效塑性應(yīng)變云圖，(b)為降雨滲流4.8 h后有效塑性應(yīng)變云圖，(c)為降雨滲流120 h后有效塑性應(yīng)變云圖.土質(zhì)邊坡有效塑性應(yīng)變主要表現(xiàn)為由坡角處向邊坡體內(nèi)延伸，即坡角處為最先發(fā)展并開始出現(xiàn)的位置，且以較慢的速度向邊坡內(nèi)部擴大，呈弧形狀，而后受到周圍土體的約束，在降雨滲流4.8 h后無明顯變化.

圖14 有效塑性應(yīng)變云圖Fig.14 Effective plastic strain nephogram

(4)位移特征

由圖15所示，其中：(a)為降雨滲流24 h后位移云圖，(b)為降雨滲流72 h后位移云圖，(c)降雨滲流120 h后位移云圖.土質(zhì)邊坡在降雨-滲流-應(yīng)力場耦合作用下的位移表現(xiàn)為由邊坡內(nèi)部向坡表面逐漸增大，并出現(xiàn)了由上坡面至坡角聚集的現(xiàn)象，在坡角處的位移最大，且最為明顯.針對降雨滲流后的邊坡位移數(shù)值進(jìn)行提取，以坡肩位移和坡角位移為例，并繪制曲線圖，如圖16所示.

圖15 位移云圖Fig.15 Displacement nephogram

圖16 坡肩和坡角位移曲線圖Fig.16 Displacement curve of slope shoulder and slope angle

由圖16可知，坡肩與坡角的位移變化可分為兩個階段，即以降雨滲流1.4 h為分界點，在降雨滲流1.4 h前呈直線上升，降雨滲流1.4 h后，坡角位移曲線呈弧線上升，而坡肩的位移曲線上下浮動較大，呈波形狀.

2.4 結(jié)合有限元強度折減法應(yīng)用

有限元強度折減法是邊坡穩(wěn)定性有限元分析領(lǐng)域的常用方法，其實質(zhì)就是通過對模型進(jìn)行單元劃分，將其劃分為微小單元，從而達(dá)到更為精確的計算，進(jìn)而對材料的黏聚力和內(nèi)摩擦角折減并不斷帶入模型迭代計算，直至不收斂則表示模型已經(jīng)達(dá)到極限，發(fā)生土體失穩(wěn)[17]；此方法能夠很好地解決有限元模擬復(fù)雜的問題，結(jié)果更為直觀，在工程中有著廣泛的應(yīng)用，下面將對土質(zhì)邊坡在降雨-滲流-應(yīng)力場耦合作用下的位移及塑性變化結(jié)果進(jìn)行分析，如圖17所示.

圖17 耦合作用下土質(zhì)邊坡潛在破壞圖Fig.17 Potential failure diagram of soil slope under coupling action

由圖17可知，其中：(a)為有限元強度折減法計算后土質(zhì)邊坡有效塑性應(yīng)變等值云圖，(b)為位移等值云圖，(c)為坡面位移曲線圖.在結(jié)合有限元強度折減法計算后，可得到土質(zhì)邊坡在降雨-滲流-應(yīng)力耦合作用下的潛在破壞面，所表現(xiàn)出的失穩(wěn)破壞類型為滑坡破壞，其主要表現(xiàn)特征是邊坡土體內(nèi)部朝向坡面，由坡角延伸至坡頂處形成滑動弧面，折減后的安全系數(shù)為1.283，最小位移為0.04 m，最大位移為1.45 m，其中以坡角處位移最為敏感.

3 安全性分析

土質(zhì)邊坡的穩(wěn)定性與其安全性密不可分，邊坡穩(wěn)定性是實現(xiàn)邊坡安全性的前提，即要實現(xiàn)邊坡安全性必須要先解決邊坡穩(wěn)定性問題.從安全學(xué)原理角度出發(fā)，將土質(zhì)邊坡發(fā)生的失穩(wěn)破壞假定為事故，并將其劃分為三個階段，即事故孕育階段、事故生長階段、事故損失階段，根據(jù)COMSOL Multiphysics軟件對降雨-滲流-應(yīng)力耦合作用下土質(zhì)邊坡穩(wěn)定性數(shù)值模擬云圖結(jié)果，可將土質(zhì)邊坡在降雨滲流前1.4 h定義為事故孕育階段，此時天然土質(zhì)邊坡為低風(fēng)險狀態(tài)，較為安全；在1.4 h至邊坡發(fā)生失穩(wěn)破壞前定義為事故生長階段，此時天然土質(zhì)邊坡為高風(fēng)險狀態(tài)，即較為不安全；在持續(xù)作用直至失穩(wěn)破壞后，轉(zhuǎn)變?yōu)槭鹿蕮p失階段，此時事故已經(jīng)產(chǎn)生，又將由不安全逐漸轉(zhuǎn)變?yōu)橼呌诎踩珷顟B(tài)，即逐漸向低風(fēng)險狀態(tài)過渡，形成新的力學(xué)平衡狀態(tài).而從有限元強度折減法對安全系數(shù)的定義出發(fā)，土質(zhì)邊坡在降雨滲流120 h后的安全系數(shù)為1.283，按照安全系數(shù)大于1時為穩(wěn)定狀態(tài)，小于1時為不穩(wěn)定狀態(tài)的定義，理論上來說此時的邊坡處于一個較為安全的狀態(tài)，但隨著降雨-滲流-應(yīng)力耦合的不斷作用，土質(zhì)邊坡持續(xù)受作用5天后的位移變形是直線上升的，發(fā)展速度快，且依據(jù)在現(xiàn)實生活中的實際情況，安全系數(shù)值只能作為一個參考，不完全符合實際，故而將降雨-滲流-應(yīng)力耦合作用下的土質(zhì)邊坡定義為不安全狀態(tài).

4 結(jié)論

以天然土質(zhì)邊坡為研究對象，運用COMSOL Multiphysics數(shù)值模擬軟件，通過構(gòu)建幾何二維平面有限元土質(zhì)邊坡模型，揭示了天然土質(zhì)邊坡在多物理場耦合作用下的演化規(guī)律，獲取了在降雨-滲流-應(yīng)力場耦合作用下天然土質(zhì)邊坡的安全性，結(jié)果表明：隨著降雨滲流的不斷進(jìn)行，在降雨滲流前1.4 h，天然土質(zhì)邊坡的位移變形呈直線上升，此時天然土質(zhì)邊坡為低風(fēng)險狀態(tài)，處于事故孕育階段，較為安全；在1.4 h至24 h為緩慢變形時間段，24 h后變形程度逐漸地加大，其中以坡角處變形最為明顯，此時天然土質(zhì)邊坡為高風(fēng)險狀態(tài)，處于事故生長階段，即較為不安全；在持續(xù)作用直至失穩(wěn)破壞后，轉(zhuǎn)變?yōu)槭鹿蕮p失階段.通過有限元強度折減法對有限元模型進(jìn)行了計算分析，在降雨-滲流-應(yīng)力場耦合作用下，天然土質(zhì)邊坡潛在破壞面為滑動破壞，其主要表現(xiàn)特征是邊坡土體內(nèi)部朝向坡面，由坡角延伸至坡頂處形成滑動弧面，最小位移為0.04 m，最大位移為1.45 m，安全系數(shù)為1.283.