淺基礎地基荷載－變形特性的分析研究

羅 強，賈 虎，欒茂田

(1.南陽師范學院 土木建筑工程學院，河南 南陽 473061;2.大連理工大學 土木工程學院 巖土工程研究所，遼寧 大連 116024)

在主應力方向的旋轉過程中，主應力方向與塑性主應變增量方向不是完全重合的，即非共軸現(xiàn)象［1-3］。傳統(tǒng)的巖土彈塑性本構模型建立在共軸理論上，不能對非共軸現(xiàn)象進行合理地反映［4-6］。國內外研究學者相繼提出和完善了許多顆粒狀材料的非共軸理論，Bardet［7］在雙軸加載條件下研究了非共軸現(xiàn)象對砂土剪切帶的角度的影響。Papamichos 和Vardoulakis［8］在應變硬化模型的基礎上引入了非共軸理論，對砂土剪切帶的形成機理進行了研究。Hashiguchi 和Tsutsumi［9］采用非共軸模型研究了不排水雙軸壓縮試驗中的剪切帶的發(fā)展情況。黃茂松等［10-11］在彈塑性本構模型中引入非共軸塑性流動理論來描述非共軸現(xiàn)象。Yang 和Yu［12-13］將角點結構非共軸理論運用到有限元數(shù)值計算中，對砂土單剪試驗進行數(shù)值模擬，研究了非共軸現(xiàn)象對試樣的應力-應變關系的影響。

在許多實際工程問題中，例如:淺基礎作用下的地基承載力問題、地震或波浪荷載對海床的作用等，主應力方向旋轉引起的非共軸現(xiàn)象是不容忽視的［5，14］。Yang 和Yu［15-16］采用非共軸本構模型對淺基礎承載力與變形問題進行了有限元數(shù)值分析，其研究結果表明非共軸模型計算得到的基礎沉降高于共軸模型的結果。然而，上述研究工作大多是在理想彈塑性情況下進行的，很少在密砂的應變軟化情況下對非共軸現(xiàn)象進行研究。

在飽和密砂條件下，針對非共軸現(xiàn)象對淺基礎地基荷載-變形特性的影響進行研究。通過有限元二次開發(fā)，將非共軸本構模型應用到密砂單剪試驗的數(shù)值模擬中，對非共軸現(xiàn)象及其對應力-應變關系的影響進行研究;采用離心模型試驗方法，對圓形淺基礎作用下飽和密砂地基的荷載-變形特性進行研究;對離心模型試驗進行數(shù)值模擬，將試驗結果與數(shù)值計算結果進行對比，對非共軸模型的計算結果的合理性進行驗證。

圖1 非共軸、共軸塑性應變增量在屈服面上的關系Fig.1 Schematic illustration of non-coaxial and coaxial plastic strain rates

1 非共軸本構模型

基于角點結構理論［17］的非共軸模型是在Drucker-Prager 彈塑性理論的圓形屈服面上增加了一個與屈服面相切的非共軸塑性應變增量，它與傳統(tǒng)塑性應變增量是相正交的，如圖1 所示。

根據(jù)角點結構非共軸彈塑性理論可知:

采用Rokonuzzaman 和Toshinori Sakai［18］所提出的密砂應變軟化模型，該模型的屈服函數(shù)(f)和塑性勢函數(shù)(g)分別采用Mohr-Coulomb 函數(shù)和Drucker-Prager 函數(shù)的形式:

式(2)中的函數(shù)μ = μ ( ξp)為該模型的硬化參數(shù)，累積塑性應變?yōu)樗苄詰冊隽俊?/p>

式(3)中的洛德角θ 表示為

為了描述密砂材料的應變軟化特性，通過函數(shù)μ ( ξp)和α ( ξp)，該模型中的機動摩擦角φm和膨脹角φm分別從其峰值減小至殘余強度。

其中，φult為內摩擦角的峰值，φr為內摩擦角的殘余強度，φo為膨脹角的峰值。

塑性硬化模量Kp的表達形式:

彈性模量:

其中，eo為材料的初始孔隙比;ν 為泊松比;常數(shù)項Go為初始模量;Pa= 98 kPa 為大氣壓力。

其中，hnc為非共軸塑性模量，它為累積塑性應變ξp的函數(shù)，

其中，hnco為初始非共軸塑性模量;b1，b2為模型系數(shù)，其值分別為－16 和0.7。

式(19)表明:當skl和的方向重合時，非共軸塑性應變增量為0，那么，式(1)轉變?yōu)榕c傳統(tǒng)彈塑性理論相對應的共軸模型。

式(19)可以重新表達為如下形式

將式(1)、(18)、(21)和(22)結合在一起，得到非共軸彈塑性本構模型的表達形式

其中，Rij= ?g/?σij，lij= ?f/?σij，為非共軸模型的彈塑性剛度矩陣。與傳統(tǒng)共軸模型的表達形式相比，式(23)的右端增加了一項非共軸項。

一般來講，應變軟化模型的流動法則主要分為三種情況:相關聯(lián)流動法則、非關聯(lián)流動法則和塑性體積應變?yōu)榱?。在大部分情況下，土體的真實應力-應變關系分布在上述三種情況的計算結果中間。當采用非關聯(lián)流動法則時，屈服函數(shù)和塑性勢函數(shù)分別采用式(2)和(4)的形式;當采用相關聯(lián)流動法則時，屈服函數(shù)和塑性勢函數(shù)均采用式(2)的形式;當塑性體積應變?yōu)榱銜r，屈服函數(shù)和塑性勢函數(shù)分別采用式(2)和(4)的形式，并且

基于有限元軟件ABAQUS 的二次開發(fā)子程序UMAT，采用顯式積分算法和自動分步法，對非共軸本構模型進行數(shù)值積分［19］。

圖2 主應力方向和塑性主應變增量方向的旋轉Fig.2 Rotation of directions of principal stress and principal plastic strain rate

2 單剪試驗的數(shù)值模擬及結果分析

2.1 單剪試驗的數(shù)值模擬

有限元模型采用四邊形平面應變單元，其類型為八節(jié)點二次縮減積分單元。在模型頂邊施加水平位移邊界條件，同時，模型的左右兩邊保持直線狀態(tài)。模型底邊的豎向和水平方向位移均被固定，同時，豎向應力σyy施加在模型的頂面并保持不變。

在單剪過程中，由于剪應力τxy的作用，模型沿水平方向將會產生應力變化Δσxx，以及豎直方向的應變εyy。模型主應力方向的旋轉是由剪應力τxy的變化引起的。主應力方向和塑性主應變增量方向的旋轉如圖2 所示。

在圖2 中，虛線為變形后的狀態(tài)，實線為初始狀態(tài);α 為主應力或塑性主應變增量方向的旋轉角度為最大主應力或最大塑性主應變增量;為最小主應力或最小塑性主應變增量。

有限元數(shù)值模擬采用相對密實度Dr =80%的密砂，材料參數(shù):內摩擦角峰值φult= 40o，干密度ρd=1.593 g/cm3，泊松比ν = 0.3 ，靜止側壓力系數(shù)Ko=0.5，初始模量Go=125，內摩擦角殘余強度φr=36o，膨脹角峰值φo=20o，初始孔隙比eo=0.62。

2.2 數(shù)值計算結果分析

2.2.1 剪應力比-剪應變關系

當作用在試樣的豎向應力σyy=135 kPa 時，在流動法則的三種情況下，由共軸和非共軸模型計算所得到的剪應力比－剪應變關系曲線如圖3 所示。圖中豎向坐標為剪應力比，即τxy/σyy。非共軸模型中的hnco/G 分別取為0.2、0.4 和0.8。

圖3 數(shù)值計算結果Fig.3 Numerical calculation results

由圖3 可以發(fā)現(xiàn):1)在剪切變形的初期，共軸與非共軸模型的計算結果沒有顯著差異;隨著剪切變形的發(fā)展，非共軸模型計算結果的增長速度滯后于共軸模型計算結果的增加速度，兩者之間的差異逐漸顯著;當剪應力比達到峰值時，兩種模型計算結果之間的差異達到最大;隨著剪切變形的深入發(fā)展，試樣的抗剪強度由峰值向殘余強度發(fā)展，剪應力比逐漸減小，非共軸與共軸模型結果之間的差異逐漸減小;當試樣的抗剪強度達到臨界狀態(tài)以后，非共軸與共軸模型計算得到的剪應力比完全一致，非共軸現(xiàn)象的影響將不存在。2)隨著hnco/G 的減小，非共軸與共軸模型計算結果之間的差異逐漸增加。3)由關聯(lián)法則變化到塑性體積應變?yōu)榱銜r，應力-應變關系的應變軟化特性逐漸減弱，非共軸與共軸模型計算結果之間的差異逐漸減小。

2.2.2 主應力方向與塑性主應變增量方向的旋轉

當采用非關聯(lián)流動法則時，主應力方向和塑性主應變增量方向的旋轉變化如圖4 所示。

圖4 主應力方向與塑性主應變增量方向的旋轉Fig.4 Rotations of the orientations of major principal stress and plastic strain rate

由圖4 可以發(fā)現(xiàn):1)在共軸模型的計算結果中，主應力方向和塑性主應變增量方向是完全重合的。2)在非共軸模型的計算結果中，在剪切變形的初期，主應力方向的增長趨勢滯后于塑性主應變增量方向的增長趨勢，前者的旋轉角度要低于后者的角度;隨著剪切變形的增加，主應力方向的增長趨勢領先于塑性主應變增量方向的增長趨勢，前者的旋轉角度高于后者的角度;在剪切變形的后期，塑性主應變增量方向與主應力方向逐漸重合，非共軸現(xiàn)象逐漸消失。3)當hnco/G 較小時，非共軸現(xiàn)象比較明顯;隨著hnco/G 的增加，非共軸現(xiàn)象逐漸減弱。

2.3 數(shù)值計算結果與試驗結果的對比

在σyy=135 kPa 的條件下，Roscoe［4］研究了密砂試樣在單剪試驗過程中的非共軸現(xiàn)象。試驗得到的剪應力比-剪應變關系以及主應力方向和塑性主應變增量方向的旋轉趨勢如圖5 所示。

將圖3 的數(shù)值計算結果與圖5 的試驗結果進行對比，可以發(fā)現(xiàn):1)在試驗和數(shù)值計算結果中，均能夠觀測到非共軸現(xiàn)象。2)當選取比較合適的hnco/G 時，非共軸模型的計算結果與試驗結果比較接近。

圖5 Roscoe 的試驗結果Fig.5 Results of Roscoe’s test

3 淺基礎地基荷載-變形特性的數(shù)值分析及離心模型試驗研究

3.1 數(shù)值分析

圓形淺基礎的直徑D=1 m，基礎埋深Df與D 的比值Df/D 分別為0.0 和0.2。采用軸對稱單元建立有限元模型，為了模擬淺基礎與地基之間載荷的傳遞和相對滑動，淺基礎與地基在接觸面上采用摩擦接觸對算法。

采用Dr=80%的密砂，材料參數(shù)與前述單剪試驗中的參數(shù)相同。數(shù)值計算采用非關聯(lián)流動法則，非共軸模型中的hnco/G 分別取0.1、0.2 和1.0。

3.1.1 主應力方向旋轉及非共軸現(xiàn)象

針對淺基礎下方第一層土體單元，研究主應力方向的旋轉與單元水平位置之間的關系，如圖6 所示，其中，S 表示土體單元的水平位置(土體單元中心與淺基礎中心之間的距離)，V 表示基礎沉降。

圖6 土體單元的主應力方向旋轉Fig.6 Orientations of major principal stress of soil elements

由圖6 可知:主應力方向隨著基礎沉降的增加而逐漸增長到極值。隨著土體單元水平位置的增加，主應力方向的旋轉角度越來越大，淺基礎邊緣下方土體單元(S=D/2)的主應力方向旋轉角度最大。

當Df/D=0.0 和hnco/G=0.1 時，對不同水平位置處土體單元的非共軸現(xiàn)象進行研究，如圖7 所示。

圖7 土體單元的非共軸現(xiàn)象Fig.7 Non-coaxial phenomenon of soil element

由圖7 可知:1)主應力方向的旋轉變化滯后于塑性主應變增量方向的旋轉變化，兩者之間的差異隨著基礎沉降的增加而逐漸減小。2)隨著土體單元水平位置的增加，非共軸現(xiàn)象逐漸明顯。

3.1.2 非共軸現(xiàn)象對地基荷載-變形特性的影響

選取淺基礎下方第一層土體單元豎向應力的平均值作為地基豎向荷載，其與基礎沉降之間的關系如圖8 所示。

圖8 豎向荷載-基礎沉降關系Fig.8 Relationship between vertical load and settlement

由圖8 可知:1)在地基變形的初期，非共軸與共軸模型計算結果之間的差異不明顯;隨著基礎沉降的增加，非共軸模型所得到的豎向荷載的增長趨勢滯后于共軸模型結果的增長趨勢，兩者之間的差異逐漸增加;當?shù)鼗Q向荷載達到極值時，非共軸模型所得到的極值要低于共軸模型的計算結果;在地基變形的后期，非共軸模型與共軸模型得到的豎向荷載均由極值逐漸降低，兩者之間的差異逐漸減小。2)非共軸現(xiàn)象減緩了豎向荷載向其極值增長的速度，并降低了豎向荷載的極值。如果采用基礎沉降來確定地基豎向荷載，那么，非共軸模型的計算結果要低于共軸模型的結果。3)隨著hnco/G 的增加，兩種模型計算結果之間的差異逐漸減小。

3.2 離心模型試驗研究

采用大連理工大學的國內首臺土工鼓式離心機GT450/1.4，對圓形淺基礎作用下飽和密砂地基的荷載-變形特性進行研究。離心機試驗設備的整體構造如圖9 所示，鼓槽尺寸為1.4 m(直徑)×0.35 m(豎向寬度)×0.27 m(徑向深度)。鼓槽的總容量為450 g-t，質量為1 300 kg，體積為0.335 m3。鼓槽最大轉速為875 rpm，外側和內側最大離心加速度分別為600g 和369g，最大允許不平衡力為67.1 kN，最大允許徑向集中荷載10 kN。在鼓槽的直徑兩端各固定一個模型箱(如圖10 所示)，模型箱的尺寸為200 mm ×200 mm ×280 mm(寬×深×高)。

圖10 模型箱安置在鼓槽內Fig.10 Installation of model box with sand

在1g 狀態(tài)下，采用砂雨法將干砂均勻地灑入模型箱中，砂樣的深度為120 mm。然后，將模型箱放入水箱中，保證模型箱外水位高于箱內透水板約30 mm。在水壓差的作用下，模型箱外水流滲透進模型箱底部。通過干砂試樣顆粒間的毛細效應，水流被均勻地吸入砂樣中，可簡稱該過程為毛細滲透過程，具體裝樣過程詳見文獻［20］。當毛細滲透過程結束后，將模型箱從水箱中取出并豎直放置，箱內砂樣的初始狀態(tài)沒有發(fā)生變化，如圖10 所示。在毛細效應的作用下，砂樣顆粒緊緊吸附在一起，砂樣的初始狀態(tài)不受擾動。此時，砂樣的飽和度是比較低的，一般只能達到60% ～70%。

在離心模型試驗中，淺基礎的Df/D 分別為0.0 和0.5。淺基礎模型的直徑Dm分別為25 mm 和30 mm，通過變化離心加速N，可以得到不同的淺基礎實際直徑D，即D=N×Dm。

通過常規(guī)室內測試試驗，得到砂土的基本參數(shù):土粒相對密度Gs=2.627;顆粒尺寸d10=0.11 mm，d50=0.17 mm，d90=0.28 mm;不均勻系數(shù)Cu=1.727;最大、最小孔隙比為emax=0.913，emin=0.583;最大、最小干密度為1.66 g/cm3和1.37 g/cm3。砂土材料的相對密實度Dr =80%，干密度為1.593 g/cm3，浮重度為9.9 kN/m3。

淺基礎模型的寬度Dm的最小尺寸為25 mm，試驗材料的d50=0.17 mm，Dm/d50的數(shù)值為147;當Dm/d50＞30 時，粒徑效應對試驗結果的影響是可以忽略的［21］。試驗中的加載速度控制為0.01 mm/s。

試驗所得到的豎向荷載-基礎沉降關系如圖11 所示，在圖11(a)中若干試驗工況重復進行。

由圖11 可知:1)當離心機試驗工況重復進行時，所得到的豎向荷載-位移關系曲線比較接近，表明試驗方法具有良好的可重復性。2)試驗所得到的豎向荷載-位移關系曲線均有明顯的拐點，具有整體剪切破壞的特點。當豎向荷載-位移關系曲線出現(xiàn)拐點以后，隨著基礎沉降的增加，豎向荷載的變化不明顯，豎向荷載-位移關系曲線近似呈現(xiàn)水平狀態(tài)。3)拐點處的豎向荷載可以被視為地基承載力。隨著基礎寬度和埋深的增加，地基承載力也逐漸增加。

圖11 豎向荷載-基礎沉降關系Fig.11 Curves of vertical load-settlement

3.3 數(shù)值分析與模型試驗的結果對比

對離心模型試驗進行有限元數(shù)值模擬，將計算結果與試驗結果進行對比，如圖12、13 所示。

圖12 試驗結果與數(shù)值計算結果對比Df/D=0.0Fig.12 Comparisons between results of centrifuge tests and numerical calculations Df/D=0.0

由圖12 和13 可知:1)離心模型試驗和數(shù)值計算所得到的豎向荷載-位移關系曲線均有明顯的拐點，在拐點之后，豎向荷載-位移關系曲線呈現(xiàn)近似水平，應變軟化現(xiàn)象并不明顯。2)在地基變形初期，例如V/D ＜0.05，數(shù)值計算與離心模型試驗所得到的豎向荷載-位移關系曲線非常接近，非共軸模型與共軸模型的計算結果之間的差異不明顯;隨著基礎沉降的增加，非共軸與共軸模型計算結果之間的差異越來越明顯，離心模型試驗結果分布在非共軸模型的計算結果中間;當豎向荷載-位移關系曲線達到拐點以后，離心模型試驗所得到的地基承載力低于共軸模型的計算結果，分布在非共軸模型的計算結果中間。3)當選取比較合理的hnco/G 時，非共軸模型的數(shù)值計算結果與離心模型試驗結果比較接近。

圖13 試驗結果與數(shù)值計算結果對比Df/D=0.5Fig.13 Comparisons between results of centrifuge tests and numerical calculations Df/D=0.5

4 結 語

通過有限元二次開發(fā)，將非共軸本構模型應用到密砂單剪試驗的數(shù)值模擬中，對非共軸現(xiàn)象及其對應力-應變關系的影響進行研究;對圓形淺基礎作用下飽和密砂地基的荷載-變形特性進行離心模型試驗研究;對離心模型試驗進行數(shù)值模擬，將試驗結果與數(shù)值計算結果進行對比，對非共軸模型的計算結果的合理性進行驗證。主要得到以下結論:

1)在單剪試驗中，非共軸模型計算得到的剪應力比的增長速度滯后于共軸模型計算結果的增加速度;當試樣的抗剪強度達到臨界狀態(tài)以后，非共軸與共軸模型計算結果完全一致，非共軸現(xiàn)象的影響將消失。

2)對淺基礎地基荷載-變形特性進行數(shù)值分析時，淺基礎下方土體單元的主應力方向旋轉和非共軸現(xiàn)象比較顯著。非共軸模型所得到的豎向荷載的增長趨勢滯后于共軸模型結果的增長趨勢，非共軸模型所得到的地基承載力要低于共軸模型的計算結果。

3)離心模型試驗方法具有良好的可重復性。試驗所得到的豎向荷載-位移關系曲線均有明顯的拐點，在拐點之后的應變軟化現(xiàn)象并不明顯。

4)非共軸現(xiàn)象對淺基礎地基荷載-變形特性具有顯著的影響;當選取比較合理的非共軸塑性模型時，非共軸本構模型的計算結果與離心模型試驗結果比較接近。

［1］Roscoe K H，Bassett R H，Cole，E R.Principal axes observed during simple shear of a sand［C］//Proceedings of the Geotechnique.1967:231-237.

［2］Roscoe K H.The influence of strains in soil mechanics［J］.Geotechnique，1970，20(2):129-170.

［3］Oda M，Konishi J.Microscopic deformation mechanism of granular material in simple shear［J］.Soils and Foundations，1974，14(4):25-38.

［4］Madsen O S.Wave-induced pore pressures and effective stresses in a porous bed［J］.Geotechnique，1978，28(4):377-393.

［5］Ishihara K，Towhata I.Sand response to cyclic rotation of principal stress directions as induced by wave loads［J］.Soils and Foundations，1983，23(4):11-26.

［6］Yamamoto T，Koning H L，Spellmeigher H.On the response of a poro-elastic bed to water waves［J］.Journal of Fluid Mechanics，1978，87(1):193-206.

［7］Bardet J P.Orientation of shear bands in frictional soils［J］.Journal of Engineering Mechanics，1991，117(7):1466-1484.

［8］Papamichos E，Vardoulakis I.Shear band formation in sand according to non-coaxial plasticity model［J］.Geotechnique，1995，45(5):649-661.

［9］Hashiguchi K，Tsutsumi S.Shear band formation analysis in soils by the subloading surface model with tangential stress rate effect［J］.International Journal of Plasticity，2003，19(10):1651-1677.

［10］黃茂松，孫海忠，錢建固.粗粒土的非共軸性及其離散元數(shù)值模擬［J］.水利學報，2010，41(2):172-181.(HUANG Mao-song，SUN Hai-zhong，QIAN Jian-gu.Non-coaxial behavior of coarse granular aggregates simulated by DEM［J］.Shuili Xuebao，2010，41(2):172-181.(in Chinese))

［11］扈 萍，黃茂松，錢建固，等.砂土非共軸特性的本構模擬［J］.巖土工程學報，2009，31(5):793-798.(HU Ping，HUANG Mao-song，QIAN Jian-gu，et al.Non-coaxial plasticity constitutive modeling of sands［J］.Chinese Journal of Geotechnical Engineering，2009，31(5):793-798.(in Chinese))

［12］Yang Y M，Yu H S.Numerical simulations of simple shear with non-coaxial models［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2006，30(1):1-19.

［13］Yang Y M，Yu H S.A non-coaxial critical state soil model and its application to simple shear simulations［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2006，30(13):1369-1390.

［14］Manzari M T，Dafalias Y F.A critical state two-surface plasticity model for sands［J］.Geotechnique，1997，47(2):255-272.

［15］Yang Y M，Yu H S.Application of a non-coaxial soil model in shallow foundations［J］.Geomechanics and Geoengineering:An International Journal，2006，1(2):139-150.

［16］Yang Y M，Yu H S.Numerical aspects of non-coaxial model implementations［J］.Computers and Geotechnics，2010，37(1):93-102.

［17］Rudnicki J W，Rice J R.Conditions for the localisation of deformation in pressure-sensitive dilatant materials［J］.Journal of Mechanics and Physics of Solids，1975，23(6):371-394.

［18］Rokonuzzaman M D，Sakai T.Model tests and 3D finite element simulations of uplift resistance of shallow rectangular anchor foundations［J］.International Journal of Geomechanics，2012，12(2):1-24.

［19］羅 強，王忠濤，欒茂田，等.非共軸本構模型在地基承載力數(shù)值計算中若干影響因素的探討［J］.巖土力學，2011，32(S1):732-737.(LUO Qiang，WANG Zhong-tao，LUAN Mao-tian，et al.Factors analysis of non-coaxial constitutive model’s application to numerical analysis of foundation bearing capacity［J］.Rock and Soil Mechanics，2011，32(S1):732-737.(in Chinese))

［20］LUO Qiang，LUAN Mao-tian，YANG Yun-ming，et al.Numerical analysis and centrifuge modeling of shallow foundations［J］.China Ocean Engineering，2014，28(2):163-180.

［21］Ovesen N K.Centrifuge testing applied to bearing capacity problems of footings on sand［J］.Geotechnique，1975，25(3):394-401.