赵 莉，李亦宁

(1.东北大学 资源与土木工程学院， 辽宁 沈阳 110819；2.曼尼托巴大学 科学学院， 温尼伯格， MB R3T 2N2)

建(构)筑物地基一般由层状土构成的。研究表明，层状土的力学性状并非由各单一土层力学性状的简单组合。目前的地基静动力特性分析大多针对简单浅埋基础下的单一土层，而对于层状土整体力学特性的研究并不多见，开展层状土的静动力力学性能分析已经成为国内外岩土工程防灾减灾领域的重要课题[1-3]。应用单一岩土材料的本构模型，地基沉降计算的方法仅考虑弹性的范围。由于缺乏理论支撑，工程经验仍在工程设计中占主导地位，如现行规范法尽管是分层总和法的改进方法[4]，但仍然是一种半经验半理论的方法，仅适用于计算柱下扩展基础、条形基础等简单浅埋基础下有限点的地基沉降，应用局限性较大。随着高层和超高层建筑的不断涌现，大型复杂基础越来越多，基础埋深越来越大，要求进行变形验算和稳定性评价的土层种类和土层数量也越来越多。如果继续采用不合理的传统计算方法，将导致其计算结果的不准确性，同时给建(构)筑物的安全性带来不可忽视的影响，因此提高计算层状地基土变形的准确性成为研究学者关注的焦点[5-7]，迫切需要一个恰当的理论分析模型能合理预估层状土地基变形及其整体稳定性，进而建立性能化的抗失稳损伤设计与控制方法。

Li等[8]针对大型基础下的地基沉降问题，创建了基底平面网格系统，可以根据实际需要划分子域，采用该方法计算大型基础下基底平面任意点的地基沉降时，无需考虑作用在地基上的荷载类型和平均附加应力系数，克服了传统分层总和法只能单点计算基底沉降的不足，基底平面上所有子域形心处的沉降量都能同时被计算和体现出来。Díaz等[9]提出了一种新的表达式，用于计算两个不同变形能力的地基土上的地基沉降。Bezvolev[10]解决了空间不均匀非线性变形土层实际分布能力的计算问题，提出了一种考虑非线性土变形的程序。Sheinin等[11]提出了考虑复杂平面板、地基上的非均匀荷载以及可直接从测量数据分配的实际地基不均匀性的方法。Gazetas[12]提出了一种非均质弹性各向同性介质上条形基础动静力分析的数值分析方法。Enayat Masoud等[13]提到土体结构和孔隙比是土体渗透性变化的主要因素，从而控制土体变形，并提出了一个概念模型。李纯等[14]、修占国等[15]针对大型基础下层状土地基的变形问题，分别采用弹性半无限空间模型、修正的剑桥黏土模型、Mohr-Coulomb模型以及 Pastor-Zienkiewicz砂土模型等，建立了K0条件下和部分侧限条件下的层状土变形计算方法，并编写了解决饱和土和非饱和土变形问题的分析程序。同Mohr-Coulomb模型相比较，Drucker-Prager模型考虑了中间主应力的影响，对Mohr-Coulomb模型的屈服面函数作了适当的修改，运用圆锥形屈服面替代Mohr-Coulomb模型的六棱锥形屈服面。结果表明，由Drucker-Prager模型计算层状砂土的变形要比Mohr-Coulomb模型更加符合实际。

本文结合Drucker-Prager弹塑性本构模型，基于网格系统分层法推导建立一种非自由场层状砂土地基变形计算方法，利用该方法不但可以实现分析大型基础下复杂地基的变形性状，还为预测和预防建(构)筑物地基在发生变形失稳时突出灾害的发生提供一种新的分析方法。

1 层状土变形计算模型

1.1 平面网格子域系统构建

构建平面网格子域系统的主要目的是：(1)计算子域内基底有效附加应力的等效集中力P；(2)在宏观层面上建立基底荷载与土层变形之间的理论关系；(3)不但克服了传统方法只能用于计算小型简单基础的地基沉降，而不能解决大型基础下的复杂地基变形等缺陷，还能有效地解决地基内的复杂应力的空间传递与叠加问题。平面网格子域系统的实现原理及建立过程如下：

设基础底板长为L，宽为B。将基础下z=0处的地基平面被划分成m×n个子域，x方向m等份，y方向n等份。等效集中力Pij(i=1,2,…,m;j=1,2,…,n)的作用点为Aij(xi,yj)，且Aij(xi,yj)位于子域的形心上。平面网格子域系统图示见图1。

图1 平面网格系统图示

1.2 有效附加应力计算公式推导

(1)

对于未降水饱和地基土，有效应力增量Δσ′=Δσ-Δuw=0，则:

(2)

(3)

对于降水饱和地基土，基底平面上第i行、第j列子域内的有效附加应力为:

(4)

式中:Δσ为土的总应力增量,kPa；Δuw为孔隙水压力增量,kPa。

(5)

1.3 变形计算方法建立

在分析土体内的附加应力时，J.布辛奈斯克曾假设附加应力以土的内摩擦角形式向下扩散，并认为这种扩散作用一直以均匀应力下传。由于基础一般为钢筋混凝土结构，其刚度远大于土体的刚度，因此，在分析地基应力和计算荷载的分布规律时，又常常把基础底板视为刚体，将地基土近似地看成半无限空间弹性体。在此基础上，J.布辛奈斯克给出了弹性半空间内任意点M(x,y,z)处的应力和位移的解析解答，采用了公式(6)来计算竖向附加应力，如下:

(6)

式中:R表示的是任意点M到集中力作用点的距离，计算公示为R2=x2+y2+z2。

将地基有效压缩层分成l层，则第k分层上，点Mijk(xi,yj,zk)的沉降量计算公式为:

(7)

(8)

(9)

(10)

(11)

(12)

1.4 Drucker-prager本构方程

Drucker-prager屈服准则(下称“D-P准则”)与Mohr-Coulomb屈服准则(下称“M-C准则”)存在外角点外接圆、内角点外接圆、相关联内切圆、非相关联内切圆、等面积圆等匹配形式。

在三维应力空间中，D-P准则可以写成:

(13)

式中:aφ、k为参数系数；I1为应力张量的第一不变量;J2为应力偏张量的第二不变量。

M-C准则公式可表示为:

(14)

式中:c为黏聚力;φ为内摩擦角；θσ为Lode角，可以用J2和应力偏张量第三不变量J3的函数表示。

当M-C屈服准则中θσ=π/6时，将D-P屈服函数(见公式(13))与M-C屈服函数(见公式(14))联立，可得:

(15)

即为平面上屈服曲线通过摩尔-库伦不等角六角形外角点的外接圆，则外角点外接D-P屈服面的屈服函数为:

(16)

对式(16)求导可得:

(17)

又因为:

(18)

(19)

对于常量C1、C2，计算可得:

(20)

(21)

则

(22)

针对D-P本构模型，采用理想塑性加载准则，流动法则不考虑硬化规律，可确定所需的弹塑性应力-应变关系，即为D-P理想弹塑性增量本构关系，可以表示为:

dσ=Dedεe=D(dε-dεp)
=(De-Dp)dε=Depdε

(23)

式中:εe、εp分别为弹性应变和塑性应变；Dep为弹塑性刚度矩阵。

将(23)式展开得到:

(24)

在完全侧限K0条件下，土体压缩模量用增量表示为竖向应力增量与竖向应变增量之比，即:

(25)

按照选取的屈服准则、加载准则、流动准则、硬化规律关系，可求得最终简化的弹塑性刚度矩阵Dep，即:

(26)

2 程序编制与算法

SSBS的计算程序基于MATLAB语言平台，是针对解决层状土变形问题而自行编写的计算分析程序。该程序通过命令的设置，实现对于所计算的网格系统中平面尺寸、测试作用点的精确位置进行识别；生成基底初始平面网格系统；按照要求对初始网格系统进行网格细分和z坐标Lagrange二元插值，以生成子域等效集中力矩阵和(x,y)坐标子矩阵；实现基底有效附加应力的空间传递与计算过程以及计算各分层土的竖向变形和总变形值等。SSBS程序不但能计算拟建建筑物的基底变形，还能计算地表沉降及邻近建筑物地基受有效附加应力扩散的影响而产生的基底附加沉降等。

3 工程算例

以辽宁省文化艺术中心、省博物馆新馆为工程实例，根据实测的数据结果显示，地质岩性自上到下的土层分布由杂填土、砾砂、圆砾、砾砂、粉质黏土、圆砾、粉质黏土、粗砂、粉质黏土、粗砂组成。地基处于无降水影响的饱和状态，基础型式采用大型构造板扩展基础形式，基础埋深-6.46 m，地下水的稳定标高为-6.60 m。基础与地基剖面示意见图2，不同地基土的物理力学参数见表1。

表1 地基土的主要物理力学参数

图2 基础和地基剖面示意

3.1 基底附加应力实测

测试区域基底平面网格系统尺寸长L=18 m，宽B=9 m，结构柱网跨距9.0 m×9.0 m。为了便于测试工作的操作和不影响主体结构的施工进度，测试区域取横向两跨，竖向一跨。基底平面网格系统划分距离设置为2.25 m，测试点分布个数设置39个，其中30#失效，18#和37#的值为56.87 kPa、42.54 kPa。基底压力传感器布设如图3所示。

图3 基底压力传感器布设

将测得的基底接触应力值减去基底上覆土层的自重力便得到基底附加应力值。各监测点的基底接触应力测试结果见图4(a)；测试结果随时间变化情况见图4(b)。

图4 各监测点基底接触应力及随时间变化情况

图5 基底附加应力二元插值网线图

3.2 基于Drucker-Prager模型的计算结果

变形计算完全采用SSBS程序[11]来完成。主程序首先调用生成平面坐标子程序、网格系统的初始插值子程序和空间多维数组子程序，生成网格系统各子域等效集中力矩阵，计算基础底部的有效附加应力，而且可以得到所需深度处的应力值。所需的计算分层可以设置若干个，可以根据需要任意划分和给定，非常灵活。待输入不同土层的物理力学参数后，主程序自动辨识各土层的物理力学属性，自动调用程序库中恰当的材料本构关系分别计算完全侧限条件下和非完全侧限条件下各层土的竖向变形及基底总变形。由于算例中除了3个较薄的粉质黏土层外，大部分都是较密实的砂土层，故本文考虑采用Drucker-Prager本构模型计算分析受载层状土的变形性能。

图6为基底竖向变形的计算结果，最大竖向变形发生点坐标为(12.75,1.25)，最大竖向变形值为14.34 mm。而基底附加应力的最大值为157.24 kPa，相应受力点坐标为(9.00,0.00)。从表2可以看出，在地基有限计算深度范围内，地基最大竖向变形点和最大荷载作用点位置可能不同，计算深度越浅差异越大；计算深度越大，差异逐渐变小，这就是砂土地基变形工程实测值往往小于理论值的缘故。因此，传统方法仅计算基础中心点和柱下地基有效计算深度内的地基变形显然是不妥当的，容易造成工程事故。

图6 基底竖向变形

表2 各自然层的最大竖向变形计算值及位置

3.3 计算模型验证

为了进一步证明本文建立的大型基础下层状土变形计算模型的合理性，利用FLAC3D分析软件建立工程算例数值模型，分别调用软件中的M-C本构模型和D-P本构模型，输入相关参数和测试数据，进行数值模拟计算。在数值模型中，采用实体单元建模，其中间部分作用基底荷载，分布尺寸设置为18 m×9 m；z轴的上边界设置为基底，下边界为地基土分层最深处负19.04 m；边界条件的设置采用的是对于模型四周以及底部限制x、y、z方向的位移，模型的x、y方向的尺寸设置为荷载分布尺寸的5倍，达到工程实际的边界要求。

图7为本文计算模型及验证方法分别得到的基底竖向变形沿深度方向的变化曲线。可以看出，在地基同一深度平面上，由本文计算方法得到的基底竖向变形最小，而由分析软件模拟得到的变形值均偏大，这是因为本文计算采用的是理想弹塑性增量Drucker-Prager本构模型、而验证分析分别选择的是粘弹-双曲线Drucker-Prager塑性本构模型和通用弹塑性Mohr-Coulomb本构模型的缘故。由此表明，本文基于平面网格子域系统的层状土变形计算模型合理、实用。

图7 变形沿深度方向的变化曲线

3.4 计算结果比较分析

表3给出了分别由本文计算方法、验证方法1和验证方法2三种途径得到的基底最大变形计算结果和误差数据。

首先，将本文方法的计算结果与验证方法1的模拟结果进行比较后发现，最大竖向变形值绝对误差0.132 mm，相对误差为0.92%。发生最大竖向变形点坐标在x方向无偏差，在y方向的偏差为0.25 m，最大变形发生点位置基本一致。

其次，将本文方法的计算结果与验证方法2的模拟结果进行比较后发现，最大竖向变形值绝对误差1.182 mm，相对误差为7.61%。发生最大竖向变形点坐标在x方向仍然无偏差，在y方向的偏差为0.50 m，最大变形发生点位置偏差不大。

表3 基底最大变形计算结果及误差比较

为了更加直观地比较表3中三种方法得到的计算结果，图8为该三种方法得到的最大竖向变形发生位置沿x方向的变形曲线和沿y方向的变形曲线。由曲线的变化趋势可以看出，三种方法的计算结果变化趋势完全相同。

图8 各算法地基总变形分布曲线

基于以往工程实测和上述比较分析结果，充分证明了本文提供的大型基础下受载层状土变形计算方法及自编程序是合理和可靠的。

4 结 论

(1) 基底受力状态复杂，基底附加应力不一定呈均匀分布状态，因此，开展地基性能评估和基础设计时，不能简单地将基底荷载简化为均匀荷载。

(2) 基底最大竖向变形并不一定与基底最大应力作用点相重合。因此，基础设计时，将基础中心点的竖向变形作为基底最大竖向变形来控制基础截面尺寸的做法是不准确的。

(3) 考虑D-P本构模型中间主应力的影响，计算模型采用圆锥形屈服面，得到的计算结果要比采用其他类型的屈服面更加符合工程实际。

(4) 采用网格系统分层法，有效解决了大型基础下层状土地基变形计算问题，克服了采用单一的岩土本构模型在弹性范围内单点逐次计算沉降的困难。