考慮時變效應的堆載引發鄰近單樁水平變形分析

閔鵬，申玉生，林作忠，董俊，雷龍

(1.西南交通大學交通隧道工程教育部重點實驗室，四川成都，610031；2.中鐵第四勘察設計院集團有限公司，湖北武漢，430063)

隨著城鎮化進程的發展，城市建設用地日益緊張。樁基被廣泛應用于高層建筑、橋梁墩臺、港口碼頭等基礎工程中。近年來，由鄰近地面堆載引發的工程事故屢見不鮮[1?3]。地面堆載將改變地層原有的應力狀態，導致既有樁基撓曲變形，進而引發一系列不利影響。因此，準確地預測地面堆載對鄰近既有樁基產生的變形具有重要的現實意義。目前，各國學者對此類似問題進行了廣泛研究，主要包括理論解析、數值模擬、現場監測和模型試驗等。

與數值模擬、現場監測和模型試驗相比，理論解析作為一種低成本而高效率的分析方法，適用于指導工程的初步設計。在目前的理論研究中，通常采用兩階段法來分析地面堆載對樁基的影響[4?6]。基于兩階段法，許多研究者從地基彈性入手，將樁基簡化為擱置在Winkler彈性地基上的連續梁[7?9]。Winkler模型將地基看作由若干互不影響的獨立彈簧組成，地基表面任一點處反力僅與該點彈簧的變形量成正比，該模型簡化了地基復雜的力學現象，但忽略了彈簧間的剪切作用和地基的連續性，不能很好地反映土體的力學特性。為了克服這一缺陷，梁發云等[10?11]采用Pasternak 雙參數地基模型來模擬樁?土間的相互作用，將樁基簡化為Euler-Bernoulli 梁，建立了樁身的撓曲控制微分方程。Pasternak 雙參數地基模型在Winkler 模型的基礎上假定地基彈簧之間存在剪切效應，在各彈簧頂部增加一層剪切層，從理論上改進了Winkler 地基不連續的缺陷。此外，大多研究將樁基簡化為Euler-Bernoulli 梁，僅考慮了彎曲變形而忽略了剪切效應的影響，導致樁基的剪切剛度被高估，不能很好地反映樁基的實際受力狀態。而Timoshenko 梁理論適用于大多數樁基，并能得到更加精確的結果[12]。目前已有究者采用Timoshenko 梁來模擬樁基[13?16]，同時考慮了樁基的彎曲效應與剪切效應。

上述的理論研究中很少考慮時間對樁?土間相互作用的影響。然而，實際工程中樁基的變形會隨著時間的變化逐漸發展，在長期的變形過程中應考慮土體的流變性對樁基的影響[16]。采用簡化黏彈性模型和彈性?黏彈性對應原理來研究土體流變性問題是目前較為普遍的研究方法[17]。張強[18]基于黏彈性地基梁理論，建立了開挖卸荷導致下方既有隧道發生隆起位移的豎向變形公式，分析了時間對變形的影響。張治國等[19?20]基于Boltzmann黏彈性模型，采用兩階段法研究了黏彈性地基中基坑開挖對鄰近樁基水平位移的影響及其隨時間的變化趨勢。但在地面堆載引發鄰近既有單樁水平變形的研究中大多未考慮時間的影響。

因此，本文作者從地基黏彈性角度出發，基于兩階段法，提出一種預測黏彈性地基中地面堆載引發鄰近既有單樁水平變形的簡化分析方法，并通過有限差分和Laplace 逆變換求得單樁水平變形黏彈性解。將本文所得理論解與數值模擬結果進行比較，并與已有解進行對比分析，驗證本文方法的正確性。最后，分析Burgers 模型參數、豎向均布堆載、水平均布堆載以及堆載?樁體水平距離對既有單樁水平變形的影響。

1 樁體附加應力黏彈性解

1.1 黏彈性半空間Boussinesq解

Boussinesq解為半無限體表面一點受法向集中力作用的經典彈性解[21]，可以用于研究地面荷載作用下在巖土體內部產生的附加應力場問題。計算模型如圖1所示。

圖1 半無限體表面受集中力示意圖Fig.1 Schematic diagram of half space with forces applied on surface

在半無限體表面O(0，0，0)處豎向集中荷載F作用下，半無限體內一點H(x，y，z)處引起的水平附加應力σx(x，y，z)為

式中：ν為泊松比；R為H點到O點的距離，R=

Burgers 模型是一種常用的線性黏彈性模型，可以較好地反映土體的流變性，既能描述加載初期的瞬時彈性變形，又能描述土體隨時間變化的特性[22?24]。本文選用Burgers 模型模擬巖土體的黏彈性本構關系，研究黏彈性半空間體表面堆載在空間內部某一點所產生附加應力的擬靜態黏彈性解。所采用的Burgers 模型由一個Maxwell 單元和一個Kelvin單元串聯組成，如圖2所示。

圖2 Burgers模型示意圖Fig.2 Schematic diagram of Burgers model

根據Burgers 黏彈性模型，在一維條件下巖土體的本構關系可以表示為

對于各向同性黏彈性材料，在三維應力狀態下，可將應力、應變張量分別分解為球張量與偏張量，應變球張量只與應力球張量有關，應變偏張量只與應力偏張量有關。因此，三維黏彈性體本構關系一般為：

式中：Sij為應力偏張量；eij為應變偏張量；σ為應力球張量；e為應變球張量；P′，Q′，P″和Q″為黏彈性本構方程的微分算子，

空間三維應力狀態下，Burgers 黏彈性模型中反映巖土體彈性和黏性的本構關系分別為：

式中：G為剪切模量；η為黏滯系數。

將巖土體的本構關系分為剪切變形和體積變化2 部分，其中，剪切變形服從Burgers 黏彈性本構關系，則其偏張量部分的表達式為

體積變化服從彈性本構關系，則其球張量部分的表達式為

式中：K為體積模量。

將式(11)和(12)代入式(5)～(8)可得黏彈性本構關系各微分算子分別為：

對式(13)～(16)進行Laplace變換可得：

黏彈性材料參數的Laplace 變換與各微分算子間的關系為：

由彈性力學關系可得：

式中：ν為泊松比；μ和λ為拉梅系數。

根據彈性?黏彈性對應原理，將(s)和(s)代替式(23)和(24)中的μ和λ，即可得到黏彈性條件下泊松比ν(s)的Lapalce變換式分別為：

將式(17)～(20)代入式(25)和(26)可得：

假定半無限黏彈性體表面原點O處的法向集中荷載F(t)=FH(t)，并對F(t)進行Laplace 變換可得：

根據彈性?黏彈性對應原理，將式(26)代入式(1)，即可得到在表面集中荷載F作用下半無限體內部任意一點H(x，y，z)處x向水平附加應力關于時間t的Laplace變換：

其中，

對式(31)～(33)進行Laplace 逆變換，可得到集中荷載F作用下基于Burgers 模型的水平附加應力在時間域上的黏彈性解：

其中，

1.2 樁體水平附加應力黏彈性解

以荷載中心為坐標原點建立如圖3所示的附加應力計算模型，圖中，H為樁長，B和L分別為堆載寬度和堆載長度，在B×L內作用法向均布荷載F。

圖3 單樁計算模型示意圖Fig.3 Schematic diagram of computing model of pile

根據式(34)，由荷載微元Fdxdy在樁體任意深度M(L0，Y0，z)處所產生的水平附加應力為

在B×L內對式(36)進行積分，即可求出地面均布荷載F引起樁體軸線上任意一點P(L0，Y0，z)處水平附加應力為

引入數值積分進行求解式(36)，為了使結果能夠達到較高的代數精度，還需對積分式進行復化處理，對荷載積分區間進行等分，要求每個矩形子區間邊長不超過1 m，對各子區間運用Gauss-Legendre數值積分公式求解，將各子區間積分值疊加即可得到樁體水平附加應力σx(z,t)。

2 樁體水平變形黏彈性解

通過兩階段分析法研究地面堆載作用下單樁水平撓曲黏彈性解。第一階段，忽略樁體的存在，考慮土體的流變性并結合傳統的Boussinesq解推導出地面堆載引發的土體附加應力黏彈性解；第二階段，將地面堆載引發的附加應力作為外荷載施加在樁體上，根據黏彈性地基梁理論，并結合有限差分法求得單樁水平撓曲的理論解。

2.1 樁?土系統計算模型

在Pasternak 地基模型的基礎上，引入能夠考慮土體流變性的線性黏彈性地基模型，如圖4所示。該模型由若干個Burgers 線性黏彈性單元體和能反映單元體間剪切相互作用的剪切層組成，其表達式為

圖4 黏彈性地基模型示意圖Fig.4 Schematic diagram of viscoelastic foundation model

其中：

p(z,t)為地基反力；ω(z,t)為樁體的水平位移；Gp(t)為剪切層剛度；A(t)為蠕變柔量；E(t)和μ(t)分別為考慮時變效應的地基彈性模量和地基泊松比；ht為Pasternak 剪切層厚度，本文取ht=11D[19]，D為樁體等效寬度。

為簡化計算，對該模型進行如下基本假定：1)樁體為考慮剪切效應的Timoshenko 矩形梁，其等效剪切剛度為κGA0，等效抗彎剛度為EI，等效寬度為D；2)樁體與周圍土體接觸緊密，不考慮樁?土間相對滑動及摩擦力；3)假設地基連續均勻，不考慮土體的成層性；4)Pasternak 剪切層不可壓縮，僅產生剪切變形；5)不考慮重力產生的影響。

取樁體上一微元段進行受力分析，如圖5所示。根據材料力學，微元段水平受力平衡方程和彎矩平衡方程分別為：

圖5 微單元受力分析Fig.5 Viscoelastic foundation model

式中：q(z,t) 為作用在樁體上的附加荷載，由式(37)計算；Q為剪切力；p(z,t) 為黏彈性地基反力；M為彎矩。

由Timoshenko 梁理論，樁體剪力Q與彎矩M表達式為：

式中：κ為截面系數，此處取κ=0.83；G為樁體剪切模量，G=E/2(1+μ0)，μ0為樁體泊松比；A0為樁體橫截面積；θ為轉角。

根據式(41)～(44)可推導得地面堆載作用下，樁體水平位移ω(z,t)、樁體彎矩M(z,t)和剪力Q(z,t)的控制微分方程分別為：

2.2 樁體水平變形求解

對式(45)進行關于時間t的Laplace變換可得：

式中：ω(z,s)和q(z,s)分別為ω(z,t)和q(z,t)的Laplace 變換；A(s)和Gp(s)分別為A(t)和Gp(t)的Laplace 變換；E(s)和μ(s)分別為地基彈性模量和地基泊松比的Laplace變換。

對式(48)采用有限差分法進行數值求解。樁體有限差分示意圖如圖6所示，將樁體劃分為n個長度相同的微元段，并在樁體頂部和底部分別增設2個虛擬節點用以計算樁體的邊界條件。

圖6 樁體有限差分示意圖Fig.6 Schematic diagram of finite difference of pile

根據有限差分原理，式(48)的差分形式為

式中：l為單元長度，l=H/n，H為樁長；下標i為節點編號，i∈(0，n)。

考慮樁體兩端自由無約束時，樁體頂部和底部彎矩和剪力均為0。結合有限差分原理，樁體彎矩和剪力分別為：

式中，當下標i=0 時，表示樁體頂部處彎矩和剪力；當下標i=n時，表示樁體底部處彎矩和剪力。

根據式(52)～(53)可以得到樁體兩端虛擬節點位移為：

將虛擬節點位移式(54)～(57)代入式(51)，可將差分方程改寫為向量和矩陣表達式：

根據式(58)即可解得在地面荷載作用下樁體的水平變形ω(z,t)，并對ω(z,t)進行Laplace 逆變換即可得到樁體在時間域上的水平變形黏彈性解ω(z,t)：

特別地，當樁體等效剪切剛度κGA0取無窮大時，ω(z,t) 將退化為歐拉梁上的解，當t取0 時，退化為Pasternak地基上的彈性解。

3 分析與驗證

3.1 收斂性研究

樁單元劃分數量是本文求解的關鍵，微元段劃分精度直接決定最終求解的精度，因此在上述公式的基礎上，編制專門的Matlab 計算程序，通過收斂性分析，得到微元段劃分的最優個數。收斂性判斷標準如下：

式中：k為當前迭代次數；k-1為前一次迭代次數；c為容許誤差，此處取為0.1 mm。

假設地面堆載范圍B×L=30 m×50 m，堆載載荷F=150 kN/m2。荷載中心距離樁體中線20 m。樁體等效寬度為1.0 m，樁長為40 m，樁體彈性模量為31.5 GPa，泊松比為0.2。Burgers 模型參數可根據實際工程的監測數據進行參數反演或對室內蠕變試驗曲線進行擬合來確定[25?26]，本文地基體積模量K=0.68 GPa，Kelvin 黏滯系數ηk=0.1 GPa·d，Kelvin 彈性模量Gk=3.0 MPa，Maxwell 黏滯系數ηm=75.0 GPa?d，Maxwell 彈性模量Gm=5.0 MPa。圖7所示為深度z=20 m 處樁體水平位移與微單元劃分數量的關系曲線。從圖7可以看出，當微元段數量達到40，樁體中點水平位移穩定收斂，位移變化可以忽略不計，因此，在同時考慮其他因素的條件下，樁體微單元長度不超過1 m即能滿足計算收斂。

圖7 水平變形隨單元劃分數的變化Fig.7 Variation of deformation with number of element

3.2 數值驗證

為了驗證所提理論解的正確性，采用FLAC3D有限差分軟件建立三維數值模型模擬地面堆載作用下鄰近單樁的水平撓曲，通過對比數值模擬與理論解的計算結果驗證本文分析方法的正確性。三維數值模型如圖8所示。為了提高數值計算的精度，對樁體及其周圍網格進行了細化處理，模型x×y×z=150 m×120 m×60 m，其余計算參數同3.1節。

圖8 數值計算模型Fig.8 Numerical model

圖9所示為深度z為5.0，15.0，25.0 和35.0 m處樁體的水平變形曲線。從圖9可以看出，地面堆載作用下樁體的水平變形有明顯的時變效應。在t=0 d時，深度z為5.0，15.0，25.0和35.0 m處，本文方法計算的水平變形分別為7.17，6.98，3.87 和1.87 mm；t=500 d時，相應位置的水平變形分別為19.98，18.05，11.45 和5.42 mm，隨時間的延長，相應位置的變形分別為12.81，11.07，7.58 和3.55 mm，依次占總水平變形的64.1%，61.3%，66.2%和68.3%。在t=0 d 時，由于地基的彈性效應，加載瞬間樁體產生變形，在t=0～120 d 時，樁體變形速率持續減小，但此階段樁體水平變形隨時間不斷累積，土體流變性的影響最為顯著，120 d 后，進入穩定變形階段，樁體的水平變形速率不再發生變化，變形趨于平穩。因此，地面堆載施加后的120 d內是對樁體水平變形影響最大的時間段。對比本文解與有限差分數值解可以看出，本文解與有限差分數值解有較好的一致性，說明本文推導的地面堆載引發鄰近單樁水平變形的黏彈性解具有較高的計算精度。

圖9 樁體水平變形?時間關系曲線Fig.9 Curves of pile horizontal deformation and time

3.3 案例驗證

與馮昌明等[6,27]提出的堆載引發樁體水平變形計算方法進行對比分析，其中案例基本工況如下：地面矩形均布堆載范圍長×寬為6 m×6 m、載荷為20 kPa，樁體軸線與堆載邊緣凈距為2 m，樁長為21 m，樁徑為0.5 m，樁體和土的彈性模量分別為40 GPa 和5 MPa，泊松比為0.3[27]。計算結果如圖10所示，從圖10可以看出：本文方法所得樁體水平變形與文獻[6,27]所得結果基本吻合，從而驗證了本文方法的正確性。

圖10 樁體水平變形曲線Fig.10 Horizontal displacement of pile

4 Burgers模型參數分析

4.1 Maxwell彈性模量的影響

在其他參數不變的情況下，取Maxwell彈性模量Gm分別為2.5，5.0 和7.5 MPa，探究Gm對地面堆載作用下樁體水平撓曲的影響。圖11所示為t=0 d 和t=120 d 時不同Gm下樁體的水平變形。從圖11可以看出：Gm對樁體水平變形的影響顯著，Gm越小，樁體水平變形越大，樁體水平變形對Gm的變化越敏感。對比不同時刻樁體水平變形曲線還可以看出：當t=120 d 時，水平變形量顯著增加，說明樁體水平變形有明顯的時變效應。

圖11 不同Gm下樁體水平變形對比Fig.11 Horizontal deformation comparison of pile with different Gm

4.2 Kelvin彈性模量的影響

在其他參數不變的情況下，取Kelvin 彈性模量Gk，分別為1.5，3.0和4.5 MPa，探究Gk對地面堆載作用下樁體水平撓曲的影響。圖12所示為t=0 d 和t=120 d 時不同Gk下樁體的水平變形曲線。從圖12可以看出，t=0 d時，不同Gk下樁體的水平變形幾乎一致，說明堆載施加瞬間樁體的水平變形不受Gk的影響。在t=120 d時，樁體的水平變形與t=0 d 時相比有明顯增加，說明土體的流變性對樁體水平變形有顯著影響，同時，t=120 d 時，Gk越小，樁體的水平變形越大，說明在Burgers 體中瞬時變形主要受Gm的影響而與Gk無關，Gk主要反映由時變效應產生的黏彈性變形。

圖12 不同Gk下樁體水平變形對比Fig.12 Horizontal deformation comparison of pile with different Gk

4.3 Maxwell黏滯系數的影響

在其他參數不變的情況下，取Maxwell黏滯系數ηm分別為15.0，30.0和75.0 GPa·d，探究ηm對地面堆載作用下樁體水平撓曲的影響。圖13所示為不同ηm下深度20.0 m 處樁體水平變形曲線。從圖13可以看出，在t=0～120 d時，不同ηm下樁體水平變形對于時間的響應規律基本一致，說明在此時間段樁體水平變形幾乎不受ηm變化的影響。當t>120 d 后，樁體進入穩定變形階段。ηm越小，樁體水平變形速率越快，導致樁體隨時間發展累積更大的水平變形，說明ηm對最終樁體水平變形值有直接影響。

圖13 不同ηm下樁體水平變形?時間關系曲線Fig.13 Curves of horizontal deformation and time of pile with different ηm

4.4 Kelvin黏滯系數的影響

在其他參數不變的情況下，取Kelvin 黏滯系數ηk分別為0.1，0.2和0.4 GPa·d，探究其對地面堆載作用下樁體水平撓曲的影響。圖14所示為不同ηk下深度20.0 m 處樁體水平變形曲線。從圖14可以看出，ηk對樁體瞬時變形幾乎沒有影響，對樁體最終水平位移也無明顯影響，但是對于流變變形進入穩定變形階段的時間有明顯影響。ηk越大，進入穩定變形階段的時間越長，當ηk=0.1 GPa·d時，樁體水平變形在t=70 d后基本穩定，而當ηk=0.4 GPa·d時，樁體水平變形在t=240 d 后才基本穩定，歷時明顯增加。因此，在實際工程中應根據ηk的實際取值制定監測方案。

圖14 不同ηk下樁體水平變形?時間關系曲線Fig.14 Curves of horizontal deformation and time of pile with different ηk

5 堆載影響分析

5.1 堆載寬度的影響

在其他參數不變的情況下，取堆載寬度B分別為20.0，50.0 和100.0 m，探究在黏彈性條件下B對樁體水平變形的影響。圖15所示為t=0 d 和t=120 d 時，不同堆載寬度下樁體的水平變形。從圖15可以看出，B對樁體水平變形有明顯的影響，B越大，樁體水平變形越大，同時樁體水平變形對B的變化也越敏感。對比t=0 d和t=70 d時樁體水平變形曲線，不難發現樁體水平變形隨著時間的發展逐漸累積，地基的流變性對樁體水平變形有明顯影響。

圖15 不同B時樁體水平變形對比Fig.15 Horizontal deformation comparison of pile of with different B

圖16所示為B為20.0 m和100.0 m，t=120 d時樁體周圍土體的水平應力分布。從圖16可以看出：隨著堆載寬度增加，土體水平應力明顯增加，樁體產生更大的水平變形。

圖16 不同B時土體的水平應力分布Fig.16 Horizontal stress distribution of soil with different B

5.2 堆載載荷的影響

在其他參數不變的情況下，取均布載荷F分別為100.0，150.0 和200.0 kN/m2，探究在黏彈性條件下F對樁體水平變形的影響。圖17所示為t=0 d和t=120 d 時不同F下樁體的水平變形。從圖17可以看出，隨著F增加，樁體產生更大的水平變形，說明F對樁體變形有顯著影響。

圖17 不同F時樁體水平變形對比Fig.17 Horizontal deformation comparison of pile with different F

F分別為100.0 kN/m2和200.0 kN/m2，t=120 d時樁體周圍土體的水平應力分布如圖18所示。從圖18可以看出：隨著F增加，土體水平應力明顯增加，導致樁體處水平附加應力增加，樁體產生更大的水平變形。

圖18 不同F時土體水平應力分布Fig.18 Horizontal stress distribution of soil with different F

5.3 堆載?樁體水平距離的影響

在其他參數不變的情況下，取堆載?樁體水平距離L0分別為20，25 和30 m，探究其對地面堆載作用下樁體水平變形的影響。圖19所示為t=0 d和t=120 d時不同L0下樁體的水平變形。從圖19可以看出，在深度0～25 m 范圍內，L0越小樁體水平變形越大，在深度大于25 m 時，樁體水平變形基本不受L0影響。對比t=0 d和t=120 d時樁體的水平位移還可以看出：樁體的水平變形有明顯的時變效應。

圖19 不同L0時樁體水平變形對比Fig.19 Horizontal deformation comparison of pile with different L0

6 結論

1)微元段劃分數量直接影響黏彈性解的求解精度，當微元段劃分數量達到40，樁體水平變形即可穩定收斂。

2)黏彈性地基上樁體水平變形有明顯的時變效應，與數值模擬結果有較好的一致性，驗證了本文方法預測地面堆載作用下樁體水平變形的準確性，通過與已有解進行對比分析，進一步驗證了本文方法的正確性。

3)Burgers 模型各參數對樁體水平變形均有較大影響。Gm對樁體瞬時變形有較大影響，Gk主要反映由時變效應產生的黏彈性變形，ηk決定變形進入穩定階段的時間跨度而對最終變形量基本沒有影響，ηm主要影響樁體穩定變形階段的變形速度。

4)地面堆載寬度和堆載載荷對樁體水平變形均有較大影響，并且樁體水平變形隨著堆載中心與樁體水平距離的減小而明顯增加。

5)本文主要研究均質地基中黏彈性條件下堆載引發鄰近既有實心單樁的水平變形響應，并未考慮地基土的成層性、開口管樁和群樁基礎等工況，因此，本文方法不能分析地基成層性、開口管樁土塞效應和群樁效應等對樁體水平變形的影響，針對這些問題，還需要進一步研究。