浅基础地基荷载－变形特性的分析研究

罗 强，贾 虎，栾茂田

(1.南阳师范学院 土木建筑工程学院，河南 南阳 473061;2.大连理工大学 土木工程学院 岩土工程研究所，辽宁 大连 116024)

在主应力方向的旋转过程中，主应力方向与塑性主应变增量方向不是完全重合的，即非共轴现象［1-3］。传统的岩土弹塑性本构模型建立在共轴理论上，不能对非共轴现象进行合理地反映［4-6］。国内外研究学者相继提出和完善了许多颗粒状材料的非共轴理论，Bardet［7］在双轴加载条件下研究了非共轴现象对砂土剪切带的角度的影响。Papamichos 和Vardoulakis［8］在应变硬化模型的基础上引入了非共轴理论，对砂土剪切带的形成机理进行了研究。Hashiguchi 和Tsutsumi［9］采用非共轴模型研究了不排水双轴压缩试验中的剪切带的发展情况。黄茂松等［10-11］在弹塑性本构模型中引入非共轴塑性流动理论来描述非共轴现象。Yang 和Yu［12-13］将角点结构非共轴理论运用到有限元数值计算中，对砂土单剪试验进行数值模拟，研究了非共轴现象对试样的应力-应变关系的影响。

在许多实际工程问题中，例如:浅基础作用下的地基承载力问题、地震或波浪荷载对海床的作用等，主应力方向旋转引起的非共轴现象是不容忽视的［5，14］。Yang 和Yu［15-16］采用非共轴本构模型对浅基础承载力与变形问题进行了有限元数值分析，其研究结果表明非共轴模型计算得到的基础沉降高于共轴模型的结果。然而，上述研究工作大多是在理想弹塑性情况下进行的，很少在密砂的应变软化情况下对非共轴现象进行研究。

在饱和密砂条件下，针对非共轴现象对浅基础地基荷载-变形特性的影响进行研究。通过有限元二次开发，将非共轴本构模型应用到密砂单剪试验的数值模拟中，对非共轴现象及其对应力-应变关系的影响进行研究;采用离心模型试验方法，对圆形浅基础作用下饱和密砂地基的荷载-变形特性进行研究;对离心模型试验进行数值模拟，将试验结果与数值计算结果进行对比，对非共轴模型的计算结果的合理性进行验证。

图1 非共轴、共轴塑性应变增量在屈服面上的关系Fig.1 Schematic illustration of non-coaxial and coaxial plastic strain rates

1 非共轴本构模型

基于角点结构理论［17］的非共轴模型是在Drucker-Prager 弹塑性理论的圆形屈服面上增加了一个与屈服面相切的非共轴塑性应变增量，它与传统塑性应变增量是相正交的，如图1 所示。

根据角点结构非共轴弹塑性理论可知:

采用Rokonuzzaman 和Toshinori Sakai［18］所提出的密砂应变软化模型，该模型的屈服函数(f)和塑性势函数(g)分别采用Mohr-Coulomb 函数和Drucker-Prager 函数的形式:

式(2)中的函数μ = μ ( ξp)为该模型的硬化参数，累积塑性应变为塑性应变增量。

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

为了描述密砂材料的应变软化特性，通过函数μ ( ξp)和α ( ξp)，该模型中的机动摩擦角φm和膨胀角φm分别从其峰值减小至残余强度。

其中，φult为内摩擦角的峰值，φr为内摩擦角的残余强度，φo为膨胀角的峰值。

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

弹性模量:

其中，eo为材料的初始孔隙比;ν 为泊松比;常数项Go为初始模量;Pa= 98 kPa 为大气压力。

其中，hnc为非共轴塑性模量，它为累积塑性应变ξp的函数，

其中，hnco为初始非共轴塑性模量;b1，b2为模型系数，其值分别为－16 和0.7。

式(19)表明:当skl和的方向重合时，非共轴塑性应变增量为0，那么，式(1)转变为与传统弹塑性理论相对应的共轴模型。

式(19)可以重新表达为如下形式

将式(1)、(18)、(21)和(22)结合在一起，得到非共轴弹塑性本构模型的表达形式

其中，Rij= ∂g/∂σij，lij= ∂f/∂σij，为非共轴模型的弹塑性刚度矩阵。与传统共轴模型的表达形式相比，式(23)的右端增加了一项非共轴项。

一般来讲，应变软化模型的流动法则主要分为三种情况:相关联流动法则、非关联流动法则和塑性体积应变为零。在大部分情况下，土体的真实应力-应变关系分布在上述三种情况的计算结果中间。当采用非关联流动法则时，屈服函数和塑性势函数分别采用式(2)和(4)的形式;当采用相关联流动法则时，屈服函数和塑性势函数均采用式(2)的形式;当塑性体积应变为零时，屈服函数和塑性势函数分别采用式(2)和(4)的形式，并且

基于有限元软件ABAQUS 的二次开发子程序UMAT，采用显式积分算法和自动分步法，对非共轴本构模型进行数值积分［19］。

图2 主应力方向和塑性主应变增量方向的旋转Fig.2 Rotation of directions of principal stress and principal plastic strain rate

2 单剪试验的数值模拟及结果分析

2.1 单剪试验的数值模拟

有限元模型采用四边形平面应变单元，其类型为八节点二次缩减积分单元。在模型顶边施加水平位移边界条件，同时，模型的左右两边保持直线状态。模型底边的竖向和水平方向位移均被固定，同时，竖向应力σyy施加在模型的顶面并保持不变。

在单剪过程中，由于剪应力τxy的作用，模型沿水平方向将会产生应力变化Δσxx，以及竖直方向的应变εyy。模型主应力方向的旋转是由剪应力τxy的变化引起的。主应力方向和塑性主应变增量方向的旋转如图2 所示。

在图2 中，虚线为变形后的状态，实线为初始状态;α 为主应力或塑性主应变增量方向的旋转角度为最大主应力或最大塑性主应变增量;为最小主应力或最小塑性主应变增量。

有限元数值模拟采用相对密实度Dr =80%的密砂，材料参数:内摩擦角峰值φult= 40o，干密度ρd=1.593 g/cm3，泊松比ν = 0.3 ，静止侧压力系数Ko=0.5，初始模量Go=125，内摩擦角残余强度φr=36o，膨胀角峰值φo=20o，初始孔隙比eo=0.62。

2.2 数值计算结果分析

2.2.1 剪应力比-剪应变关系

当作用在试样的竖向应力σyy=135 kPa 时，在流动法则的三种情况下，由共轴和非共轴模型计算所得到的剪应力比－剪应变关系曲线如图3 所示。图中竖向坐标为剪应力比，即τxy/σyy。非共轴模型中的hnco/G 分别取为0.2、0.4 和0.8。

图3 数值计算结果Fig.3 Numerical calculation results

由图3 可以发现:1)在剪切变形的初期，共轴与非共轴模型的计算结果没有显著差异;随着剪切变形的发展，非共轴模型计算结果的增长速度滞后于共轴模型计算结果的增加速度，两者之间的差异逐渐显著;当剪应力比达到峰值时，两种模型计算结果之间的差异达到最大;随着剪切变形的深入发展，试样的抗剪强度由峰值向残余强度发展，剪应力比逐渐减小，非共轴与共轴模型结果之间的差异逐渐减小;当试样的抗剪强度达到临界状态以后，非共轴与共轴模型计算得到的剪应力比完全一致，非共轴现象的影响将不存在。2)随着hnco/G 的减小，非共轴与共轴模型计算结果之间的差异逐渐增加。3)由关联法则变化到塑性体积应变为零时，应力-应变关系的应变软化特性逐渐减弱，非共轴与共轴模型计算结果之间的差异逐渐减小。

2.2.2 主应力方向与塑性主应变增量方向的旋转

当采用非关联流动法则时，主应力方向和塑性主应变增量方向的旋转变化如图4 所示。

图4 主应力方向与塑性主应变增量方向的旋转Fig.4 Rotations of the orientations of major principal stress and plastic strain rate

由图4 可以发现:1)在共轴模型的计算结果中，主应力方向和塑性主应变增量方向是完全重合的。2)在非共轴模型的计算结果中，在剪切变形的初期，主应力方向的增长趋势滞后于塑性主应变增量方向的增长趋势，前者的旋转角度要低于后者的角度;随着剪切变形的增加，主应力方向的增长趋势领先于塑性主应变增量方向的增长趋势，前者的旋转角度高于后者的角度;在剪切变形的后期，塑性主应变增量方向与主应力方向逐渐重合，非共轴现象逐渐消失。3)当hnco/G 较小时，非共轴现象比较明显;随着hnco/G 的增加，非共轴现象逐渐减弱。

2.3 数值计算结果与试验结果的对比

在σyy=135 kPa 的条件下，Roscoe［4］研究了密砂试样在单剪试验过程中的非共轴现象。试验得到的剪应力比-剪应变关系以及主应力方向和塑性主应变增量方向的旋转趋势如图5 所示。

将图3 的数值计算结果与图5 的试验结果进行对比，可以发现:1)在试验和数值计算结果中，均能够观测到非共轴现象。2)当选取比较合适的hnco/G 时，非共轴模型的计算结果与试验结果比较接近。

图5 Roscoe 的试验结果Fig.5 Results of Roscoe’s test

3 浅基础地基荷载-变形特性的数值分析及离心模型试验研究

3.1 数值分析

圆形浅基础的直径D=1 m，基础埋深Df与D 的比值Df/D 分别为0.0 和0.2。采用轴对称单元建立有限元模型，为了模拟浅基础与地基之间载荷的传递和相对滑动，浅基础与地基在接触面上采用摩擦接触对算法。

采用Dr=80%的密砂，材料参数与前述单剪试验中的参数相同。数值计算采用非关联流动法则，非共轴模型中的hnco/G 分别取0.1、0.2 和1.0。

3.1.1 主应力方向旋转及非共轴现象

针对浅基础下方第一层土体单元，研究主应力方向的旋转与单元水平位置之间的关系，如图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 时，对不同水平位置处土体单元的非共轴现象进行研究，如图7 所示。

图7 土体单元的非共轴现象Fig.7 Non-coaxial phenomenon of soil element

由图7 可知:1)主应力方向的旋转变化滞后于塑性主应变增量方向的旋转变化，两者之间的差异随着基础沉降的增加而逐渐减小。2)随着土体单元水平位置的增加，非共轴现象逐渐明显。

3.1.2 非共轴现象对地基荷载-变形特性的影响

选取浅基础下方第一层土体单元竖向应力的平均值作为地基竖向荷载，其与基础沉降之间的关系如图8 所示。

图8 竖向荷载-基础沉降关系Fig.8 Relationship between vertical load and settlement

由图8 可知:1)在地基变形的初期，非共轴与共轴模型计算结果之间的差异不明显;随着基础沉降的增加，非共轴模型所得到的竖向荷载的增长趋势滞后于共轴模型结果的增长趋势，两者之间的差异逐渐增加;当地基竖向荷载达到极值时，非共轴模型所得到的极值要低于共轴模型的计算结果;在地基变形的后期，非共轴模型与共轴模型得到的竖向荷载均由极值逐渐降低，两者之间的差异逐渐减小。2)非共轴现象减缓了竖向荷载向其极值增长的速度，并降低了竖向荷载的极值。如果采用基础沉降来确定地基竖向荷载，那么，非共轴模型的计算结果要低于共轴模型的结果。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 状态下，采用砂雨法将干砂均匀地洒入模型箱中，砂样的深度为120 mm。然后，将模型箱放入水箱中，保证模型箱外水位高于箱内透水板约30 mm。在水压差的作用下，模型箱外水流渗透进模型箱底部。通过干砂试样颗粒间的毛细效应，水流被均匀地吸入砂样中，可简称该过程为毛细渗透过程，具体装样过程详见文献［20］。当毛细渗透过程结束后，将模型箱从水箱中取出并竖直放置，箱内砂样的初始状态没有发生变化，如图10 所示。在毛细效应的作用下，砂样颗粒紧紧吸附在一起，砂样的初始状态不受扰动。此时，砂样的饱和度是比较低的，一般只能达到60% ～70%。

在离心模型试验中，浅基础的Df/D 分别为0.0 和0.5。浅基础模型的直径Dm分别为25 mm 和30 mm，通过变化离心加速N，可以得到不同的浅基础实际直径D，即D=N×Dm。

通过常规室内测试试验，得到砂土的基本参数:土粒相对密度Gs=2.627;颗粒尺寸d10=0.11 mm，d50=0.17 mm，d90=0.28 mm;不均匀系数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的数值为147;当Dm/d50＞30 时，粒径效应对试验结果的影响是可以忽略的［21］。试验中的加载速度控制为0.01 mm/s。

试验所得到的竖向荷载-基础沉降关系如图11 所示，在图11(a)中若干试验工况重复进行。

由图11 可知:1)当离心机试验工况重复进行时，所得到的竖向荷载-位移关系曲线比较接近，表明试验方法具有良好的可重复性。2)试验所得到的竖向荷载-位移关系曲线均有明显的拐点，具有整体剪切破坏的特点。当竖向荷载-位移关系曲线出现拐点以后，随着基础沉降的增加，竖向荷载的变化不明显，竖向荷载-位移关系曲线近似呈现水平状态。3)拐点处的竖向荷载可以被视为地基承载力。随着基础宽度和埋深的增加，地基承载力也逐渐增加。

图11 竖向荷载-基础沉降关系Fig.11 Curves of vertical load-settlement

3.3 数值分析与模型试验的结果对比

对离心模型试验进行有限元数值模拟，将计算结果与试验结果进行对比，如图12、13 所示。

图12 试验结果与数值计算结果对比Df/D=0.0Fig.12 Comparisons between results of centrifuge tests and numerical calculations Df/D=0.0

由图12 和13 可知:1)离心模型试验和数值计算所得到的竖向荷载-位移关系曲线均有明显的拐点，在拐点之后，竖向荷载-位移关系曲线呈现近似水平，应变软化现象并不明显。2)在地基变形初期，例如V/D ＜0.05，数值计算与离心模型试验所得到的竖向荷载-位移关系曲线非常接近，非共轴模型与共轴模型的计算结果之间的差异不明显;随着基础沉降的增加，非共轴与共轴模型计算结果之间的差异越来越明显，离心模型试验结果分布在非共轴模型的计算结果中间;当竖向荷载-位移关系曲线达到拐点以后，离心模型试验所得到的地基承载力低于共轴模型的计算结果，分布在非共轴模型的计算结果中间。3)当选取比较合理的hnco/G 时，非共轴模型的数值计算结果与离心模型试验结果比较接近。

图13 试验结果与数值计算结果对比Df/D=0.5Fig.13 Comparisons between results of centrifuge tests and numerical calculations Df/D=0.5

4 结 语

通过有限元二次开发，将非共轴本构模型应用到密砂单剪试验的数值模拟中，对非共轴现象及其对应力-应变关系的影响进行研究;对圆形浅基础作用下饱和密砂地基的荷载-变形特性进行离心模型试验研究;对离心模型试验进行数值模拟，将试验结果与数值计算结果进行对比，对非共轴模型的计算结果的合理性进行验证。主要得到以下结论:

1)在单剪试验中，非共轴模型计算得到的剪应力比的增长速度滞后于共轴模型计算结果的增加速度;当试样的抗剪强度达到临界状态以后，非共轴与共轴模型计算结果完全一致，非共轴现象的影响将消失。

2)对浅基础地基荷载-变形特性进行数值分析时，浅基础下方土体单元的主应力方向旋转和非共轴现象比较显著。非共轴模型所得到的竖向荷载的增长趋势滞后于共轴模型结果的增长趋势，非共轴模型所得到的地基承载力要低于共轴模型的计算结果。

3)离心模型试验方法具有良好的可重复性。试验所得到的竖向荷载-位移关系曲线均有明显的拐点，在拐点之后的应变软化现象并不明显。

4)非共轴现象对浅基础地基荷载-变形特性具有显著的影响;当选取比较合理的非共轴塑性模型时，非共轴本构模型的计算结果与离心模型试验结果比较接近。

［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］黄茂松，孙海忠，钱建固.粗粒土的非共轴性及其离散元数值模拟［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］罗 强，王忠涛，栾茂田，等.非共轴本构模型在地基承载力数值计算中若干影响因素的探讨［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.