黏弹性地基上厚筏基础蠕变沉降的耦合边界元法.pdf

第 34 卷 第 1 期 岩 土 工 程 学 报 Vol.34 No.1 2012 年 .1 月 Chinese Journal of Geotechnical Engineering Jan. 2012 黏弹性地基上厚筏基础蠕变沉降的耦合边界元法 闫富有 1，刘忠玉1，殷伟希2 1. 郑州大学土木工程学院，河南 郑州 450001；2. 郑州大学综合设计研究院，河南 郑州 450002 摘 要 采用 Boltzmann 剪切松弛模型和体积弹性模型组合的黏弹性半空间地基模型， 考虑筏板的横向剪切变形效应和 柱荷载随时间的变化，在时间增量段对基底反力进行 Newton 二次插值，建立了时间域内黏弹性地基与厚筏基础相互作 用分析的耦合边界元数值方法，并导出了相应的数值计算公式，为分析施工过程以及施工结束后筏板基础随时间持续 发展的沉降和差异沉降以及弯矩和剪力提供了一种有效的数值计算方法。计算表明，在建筑物施工结束后，随着地基 的蠕变，基底反力由基础中间向四周转移，但这种变化并不显著。筏板的沉降、差异沉降和最大弯矩则随地基的蠕变 在一个较长时间内显著增加，最终接近于一稳定值。对于蠕变特性明显的地基，设计时应予以重视。 关键词筏板基础；蠕变沉降；黏弹性半空间；边界元法 中图分类号TU471 文献标识码A 文章编号1000–4548201201–0094–08 作者简介闫富有1963– ，男，河南许昌人，硕士，副教授，主要从事岩土力学、地基基础和地下工程等方面的教 学与研究工作。E-mail yfy。 Coupled boundary element for creep settlement of thick raft foundation on viscoelastic ground YAN Fu-you1, LIU Zhong-yu1, YIN Wei-xi2 1. School of Civil Engineering, Zhengzhou University, Zhengzhou 450001, China; 2. Multi-Functional Design and Research Academy of Zhengzhou University, Zhengzhou 450002, China Abstract A coupled boundary element for analyzing the interactive behaviour between thick raft foundation and viscoelastic half-space soil medium below the foundation in the time domain is presented. The viscoelastic rheological model for soils, combining a standard Boltzmann model in deviatoric component with an elastic model in volumetric component, is employed, and the effect of transverse shear deation in the raft foundation is considered as it bends under varying columnn loading during construction. The boundary element is employed to simulate the raft foundation and the ground creep respectively, and the relevant numerical ulae are derived, in which the Newton quadratic interpolation for subgrade reaction in the time interval is introduced. The present , as well as the computing pragrom, can be used to calculate settlements and differential settlements as well as bending moments and shear forces in the raft. It is shown that the values of subgrade reaction will slightly increase on the edges of the raft and decrease in the centric domain with soil creep after the completion of construction, but the values of variation are not significant and may be ignored reasonably. Howerver, the values of the raft settlement and differential settlement and the maximum bending moments increase pronouncedly with soil creep up to a few years, and tend towards their stable values. And thus the behaviors should be paid attention to in foundation designs for some soil strata with high rheology. Key words raft foundation; creep settlement; viscoelastic half-space; boundary element 0 引 言 一些膨胀土、黄土、饱和砂土以及黏性土在荷载 长期作用下具有明显的流变特性，其本构关系可用黏 弹性模型比较好地拟合[1-2]， 在该类地基上的基础则随 土体的流变表现为长时间的持续沉降。基础的这种沉 降称为蠕变沉降，应采用土的流变模型进行分析[3]。 Viladkar 应用黏弹性有限元法分析了黏弹性地基上独 立基础的次固结沉降行为，并基于 Skempton 孔隙水 压力理论计算了一维竖向固结沉降，表明土体的黏弹 性变形与固结沉降同样不可忽略[4]。Bahar 等基于 Kelvin-Voigt 模型，建立了流变分析的有限元数值格 式，并预测了核电站基础的长期沉降变形[5]。Justo 等 ─────── 收稿日期2010–12–03 第 1 期 闫富有，等. 黏弹性地基上厚筏基础蠕变沉降的耦合边界元法 95 通过现场试验测定了土的流变参数，针对 Boltzmann 流变模型讨论了坝基的蠕变沉降行为[6]。考虑土体流 变特性的地基与筏型基础相互作用分析，主要用于基 础的动力分析[7]，多基于 Winkler 地基模式的黏弹性 地基模型[8]。基于黏弹性半空间地基模式下地基基础 的相互作用分析是一个比基于 Winkler 地基模式的相 互作用分析更为复杂的问题，因为筏板无给定位移等 简单支撑条件，而是随着地基的蠕变产生横向位移， 其解不是简单封闭的形式，必须在每一离散时刻将结 构控制方程与地基的柔度或刚度方程联合求解，从而 产生了数值计算的收敛性问题[9]。虽然有限元法是地 基基础相互作用分析的常用方法，但对于实际工程中 作用有横向荷载和弯矩的大尺寸筏板基础以及半空间 黏弹性地基模型， 其网格划分与后处理工作量非常大， 尤其是在柱作用区域附近对分布力或弯矩的处理并非 令人满意[10]。用边界元法求解厚板的弯曲问题，具有 较高的计算精度且避免了剪力自锁现象[11]；用边界元 法计算黏弹性半空间问题，在处理黏性参数随时间变 化以及具有与时间相关的复杂边界条件时显示出很大 的灵活性[12-14]，且避免了区域离散问题。将二者的积 分方程耦合求解，是进行地基基础相互作用分析的一 种有效方法。 本文针对土体的 Boltzmann 剪切松弛模型和体积 弹性模型组合的黏弹性半空间模型，基于黏弹性边界 元法的时间增量法的研究成果[12-14]，建立半空间问题 增量叠加形式的黏弹性边界积分方程；把 Reissner 板 的边界元法[11,15-16]以及 Winkler 弹性地基上厚板问题 的边界元法[10]应用于黏弹性半空间问题，建立了筏板 基础与黏弹性地基相互作用分析的数值方法，为研究 筏板基础的蠕变沉降提供一种有效方法。 1 筏板的积分方程 考虑由边界所围区域Ω、 厚度为h的筏板基础， 其中面为 3 x0 的x1, x2平面，如图 1 所示。假设其自 重密度q在区域Ω内均匀分布，柱或者剪力墙荷载 Fjt j1, 2 分别表示x1和x2方向的弯矩，j3 表示竖 向力作用的小区域为Ωc。 施工期间Fjt随时间t逐渐 增大，施工结束后保持恒定。筏板底面作用有基底反 力Rt，不考虑筏板与地基接触面上的摩擦力。 根据 Reissner 板理论，仅在自重q作用下板的边 界积分方程为[15] 3 , d , d [ , ijjijj ijji cx uxTx y uyy Ux y tyyUx y         , , ] d i MUx y q yy    ， 1 式中，Tijx, y和U ijx, y分别为面力和广义位移基本 解， ,i U x, y为基本解的导数，其表达式见文献[15， 16]。x和y分别为源点和场点，ti为边界面力，tαMαβnβ， t3Q n ，n为边界Γ的外法线方向余弦，M和Q 分别为板的弯矩和剪力。uj是板的转角（j1，2）或 挠度（j3） ，cijx为与边界几何形状相关的系数。本 文中，拉丁字母下标变化为 1，2，3，希腊字母下标 变化为 1， 2。Mν/[1- 2 ]， 对于 Midlin 板，M0， 和分别为板的特征常数和泊松比， 2 10/h。 图 1 筏板基础模型 Fig. 1 Raft model resting on ground 由于假设板的自重q在Ω内均匀分布，式（1） 中的区域积分可转化为边界积分[16]。一般情况下，筏 板基础的柱荷载均作用在板的内部，故假设筏板为自 由边界条件，即在边界上tj0，方程（1）右端第一项 等于 0。考虑柱荷载后，方程（1）为[10] c , ,3 column , d [ , , ]d [ , , ] d ijjijj ii ikikk cx uxTx y uyy qx yMUx y ny Ux yMUx yF yy                 3, [ , , ] d ii Ux yMUx y R yy      ， 2 式中，δkα为 Kronecker delta 符号，Λi,α为相关的基本 解张量[15-16]，右端第二项表示积分在所有柱的作用域 Ωc内进行。 式（2）的解并不唯一，需建立补充方程。由于含 有基底反力项的积分需要在整个筏板作用区域内进 行，把区域Ω划分为一系列的内部网格，当点x位于 这些网格内时，令i3，式（2）为 c 33 3,3 33 ,3 column , d [ , , ]d [ , , ] d jj kkk u xTx y uyy qx yMUx y ny Ux yMUx yF yy                 333 , [ , , ] d Ux yMUx y R yy      。 3 由于考虑施工期柱荷载的变化，式（2） 、 （3）中 的未知量均为时间t的函数。 2 黏弹性地基的积分方程 由于 Boltzmann 黏弹性模型可以比较好的反映土 96 岩 土 工 程 学 报 2012 年 的蠕变特性且相关参量易于测试[2]，本文采用土的应 力与应变偏量符合 Boltzmann 模型、体积应变与平均 应力为线弹性的组合模型，如图 2 所示。 图 2 土体黏弹性模型 Fig. 2 Viscoelastic model 对于图 2 所示的组合黏弹性模型，应用弹性与黏 弹性对应原理[13]，经推导，在时间t域内相应的黏弹 性半空间问题 Boussinesq 解为  tata bbb r tu 21 ee π4 1 21033   ， 4 式中，r为计算点至基底附加压力Rt作用点的水平距 离，r|y-x|。 01 0 0100101 1 1 2 1 2 00101 1 1 1 00101 2 001 311 3 1 3 3 3 3 GG b GGK GGG G b G G b K GGG G G a K GGG G a KG                          ， ， ， ， 。 5 在基底附加压力Rt作用下，地基表面的沉降量 为 33333 0 1 [ 0d ]d 4π t R u tut Rut r         ，6 式中，R0为t0时刻的基底反力。 把式（4）代入式（6） ，整理得 301122 1 [ , ,]d 4π u tb R tbt abt a r    ，7 式中， 0 , e0ed t ata t R t aR          ， 8 其中，a为a1或a2。 把时间区域离散，令∆titi-ti-1，i1,2,，t00， 函数, at可表示为递推形式 1 11 1 , e , ed i ii i t a ta t ii t R tat a              。9 在时间段[ti-1, ti, ti1]内，对Rt进行Newton二次 插值，式（9）为 1 1211 , e , , i a t iiii tat ag ta R t      11011 , , iiii g ta R tg ta R t   。 10 其初值0, 0 Rat，函数 0 g， 1 g和 2 g分别为 1 01 11 11 11 1 21 111 2 , [ 12] 12 , [ 12] 2 , [ 2] i i iiii i i iii ii i iiii t g ta a ttta t t g ta a tta t tt g ta attta t                                ， ， ， 11 式中， 1 1e i a t     ， 0 g ， 1 g 和 2 g 的初值分别为 0, 10 atg， 111 /,taatg，,, 1112 atgatg。 把式（10）代入式（6） ，得到 1i t时刻基底反力与 竖向沉降量的关系     3101211 22121 11112112 101120121 1 , 4π , ,, ,, ii ii iii iii u tbb g ta r b g taR t b g tab g taR t b g tab g taR t            1121 1122 e ,e ,d ii atat ii bt abt a     。 12 由于 1i t和 i t时刻的基底反力和沉降已经求出，式 （12）便建立了 1i t时刻基底反力与沉降的关系。 3 数值方法 3.1 方程离散求解 把筏板边界Γ离散为NB个结点，mB个线性单元， 假设ui在边界单元内线性变化；把区域Ω划分为Ni 个四边形或三角形网格，假设基底反力Rt在每个网 格内均匀分布。虽然基底反力在内部网格相接处不连 续，但并不影响计算结果，且与现有的地基与基础相 互作用分析方法相一致。计算表明，由于假定了基底 反力Rt在每个网格内均匀分布，消除了边界基底压 力计算值过大的现象。 假设筏板的自重密度q在t0时骤然施加上去； 柱荷载在施工期tc内线性增加，在施工结束后保持常 数。在式（2） 、 （3）中，选取点x为内部网格的中点， 式（2） 、 （3） 、 （12）离散后用矩阵分别表示为 1B11111 2B112121 [ ]{}[]{ }{} []{}{}[]{ }{} iii iiniii TuURQ TuuURQ       ， ， 13 1111 {}[ ] { }{ } iniiii uuRC   。 14 式中，[T1]和[U1]、[T2]和[U2]分别对应于式（2） 、 （3） 中的边界积分项和反力R的区域积分项所组成的系数 矩阵；{uB}和{uin}分别是边界结点的广义位移ui和内 第 1 期 闫富有，等. 黏弹性地基上厚筏基础蠕变沉降的耦合边界元法 97 部网格中点的位移u3所组成的结点位移列阵，{R}是 由内部网格的地基反力排列而成的列阵，{Q1}和{Q2} 分别对应式（2） 、 （3）中均布荷载和柱荷载积分后所 形成的列阵。式（12）中的Rti1项积分后所形成的 系数矩阵为[u′]， 其他项积分后所形成的已知列向量为 {C}。下标i1表示各参量为ti1时刻的值。 式（13） 、 （14）中的未知量{uB}、{uin}和{R}的个 数共计3NB2Ni，方程的个数分别为Ni，3NBNi和 Ni，其解答是封闭的。式（13）的两式组合为   B111 1 222 1 {}[][0][]{} {}[][I][]{} i in i uTUQ R uTUQ       ，15 式中，[I]和[0]分别为单位矩阵和零矩阵。式（15）进 一步表示为 111 [ ]{ }[ ]{ }{ } iii TuURQ   。 16 把式（14）代入式（16）得 11 11in1111 [ ]{ }[ ][ ] {}{ }[ ][ ] { } iiiiii TuU uuQU uC   。 17 对于式（17） ，把{uin}i1的系数按变量排列顺序 叠加到[T]中，同样把右端第二项叠加到{Q}i1中。最 后形成一线性方程组 111 [ ] { }{ } iii TuQ   。 18 求解时，首先形成筏板积分方程的系数矩阵[T] 和[U]。在ti时刻求解结束后，由式（11）计算出 0 g， 1 g和 2 g，并由式（10）递推得到,ati；然后给定一 时间增量 1  i t， 11  iii ttt，根据柱荷载随时间的 变化关系， 得到此时的柱荷载， 由式 （12） 求得 1i tR 的系数矩阵 1 [ ]iu  ，并根据已知的 i tR， 1i tR和 ,ati求得列向量 1 ][ i C；在形成式（13） 、 （14）的系 数矩阵后，进而求得所有边界未知量。 3.2 筏板内力计算 对于Reissner板，其广义应力与广义位移的关系 为[11]     MquuuDM           ,,, 1 2 2 1 ，19     , 3 2 2 1 uuDQ   ， 20 式中，D为板的刚度，]1 12/[ 3 EhD。 令式（2）中的x为板的内部点，并代入式（19） 、 （20） ，得广义应力的计算公式 c 3 column 3 , d [ , , ]d [ , , ] d [ , , ] d iijj ii ikikk ii xTx y uyy qYx yMUx y ny Ux yMVx yF yy Ux yMVx y R yy                         [ ] i M qR y   ， 21 式中， 11x M， 12xy M， 22y M， 13x Q， 23y Q， 张量 , ij Tx y  ， , ij Ux y  ， , i Yx y   和 , i Vx y  分别为相关的广义应力基本解张量[15-16]。 在求得边界结点与内部网格中点的位移和基底反 力后，应用上式便可求得筏板内部任意一点的弯矩和 剪力。 3.3 奇异积分 在式（2） 、 （3）的边界积分项中，当点x和y重 合时， , ij Tx y具有强奇异性，可采用刚体位移法求得 1 [ ]T的主对角子矩阵[11]。当点x和y分别是同一单元 的两个结点时， , ij Tx y在点y的积分具有弱奇异性， 可求出其解析式。区域积分项的奇异性可分为以下二 种情况 （1）点x位于所积分的内部网格之内。取点x 为圆心、半径为的圆，在该圆内进行积分。当趋 于0时，该积分值也趋于0。说明该积分在点x无跳 跃，可把区域积分转化为该内部网格的边界积分。 （2） 点x位于所积分的内部网格的边界， 如筏板 内力计算式（21）中基底反力项的积分。在土体积分 方程12中，1/r的域内积分也归结为这种情况。 此时， 奇异积分可分为如图3（a） 、3（b）和3（c）3种情 况。分别自源点x把单元1234分割为2个、3个和4 个三角形子域，在该区域内的积分等于在所有三角形 子域内的积分之和。由于式（21）中 3 , i Ux y  的积分 具有弱奇异性， , i Vx y  具有强奇异性，下面仅讨论 , i Vx y  的积分。 图 3 区域单元积分奇异性的三种情况 Fig. 3 Three cases of singularity , i Vx y  的表达式为[16] ,, 2 3 1 2 2π 0 Vrr r r Vr            ， 。 22 取图 3 中任意一个三角形子域进行分析。假设三 角形子域的顶点 123 按逆时针顺序排列，顶点 1 为奇 异点。在坐标系, 21 xx中把三角形区域映射为坐标系 , 21 中单位边长的正方形区域[17]，如图 4 所示。 其变换关系为 3 21 2 21 1 1 1 1  xxxx，23 式中， j x的下标表示坐标，上标 j 表示顶点号。 把 Vr  在三角形子域123中的积分转化为沿 98 岩 土 工 程 学 报 2012 年 三角形周界的积分，即 * 12312 23 31 d dVrVr nl         ， 24 ]2[ π4 1 ,,,,, *    rrrrr r rV  ， 25 图 4 三角形单元坐标映射 Fig. 4 Mapping of triangle element 式中，1/1 ，  n为边的外法线方向余弦。 积分沿三角形周界正向12，23和31边进行。 由图3可知，在相邻两个三角形的共用边，积分 方向相反，因而在该边上的积分相互抵消，故式（24） 中仅需要计算在边23上的线积分即可。 由于在边3 2   上1 1 ，则 2 32 ddd kkkk xln ， 26 式中，1,2,1,2k。0 2211  ee，1 12 e，1 21 e。 把式（25） 、 （26）代入式（24） ，得 *32 12312 31 21322132 0101 213201 0101 1 d d 4π [ 2 kk VrVr nle x I xI xI xI x I xI xxIxI                 23 23 ]xIxI   ， 27 式中， 0212121 1322121213221212132 2213232322132323221 3323232 xxxx xxxxxxxxxx xxxxxxxxxx xxxx               ， ， ， ， 28 1 22 2 0 01222 1 22 22 0 01222 d d i i i i I aaa I aaa                   ， 。 29 式（29）中的积分可解析求得。 4 不同算法结果比较 筏板平面尺寸如图 5 所示[18]，板厚 0.45 m，混凝 土强度等级为 C25，柱的几何尺寸和荷载大小见表 1。 筏板的弹性模量和泊松比分别为 E2.8104 MPa， 2 . 0，横向自重荷载 q11.25 kPa 在 t0 时骤然施 加上去。 黏弹性地基模型参数为 G02.2 MPa， G13.5 MPa，K02.86 MPa，η11510 MPad。施工期 tc180 d， 假设在施工期内柱荷载线性变化。 分别采用本文算 法所编制的 C程序和 ABAQUS 有限元软件计算。 边界元法计算所采用的内部网格尺寸为 2 m2 m，共 采用了 91 个内部网格和 40 个边界单元。有限元计算 时，土体计算区域自筏板边界向外取 40 m，计算深度 取 30 m，筏板离散为 1 m1 m 的 S4R 单元，土体离 散为 2 m2 m2 m 的 C3D8R 单元，筏板与土体之间 的接触采用无摩擦接触，柱荷载按集中荷载考虑。 图 5 筏板几何尺寸 Fig. 5 Geometry of raft foundation 表 1 柱几何参数与荷载 Table 1 Geometric parameters and loads of columns 柱编号 尺寸/m 荷载/kN C1 0.50.5 1200 C2 0.50.5 1800 C3 0.50.5 2400 C4 0.50.5 3600 图 6 中心轴线 AB 上筏板的竖向位移结果比较 Fig. 6 Comparison of vertical displacement results along line AB 图 7 轴线 A′B′上筏板的竖向位移结果比较 Fig. 7 Comparison of vertical displacement results along line A′B′ 第 1 期 闫富有，等. 黏弹性地基上厚筏基础蠕变沉降的耦合边界元法 99 在施工结束（ttc）和 1 a 后（ttc1y） ，筏板中 心轴线 AB 和边缘轴线 A′B′（柱下）上竖向位移的两 种算法计算结果比较如图 6，7 所示（图中，L 为筏板 的纵向长度） 。由图 6，7 可知，本文算法和有限元法 的竖向位移计算结果略有差异，即在柱下轴线上有限 元法的结果较边界元法的结果偏大，但其蠕变规律相 同。其差异的主要原因在于有限元法对区域的离散增 加了筏板的柔度，从而导致柱下筏板的竖向位移计算 值偏大[10]。文献[18]分别采用有限元法和有限差分法 的计算结果对比和文献[10]针对 Winkler 弹性地基上 的筏板基础分别采用边界元法和有限元法的计算结果 对比均说明了这一点。 施工结束和一年后筏板轴线 AB 上弯矩的计算结果对比如图 8，9 所示。由图可知，二 者的计算结果比较吻合，表明本文算法是正确的。 图 8 施工结束时轴线 AB 上弯矩结果比较 Fig. 8 Comparison between bending moments by BEM and results by FEM along line AB after completion of construction 图 9 施工结束 1 a 后轴线 AB 上弯矩结果比较 Fig. 9 Comparison between bending moments by BEM and results by FEM along line AB a year after completion of construction 5 算例分析 本节采用 4 中的模型进行分析。在施工期内，沿 筏板纵向中线 AB 上基底反力的分布如图 10 所示。 在 施工结束后，轴线 A，B 两点的基底反力随时间的变 化关系如图 11 所示。 由图可知， 由于假设柱荷载呈线 性增加，基底反力几乎呈线性增大。此后随着地基的 蠕变，筏板中心的基底反力逐渐减少，靠近边缘的基 底反力逐渐增大。说明随着地基的蠕变，基底反力由 基础中间向四周转移，但这种变化并不显著。因此， 可以采用恒定的基底反力计算次固结沉降，这与现有 的次固结沉降的简单计算方法相一致[10]。 图 10 沿轴线 AB 上筏板的基底反力分布 Fig. 10 Results of foundation reaction along line AB 图 11 施工结束后 A，B 点的基底反力随时间的变化 Fig. 11 Results of foundation reaction of points A B vs. time after completion of construction 不同时刻沿筏板中心轴线 AB 上的竖向位移的变 化如图 12 所示。在施工期间，挠度增加较快。施工结 束后，挠度随地基的蠕变长时间持续增加（图 13） ， 板的最大差异沉降持续增加，约经历 4 a 后达到一个 稳定值，稳定时的最大差异沉降较施工结束时的最大 差异沉降增加了 22。 图 12 不同时刻筏板轴线 AB 上的竖向位移 Fig. 12 Results of vertical displacement along line AB vs. time 随着施工荷载的增加，筏板沿长边轴线上的弯矩 100 岩 土 工 程 学 报 2012 年 在不同时刻的分布如图 14 所示。 由于采用了半空间模 型， 距筏板左端 x7 m 和筏板中心 B 点的负弯矩相差 不大。与挠度变化类似，筏板中心的弯矩随地基的蠕 变在一个较长时间内持续增加（图 15） ，达到稳定时 的最大弯矩较施工结束时最大弯矩增加了约 4。 黏弹性地基的计算特点是计算时间较长。计算过 程表明， 由于对基底反力在时间步长内采用了 Newton 二次插值，在时间步长 1  i t的选择上可适当放宽，经 比较发现，时间步长分别取为 1 d 和 10 d 时对计算结 果几乎没有影响。 虽然本例中迭代计算周期为 1800 d， 所用的计算时间很短，其收敛效果很好。 图 13 施工结束后轴线 A，B 两点的竖向位移随时间的变化 Fig. 13 Results of vertical displacement of points A B vs. time ..after completion of construction 图 14 筏板沿中线 AB 上弯矩的分布 Fig. 14 Results of bending moment along line AB 图 15 施工结束后轴线 AB 上最大弯矩随时间的变化 Fig. 15 Results of maximum bending moment of line AB vs. time .after completion of construction 内部网格尺寸的设定对筏板和地基的控制方程耦 合求解能否收敛具有重要影响。以往对筏板与地基相 互作用的数值分析，应用相同模型但不同数值方法的 计算结果相差较大[18]。采用本文方法对多个不同模型 的计算结果表明，当内部网格的边长为 1～3 m 时， 对计算结果的影响很小，说明边界元法具有较高的计 算精度。如果内部网格划分过细，筏板边缘地基反力 的计算值明显较实际值偏大，甚至得到错误的地基反 力，这与筏板边缘土的应力集中和局部破坏有关[9]。 6 结 论 （1）针对应力与应变偏量符合 Boltzmann 模型、 体积应变与平均应力为线弹性的黏弹性地基组合模 型，建立了黏弹性半空间地基上厚筏基础相互作用分 析的边界元直接解法， 并导出了相应的数值计算公式， 为分析厚筏基础施工过程以及施工结束后因地基蠕变 所产生的次固结沉降提供了一种有效的数值方法。 （2） 在时间域内， 黏弹性地基和厚筏基础的边界 元耦合算法具有较高的计算精度和效率。由于对基底 反力在时间步长内采用了 Newton 二次插值，在时间 步长的选择上可适当放宽，在较长的时间周期内进行 地基蠕变分析显得方便和简洁。对