巖溶區(qū)橋梁雙樁基礎(chǔ)有限元極限分析*

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

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

0 引 言

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

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

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

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

1.1 單元離散

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

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

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

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

1.2 數(shù)學(xué)規(guī)劃模型建立

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

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

MinimizeQ=σTBu-cTu，

(1)

Subject toAu=b，

(2)

(3)

(4)

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

(5)

(6)

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

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

MaximizeQ=cTx，

(7)

Subject toAx=b，

(8)

fi(x)≤0，j∈J

(9)

x∈Rn

(10)

式中，x為全局節(jié)點(diǎn)應(yīng)力列向量；fj(x)為屈服準(zhǔn)則或其他條件產(chǎn)生的不等式約束，其他符號意義同前。

1.3 有限元極限分析的計算機(jī)實現(xiàn)

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

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

1.4 基于Hoek-Brown準(zhǔn)則的有限元極限分析

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

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

(11)

式中，σ′1和σ′3分別為最大、最小主應(yīng)力；σc為巖石單軸飽和抗壓強(qiáng)度；mb，s和a為與地質(zhì)強(qiáng)度指標(biāo)GSI有關(guān)的參數(shù)，其表達(dá)式為，

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為與巖石完整程度有關(guān)的參數(shù)；D為擾動系數(shù)，沒有擾動取0，完全擾動取1，本文中D=0。

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

在優(yōu)化模型的求解過程中，除了需要計算當(dāng)前迭代點(diǎn)的屈服函數(shù)值，還必須獲得屈服函數(shù)對應(yīng)力變量的一階、二階導(dǎo)數(shù)，即屈服函數(shù)的梯度向量和Hessian矩陣，具體過程如下：

為了便于推導(dǎo)，將Hoek-Brown準(zhǔn)則轉(zhuǎn)化為應(yīng)力不變量表達(dá)的形式，

(15)

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

g(θ)=-2cos(θ)

(16)

(17)

(18)

(19)

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

(20)

ε是一個與材料抗拉強(qiáng)度相關(guān)的常數(shù)，其值可由式(21)確定。

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

(21)

其中，μ=10-1，δ=10-6，ρ為屈服面與縱軸對應(yīng)的交點(diǎn)。

則雙曲型近似Hoek-Brown準(zhǔn)則可表示為：

(22)

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

2 問題描述

2.1 計算模型及假定

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

圖 4 計算模型Fig. 4 Calculation model

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

圖 5 網(wǎng)格劃分示意圖Fig. 5 Sketch of meshing

2.2 網(wǎng)格劃分

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

(23)

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

(24)

ν=-0.002GSI-0.003mi+0.457

(25)

2.3 評價指標(biāo)

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

Nσ=qu/σc

(26)

表 1 誤差分析Table 1 Error analysis

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

3 結(jié)果與討論

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

3.1 上覆土層的影響

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

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

3.2 嵌巖深度的影響

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

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

3.3 溶洞半徑R的影響

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

圖 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條件下雙樁與單樁承載力系數(shù)對比Fig. 10 Comparison of bearing capacity factor for double-piles and single pile with different X

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

3.5 垂直距離Y的影響

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

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

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

3.6 樁間距L的影響

在X確定的情況下，樁間距L的影響主要體現(xiàn)在右樁與溶洞的位置關(guān)系，當(dāng)右樁離溶洞較遠(yuǎn)時，雙樁基礎(chǔ)的承載力也會較高。雙樁基礎(chǔ)與溶洞的位置關(guān)系在前面已論述，在此不再贅述。

3.7 地質(zhì)強(qiáng)度指標(biāo)GSI的影響

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

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

3.8 參數(shù)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 破壞模式

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

圖 16給出了不同X條件下溶洞的破壞模式，其中，R=2id、Y=2id，當(dāng)-1.5d≤X≤0，發(fā)生整體剪切破壞； 隨著X的增大，破壞模式逐漸變?yōu)橛勺髽犊刂频臎_切破壞，同時右樁底部發(fā)生地基破壞模式。

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

5 結(jié)果驗證

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

6 結(jié) 論

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

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

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

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

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