李纯, 修占国, 王斐笠
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2016-10-24
基金项目：国家自然科学基金重点资助项目(90915005)；辽宁省自然科学基金资助项目(2013010311)。
作者简介：李纯(1971-)，男，安徽潜山人，东北大学副教授。

摘要：基于饱和土的相关理论, 推导了大型基础下层状地基的沉降计算模型.考虑不同的土性, 根据修正剑桥模型和摩尔-库伦模型, 推导了K₀条件下且能够考虑其塑性变形的压缩模量.基于平面网格子域法结合分层总和法建立了一个完整的、更加贴近工程实际的大型基础下含软弱夹层的层状砂土地基沉降计算模型.研究结果表明:该模型的计算结果与模拟结果很贴近, 且能够突出体现出软弱夹层的变形特性.研究结论可为含软弱夹层的大型基础地基沉降计算提供直接参考依据.
关键词：层状砂土地基软弱夹层平面网格子域本构关系沉降
Calculation on Deformation of Layered Sandy with Weak Interlayer Under Large Foundation
LI Chun, XIU Zhan-guo, WANG Fei-li
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: XIU Zhan-guo, E-mail: xiuzhanguo109@126.com
Abstract: A settlement calculation models of layered subgrade under large foundation was established based on theories of saturated soil. Taking the different soil properties and the modified Cam-Clay and Mohr-Coulomb models into account, the compressive modulus were derived by considering the plastic deformations under the condition of K₀. Based on the subfield of planar mesh and the layer-wise summation method, a settlement calculation model of layered sandy foundation base with weak interlayer under large foundation was established, which is more close to the engineering practice. The results showed that the calculation results by the suggested model are comparatively close to the numerical simulation results, and can clearly reflect the characteristic of deformation in weak interlayer. The suggested calculation model can provide basis to the settlement calculation of layered sandy foundation base with weak interlayer under large foundation.
Key Words: layered sandy foundation baseweak interlayersubfield of planar meshconstitutive relationsettlement
当前, 用于计算地基变形的方法很多, 分层总和法是计算地基最终沉降量最常用的方法.现行规范法在分层总和法的基础上进行了改进, 使沉降计算结果在较大程度上更加接近实测值, 基本上可解决基础形状简单、尺寸不大的中小型基础沉降问题.然而, 对于大型、复杂、重要的基础, 如筏型基础、构造板扩展基础等, 误认为最大荷载作用点就是最大沉降发生点, 将最大荷载作用点作为最大变形控制点存在安全隐患^[1].
大型基础的变形(沉降)问题一直是国内外学者关注的热点, 其影响因素繁多, 合理的计算过程过于复杂, 没有得到广泛的工程应用.文献[2]提出了一种基于两种不同土性的大型基础沉降的简便计算方法.文献[3]提出了一种层状地基土的弹塑性计算方法, 并编写了相应的计算程序.文献[4]针对不均匀多层地基土提出了一种解析的数值计算方法.文献[5]从定性的角度分析了在连续扩展基础下土体纵向变化对地基变形的影响.文献[6]针对大型基础下层状地基变形实测, 提出一种基于网格系统分层法结合分层总和法的地基沉降计算方法.
层状砂土地基是大型基础下地基土的典型类型, 基底下成层土分布极其不均匀, 土性相差悬殊.层状砂土地基中往往含有软弱夹层, 软弱夹层的存在将极大程度上减弱地基的整体刚度, 产生较大变形^[7-8].因此, 在计算中考虑不同土性，选用适用于该土性的本构模型是准确计算的保证.本文将依据理论推导和工程算例, 建立一种完整的、更接近工程实际的地基沉降计算模型.
1 数学计算模型1.1 基于平面网格子域法的有效附加应力基底下有效附加应力由有效应力增量和附加应力共同组成, 有效应力的产生主要是由于施工过程中地基土的容重和孔隙水压力变化产生的.本文针对大型基础下未降水饱和地基土建立计算模型, 土的总应力增量Δσ和孔隙水压力增量Δu均为零.因此, 有效应力增量也为零.也就是说, 只有建筑物施加给地基的荷载在变化, 有效附加应力在数值上等于最大附加应力增量.图 1为未降水基坑饱和地基土的有关情况.
图 1(Fig. 1)

图 1 未降水基坑全饱和地基土示意图Fig.1 Sketch map for saturation soils before dewatering

假定地基土为均质单一土层, 选择M₁, M₂和M₃特征点分析, 其中M₁代表建筑物室外自然地面处土体单元, M₂代表水位线处的土体单元, M₃代表水位线下h₂深度处的土体单元.M₁M₂区段为非饱和土, 其天然容重为γ; M₂M₃区段为饱和土, 土体容重为γ_sat.M₁, M₂和M₃处的有效应力分别为0, γh₁和γh₁+γ_sath₂－γ_wh₂, 由于土体容重和饱和度均未发生改变, 有效应力增量为零, 有效附加应力(σ₀′)在数值上等于最大有效附加应力增量.M₁, M₂和M₃处的有效附加应力分别为0, Δσ_0max^M₂, Δσ_0max^M₃.其中, Δσ_0max^M₂, Δσ_0max^M₃分别为基础底板与地基接触面处的最大附加应力增量和基础底板以下h₂深度处的最大附加应力增量.
基于布辛奈斯克给出的弹性半空间任一点M(x, y, z)处的应力(σ_x, σ_y, σ_z, τ_xy, τ_yz, τ_zx), 计算地基土竖向附加应力σ_z的公式为

(1)

式中, R为M点与集中力P作用点O之间的距离.
在地基中，基底有效附加应力(σ′₀)与竖直方向成φ角(内摩擦角)连续均匀下传.在基底平面上, 假定基础底板长L, 宽B.将基础下地基划分为m×n个子域, 作用于每个子域的有效附加应力可视为均布分布, 可用作用于中心处的1个等效集中力P_ij(i=1, 2, …, m; j=1, 2, …, n)代替.设基底平面网格系统上存在网格中心点A_ij(i=1, 2, …, m; j=1, 2, …, n), A_ij与等效集中力P_ij的作用点重合, A_ij点下深度为h₃(h₃为计算地基变形时的总分层数)标高处的水平面上存在竖向对应网格中心点A′_ij, 点A′_ij与等效集中力P′_ij的作用点重合.这样的等效集中力分量共计m×n个, 所有等效集中力分量的叠加结果则为点A′_ij的等效集中力.重复上述过程, 可得到新的目标计算面上各相应点的等效集中力.此过程可持续到地基变形有效压缩层深度为止.子域内等效集中应力分布情况如图 2所示.
图 2(Fig. 2)

图 2 等效集中力的分布Fig.2 Distribution of equivalent concentrated force

1.2 沉降计算模型的理论推导地基表面任意点A_ij以下深度z处的分层平面上相应各点的有效附加应力σ′_0z为所有P_ij单独作用时, 在该分层面上对应点处所产生的各附加应力σ′_0z^ij之和, 即

(2)

式中:P_ij=σ′₀^ijab, a=L/n, b=B/m, σ′₀^ij为基底平面上第i行、第j列子域内的有效附加应力.
设地基有效压缩层深度共分为l层, 则第k(k= 1, 2, …, l)分层上点M_ijk(x_i, y_j, z_k)所对应的最终沉降值S的表达式为

(3)

式中:

(4)

(5)

(6)

(7)

2 地基土的压缩模量分层总和法和现行规范法统一使用K₀加荷条件下的压缩模量^[9].多年来许多学者对此进行了大量的研究^[10], 普遍认为问题的关键是确定合理的、通用的本构关系.故本文将针对不同土性讨论不同的本构模型, 而不采用统一的本构关系推导压缩模量.
2.1 弹塑性本构模型通析经典弹塑性理论将应变ε分为弹性应变ε^e和塑性应变ε^p, 弹性部分按胡克定律计算, 塑性部分按塑性理论计算.采用增量法表示:

(8)

应力增量和应变增量表示为

(9)

式中, C^ep是弹塑性刚度矩阵, C^ep=C^e-C^p.
对于各向同性的硬化材料, 已知屈服函数f, 塑性势函数g分别为f(σ_ij, H)=0和g(σ_ij, C)=0.其中, H为硬化参量, H=H(ε^p); C为常数.
弹性应变增量dε^e由广义胡克定律确定:

(10)

式中, C^e为弹性刚度矩阵.
塑性应变增量dε^p由塑性位势确定:

(11)

式中, Λ表示塑性应变增量大小的一个标量, 称为塑性标量因子,

(12)

由此可得弹塑性刚度矩阵C^ep为

(13)

(14)

(15)

(16)

将式(10)展开为

(17)

在K₀条件下(土体变形只通过z方向体现), 土体压缩模量用增量表示为竖向应力增量与竖向应变增量之比, 即

(18)

式(18)虽是弹性力学公式, 但通过弹塑性刚度矩阵C^ep结合分层总和法的基本假定进行了修正, 使其能够考虑土体的塑性变形.由式(18)可知, 只要求得相应的本构关系的弹塑性刚度矩阵C^ep, 便可得到修正的压缩模量值.
2.2 基于修正剑桥模型的弹塑性刚度矩阵修正剑桥模型可以反映正常固结黏土的基本变形特征, 其模型参数少, 确定方法简单, 在工程中得到广泛应用.
采用屈服函数f与塑性势函数g相等的相关联流动法则, 则修正剑桥模型的屈服函数为

(19)

式中:λ, κ, M是土性基本参数; e₀是初始平均应力为p₀时的初始孔隙比; ε_v^p为塑性体积变形; p, q分别为平均正应力和广义剪应力.
依据修正剑桥模型的塑性势面g, 屈服面f=g和硬化参量H=ε_v^p, 求得

(20)

(21)

(22)

式中:η=q/p; M为临界状态比, 与土的强度指标内摩擦角φ值有关,

(23)

依据屈服函数f的正负判定土体所处的状态.当f＜0时, 依据式(11)求解竖向变形值; 当f≥0时, 依据弹塑性刚度矩阵C^ep求解能够考虑塑性变形的修正压缩模量E_s, 然后求解竖向变形值.
2.3 基于摩尔-库伦模型的弹塑性刚度矩阵对于砂土, 本文采用摩尔-库伦本构模型, 其模型参数较少且能够采用常规的土工试验获得.
采用屈服函数f与塑性势函数g相等的相关联流动法则, 则摩尔-库伦模型的屈服函数为

(24)

式中:I₁, J₂为应力张量的第一不变量和应力偏张量的第二不变量; c, φ为黏聚力和内摩擦角; θ_σ为Lode角, 可表示为应力偏张量第二不变量J₂和应力偏张量第三不变量J₃的函数:

(25)

屈服函数f对应力的导数为

(26)

式中:

(27)

(28)

(29)

同样, 可计算摩尔-库伦的弹塑性刚度矩阵C^ep, 以便获得能够考虑其塑性变形的压缩模量值.依据屈服函数f的正负状态, 求解其不同状态下地基土的竖向变形值.
本文旨在建立一个完整的、更符合工程实际的大型基础下含软弱夹层的层状砂土地基沉降计算模型, 该模型能够依据不同的土性选用适宜的本构模型.根据弹塑性刚度矩阵C^ep, 推导在K₀条件下能够考虑土体塑性变形的压缩模量.利用平面网格子域法结合分层总和法, 更加准确地计算地基土的沉降.
3 工程算例辽宁省文化艺术中心采用构造板柱下扩展基础.依据岩土工程勘察报告, 地基土在钻探深度范围内自上而下共涉及7个自然分层, 含2个粉质黏土软弱夹层.基础埋深-6.46 m, 地下水的稳定标高为-6.60 m, 地基处于无降水影响的饱和状态.其土层分布与计算参数如表 1所示.
表 1(Table 1)

表 1 地基土的主要力学参数Table 1 Main mechanical parameters of foundation soils 土层名称 土层厚度 分层厚度 变形模量 泊松比 μ 黏聚力标准值 内摩擦角标准值 m m MPa kPa (°) 砾砂 3.54 0.5 28.8 0.3 3.0 36.5 圆砾 3.50 0.5 35.7 0.2 2.0 36.4 砾砂 4.00 1.0 32.7 0.3 2.5 37.0 粉质黏土 0.40 0.4 5.1 0.4 29.3 4.8 圆砾 3.40 1.1 32.9 0.2 2.0 38.0 粉质黏土 0.30 0.3 6.4 0.4 32.1 9.0 粗砂 2.20 1.2 27.0 0.3 1.3 35.0

表 1 地基土的主要力学参数 Table 1 Main mechanical parameters of foundation soils

基底平面网格系统尺寸L=18 m, 宽B=9 m, 将该平面划分为72×36个子域, 对于每一个子域, a=b=250 mm.等效集中力P_ij=0.062 5σ′₀^ij, 作用点为子域的形心.利用工程实测的地基反力结合二元插值得到的有效附加应力空间分布如图 3所示.
图 3(Fig. 3)

图 3 基底附加应力的空间曲面图Fig.3 Spatial distribution of subsidiary stresses under foundation

由图 3可知, 构造板扩展基础下的基底附加应力呈空间曲面分布, 且随深度的增加最大荷载作用点逐渐向基底中心点移动.与小型基础相比更加复杂, 故其变形计算更为复杂, 传统的计算方法并不能适用大型基础下的地基土变形计算.
基于平面子域网格划分和分层总和法的思想, 结合Matlab程序, 利用上述建立的数学模型计算得到基底最终沉降结果, 如表 2所示.
表 2(Table 2)

表 2 各自然层的最大竖向变形及坐标Table 2 Maximum vertical deformation and its coordinates of each natural layer 变量 第1自然层 h=3.54 m 第2自然层 h=3.50 m 第3自然层 h=4 m 第4自然层 h=0.40 m 第5自然层 h=3.40 m 第6自然层 h=0.30 m 第7自然层 h=2.20 m 最大变形值 10.8454 5.0090 3.0025 0.9457 2.4479 0.3482 0.9494 mm 最大变形点 (x, y) (13.00, 1.00) (11.755, 2.00) (10.00, 4.00) (10.00, 4.00) (9.75, 4.00) (9.50, 4.25) (9.50, 4.25) 最大荷载点 (x, y) (12.25, 1.50) (11.00, 3.00) (10.00, 4.00) (10.00, 4.00) (9.75, 4.00) (9.50, 4.25) (9.50, 4.25)

表 2 各自然层的最大竖向变形及坐标 Table 2 Maximum vertical deformation and its coordinates of each natural layer

由计算结果可知, 基底最大附加应力的作用点并不是固定不变的, 随着深度的增加逐渐向地基中心点移动, 这是附加应力随深度的传递和叠加造成的; 最大沉降点与最大荷载作用点并不完全相同, 但随深度的增加最终趋于重合.
地基最终沉降计算结果如图 4所示.最大总变形计算值为22.614 0 mm, 发生点坐标为(12.75 m, 1.25 m), 而基底平面最大附加应力作用点坐标为(9.391 3 m, 0.000 0 m).基础中心点(9.00 m, 4.50 m)处对应的沉降计算值为16.240 0 mm.结果表明, 含软弱夹层的大型层状砂土地基, 地基最大沉降并不一定发生在基础中心处, 且基底最大荷载作用点沉降值也不一定是最大沉降值.因此, 仅计算基础中心点位置的地基沉降是不够的, 中心点与最大值点的误差为28.19%, 存在一定程度的安全隐患.
图 4(Fig. 4)

图 4 地基总竖向变形Fig.4 Total vertical deformations of foundation (a)—空间曲面图；(b)—等值线图.

基于网格划分与分层总和法，对大型基础下计算最大沉降量和寻找最大沉降发生点时, 不存在盲目性, 计算结果直观.
为了进一步证明本文建立的大型基础下含软弱夹层的层状砂土地基沉降计算模型的合理性, 将本文的工程算例分别采用ABAQUS和FLAC^3D进行数值模拟计算, 计算结果如图 5所示.
图 5(Fig. 5)

图 5 各层地基土最大竖向应变曲线图Fig.5 Maximum vertical strain of each layer soil

由图 5可知, 在砂土层中, 本文计算模型的计算结果与ABAQUS和FLAC^3D模拟结果比较接近, 但FLAC^3D模拟黏土层的效果较差.由于ABAQUS中自带修正剑桥模型, 其模拟计算结果与本文模型在黏土层(第4, 6自然层)中的计算结果十分相近, 表明了本文针对不同土性采用不同本构模型的合理性.
4 结论1) 在地基未降水条件下, 基于平面子域网格划分技术, 推导了大型基础下层状地基的沉降计算模型.
2) 基于对基底不同土性的考虑, 分别推导了修正剑桥模型和摩尔-库伦模型的压缩模量的计算方法.
3) 针对大型基础下含有软弱夹层的层状砂土地基, 建立一种完整的、更贴近工程实际的、能够考虑不同土性的地基沉降计算模型.通过工程算例验证了该模型的合理性.
4) 计算实例表明, 地基最大沉降并不一定发生在基础中心处, 且基底最大荷载作用点沉降值也不一定是沉降最大值.按照传统方法计算, 存在一定的安全隐患.
参考文献

[1]	Merlos J, Romo M P. Fluctuant bearing capacity of shallow foundations during earthquakes[J].Soil Dynamics and Earthquake Engineering, 2006, 26: 103–114.DOI:10.1016/j.soildyn.2004.11.029
[2]	Díaz E, Tomás R. A simple method to predict elastic settlements in foundations resting on two soils of differing deformability[J].European Journal of Environmental and Civil Engineering, 2016, 20(3): 263–281.DOI:10.1080/19648189.2015.1035806
[3]	Bezvolev S G. Method of accounting for the deformability of an inhomogeneous elastoplastic bed in analyzing foundation slabs[J].Soil Mechanics and Foundation Engineering, 2002, 39: 162–170.DOI:10.1023/A:1021665429148
[4]	Gazetas G. Static and dynamic displacements of foundations on heterogeneous multilayered soils[J].Géotechnique, 1980, 30: 159–177.DOI:10.1680/geot.1980.30.2.159
[5]	Denis A, Elachachi S M. Effects of longitudinal variability of soil on a continuous spread footing[J].Engineering Geology, 2011, 122: 179–190.DOI:10.1016/j.enggeo.2011.05.015
[6]	李纯, 朱浮生, 张学元. 层状砂土地基静力变形形状测试[J].东北大学学报(自然科学版), 2012, 33(7): 1038–1042. ( Li Chun, Zhu Fu-sheng, Zhang Xue-yuan. Deformation properties of layered sandy foundation soil tested with static additional load[J].Journal of Northeastern University(Natural Science), 2012, 33(7): 1038–1042.)
[7]	Karstunen M, Leoni M. Geotechnics of soft soils:focus on ground improvement[M]. London: Taylor & Francis Group, 2009.
[8]	Peduto D, Nicodemo G, Maccabiani J, et al. Investigating the behavior of buildings with different foundation types on soft soils:two case studies in the Netherlands[J].Procedia Engineering, 2016, 158: 529–534.DOI:10.1016/j.proeng.2016.08.484
[9]	Bezvolev S G, Fedorovskii V G, Aleksandrovich V F. Improvement of foundation settlement calculations by the layerwise summation method[J].Hydrotechnical Construction, 1991, 25(10): 638–644.DOI:10.1007/BF01665652
[10]	Han X L, Xiao C A, Li J X, et al. Spring constitutive model of rigid pile composite foundation and application in design of raft foundation[J].Journal of Central South University, 2013, 20(4): 1079–1083.DOI:10.1007/s11771-013-1587-7