飽和土不排水柱孔擴張計算的通用半解析解

劉貫飛，雷勝友

(1.長安大學 公路學院, 陜西 西安 710064； 2.中國十九冶集團有限公司, 四川 成都 610031)

在巖土工程中，孔的擴張理論是重要研究課題，涉及如隧道、井筒、沉樁等工程問題[1-2]，目前已有許多學者對此展開研究。早期的孔擴張理論大多基于以Mohr-Coulomb 屈服準則為代表的彈-理想塑性模型，如Vesi?[3]、Carter 等[4]提出的孔擴張模型。由于推導方法和土體本構模型的限制，此時的孔擴張解僅可得到孔周的應力分布和基于不同屈服準則時的孔周彈塑性圈半徑，不能反映土體應力歷史和硬化行為的影響。為合理描述飽和土在圓孔擴張過程中的變形特點，Collins 等[5]采用原始劍橋模型（OCC）和修正劍橋模型（MCC），在假定孔周土有限應變的條件下，推導了各向同性硬化材料中的圓孔不排水擴張解，之后許多學者[6-7]以該解為驗證依據(jù)。然而，由于該解答中對土的豎向應力進行了簡化處理，不能真實反映飽和土中的柱孔擴張效應。Chen 等[8-9]基于MCC 模型使用拉格朗日分析法建立了K0固結飽和土中的不排水和排水柱孔擴張半解析解，該解答沒有簡化土體豎向應力。借鑒Chen 等[8-9]的研究方法，李鏡培等[10]使用MCC 模型建立了飽和土中的不排水球孔擴張半解析解；李林等[11]、Sivasithamparam 等[12]考慮天然沉積飽和土的各向異性特點，采用改進后的MCC 模型推導了土中的不排水柱孔擴張解答；Zhai 等[13]使用改進的MCC 模型建立了結構性土中的不排水柱孔擴張解答。此外，Zhou 等[14]基于MCC 模型建立了同時適用于兩種不同排水條件的飽和土中柱孔擴張統(tǒng)一相似解；武孝天等[15]基于CSUH 模型建立了適用于飽和黏土和砂土中的不排水柱孔擴張解。

以上基于MCC 模型或改進型MCC 模型建立的土中柱孔擴張解答，可以較好反映飽和土在徑向受壓變形過程中的應力變化和硬化特點，具有較大的理論意義。然而，由于臨界狀態(tài)土力學的復雜性，MCC 模型對正常固結和輕度超固結飽和黏土具有較好的適用性，但不適用于嚴重超固結土和顆粒材料[16-18]；另外，該類模型假設土的屈服函數(shù)形狀及臨界狀態(tài)平均有效應力與硬化參數(shù)之間的比值均固定不變，這與實際并不相符。為此，本文使用統(tǒng)一狀態(tài)參數(shù)模型（CASM）[18]和Rowe 剪脹方程[19]來描述土的彈塑性行為，結合大變形理論，通過引入輔助變量，使用拉格朗日分析法建立飽和土中的不排水柱孔擴張半解析解。

1 問題描述和本構模型介紹

1.1 柱孔擴張問題描述

無限均質飽和土中的不排水柱孔擴張如圖1 所示。飽和黏土具有原始的豎向應力σ'v0、水平應力σ'h0=K0σ'v0、孔隙水壓力u0及體積v0。柱孔初始孔徑為a0，在均勻分布的內(nèi)壓力σa作用下發(fā)生擴張，隨著σa的不斷增大，孔周土隨土質點與孔壁之間距離的增大逐漸發(fā)生屈服。擴孔結束后，孔周可分為已經(jīng)屈服的臨界狀態(tài)區(qū)和彈塑性區(qū)，以及未屈服的彈性區(qū)，彈塑性區(qū)的半徑為rp。假設此時的柱孔孔徑為a，孔周某一土質點則從其初始位置rx0移動到了rx，應力狀態(tài)從(σ'h0,σ'h0,σ'v0,u0)變?yōu)?σ'r,σ'θ,σ'z,u)。在上述柱孔擴張過程中，使用統(tǒng)一狀態(tài)參數(shù)模型（CASM）來描述飽和土屈服后的彈塑性變形。

圖1 不排水柱孔擴張示意Fig.1 Schematic diagram of undrained cylindrical cavity expansion

如圖1 所示，可將飽和土中的柱孔擴張視為平面應變過程，在柱坐標系中建立土的應力平衡方程：

式中：σ'r、σ'θ分別為飽和土的徑向、環(huán)向有效應力；u為土中孔隙水壓力；r為土質點的徑向位置。

1.2 統(tǒng)一狀態(tài)參數(shù)模型（CASM）

基于土的狀態(tài)參數(shù)概念和臨界狀態(tài)土力學理論，Yu[18]提出了能靈活描述黏性土和砂土受力變形行為的統(tǒng)一狀態(tài)參數(shù)模型（CASM），近年來被許多學者用于分析飽和土中的旁壓試驗、靜力觸探和樁基沉降等工程問題[19-22]。CASM 模型的屈服面函數(shù)為：

式中：F(p',q,p'm)為土的屈服面函數(shù)；p'm為硬化參數(shù)；M為破壞應力比；如圖2 所示，r*為間距比參數(shù)，表示硬化參數(shù)p'm與臨界狀態(tài)平均有效應力p'cs之間的比值；n為應力-狀態(tài)參數(shù)，控制土的屈服面形狀；p'和q分別為土的平均有效應力和偏應力。

圖2 r* = 3 時不同n 值情況下的 CASM 模型屈服面形狀Fig.2 Yield surface shape of CASM with different n values when r* = 3

在CASM 模型中，Yu 建議用基于Rowe 剪脹方程[23]的塑性勢函數(shù)描述土的塑性變形方向，表達式為：

式中：G(p',q,β)為土的塑性勢函數(shù)；β表示塑性勢函數(shù)的大小，可由p'、q值計算得到。根據(jù)式（3）可得：

式中：H為塑性模量；和dσ'j分別為土的塑性應變增量和有效應力增量，i或j=r、θ、z。

1.3 材料的應力應變關系

飽和土的彈性變形服從廣義胡克定律，應力應變關系為：

將式（4）與式（6）相加，整理后可得到CASM 模型的彈塑性應力應變關系為：

式中：dεi為土的應變增量，i=r、θ、z；dσ'i（或dσ'j）為土的應力增量，i或j=r、θ、z。

2 不排水柱孔擴張半解析解

2.1 柱孔擴張彈性解

不排水柱孔擴張時，土的體應變dεv= dεr+ dεθ+ dεz和豎向應變dεz均為0。并且飽和土的純彈性變形不產(chǎn)生超孔隙水壓力（Δu= 0）。根據(jù)式（1）、（6）可解得彈性區(qū)的應力解[8]為：

式中：σ'rp和rp分別為彈塑性交界面處土的徑向有效應力和彈塑性區(qū)的半徑；ure為彈性區(qū)土的位移。

2.2 柱孔擴張彈塑性解

2.2.1 土體應力計算式 在求解彈塑性區(qū)內(nèi)土體應力時，Chen 等[8]的方法需要在孔周彈塑性區(qū)內(nèi)每個計算點上，根據(jù)土的總應變量分別進行差分求解，使得計算效率不高，且土中孔隙水壓力的求解必須等待上述所有計算點上的土體應力全部求出以后，才能用式（1）進行積分計算，不能直觀反映土中孔隙水壓力在孔徑擴張時的動態(tài)變化。為此，本文引入輔助變量ξ，通過一系列坐標轉換，在拉格朗日坐標系中建立以ξ為自變量，適用于孔周任一土質點上有效應力和孔隙水壓力計算的一階常微分方程組。與Chen 等[8]的方法相比，本文方法更高效，所得結果也可直觀反映孔周土應力和孔隙水壓力在柱孔擴張時的變化。

在孔周彈塑性區(qū)內(nèi)，土體發(fā)生大變形，用對數(shù)應變表示土的徑向應變εr和環(huán)向應變εθ，表達式[8]為：

式中：ξ為輔助變量，ξ=urp/r= (r-r0)/r，urp為孔周彈塑性區(qū)內(nèi)某點土的徑向位移。在式（12）中，以ξ為自變量對εr、εθ取微分，然后將所得dεr和dεθ的表達式及dεz= 0 一起代入式（7），可得：

式中：c11=1-μ2+Ehμbzaθ+Ehμazbθ+Ehazbz+Ehaθbθ；c12=-Ehar(bθ+μbz)+μ(1+μ-Ehbθaz+Ehazbz)；c21=-Ehar(aθ+μaz)+μ(1+μ-Ehaθbz+Ehazbz) ；c22=1-μ2+Eharbr+Ehμarbz+Ehμbraz+Ehazbz；c31=-Ehbr×(μaθ+az)+μ(1+μ+Ehaθbθ-Ehbθaz)；c32=μ+μ2+Ehμarbr-Ehbθaz-Ehμ]arbθ-Ehμbraz；Ψ =(1+μ)×[1-μ-2μ2+Eh(1-μ)(arbr+aθbθ+azbz)+Ehμ(aθbr+bθar+azbθ+bzaθ+azbr+bzar)/E；ai=?F/?σ′i，bi=?G/?σ′i，i=r、θ、z。

在式子ξ=urp/r中令ξ對r求導，結合式（12），可得：

durp/dr=1-exp(εr)=1-1/(1-ξ)式中： 為土質點位移urp對r的空間導數(shù)；dr為孔周任意相鄰兩個土質點之間的徑向距離。根據(jù)式（1）、式（13）和式（14），可以解得：

式（13）和式（16）分別為彈塑性區(qū)內(nèi)基于拉格朗日描述的土體有效應力和孔隙水壓力求解方程式，初始條件為σ'r(ξp)、σ'θ(ξp)、σ'z(ξp)和u(ξp)，ξp對應土質點剛開始發(fā)生塑性變形時的ξ；超孔隙水壓力Δu=u(ξ) –u0。

2.2.2 邊界條件 正常固結土的初始狀態(tài)為式（13）和式（15）的初始條件。對于超固結土，根據(jù)式（2）、式（10），可得：

式中：q(ξp)=M[(ln(Rp′m0)-ln(p′0))/lnr?]1/n；R=p'C/p'm0為表示土體超固結程度的參數(shù)[8]；p'C為土體歷史上的最大硬化程度；p'm0為土體在當前應力條件下對應虛擬的硬化程度。

2.2.3 計算結果的轉換 由式（13）和式（16）求解得到的σ'r、σ'θ、σ'z和u均是關于輔助變量ξ的拉格朗日描述。為便于分析，須將計算結果變換為關于徑向位置r的歐拉描述。

從ξp到ξ對式（14）進行積分，可以得到徑向位置r與輔助變量ξ之間的關系式為：

從ξp到ξa（孔壁處的ξ值）對式（14）進行積分，可以得到擴孔結束后的彈塑性區(qū)半徑rp為：

3 算例分析

3.1 參數(shù)分析

以飽和倫敦黏土為例，分析應力-狀態(tài)參數(shù)n和間距比r*對飽和土不排水柱孔擴張結果的影響。根據(jù)文獻[18]，倫敦黏土的臨界狀態(tài)參數(shù)為：Γ= 2.759，λ= 0.161，κ= 0.062，μ= 0.3，M= 0.888，φ'cs= 22.75°，正常固結土和超固結土的靜止側壓力系數(shù)分別為 1-sinφ和為用土體原位豎向應力定義的土的超固結比）。土的初始體積v0按下式計算：

式中：Γ為p'=1.00 kPa 時土的體積；λ為土的等向壓縮系數(shù)；φ'cs為土的臨界內(nèi)摩擦角。

3.1.1 應力-狀態(tài)參數(shù)n針對R= 1.0 和4.0 的兩種倫敦黏土，計算參數(shù)如表1 所示，土中不排水柱孔擴張結果分別如圖3、4 所示。

表1 r＊ = 3.00 時不同n 值的倫敦黏土初始參數(shù)Tab.1 Initial parameters of London clay with different n values when r* = 3.00

圖3 不同n 值時正常固結倫敦黏土（R = 1.0）中的柱孔擴張結果Fig.3 Results of cavity expansion in normally consolidated London clay (R = 1.0) with different n values

圖4 不同n 值時超固結倫敦黏土（R = 4.0）中的柱孔擴張結果Fig.4 Results of cavity expansion in overconsolidated London clay (R = 4.0) with different n values

如圖3(a)、4(a)所示，當n值不同時，土的屈服面形狀和相應的不排水應力路徑也不同；4 種n值時的正常固結土在柱孔擴張過程中均表現(xiàn)為硬化，而超固結土則均表現(xiàn)為軟化；結合表1 可知，由于同樣初始應力條件下兩種固結程度飽和土的初始體積v0均隨著參數(shù)n的增大而增大，相應不排水應力路徑終點處的p'、q值將隨n的增大而減小。如圖3(b)、4(b)所示，參數(shù)n的取值對孔周σ'r和Δu的分布有較大影響：孔壁附近的σ'r均隨n的增大而減小，并且孔周圍臨界狀態(tài)區(qū)的半徑也隨n的增大而減?。徽９探Y土中擴孔所得Δu隨n的增大而增大，而超固結土中的Δu在孔壁附近（1.0≤r/a≤1.5）隨n的增大而減小，在較遠區(qū)域（1.5

3.1.2 間距比r*針對R= 1.0 和4.0 的兩種倫敦黏土，計算參數(shù)如表2 所示，土中不排水柱孔擴張結果分別如圖5、6 所示。

表2 n = 2.00 時不同r＊情況下的倫敦黏土初始參數(shù)Tab.2 Initial parameters of London clay with different r* values when n = 2.00

圖5 不同r*值時正常固結倫敦黏土（R = 1.0）中的柱孔擴張結果Fig.5 Results of cavity expansion in normally consolidated London clay (R = 1.0) with different r* values

圖6 不同r*值時超固結倫敦黏土（R = 4.0）中的柱孔擴張結果Fig.6 Results of cavity expansion in overconsolidated London clay (R = 4.0) with different r* values

如圖5(a)、6(a) 所示，參數(shù)r*的不同取值使飽和土具有不同的屈服面形狀和不排水應力路徑；結合表2 可知，同樣初始應力條件下土的初始體積v0隨r*的增大而增大，因此飽和土不排水應力路徑終點處的p'、q值隨r*的增大而減小；4 種r*值情況下的正常固結土在柱孔擴張過程中均發(fā)生硬化，而超固結土則在r*= 2.00、3.00、4.00 時發(fā)生軟化，在r*= 5.00 時發(fā)生硬化，這是因為r*值的增加使土的初始狀態(tài)從臨界狀態(tài)線的下方移動到了上方（表2 中狀態(tài)參數(shù)ψ0從負值變?yōu)檎担?，進而使土的硬化行為發(fā)生了改變。如圖5(b)和6(b)所示，參數(shù)r*的取值對孔周的σ'r和Δu分布有較大影響：孔壁附近的σ'r值和臨界狀態(tài)區(qū)的半徑隨r*的增大而減小，Δu隨n的增大而增大。

3.2 不同模型結果比較

選取BBC（Boston blue clay）黏土作為案例進行柱孔擴張計算，并將所得結果與文獻[8]中基于修正劍橋模型（MCC）的結果進行對比，用以驗證本文方法的合理性與創(chuàng)新性。BBC 黏土的基本力學參數(shù)為Γ=2.759，λ= 0.150，κ= 0.030，ν= 0.278，M= 1.200，不同R值情況下的BBC 黏土其余計算參數(shù)如表3 所示。CASM 模型中的屈服面參數(shù)取n= 1.65、r*= 2.00，此時用CASM 模型計算得到4 種BBC 黏土的初始孔隙比與使用MCC 模型時的計算結果相同。另外，本文在使用CASM 模型進行計算時，還分別使用了兩種不同的剪脹關系，Rowe 剪脹關系和MCC 剪脹關系（塑性勢函數(shù)與MCC 模型的屈服函數(shù)相同），進而可以分析土的剪脹關系改變對柱孔擴張結果的影響。為了對比需要，本節(jié)計算得到的所有土體應力均使用土的不排水剪切強度su進行量綱歸一化。

表3 BBC 黏土的計算參數(shù)Tab.3 Calculation parameters for BBC clay

圖7、8 分別給出了在4 種BBC 黏土中不排水柱孔擴張時，歸一化擴孔壓力σa/su和孔壁處超孔隙水壓力Δu/su隨擴孔孔徑a/a0的變化?？梢园l(fā)現(xiàn)，無論使用哪一種剪脹關系，基于CASM 模型的計算結果都與基于MCC 模型的基本一致，說明若選擇合適的屈服面參數(shù)n、r*，本文解答可近似計算修正劍橋類黏土在不排水柱孔擴張時的應力變化。

圖7 歸一化的擴孔壓力隨孔徑的變化Fig.7 Variation of normalised internal cavity pressure with cavity radius

圖8 歸一化的孔壁超孔隙水壓力隨孔徑的變化Fig.8 Variation of normalised excess pore pressure at cavity wall with cavity radius

圖9 給出了4 種BBC 黏土中不排水柱孔擴張a/a0= 2.0 時，孔周圍σ'r、σ'θ和σ'z的分布。可以看到，3 種方法計算所得不同R值飽和土的σ'r、σ'θ和σ'z在孔壁附近均為定值，說明土體達到了臨界狀態(tài)。對于正常固結土（R=1.0）和中等超固結土（R= 1.2、3.0），當使用MCC 剪脹關系時，本文方法所得σ'r/su、σ'θ/su和σ'z/su與文獻[8]的在孔周3 個區(qū)域內(nèi)均吻合；當使用Rowe 剪脹關系時，所得結果與文獻[8]的在臨界區(qū)和彈性區(qū)內(nèi)吻合，在彈塑性區(qū)存在一些小的差異。對于嚴重超固結土（R= 10.0），文中分別使用兩種剪脹關系所得孔周土體應力分布在臨界狀態(tài)區(qū)和彈性區(qū)內(nèi)與文獻[8]的結果吻合，在彈塑性區(qū)內(nèi)存在差異，并且本文所得孔周的彈塑性區(qū)半徑明顯小于文獻[8]的，這主要是因為兩種本構模型中屈服面函數(shù)的形狀不同所致。圖10 給出了n= 1.65、r*= 2.00 時CASM 模型中的屈服面與同樣硬化參數(shù)時MCC 模型中的對比?？梢园l(fā)現(xiàn)，CASM 模型的屈服面在臨界狀態(tài)線的下方與MCC 模型吻合較好，但在臨界狀態(tài)線的上方低于MCC 模型，這使得同樣初始條件的嚴重超固結土（R= 3.0、10.0）在CASM 模型中初始屈服時的偏應力q要小些，因此更為合理[16]。

圖9 a/a0 = 2.0 時孔周圍徑向、環(huán)向和豎向有效應力分布Fig.9 Radial, tangential and vertical effective stress distributions around the cavity with a/a0 = 2.0

圖10 CASM 模型和MCC 模型的屈服面形狀對比Fig.10 Contrast of yield surface shapes in CASM and MCC

3.3 實例應用

使用Lehane[24]在弱超固結Bothkennar 黏土中的壓樁試驗對本文解答的可靠性進行驗證。場地內(nèi)從地下2 m 到7 m，土的超固結比Roc從2.00 降為1.46，靜止側壓力系數(shù)K0,OC從0.650 降為0.500，有效豎向壓力σ'z0從25.50 kPa 增加到54.80 kPa，均近似為線性增長或減小；土的有效內(nèi)摩擦角均值為35.50o，計算可知Mc= 1.430；土的壓縮系數(shù)Cc= 0.5，回彈系數(shù)Cs= 0.033，換算得到λ= 0.217、κ= 0.015；地下水位的埋深為2 m。如圖11(a) 所示，按照Yu[18]中給出的方法，使用正常固結Bothkennar 黏土的不排水剪切應力路徑來近似確定n、r*，得到n=2.00、r*= 2.50。

圖11 沉樁實例計算Fig.11 Calculation example for pile sinking

圖11(b)、11(c)分別給出了使用本文方法和MCC 模型計算所得Bothkennar 黏土中壓樁時樁側的孔隙水壓力、徑向總應力的分布?？梢园l(fā)現(xiàn)，本文方法所得孔隙水壓力、徑向總應力結果和使用MCC 模型時所得兩種結果的趨勢均與實測一致，但本文所得樁側徑向總壓力更接近于實測值。

4 結 論

本文使用統(tǒng)一狀態(tài)參數(shù)模型（CASM）和Rowe 剪脹關系建立了飽和土不排水柱孔擴張的通用解答，通過算例討論了間距比r*、應力-狀態(tài)參數(shù)n和土體剪脹關系等對柱孔擴張結果的影響，分析了本文方法的優(yōu)越性并結合算例驗證了其實用性，得出了如下結論：

（1）改變參數(shù)n、r*的取值可以使土的屈服面形狀發(fā)生改變，進而使CASM 模型可用于計算不同類型飽和重塑黏土中的不排水柱孔擴張力學過程。

（2）改變參數(shù)n、r*的大小對土中柱孔擴張結果有較大影響，其中孔壁處的徑向有效應力隨n或r*的增大而減小，而超孔隙水壓力隨n或r*的增大而增大。

（3）選擇合適的n、r*值和土體剪脹關系，CASM 模型可近似替代修正劍橋模型計算飽和重塑黏土中的不排水柱孔擴張力學過程，并且對嚴重超固結土的計算結果更合理。

（4）本文解答可用于估算飽和黏土中壓樁時樁側的孔隙水壓力和徑向總應力，并且參數(shù)n、r*的取值可用正常固結土的不排水剪切應力路徑來確定。