基于有限元極限分析的下伏溶洞路堤承載力計算

陳言章，趙明華

（1.廣東省建筑設計研究院有限公司，廣東 廣州 510010；2.湖南大學 巖土工程研究所，湖南 長沙 410082）

0 引 言

我國巖溶地貌分布廣泛，大量在建的山區公路不可避免地會穿越巖溶發育區，由于巖溶地質的復雜多樣性，公路路堤容易發生溶洞的塌陷破壞。因此，如何確定下伏溶洞路堤承載力具有重要的工程價值。

目前，眾多學者對該課題開展了相關研究，可歸結為試驗研究、理論研究和數值分析3個方面。試驗研究方面，Al-Tabbaa等[1]基于小比例模型試驗對下伏空洞地基的穩定性進行了研究；Kiyosumi等[2]對堅硬地層中含多個矩形空洞的條形基礎承載力開展了離心機試驗；劉庭金等[3]基于室內模型試驗研究了含矩形空洞地基的破壞過程。當通過試驗手段獲得數據后，有必要對模型做適當的假定和簡化，進一步對數據進行分析，提出符合試驗規律的理論方法，并結合數值分析方法進行驗證。理論研究方面，Wang等[4]假定作用在圓形空洞上條形基礎的破壞模式，采用極限分析的方法得到了基礎的極限承載力公式；劉輝等[5]利用極限分析上限法，建立了與假定破壞模式對應的速度場，求得了空洞上方條形基礎極限承載力；劉之葵等[6]基于彈性理論求得溶洞周圍巖體的應力狀態，并結合格里菲斯強度準則，分析了含溶洞巖石地基的穩定性。數值分析方面，Azam等[7]利用二維有限元軟件對含空洞時的地基承載力進行了數值計算，討論了空洞頂板厚度和空洞位置等因素對地基承載力的影響；彭芳樂等[8]利用PLAXIS分析了單溶洞的存在對基礎承載力和沉降的影響，對其發生機制進行了研究；在此基礎上，Kiyosumi等[9]進一步考慮了條形基礎作用在含多個溶洞地層時的極限承載力。陽軍生等[10]采用ABAQUS軟件對巖溶地基圓形基礎作用下溶洞頂板的穩定性進行了分析，探討了極限承載力與頂板跨度、頂板厚度等影響因素的關系。以上有限元方法雖取得了一定的成果，但采用傳統有限元方法研究承載力問題尚存在計算結果容易發散、效率較低等問題，有限元極限分析法通過降低強度或增加荷載使巖體最終達到極限破壞狀態，自動生成破壞面，同時可得到極限荷載，克服了傳統有限元方法需根據位移-荷載曲線確定極限承載力的不足。其基本思想是利用有限元將應力場離散化，然后在離散的應力場內按照上、下限定理的相關要求構建相應的數學規劃模型，最后選用合適的數學規劃算法求解該模型，搜索出極限應力場和上、下限荷載[11-12]。該方法不需要人為構造靜力容許的應力場、機動許可的速度場，對巖溶區路堤極限承載力的研究較為適合。

鑒于此，本文采用修正的Hoek-Brown準則對下伏溶洞路堤承載力進行計算，重點分析溶洞形狀、分布形態、溶洞位置、巖石物理力學參數對路堤承載力的影響，并給出極限破壞模式，最后，將作用在巖體上的條形基礎承載力與本文計算結果進行對比，以驗證程序的正確性及路堤承載力結果的可靠性。

1 計算原理與方法

本文將采用有限元上、下限分析方法研究下伏溶洞路堤承載力問題，下面簡單介紹其基本原理和計算程序求解過程。

1.1 單元離散

如圖1所示，采用三角形單元對計算域進行離散，其中，下限單元每個節點i有3個未知應力分量（σxi,σyi,τxyi），每個單元共2個未知體積力分量，即he=（hx,hy）；上限單元每個節點j有2個未知速度分量（uj,vj），每個單元共3個未知應力分量，即σe=（σx,σy,τxy）。

圖1 三角形單元示意圖Fig.1 Sketch of triangular element

1.2 非線性規劃模型的建立

單元離散后，根據上、下限定理，建立節點應力和節點速度的約束方程，確定合理的優化目標函數。下面將分別給出上、下限分析數學規劃模型的具體形式。

上限分析數學規劃模型具體形式為[11]：

下限分析數學規劃模型具體形式為[12]：

式中：x為全局節點應力列向量；fj(x)為屈服準則或其他條件產生的不等式約束，其他符號意義同前。

1.3 有限元上、下限分析的計算機實現

有限元上、下限分析的求解過程如圖2所示，本文以 MATLAB為平臺編制相關計算機程序，計算網格的劃分，對優化模型建構和求解，并調用Tecplot360軟件實現數據信息的可視化。具體過程的論述可參考文獻[13-15]。

圖2 有限元上、下限分析的計算機實現Fig.2 Computer implementation of upper and low bound finite element method

2 問題的描述

2.1 基本假定

廖麗萍等[16]對地基中橢球型空洞穩定性進行了分析，分析結果表明：在遠場應力狀態相同的條件下，橢球洞比橢圓孔更穩定。戴自航等[17]基于ABAQUS對室內試驗進行了數值模擬，研究表明：相同條件下，三維橢球狀溶洞的穩定性要大于二維橢圓形溶洞的穩定性，而后者的穩定性又高于同等條件下矩形溶洞的穩定性。鑒于此，本文將公路路堤下伏溶洞頂板的承載力問題簡化為二維平面模型，計算模型如圖3所示，并假定：（1）路堤荷載視為均布荷載，土層按層狀分布；（2）溶洞截面形狀為矩形，周邊巖層為均質材料，且滿足修正的Hoek-Brown準則。

圖3 計算模型Fig.3 Calculation model

圖3中：q為路堤荷載集度；d為路堤荷載的寬度；h1、hk、… hn分別為第1、k、…n層土的高度；γ為巖石的重度；γ1、γk、…γn分別為第1、k、…n層土的重度；X為路堤荷載中心與溶洞中心的水平距離；Y為巖層頂面與溶洞中心的垂直距離；W、H分別為溶洞跨度和高度；θ為溶洞的旋轉角度，以逆時針方向為正。

2.2 數值分析模型

本文將路堤荷載作為極限荷載qu，上覆土層以均布荷載qs的形式作用于巖層上，其表達式如式（11）所示。

數值模型與網格劃分如圖4所示，路堤荷載寬度d=10 m，為減小邊界條件的影響，取分析域寬為4 d，高為3 d，約束模型左右邊界的法向位移，模型底部則進行完全約束，網格采用自適應劃分技術[14]，單元總數為14 000個，以剪切耗散作為控制變量，初始單元數取1 000個，對網格進行3次自適應迭代。由圖4可見，在能量耗散劇烈的區域，單元劃分的較多，而在能量耗散較小的區域單元分布較為稀疏，在單元總數一定的條件下，使得計算域內塑性變形越劇烈的區域可獲得越大的變形自由度，從而有效地降低數值離散誤差。

圖4 數值模型及網格劃分Fig.4 Numerical model and mesh

2.3 材料參數及評價指標

為了描述巖石固有的非線性破壞特點，采用修正的 Hoek-Brown準則[18]，其表達式如式（12）所示。

式中：mi為與巖石完整程度有關的常數；D為擾動參數，現場無擾動巖體為0，而完全擾動巖體為1，本文取D=0。

由于Hoek-Brown準則是非線性準則，在屈服函數迭代計算過程中，所用的方法不同于其它線性屈服準則，在優化模型的求解過程中，除了需要計算當前迭代點的屈服函數值，還必須獲得屈服函數對應力變量的一階、二階導數，即屈服函數的梯度向量和Hessian矩陣，具體過程可參考文獻[15]。

采用Hoek-Brown準則的關鍵在于確定巖體參數GSI和巖石參數mi的范圍。對于巖體參數GSI，Sonmez等[19]考慮不連續面的分布率、風化程度、粗糙度和填充物性質等影響，提出了詳細的GSI定量評價方法。當GSI取值在40～90之間時，能基本滿足實際工程的要求。對于巖石參數 mi，Hoek等[20]結合室內試驗和工程經驗，提出了較為全面的mi取值方法。鑒于巖溶主要成分為碳酸鹽，其巖石參數mi的取值見表1，根據表1可知巖溶區的巖石參數mi變化范圍為6～15。

表1 碳酸鹽類巖石參數mi的取值Table 1 Values of parameter mifor carbonate rock

Hoek 等[21]和 Vá sá rhelyi[22]基于參數 GSI和 mi分別提出了計算巖石彈性模量E和泊松比v的經驗公式，其表達式為：

3 路堤極限承載力

影響路堤極限承載力的因素主要有：（1）上覆土層荷載qs；（2）溶洞寬度W；（3）溶洞高度H；（4）溶洞旋轉角度 θ；（5）路堤荷載中心與溶洞中心的水平距離X；（6）巖層頂面與溶洞中心的垂直距離Y；（7）巖石物理力學參數：GSI，mi等。由于本文采用了網格自適應劃分技術，所得下限解和上限解相對誤差在10%以內，為了使圖表便于分析，本文取下限解和上限解的平均值作為路堤的極限承載力。下面將分別討論各因素對路堤極限承載力的影響。

3.1 上覆土層荷載qs對極限承載力qu的影響

取H=0.1 d、θ=0°、X=0、Y=0.2 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。qs與 qu的關系如圖5所示。

由圖5可知，當W≤0.1 d時，qu隨著qs增大而增大，當qs＞100 kPa后，qu的增長幅度逐漸趨于緩慢；當W≥0.2 d時，qu基本保持不變。因此，當溶洞跨度較大時，上覆土層荷載qs對路堤極限承載力影響不大，在實際設計時可不考慮。值得注意的是，在無溶洞時，qu與qs大致呈線性關系。

圖5 qs對qu的影響Fig.5 Effect of qs on qu

3.2 溶洞跨度W對極限承載力的影響

取 qs=100 kPa、H=0.1 d、θ=0°、X=0、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。W與 qu的關系如圖6所示。

圖6 W對qu的影響Fig.6 Effect of W on qu

由圖6可知，當Y＜0.4 d時，qu隨著W的增大而減小，并趨向于0，qu減小的幅度也趨于緩慢；特別地，當Y=0.1 d時，qu接近于0，為保證溶洞頂板的穩定性，應保證溶洞頂板的厚度不能過??；當Y≥0.4 d時，qu隨著W的增大大致呈線性減小。

3.3 溶洞高度H對極限承載力的影響

取qs=100 kPa、Y-H/2（即溶洞頂板厚度）=0.1 d、θ=0°、X=0、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。H與qu的關系如圖7所示。

由圖7可知，當W≤0.1 d時，qu隨著H的增大呈非線性減??；當W≥0.2 d時，qu隨著H的增大基本保持不變，說明當溶洞跨度較大時，溶洞高度對路堤極限承載力結果影響不大，該結論與文獻[10]所得結論一致。

圖7 H對qu的影響Fig.7 Effect of H on qu

3.4 旋轉角度θ對極限承載力的影響

取 qs=100 kPa、X=0、Y=0.2 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。θ與 qu的關系如圖8 所示。

圖8 θ對qu的影響Fig.8 Effect of θ on qu

由圖8可知，當W≤0.1 d時，θ<30°，qu隨θ的增大而減?。沪龋?0°，qu隨θ的增大，先減小后增大；θ=30°，qu取最小值。當W≥0.15 d時，θ≤60°，qu基本不變；θ＞60°，qu增大或減小，要根據溶洞跨度確定。由以上分析可知，當溶洞跨度較大時，且θ≤60°時，qu隨θ的變化幅度并不大，在實際工程設計時可不考慮θ的影響。

3.5 水平距離X對極限承載力的影響

取qs=100 kPa、W=0.2 d、H=0.1 d、θ=0°、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。X與qu的關系如圖9所示。

由圖9可知，由于路堤荷載寬度遠大于溶洞跨度，因此，qu隨著X的增大先保持不變而后逐漸增大，且Y值越小，增大的幅度越大。根據圖5，沒有溶洞存在時，qu=35 MPa，當X≥0.9 d時，qu可達到無溶洞條件下 qu的 80%以上。由以上分析可知，實際工程中，軸對稱情況為最不利情況，但Y值較小時，X值的改變對qu的影響比較明顯，在設計時應當考慮其有利的影響。

圖9 X對qu的影響Fig.9 Effect of X on qu

3.6 垂直距離Y對極限承載力的影響

取 qs=200 kPa、hr=1 d、W=2 d、H=2 d、θ=0°、GSI=50、mi=7、γ=27 kN/m3、σci=40 MPa。Y 與 qu的關系如圖10所示。

圖10 Y對qu的影響Fig.10 Effect of Y on qu

由圖10可知，當W≤0.3 d時，qu隨Y的增大而增加，Y＜0.35 d時增加速度較快，Y＞0.35 d時增加速度逐漸變緩。當W≥0.4 d時，qu隨Y的增大而增加，Y＜0.25 d時增加速度較快，Y＞0.25 d時近似線性增加。由以上分析可知巖層頂面與溶洞中心的垂直距離對路堤承載力的影響較為明顯，在實際工程中，應注意保證溶洞頂板具有足夠的安全厚度。

3.7 巖石物理力學參數對極限承載力的影響

取 qs=100 kPa、W=0.2 d、H=0.1 d、θ=0°、X=0、Y=0.2 d。巖石各物理力學參數對qu的影響如表2、表3和圖11所示。

由表2、表3可知，有溶洞與無溶洞條件下，γ對qu的影響均不大，qu隨σci的增大而增大。因此，在設計時可不考慮 γ的影響。由圖 11可知，qu隨GSI的增大而增大，且增大的速度隨GSI增大而增大，此外，mi越大，qu越大。

表2 不同 γ條件下 qu的大小（GSI=50、mi=10、σci=40 MPa）Table 2 Ultimate bearing capacity quwith different γ(GSI=50、mi=10、σci=40 MPa)

表3 不同σci條件下qu的大?。℅SI=50、mi=10、γ=25 kN/m3）Table 3 Ultimate bearing capacity quwith different σci(GSI=50、mi=10、γ=25 kN/m3)

圖11 GSI對 qu的影響 (γ=25 kN/m3、σci=40 MPa)Fig.11 Effect of GSI on qu(γ=25 kN/m3、σci=40 MPa)

4 極限破壞模式

溶洞的存在是影響巖層極限破壞模式的最主要因素，因此下面將主要從溶洞各參數對極限破壞模式的影響展開討論。

4.1 溶洞寬度W對破壞模式的影響

取qs=100 kPa、H=0.1 d、θ=0°、X=0、Y=0.15 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同W條件下的極限破壞模式，如圖12所示。當W≤0.1 d時，溶洞側壁發生破壞，頂板發生冒落破壞；當W=0.2 d時，溶洞頂板發生沖切破壞，同時伴有冒落破壞；當W≥0.4 d時，溶洞頂板發生沖切破壞?？梢?，隨著溶洞寬度的增加，溶洞的破壞模式從側壁發生破壞向頂板發生沖切破壞轉變，轉變過程中伴有頂板的冒落破壞。

圖12 不同W條件下極限破壞模式Fig.12 Failure mechanisms with different values of W

4.2 溶洞跨度H對破壞模式的影響

取qs=100 kPa、Y-H/2（即溶洞頂板厚度）=0.2 d、θ=0°、X=0、W=0.2 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同H條件下的極限破壞模式，如圖13所示。溶洞頂板發生沖切破壞，H對破壞模式的類型并無影響，這也印證了H對路堤極限承載力結果影響不大的結論。

圖13 不同H條件下極限破壞模式Fig.13 Failure mechanisms with different values of H

4.3 旋轉角度θ對破壞模式的影響

取qs=100 kPa、X=0、Y=0.3 d、W=0.3 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同θ條件下的極限破壞模式，如圖14所示。破壞面由側壁延伸至巖層頂面，且隨著θ的變化而旋轉。

圖14 不同θ條件下極限破壞模式Fig.14 Failure mechanisms with different values of θ

4.4 水平距離X對破壞模式的影響

取qs=100 kPa、Y=0.15 d、W=0.2 d、θ=0°、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同X條件下的極限破壞模式，如圖15所示。當X≤6時，溶洞頂板發生沖切破壞；當X≥8時，由于路堤荷載遠大于上覆土層荷載，溶洞頂板和靠近路堤荷載的側壁將發生破壞。

圖15 不同X條件下極限破壞模式Fig.15 Failure mechanisms with different values of X

4.5 垂直距離Y對破壞模式的影響

取 qs=100 kPa、X=0、W=0.2 d、θ=0°、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同Y條件下的極限破壞模式，如圖16所示。當Y≤1.5 d時，溶洞頂板發生沖切破壞；當2 d≤Y＜3 d時，溶洞側壁發生破壞，破壞面延伸到巖層頂面；當Y≥3 d時，溶洞側壁發生破壞。可見，隨著巖層頂面至溶洞中心垂直距離的增加，溶洞的破壞模式從頂板發生沖切破壞向側壁發生破壞轉變。

圖16 不同Y條件下極限破壞模式Fig.16 Failure mechanisms with different values of Y

5 結果驗證

Serrano等[23]假定巖體的重度 γ=0，基于修正的Hoek-Brown準則提出了一種計算條形基礎極限承載力的方法。在此基礎上，Merifield等[24]利用極限理論對巖石地基的極限承載力進行了數值模擬。為了驗證本文方法的正確性，考慮無溶洞條件下條形基礎作用在巖層上的極限承載力，將本文計算結果與文獻[23-24]的研究成果進行對比，對比結果如表4所示。

由表4可知，本文所得條形基礎極限承載力與文獻[23-24]的結果基本一致，誤差在 5%以內，因此本文所提方法較為可靠。

表4 無溶洞條件條形基礎極限承載力Table 4 Ultimate bearing capacity of strip footing with no cave

6 結 論

（1）根據極限分析上、下限定理，利用MATLAB平臺編制了有限元極限分析計算程序，并基于Hoek-Brown準則計算了下伏溶洞路堤極限承載力。

（2）巖層上覆荷載、溶洞高度、巖石重度對路堤極限承載力影響不大；當溶洞旋轉角度小于60°時，可不考慮旋轉角度的影響。

（3）路堤極限承載力隨溶洞跨度W的增大而減小，隨路堤荷載中心與溶洞中心水平距離X、巖層頂面距溶洞中心垂直距離Y、巖石單軸抗壓強度σci地質強度指標GSI的增大而增大。

（4）極限破壞模式主要有溶洞頂板的沖切破壞，溶洞側壁發生破壞，溶洞頂板沖切和側壁的聯合破壞，溶洞頂板冒落和側壁的聯合破壞等。

（5）通過計算無溶洞時路堤的極限承載力，與已有研究成果進行了對比，驗證了本文所提方法的正確性。