線性增長樁側阻力半無限體多軸向應力解研究

目的 為獲得樁側線性增長分布力作用下土體的徑向及環向應力理論解析解。 方法 通過彈性力學半無限空間軸對稱問題微元體的幾何方程、本構關系方程及使用積分方法。結果在半無限空間體內部一集中豎向荷載下的cerruti位移解和Mindlin豎向正應力解的基礎上推導出了單樁在三角形分布樁側阻力荷載條件下樁外任何一點土體的環向正應力、徑向正應力及徑向剪應力等的理論解計算公式，并給出了計算機運行通過的它們各自的VisualBasic程序代碼清單和計算案例。結論得出的計算公式及計算機計算程序將為全面考慮明德林效應的樁基極限承載力效應計算和沉降效應計算等提供必要的理論支撐。

樁基礎； 多軸向應力解； 半空間無限彈性體； 明德林效應

TU473.1+1A

［定稿日期］2023-04-04

［基金項目］國家自然科學基金資助項目（項目編號：50978171）

［作者簡介］宋寶峰（1974—），男，碩士，高級實驗師，研究方向為建筑結構、建筑施工。

0 引言

美國哥侖比亞大學Raymond D.Mindlin［1］教授在1936年推導出了半無限空間體內部一集中力下的土體各向位移及多軸向應力解公式，但其計算公式復雜臃長，且沒有給出樁側側阻力對巖土產生的其他各類效應的計算公式；2010年我國的學者宋東輝教授等［2］利用對稱性將半無限空間體的Mndlin問題轉化為半無限空間體的Cerruti問題給出樁端一豎向集中力下的土體各向位移相對簡化的理論計算公式。本文在此Cerruti簡化位移公式的基礎上通過積分方法得出了樁周三角形分布荷載下的土體各徑向及環向應力解析解的具體計算公式。在樁基工程的沉降近似計算時，我國JGJ 94-2008［3］《建筑樁基技術規范》僅考慮了土體在樁側、樁端阻力作用下的Mindlin豎向應力下的貢獻，并給出了Mindlin豎向應力的計算公式及表格，實際上土體在樁側、樁端阻力作用下的Mindlin徑向及環向應力同樣會對土體的沉降效應產生貢獻，特別是在考慮群樁效應計算樁側及樁端極限承載力時產生的貢獻更不容忽視，所以對樁在復雜受力條件下土體的徑向及環向應力解析解的研究具有重要意義。

本問題屬于空間軸對稱問題，設有一無限半空間彈性體（圖1），不考慮土體自重，E為土體的彈性模量，G為土體的剪切模量，距離地表面深度為a處有一樁端豎向集中力P0;在距離縱坐標軸為r、距離樁頂平面深度為z處有一土粒微元體，ur為其徑向位移，w為其豎向位移，σz為其豎向應力，σr為其徑向應力，σθ為其環向應力，τzr為其豎向徑向剪應力。

簡化后的Cerruti樁端集中力下的半無限空間土體任一點的位移公式［2］見式（1）、式（2）。

ur=P016πG（1-μ）［r（z-a）M3+r（z+a）N3］（1）

w=P016πG（1-μ）M［（3-4μ）+（z-a）2M2］+

P016πG（1-μ）N［（3-4μ）+（z+a）2N2］（2）

式中：μ為土體的泊松比；M=r2+（z-a）2；N=r2+（z+a）2，G=E2（1+μ）。

1 樁側線性增長三角形分布荷載下四向應力σz、σr、σθ、τzr解析解的求解

1.1 本問題土體微元體的各向位移計算

本問題屬于空間軸對稱問題，設有一半無限空間彈性體（圖2），不考慮土體自重，E為土體的彈性模量，G為土體的剪切模量，在樁長l范圍內樁側土受正三角形分布側阻力荷載q（a）的作用，總合力為Q=Qs;在距離縱坐標軸為r、距離樁頂平面深度為z處有一土粒微元體，ur為其徑向位移，w為其豎向位移，σz為其豎向應力，σr為其徑向應力，σθ為其環向應力，τzr為其豎向徑向剪應力。

可以利用上述樁端集中荷載下的土體任意點的位移公式［2］，通過樁長范圍內的側阻力位移貢獻積分求出本情況的土體任意點的位移公式，因為樁側阻力成線性分布，故有：q（a）=2Qsal2，這里a為地表面至積分計算點的豎向距離，有式（3）、式（4）。

ur=∫l0q（a）16πG（1-μ）r（z-a）［r2+（z-a）2］32+r（z+a）［r2+（z+a）2］32da=rQs8πG（1-μ）l2l（l-z）2+r2-l（l+z）2+r2+ln（l+z）2+r2+（l+z）r2+z2+z-ln（l-z）2+r2+（l-z）r2+z2-z（3）

w=∫l0q（a）16πG（1-μ）3-4μr2+（z-a）2+（z-a）23r2+（z-a）2+3-4μr2+（z+a）2+（z+a）23r2+（z+a）2da

=Qs（3-4μ）8πG（1-μ）l2（l-z）2+r2+（l+z）2+r2-2r2+z2+zln［（l-z）+（l-z）2+r2r2+z2-z］-

zln［（l+z）+（l+z）2+r2r2+z2+z］+Qs8πG（1-μ）l2（l-z）2+r2+（l+z）2+r2-4r2+z2+

r2r2+（l-z）2+r2r2+（l+z）2+zln［（l-z）+（l-z）2+r2r2+z2-z］-

zln［（l+z）+（l+z）2+r2r2+z2+z］-z（l-z）（l-z）2+r2+z（l+z）（l+z）2+r2（4）

1.2 根據空間軸對稱問題的幾何方程求出微元體各向應變εr，εθ，εzr

應用上述位移解的結論，由彈性力學空間軸對稱問題［4］的結論本文進一步推導出式（5）～式（7）。

εθ=urr=Qs8πG（1-μ）l2l（l-z）2+r2-l（l+z）2+r2+

ln［（l+z）2+r2+（l+z）r2+z2+z］-ln［（l-z）2+r2+（l-z）r2+z2-z］（5）

εr=urr=εθ+rQs8πG（1-μ）l2rl3（l+z）2+r2-rl3（l-z）2+r2+r［（l+z）2+r2］-0.5（l+z）2+r2+（l+z）-r［（l-z）2+r2］-0.5（l-z）2+r2+（l-z）+r（r2+z2）-0.5（r2+z2）0.5-z-r（r2+z2）-0.5（r2+z2）0.5+z（6）

γzr=urz+wr=rQs8πG（1-μ）l2l（z+l）［（l+z）2+r2］-1.5-l（z-l）［（l-z）2+r2］-1.5+dyc

+Qs（3-4μ）8πG（1-μ）l2r（l-z）2+r2+r（l+z）2+r2-2r（r2+z2）-0.5+syc

+Qs8πG（1-μ）l2dhj-4r（r2+z2）-0.5+syc+jcm-rz（z+l）［（l+z）2+r2］-1.5-rz（z-l）［（l-z）2+r2］-1.5（7）

注意式（7）中：

dyc=（z+l）［（l+z）2+r2］-0.5+1（l+z）2+r2+（l+z）-（z-l）［（l-z）2+r2］-0.5-1（l-z）2+r2+（l-z）+z（r2+z2）-0.5-1r2+z2-z-z（r2+z2）-0.5+1r2+z2+z

syc=zr［（l-z）2+r2］-0.5（l-z）+（l-z）2+r2-zr［（l+z）2+r2］-0.5（l+z）+（l+z）2+r2+rz（r2+z2）-0.5z2+r2+z-rz（r2+z2）-0.5z2+r2-z

dhj=r［（l-z）2+r2］-0.5+r［（l+z）2+r2］-0.5

jcm=2r（l+z）2+r2-r33r2+（l+z）2+2r（l-z）2+r2-r33（l-z）2+r2

這里εr為土粒微元體的徑向應變，εθ為其環向應變，γzr為其豎向徑向剪應變。

1.3 根據空間軸對稱問題的物理方程求出微元體各向應力σr，σθ，σzr

σz可以在一些文獻資料中查找，是本問題在半無限彈性體內一向下的集中力土體豎向正應力Mindlin［1］解直接積分得出的結果，例如可從GB 50007-2011《建筑地基基礎設計規范》附錄R［5］中查找，本文這里就不再重復推導。根據空間軸對稱問題的物理方程［4］：

εr=1E［σr-μ（σθ+σZ）］，εθ=1E［σθ-μ（σr+σZ）］

εz=1E［σz-μ（σθ+σr）］，γzr=τzrG

可以求出微元體各向應力σr，σθ，σzr，見式（8）～式（10）。

σθ=Eεθ+Eμεr+μ（1+μ）σz1-μ2（8）

σr=Eεr+Eμεθ+μ（1+μ）σz1-μ2（9）

τzr=Gγzr（10）

式中的σz，εr，εθ，γzr由式（5）～式（7）可知是已知量，直接代入上式即可。但計算公式很長，不宜手算，建議使用計算機程序處理。

2 程序化計算

根據上述多軸向應力解計算公式編程，給出已經運行通過的VisualBasic程序代碼清單，物理及幾何參量與程序賦值變量對照見表1，Co mmand3命令執行運算（包括σz的計算），Co mmand4命令執行退出計算（表1）。

Dim A As Double， B As Double， F As Double， m As Double， n As Double， Is1 As Double， Is2 As Double， Is3 As Double， r As Double， z As Double， v As Double， l As Double， Q As Double， xgm As Double， r1 As Double， E As Double， G As Double， er As Double， e0 As Double， xgmr As Double， xgm0 As Double， tzr As Double， rzr As Double

Private Sub Co mmand3_Click（）

Dim zlf As Double， zlz As Double， rzlf As Double， rzlz As Double， rz2 As Double

Dim v22 As Double， nv4 As Double， v272 As Double， m6 As Double， m12 As Double， AB As Double

Dim jzt As Double， dhj As Double， syc As Double， urz As Double， wr As Double， dyc As Double

l=Val（Text1.Text）：r=Val（Text2.Text）：z=Val（Text3.Text）：v=Val（Text4.Text）：Q=Val（Text5.Text）

m=z/l：n=r/l

A=（n^2+（m-1）^2）^0.5：B=（n^2+（m+1）^2）^0.5：F=（n^2+m^2）^0.5

v22=（2*（2-v）*（4*m+1）-2*（1-2*v）*（1+m）*m^2/n^2）/B

nv4=（4*v*n^2*m+4*m^3-15*n^2*m-2*（5+2*v）*（m/n）^2*（m+1）^3+（m+1）^3）/B^3

v272=（2*（7-2*v）*m*n^2-6*m^3+2*（5+2*v）*（m/n）^2*m^3）/F^3

m6=（6*m*n^2*（n^2-m^2）+12*（m/n）^2*（m+1）^5）/B^5

m12=（12*（m/n）^2*m^5+6*m*n^2*（n^2-m^2））/F^5

AB=2*（2-v）*Log（（A+m-1）*（B+m+1）/（F+m）^2）

Is3=（2*（2-v）/A-v22-（2*（1-2*v）*m^3/n^2-8*（2-v）*m）/F-（m*n^2+（m-1）^3）/A^3-nv4-v272-m6+m12+AB）/（4*3.1415926*（1-v））

xgm=-Q*Is3/l^2

Label7.Caption=xgm

E=1：G=E/2/（1+v）

zlf=z-l：zlz=z+l：rzlf=Sqr（r^2+zlf^2）：rzlz=Sqr（r^2+zlz^2）：rz2=（r^2+z^2）^0.5

e0=Q*（l/rzlf-l/rzlz+Log（（rzlz+zlz）/（rz2+z））-Log（（rzlf-zlf）/（rz2-z）））/（8*3.1415926*G*（1-v）*l^2）

er=e0+r^2*Q*（-l/rzlf^3+l/rzlz^3+1/rzlz/（rzlz+zlz）-1/rz2/（rz2+z）+1/rz2/（rz2-z）-1/rzlf/（rzlf-zlf））/（8*3.1415926*G*（1-v）*l^2）

xgm0=（E*e0+E*v*er+v*（1+v）*xgm）/（1-v^2）

xgmr=（E*er+E*v*e0+v*（1+v）*xgm）/（1-v^2）

Label11.Caption=xgmr

Label12.Caption=xgm0

jzt=Q/（8*3.1415926*G*（1-v）*l^2）

dhj=r/rzlf+r/rzlz

syc=z*r/rzlf/（rzlf-zlf）-z*r/rzlz/（rzlz+zlz）+r*z/rz2/（rz2+z）-r*z/rz2/（rz2-z）

dyc=（zlz/rzlz+1）/（rzlz+zlz）-（zlf/rzlf-1）/（rzlf-zlf）+（z/rz2-1）/（rz2-z）-（z/rz2+1）/（rz2+z）

urz=r*jzt*（l*zlz/rzlz^3-l*zlf/rzlf^3+dyc）

wr=jzt*（3-4*v）*（dhj-2*r/rz2+syc）+jzt*（dhj-4*r/rz2+syc-r*z*zlz/rzlz^3-r*z*zlf/rzlf^3+（2*r*rzlz-r^3/rzlz）/（rzlz^2）+（2*r*rzlf-r^3/rzlf）/（rzlf^2））

rzr=urz+wr

tzr=G*rzr

Label14.Caption=tzr

End Sub

Private Sub Co mmand4_Click（）

Unload Me

End Sub

3 工程案例

設有一工程樁，樁長12 m，樁側土泊松比為0.35，三角形分布樁側摩阻力的合力為1 500 kN，求距離樁中心0.9 m，距離地表面6 m深度土體微元體的各向應力。

將上述數據輸入上述計算機程序，程序計算輸出結果：

σz=7.46863824175291 kPa

σr=0.421561243137828 kPa

σθ=-2.03487583200403 kPa

τzr=-21.6280564355353 kPa

本案例的明德林精確解：

σz=7.4686592470 kPa;σr=0.420367522097909 kPa

σθ=-2.031833913392 kPa;τzr=-21.7538514615124 kPa

4 結論

（1）由于本文基于Cerruti方法：將半無限體內一集中力的反應等價于兩個集中力共線對稱布置的全無限空間體的反應，也就是力學問題的結構正對稱、外力反對稱問題，即在兩個半無限體交界面上徑向正應力正對稱、豎向正應力反對稱、剪應力反對稱，故交界面上的豎向正應力只能為零，符合半無限體邊界條件，交界面上的剪應力可不為零，不符合半無限體邊界條件。因此本文得出的單樁在樁側線性增長阻力條件下樁外任何一點土體的環向正應力、徑向正應力解是合理的解析解；豎向正應力公式直接采用的是本問題的R.D.Mindlin積分解；本文得出的徑豎向剪應力公式對于分布力共線對稱布置的全無限空間體是合理的解析解，對于本問題的半無限空間體的水平界面的一定鄰域內是非合理的，應用于樁基問題的計算時這個鄰域范圍內的應力解應舍去，進一步分析會發現，距離水平界面及樁身較遠的區域本文給出的徑豎向剪應力的誤差會越來越小。

（2）本文得出的單樁在線性增長三角形分布樁側阻力條件下樁外任何一點土體的環向正應力、徑向正應力及豎徑向剪應力的理論解計算公式可為樁基礎的極限承載力計算提供必要的理論支撐，本文得出的同類問題的土體的豎向位移、水平位移計算公式亦可為樁基礎的沉降和變形計算提供必要的理論支撐。

（3）本文給出的VB機算程序可以使有關計算方便、快捷，可以直接進行快速工程應用。

（4）本文雖然只給出了一種力邊界條件下的單樁周圍巖土多軸向應力和位移的計算公式，但這些公式也可與其他力邊界條件同時作用下的情況進行疊加計算，綜合應用于實際復雜的樁基工況，例如群樁效應、成樁擠土穿刺效應等問題的研究。

（5）本文在建立有關力學模型時，把樁側阻力假設成柔性的線荷載，未考慮樁徑和樁體自身剛度對單樁周圍巖土效應的影響，這樣會使樁土界面附近的巖土顆粒各向應力的計算結果出現誤差，關于這方面的研究有待進一步的深入。

參考文獻

［1］ Raymond D. Mindlin. Force at a Point in the Interior of a Semi-infinite Solid［J］.Physics，1936，7（5）：195-202.

［2］ 宋東輝，徐晶.半無限空間體Mindlin問題的新解法及應用［J］.水電能源科學，2010，28（9）：101-103.

［3］ 中華人民共和國住房與城鄉建設部.建筑樁基技術規范： JGJ 94-2008［S］. 北京：中國建筑出版社，2013.

［4］ 薛守義.彈塑性力學［M］. 北京：中國建材工業出版社，2005.

［5］ 中華人民共和國住房與城鄉建設部.建筑地基基礎設計規范： GB 50007-2011［S］. 北京：中國建筑出版社，2012.

［6］ 中華人民共和國住房與城鄉建設部.建筑地基處理技術規范： JGJ 79-2012［S］. 北京：中國建筑出版社，2013.

［7］ 彭衛平.巖土工程參數分析評價［J］.中國勘察設計，2005（4）：38-40.

［8］ 葉建忠，周健.關于樁端阻力問題的分析與研究現狀［J］.建筑科學，2006，22（2）：64-68.