巖溶區橋梁雙樁基礎有限元極限分析*

肖 堯 趙明華 張 銳 徐卓君 趙 衡

(①湖南大學巖土工程研究所 長沙 410082)(②中南大學土木工程學院 長沙 410075)

0 引 言

在巖溶區進行橋梁設計與施工時，樁基的極限承載力計算關系到工程設計的安全與經濟。目前，眾多學者對該問題進行了研究，研究手段主要包括試驗研究、理論研究和數值分析等。試驗研究方面，王革力(2002)、劉鐵雄(2003)將溶洞視為梁板模型，探討了溶洞頂板厚度、跨度對基樁極限承載力的影響； 而張慧樂等(2013)則進一步研究了溶洞位置、大小和形狀等因素對基樁極限承載力的影響，江松等(2017)對巖溶樁基動力效應進行了試驗研究。理論研究方面，劉之葵等(2003)基于彈性力學平面應變假定，通過判斷圓形溶洞邊界上關鍵點是否發生破壞，進而反算得到基樁極限承載力； 趙明華等(2004， 2016a)提出了溶洞頂板抗沖切、抗剪切、抗彎拉的驗算方法，并以三者最小值作為基樁極限承載力； 韓紅艷等(2012)對巖溶路基溶洞頂板穩定性進行了研究。此外，雷勇等(2014)、趙明華等(2017， 2018)采用極限分析的方法得到了基樁承載力計算公式。數值分析方面，黎斌等(2002)、陽軍生等(2005)和黃明等(2017)采用有限元的方法，對巖溶區基樁極限承載力進行了系統研究，得出了一些有益的結論。以上研究取得了一系列的成果，但都是針對單樁方面的研究，對巖溶區橋梁雙樁基礎極限承載力的研究似尚未見報道。在實際工程中，雙樁基礎比單樁的情況更為常見，且受力更為復雜，其與單樁承載機理主要區別在于，雙樁基礎可能存在荷載疊加的效應，從而使整體承載力降低。因此，有必要針對該問題作進一步研究。

本文的目的在于，結合有限元和極限分析的優點，通過MATLAB編程求解巖溶區橋梁雙樁基礎極限承載力。該方法相較傳統有限元和傳統極限分析的優勢在于：(1)不需要人為構造可靜應力場(下限解)及機動許可的運動場(上限解)； (2)將直接給出極限承載力的下限和上限，不需要通過荷載-位移曲線來確定極限承載力，有效地克服了讀數時的誤差； (3)結果容易收斂，計算效率較高。為此，下文將首先對有限元極限分析的基本原理作簡要介紹； 其次，為了描述巖體非線性特點，采用Hoek-Brown準則，并簡要介紹將Hoek-Brown準則嵌入計算程序的方法； 然后，詳細討論溶洞大小、位置、樁距等因素對極限承載力的影響； 最后，計算無溶洞條件下單樁極限承載力，并與已有成果進行對比，驗證本文方法的正確性。

1 有限元極限分析方法簡介

有限元極限分析方法被廣泛應用于巖土工程穩定性分析中(Sloan， 2013)，該方法根據極限分析基本理論，結合有限元的方法，利用計算機求解得到嚴格的上、下限解。該方法不需要人為構造可靜應力場、機動許可的速度場，對巖溶區橋梁雙樁基礎承載力的研究尤為適合。為此，下文將對其基本原理作簡要介紹。

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

圖 2 基于退化單元的速度間斷線示意圖Fig. 2 Sketch of velocity discontinuity based on degenerated triangular elements

為了克服線性三角形單元難以模擬復雜速度場分布的缺點，本文采用Krabbenh?ft et al.(2005)提出的基于“退化”單元的速度間斷線設置方法。如圖 2所示，兩個厚度為0的單元構成了“退化”單元的速度間斷線，其有一對節點坐標重合。

1.2 數學規劃模型建立

單元離散后，根據上、下限定理，建立節點應力和節點速度的約束方程，以總的內能耗散或外力荷載作為目標函數，并得到相應的數學規劃模型。該過程在文獻(Lyamin et al.， 2002, 2005)中已有詳細介紹，本文不再贅述，下面將直接給出上、下限分析數學規劃模型的具體形式。

上限分析數學規劃模型具體形式為(Lyamin et al.， 2002)：

MinimizeQ=σTBu-cTu，

(1)

Subject toAu=b，

(2)

(3)

(4)

fi(σ)≤0，j∈Jσ，

(5)

(6)

u∈Rnu，σ∈Rnσ，λ∈RE

下限分析數學規劃模型具體形式為(Lyamin et al.， 2005)：

MaximizeQ=cTx，

(7)

Subject toAx=b，

(8)

fi(x)≤0，j∈J

(9)

x∈Rn

(10)

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

1.3 有限元極限分析的計算機實現

有限元上、下限分析的求解過程如圖 3所示，本文以MATLAB為平臺編制了相關計算機程序，計算網格的劃分，對優化模型建構和求解，并調用Tecplot360軟件實現數據信息可視化。具體過程的論述可參考文獻(張銳， 2015； 趙明華等， 2015， 2016b)。

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

1.4 基于Hoek-Brown準則的有限元極限分析

Hoek-Brown準則能較好地描述巖體非線性特征，在實際工程中得到了廣泛應用，其最新的表達形式如下(Hoek et al.， 2007)：

σ′1=σ′3+σc(mbσ′3/σc+s)a

(11)

式中，σ′1和σ′3分別為最大、最小主應力；σc為巖石單軸飽和抗壓強度；mb，s和a為與地質強度指標GSI有關的參數，其表達式為，

mb=miexp[(GSI-100)/(28-14D)]

(12)

s=exp[ (GSI-100)/(9-3D)]

(13)

a=1/2+(e-GSI152212/e-20/3)/6

(14)

式中，mi為與巖石完整程度有關的參數；D為擾動系數，沒有擾動取0，完全擾動取1，本文中D=0。

針對本文所研究的問題，也采用Hoek-Brown準則，由于Hoek-Brown準則不同于其他線性屈服準則(如Mohr-Coulomb準則、Tresca準則等)，因此有必要對屈服函數迭代計算過程中所用的方法進行介紹。

在優化模型的求解過程中，除了需要計算當前迭代點的屈服函數值，還必須獲得屈服函數對應力變量的一階、二階導數，即屈服函數的梯度向量和Hessian矩陣，具體過程如下：

為了便于推導，將Hoek-Brown準則轉化為應力不變量表達的形式，

(15)

其中，I1為第一應力不變量，J2為第二偏應力不變量，θ為第三偏應力不變量J3相關的應力Lode角，參數β、χ及函數g(θ)、h(θ)的表達式如下：

g(θ)=-2cos(θ)

(16)

(17)

(18)

(19)

由于Hoek-Brown準則的應力空間包絡面上存在奇異點，導致不能被求導，為了讓Hoek-Brown準則適用于非線性規劃算法，需對其進行“光滑化”處理。采用一個微小項ε對屈服函數進行“雙曲型近似”處理，即采用式(20)對J2進行代替。

(20)

ε是一個與材料抗拉強度相關的常數，其值可由式(21)確定。

ε=min(δ，μρ|ρg(0)+(ρh(0)+χ)α=0)

(21)

其中，μ=10-1，δ=10-6，ρ為屈服面與縱軸對應的交點。

則雙曲型近似Hoek-Brown準則可表示為：

(22)

根據式(22)求一階、二階偏導便可求得屈服函數的梯度向量和Hessian矩陣，由于篇幅所限，具體過程請參考文獻(張銳， 2015)，本文不再贅述。

2 問題描述

2.1 計算模型及假定

Serrano et al.(2014)對嵌巖樁樁端極限承載力的研究表明：三維模型比二維模型承載力大約高1.3倍。此外，廖麗萍等(2010)研究表明：橢球形溶洞比橢圓形溶洞更穩定。以上研究表明，將巖溶區樁基承載力問題簡化為平面模型是偏于安全，是有利于工程的。為了建立模型的同時突出關鍵因素的影響，本文將問題簡化為二維平面模型，如圖 4所示，并假定：(1)上部荷載全部由樁端持力巖層承擔； (2)雙樁基礎由橫系梁連接，上部荷載由雙樁基礎共同承擔，且不考慮樁體自身變形的影響； (3)溶洞截面形狀為圓形，周邊巖層為均質材料，且符合修正的Hoek-Brown準則。

圖 4 計算模型Fig. 4 Calculation model

圖 4中：q為上部荷載，d為樁徑，h1、hk、…hn分別為第1、k、…n層土的高度，γ為巖石的重度，γ1、γk、…γn分別為第1、k、…n層土的重度，hr為嵌巖深度，X、Y分別為溶洞中心與左樁樁端中心的水平距離和垂直距離，R為溶洞半徑，L為左樁與右樁中軸線的水平距離。

圖 5 網格劃分示意圖Fig. 5 Sketch of meshing

2.2 網格劃分

數值模型與網格劃分如圖 5所示，分析域寬為40id，高為25id，模型左右邊界及樁側進行法向約束，底部進行完全約束，網格采用自適應劃分技術(張銳， 2015； 趙明華等， 2016)，單元總數為6000個，進行5次網格優化，初始單元數取1000，均布極限荷載qu作用在橫梁上，上覆土層以均布荷載qs的形式作用于巖層上，其表達式如式(23)所示。

(23)

樁側與巖層的接觸面光滑，樁端與巖層接觸面完全粗糙。樁與橫梁作為整體，視為無重度的剛體，巖層為均質材料，符合修正的Hoek-Brown準則，巖體的彈性模量E和泊松比ν分別根據式(24)、式(25)確定(Hoek et al.， 2006； Vasrhelyi， 2009)。

(24)

ν=-0.002GSI-0.003mi+0.457

(25)

2.3 評價指標

為了評價溶洞對雙樁基礎極限承載力的影響，定義一個無量綱參數Nσ，表達式如下：

Nσ=qu/σc

(26)

表 1 誤差分析Table 1 Error analysis

由表 1 可知，上限解與下限解的誤差在10%以內，以上限解與下限解的平均值作為設計指標，則其與真實解的相對誤差在5%以內，為了后續便于分析，本文參數Nσ取上、下限解的平均值。

3 結果與討論

巖溶區橋梁雙樁基礎整體的極限承載力影響因素主要包括：(1)上覆土層荷載qs； (2)嵌巖深度hr； (3)溶洞半徑R； (4)左樁樁端中心與溶洞中心的水平距離X； (5)左樁樁端中心與溶洞中心的垂直距離Y； (6)樁的間距L； (7)巖石物理力學參數，GSI，mi，γ。下面將分別討論各因素對嵌巖樁樁端極限承載力的影響

3.1 上覆土層的影響

圖 6 qs對Nσ的影響Fig. 6 Effect of qson Nσ

由圖 6可知，Nσ隨著qs的增大而增大，且X越大，增長的幅度越大。當X=0時，qs對Nσ的影響基本可忽略不計。

3.2 嵌巖深度的影響

圖 7 hr對Nσ的影響Fig. 7 Effect of hron Nσ

由圖 7可知，嵌巖深度越大，橋梁雙樁基礎的承載力會得到一定提高，文獻(趙明華等， 2016)也得出了相類似的結論。

3.3 溶洞半徑R的影響

由圖 8可知，Nσ隨著溶洞半徑R的增大而不斷減小，且減小的幅度變緩，直至趨于某一穩定值。在文獻(張慧樂等， 2013)中也得出了相類似的結論。

圖 8 R對Nσ的影響Fig. 8 Effect of R on Nσ

3.4 水平距離X的影響

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

圖 10 不同X條件下雙樁與單樁承載力系數對比Fig. 10 Comparison of bearing capacity factor for double-piles and single pile with different X

為了分析雙樁與單樁承載性能的區別，將雙樁所得結果平均分配給左樁和右樁，保持所有條件不變，取Y=2id，計算得到僅有左樁、右樁時的極限承載力系數，對比結果如圖 10所示。根據圖 10，當X≥0時，雙樁基礎的承載力大致等于僅有左樁、右樁時承載力的平均值； 當X=-1.5id時，雙樁基礎承載力大于僅有單樁時的承載力，主要原因在于發生了圖 16a所示的整體剪切破壞。

3.5 垂直距離Y的影響

與圖 10類似的，取X=0，得到不同Y條件下雙樁與單樁承載力系數對比結果，如圖 12所示。從圖 12中可以看出，Y/d≥5時，雙樁基礎承載力系數大約為僅有單樁的90%。產生該現象的原因可能是雙樁基礎存在荷載疊加的效應，因此雙樁基礎相對于單樁的承載力降低了。綜合圖 10、圖12的結果，建議工程設計時，雙樁基礎的承載力計算由離溶洞最近的單樁承載力控制，并乘以0.9的折減系數。

圖 11 Y對Nσ的影響Fig. 11 Effect of Y on Nσ

圖 12 不同Y條件下雙樁與單樁承載力系數對比Fig. 12 Comparison of bearing capacity factor for double-piles and single pile with different Y

3.6 樁間距L的影響

在X確定的情況下，樁間距L的影響主要體現在右樁與溶洞的位置關系，當右樁離溶洞較遠時，雙樁基礎的承載力也會較高。雙樁基礎與溶洞的位置關系在前面已論述，在此不再贅述。

3.7 地質強度指標GSI的影響

圖 13 GSI對Nσ的影響Fig. 13 Effect of GSI on Nσ

表 2 γ對Nσ的影響Table 2 Effect of γ on Nσ

3.8 參數mi的影響

圖 14 qs對Nσ的影響Fig. 14 Effect of qs on Nσ

3.9 巖體自重γ的影響

圖 15 不同R的溶洞極限破壞模式Fig. 15 Failure patterns of the cave with different R

圖 16 不同X的溶洞極限破壞模式Fig. 16 Failure patterns of the cave with different X

圖 17 不同Y的溶洞極限破壞模式Fig. 17 Failure patterns of the cave with different Y

4 破壞模式

根據計算的結果，溶洞大小和位置是影響極限破壞模式的關鍵因素，將不同R、X、Y條件下溶洞的極限破壞模式總結出來，如圖 15～圖 17所示。圖 15中，X=0、Y=2id，當R=1id時，左樁下方的溶洞頂板發生沖切破壞，同時右樁出現地基破壞模式； 隨著R的增大，破壞模式由左樁控制，右樁對破壞模式影響不大，圖 15c、圖15d所示； 當R≥4時，發生整體剪切破壞，破壞面貫穿左、右樁的外側。

圖 16給出了不同X條件下溶洞的破壞模式，其中，R=2id、Y=2id，當-1.5d≤X≤0，發生整體剪切破壞； 隨著X的增大，破壞模式逐漸變為由左樁控制的沖切破壞，同時右樁底部發生地基破壞模式。

圖 17給出了不同Y條件下溶洞的破壞模式，其中，R=2id、X=0，溶洞與樁端豎直距離似乎對破壞模式影響不大，以整體剪切破壞為主，破壞面由左、右兩樁樁側延伸至溶洞表面，這可以用來解釋Nσ與Y大致呈線性關系。此外，值得注意的是，當Y=1id時，溶洞頂部會出現冒落區。

5 結果驗證

表 3 無溶洞條件下單樁樁端承載力系數NσTable 3 Bearing capacity factor Nσof pile tip without void

6 結 論

(1)考慮實際工程中巖溶區橋梁雙樁基礎的承載特點，根據極限分析上、下限定理，利用MATLAB平臺編制了有限元極限分析計算程序，用Hoek-Brown準則來描述巖體的非線性特征，并通過“雙曲線近似”處理，將其嵌入計算程序中。

(2)上覆土層荷載和嵌巖深度越大，雙樁基礎的承載力越大； 溶洞半徑越大，極限承載力越小，逐漸趨向于0； 極限承載力隨GSI的增大非線性增大，與mi大致呈線性關系； 當樁與溶洞距離較遠時，巖石重度越大，承載力越高，但當樁與溶洞距離較近時，巖石重度對承載力影響可忽略不計。

(3)左樁與溶洞的水平距離X增大，承載力先增大后減小，在X=0時取最小值； 承載力隨垂直距離Y的增大而大致線性增大。通過與單樁承載力的比較，在工程設計時建議，雙樁基礎的承載力計算由離溶洞最近的單樁承載力控制，并乘以0.9的折減系數。

(4)巖溶區橋梁雙樁基礎的極限破壞模式主要有3種： ①左樁下方的溶洞頂板發生沖切破壞，同時右樁出現地基破壞模式； ②由左樁控制的沖切破壞模式； ③整體剪切破壞，破壞面由左、右樁的外側貫穿至溶洞表面。

(5)本文所有計算結果都進行了無量綱化處理，并通過計算無溶洞時樁端極限承載力，與理論方法、FLAC計算結果對比，驗證了本文所提方法的正確性，可為同類工程提供參考。