三维地基极限承载力的数值极限分析.pdf

第 29 卷 第 12 期 岩 土 工 程 学 报 Vol.29 No.12 2007 年 12 月 Chinese Journal of Geotechnical Engineering Dec., 2007 三维地基极限承载力的数值极限分析 黄齐武 1，贾苍琴2 1．同济大学地下建筑与工程系，上海 200092；2.中国地质大学（北京）工程技术学院，北京 100083 摘 要根据塑性下限定理，借助有限元技术，建立了极限承载力三维情况下的锥形规划形式。采用原-对偶内点法求 解非线性规划问题，基于相应的理论计算框架，构造满足条件的静态容许应力场。该方法克服了屈服函数的线性化或 光滑近似；利用对偶性，在求解下限解的同时获得对偶命题（上限）的解。与矩形基础地基极限承载力解析结果比较， 精度和准确性有了较大提高；对刚、柔性基础下的地基极限承载力的特性也作了探讨。数值极限分析为计算矩形地基 的极限承载力提供了一种简便有效的方法。 关键词承载力；极限分析；有限元；锥形规划；刚性和柔性基础 中图分类号TU473.1 文献标识码A 文章编号1000–4548200712–1809–06 作者简介黄齐武1979– ，男，湖北仙桃人，博士研究生，从事岩土稳定性问题的数值计算方面的研究。E-mail dresearcher。 3D analysis of ultimate bearing capacity by numerical limit analysis HUANG Qi-wu1，JIA Cang-qin2 1. Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China; 2. School of Engineering and Technology, China University of Geosciences, Beijing 100083, China Abstract Based on the lower bound theorem of classical plasticity, the 3D limit analysis of rigid-perfect plastic structures was pered as a discrete nonlinear mathematical second-order cone programming problem by means of finite element technique. A statically admissible stress field was constructed derived from the solution of the lower bound nonlinear programming problem, which was obtained by primal-dual interior point . The theoretical calculation framework was also presented. Both the linearization and the approximation process for 3D yield function could be removed by adopting nonlinear programming. Based on the duality with the lower bound and the upper bound, the solutions of the primal and the dual problem could be solved at the same time. The characteristics of rigid and flexible foundations were also discussed. It was shown by comparison with analytical solutions that the numerical analysis strategy was vastly superior in terms of the correctness and effectiveness of the procedure. A simple and effective computing for the ultimate bearing capacity of foundations was provided by the 3D limit analysis. Key words ultimate bearing capacity; limit analysis; finite element ; second-order cone programming; rigid and flexible foundation 0 前 言 浅基础的地基极限承载力是基础工程的经典问题 之一。现有的计算方法均是在 Prandtl（1920）和 Reissner（1924）解的基础上，结合基础形状、荷载作 用方向和基础埋深等因素作修正考虑的，形式为 uc cc cccqqq qqq qcN s d i b gqN s d i b g 1 2 BN s d i b gγ γγγ γγγ 。 1 式中 qu为地基的极限承载力；c，q，γ 分别表示黏 聚力、荷载（γD）和土体重度；sc，sq，sγ为基础形 状系数；dc，dq，dγ为基础埋深范围的土体剪阻或修 正系数；ic，iq，iγ为荷载倾斜系数；bc，bq，bγ为基 础倾斜系数；gc，gq，gγ为地面倾斜系数；D 为基础 埋置深度。尽管所有可能因素都有涉及，但这些修正 系数大部分是通过经验方式得出的，在实际运用中往 往出现过于保守的情况[1]。 与条形基础的平面应变机理不同，矩形基础的地 基极限承载力的确定涉及复杂的土体–基础三维相互 作用研究。借助现代数值技术和极限分析定理，采用 数值极限法既可以直接通过上、 下限界定真解的范围， ─────── 收稿日期2006–11–15 1810 岩 土 工 程 学 报 2007 年 同时也避开复杂的土体本构关系；与传统的数值计算 相比，该法更具直接性和目的性。 应用数值极限法的具体思路主要有对土体进行 离散，并构造相应的单元应力模式；根据下限定理， 建立满足平衡方程、屈服条件和应力边界条件的数学 模型；最后引入数学规划的锥形规划（second-order cone programming）求问题的下限解。由于算法基于 原–对偶定理，因而在求解下限问题的同时，其上限 问题也得以解决。本文以矩形基础为例，采用数值极 （下）限法，探讨三维情况下的矩形地基受力特性， 同时也详细阐述极限下限法的数值化的实现。 1 数值下限法的一般原理及其实现 1.1 极限分析的下限定理 下限定理也称为静力相容应力场条件，是以平衡 和屈服条件为基础的，其实质内容为满足下列条件的 静力场 （1）必须是处处满足平衡微分方程的完整的应 力分布或应力场。 （2）必须满足应力边界条件。 （3）必须处处不违背屈服条件。 下限定理的关键在于应力间断场的导出， 或者说， 应力间断场的导出对下限解特别有用。借鉴有限元离 散技术， 一方面可以构造满足相关条件的静力容许场； 另一方面可以考虑复杂的几何和荷载作用形式，避免 繁琐的应力杆或“土柱” （stress legs）计算。 1.2 数值极限方法的基本原理 实际上，传统下限理论应用的局限性在于满足相 应条件的应力场的构造；借助有限单元的划分，构造 单元节点应力为主要未知量的有限元形式，可以获得 满意的结果。与常规位移型有限元不同，极限分析下 限有限元法的主要未知量为节点应力，而且相邻单元 间存在应力间断面。本节详细介绍数值极限法与常规 有限元法的不同之处。 区域离散采用如图1所示的四面体单元，与平面 应变情况相似，单元仍主要采用线性化形函数，原因 有二一方面，计算精度足够满足分析要求；另一方 面，屈服条件验算只需考虑单元节点即可，有助于问 题简化。其中，每个单元节点有六个应力分量 11 σ， 12 σ， 22 σ， 23 σ， 33 σ， 31 σ；节点的正应力和剪应力 分别为 n σ， n τ。 （1）单元间的应力间断条件 前面已经指出，极限分析有限元网格相邻单元间 存在应力不连续面，相邻单元的公共面即为应力不连 续面，如图2所示。由平衡条件要求，相邻单元沿应 力间断面允许其切向力不连续，但必须保证相应的剪 应力和正应力是连续的。图中，两单元公共面上相应 节点的正应力和剪应力分别相等，即， nn ll σσ ′ ， nn ll ττ ′ ，li，k，m 。 2 图 1 三维计算中所用的四面体单元 Fig. 1 Arrangement of elements used in three-dimensional .calculation 图 2 相邻单元的应力间断面 Fig. 2 Statically admissible stress discontinuity among adjacent ..tetrahedrons （2）屈服条件 由于单元采用线性应力分布，故只要单元节点应 力不违反屈服条件，则应力场中的任意一点都不会违 反屈服条件。应力不但要满足边界条件，而且不能违 反屈服条件，即满足 0 e ij fσ≤ e1，2，⋯，N ， 3 式中，f为一个凸的屈服函数，N为计算模型中所有节 点的数目。本文主要以Mohr-Coulomb屈服准则为基 础，其表达形式为 12 /3sin[cos1/3sinsin ] ijij fIJSσϕθϕθ− cos0cϕ≤ ， 4 式中，c，φ分别为土体的黏聚力和内摩擦角，θ为应 力Lode角，Sij为应力偏量，I1为应力张量第一不变量， J2为应力偏张量第二不变量。 （3）目标函数 地基极限承载力的下限求解问题最终就是数值优 化问题，即要寻找一个静力容许应力场，满足平衡方 程、边界条件和不违反屈服条件，并使相应的静力最 第 12 期 黄齐武，等. 三维地基极限承载力的数值极限分析 1811 大。其数学形式为 u Maximize d Subject to ijj S qnS σ σ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪⎩ ∫ ， 应力边界条件 ， 应力平衡条件 ， 应力连续条件 ， 屈服条件 ， 5 式中， ij σ为基础作用面上的应力，nj为法向矢量，Sσ 为外荷作用（或基础）边界。 上述模型是一个具有不等式约束（屈服条件）的 优化问题，必须采用数学规划方法来解决。由于内点 法（interior point ）的迭代次数与问题的规模 无关， 因而目前该法是求解优化问题的主要方法之一。 2 优化算法 对中、小规模极限承载力优化问题，一般采用最 陡边有效集法（the steepest edge active set algorithm） [2] 或Lagrange增项法求解。 但随着计算单元或优化规模 的扩大， 计算规模会以几何级数增加。 值得指出的是， Sloan所用的最陡边有效集法，采用线性化近似屈服 函数的做法， 在三维情况下会极大地增加计算的成本， 因而对于复杂的三维问题有必要寻求非线性规划方 法。 实际上，极限有限元（包括上限和下限）优化问 题关键在于屈服条件的处理。 最早由Sloan[2]采用线性 化屈服函数的方式实现该问题的求解，并得到广泛运 用。随后分别有Lyamin[3]的两阶段拟牛顿算法 （two-stage quasi-Newton algorithm） 、杨洪杰[4]的序列 二次规划法（sequential quadratic programming） 。从几 何角度而言，规划问题的可行域是凸多面体，其顶点 对应基本可行解，而其约束条件为锥形（见图3） ；根 据几何相似性，采用二阶锥形规划（second-order cone programming） 可以避免可行域凸点的平滑近似或几何 转化，提高计算精度和求解速度[5]。 图 3 Mohr-Coulomb 屈服准则空间示意图 Fig. 3 Schematic diagram for 3D Mohr-Coulomb yield surface 由于极限分析上限和下限定理是互为对偶的力学 命题[6-7]，形式如图4，因而许多学者通过数学规划的 对偶理论来考虑相关问题的求解。其中，文献[6 ，8]中 采用原–对偶问题相结合的Lagrangian函数法，对平 面应变承载力问题获得满意结果。 图 4 上限与下限定理的相互关系 Fig. 4 Relationship between the upper and the lower bound .theorem 鉴于以上分析，本文算法主要基于原–对偶内点 法，其求解过程搜索方向采用改进的预示校正法 （predictor-corrector ） 。 其主要原理是引入非负 松弛变量si将屈服函数构造为等式； 再引入Lagrangian 乘子τ，κ和障碍函数构造Lagrangian函数L。利用 原–对偶定理，将上、下限问题结合起来，加快优化 速度，依据对偶间隙和补偿条件作为收敛判定依据， 其主要思想见图5； 预示校正方法仿效Mehrotra算法， 限于篇幅，不在此赘述[9]。其中，图5中所示符号含 义分别有A，b 分别为约束矩阵和下限系数矩阵；x 为单元应力矩阵；K为屈服准则包含区域，一般为外 凸区域；c 为上限系数矩阵；s 为松弛变量矩阵；s， k为极限乘子。 原–对偶内点法（包括预示校正法）的主要步骤 为 （1）确定允许误差ε和严格初始点（x0，s0，τ0， κ0） ∈K（初始k0） 的计算， 使其满足0 0000 κτsx。 （2）优化判定若xkTskε，停止；否则，到步 骤（3） 。 （3）产生预示方向，并确定步长αp、障碍参数 p]1/[ 0000 0 ksxκτ和中心参数σ。 （4）计算校正方向，并确定步长 c α。 （5） 结合前两步所得步长和障碍参数， 确定实际 计算步长 max α，产生新的迭代点；令k k1，转到步 骤（2） 。 上述步骤只是对原–对偶内点法的简单地阐述， 事实上，该算法是相当复杂的过程，比如满足优化 KKT（Karush-Kuhn-Tucker）条件的Hessian矩阵为半 1812 岩 土 工 程 学 报 2007 年 图 5 原–对偶理想弹塑性极限分析有限元理论框架 Fig. 5 Theoretical framework of primal-dual elastic-perfecty plastic finite element 正定阵或奇异阵，以及矩阵稀疏等问题都有涉及。 一般而言，尽管本文采用的是Mohr-Coulomb屈 服准则，对于几何形如锥形的屈服准则均适用。由上 所述可知，采用原–对偶内点法在求解原问题（下限 解）的同时，也能获得其对偶问题（上限解）的解。 这也是该方法的优点之所在。 3 三维地基的极限承载力 众所周知，在二维情况下，地基的极限承载力受 基础刚柔度的影响。对柔性基础情况，采用均布荷载 直接作用在土体单元上；而在刚性基础情况下，采用 刚度远大于土体的单元来模拟基础，荷载作用在这些 单元上， 其计算模型如图6所示。 由于对称性仅取1/4 部分作为研究对象，对称轴所在面上切向应力为零， 基础作用范围受法向均布力作用，其余各面法向正应 力和切向应力均为零。 图 6 矩形基础的典型计算网格 Fig. 6 Typical lower bound mesh for rectangular foundations 以L/B 1和L/B 2两种情况来详细探究三维地 基极限承载力的相关特性。 图7和图8分别是L/B 1 情况下地基变形和应变速度矢量图。由两组图可以知 道，在相同条件下，刚性基础底端出现刚性区域，即 基础底端区域网格未变形或应变速度为零；柔性基础 下土体不存在未变形区域，而且应变速度矢量较刚性 图 7 柔性和刚性矩形基础的破坏模式 Fig. 7 Failure mechanism of flexible and rigid rectangular foundations 1813 岩 土 工 程 学 报 2007 年 基础大。另外，刚性基础下土体的影响范围较柔性基 础的广，可以达到模型长度的3/4，柔性基础则只能 达到长度的2/3处。在某种程度上而言，地表土体的 剪阻 （shear resistance） 有利于地基极限承载力的增加； 或者说，在相同条件下，刚性基础下的地基极限承载 力大于柔性基础，主要是因为有更多的地表土体剪阻 共同抵抗外力的作用。 这与Leshchinsky[10]基于试验的 观点是一致的。 图 8 柔性和刚性基础下地基的计算应变速度场 Fig. 8 Calculated strain velocity field for flexible and rigid .foundations 对于方形基础的地基极限承载力与摩擦角之间的 函数关系如图9和图10所示。 在相同条件下， 地基的 极限承载力随摩擦角的增大而增大。在图中，比较方 形基础各种计算结果， 包括Garnier[11]的数值解、Narita 和Yamaguchi[12]的极限平衡法解以及Salenon和 Matar[13]的圆形基础解。 根据极限分析原理， 真解包含 在上限解与下限解区间内；在理论上而言，由于极限 平衡法仅在破坏处满足平衡相关条件，既不是上限解 也不是下限解，更不是真解，因而其准确性有待进一 步验证。由图可知，Narita和Yamaguchi的条分法结 果落在上下限范围内，说明具有一定的准确性；尽管 Garnier采用10节点四面体高阶单元数值计算，但比 较而言，本文所得的结果更为合理，或界定区域范围 更准确。根据Vesic的承载力理论，等面积圆形基础 和方形基础的地基极限承载力相等，因此，Salenon 和Mater依据内切圆（inscribed circle）和外接圆 （circumscribed circle）基础来界定方形基础的极限承 载力， 通常认为该法结果即为方形基础的下限和上限。 比较发现，本文方法与该法结果两者吻合较好，规律 基本一致，本文数值计算结果界定范围较该法范围最 大缩小近21（φ20。图11是摩擦角φ10和30 的破坏变形图。随着摩擦角的增大，地基的变形范围 从深度和广度上都明显增大。 图 9 方形柔性基础下地基极限承载力结果比较 Fig. 9 Comparison of the ultimate bearing capacity of a flexible .square foundation 图 10 方形刚性基础地基极限承载力结果比较 Fig. 10 Comparison of the ultimate bearing capacity of a rigid .square foundation 图 11 刚性基础下地基的破坏模式 Fig. 11 Failure mechanism of rigid foundations 图12是L/B 2矩形基础的计算结果比较。整体上 讲，Salenon和Mater的解法涵括了所有的结果，但界 1814 岩 土 工 程 学 报 2007 年 定的范围随摩擦角的增大而增大。本文数值解法界定 的范围较之有极大的缩小，利于真解的确定。Garnier 的数值解在φ≥20时，高于本文的上限解，即高估了 地基的极限承载力。而Narita和Yamaguchi的解在φ 30左右低于本文的下限解，保守地估计了地基的极 限承载力。 图 12 矩形基础（L2B）的地基极限承载力比较 Fig. 12 Comparison of the ultimate bearing capacity of a .rectangular foundation L2B 4 结 语 受尺寸效应影响， 三维地基极限承载力仍是岩土工 程有待研究的难点问题之一。 采用塑性力学中的上、 下 限定理，分别求得极限荷载的上限解、下限解，进而估 计极限荷载的变化范围； 从而跳出对复杂的本构关系和 结构破坏时呈现的体胀、 软化等特性探究的窠臼， 能大 大简化问题的求解，是一种新的理论计算手段。 借助有限元离散技术，基于锥形数学规划方法， 采用原–对偶内点法，解决了矩形基础三维极限承载 力问题，得到较好的结果。与已有方法相比，本文提 出的数值极限法能大大缩小上、下限的界定范围，为 进一步确定真解打下良好基础。在相同条件下，刚、 柔性基础的地基极限承载力有所不同，主要是源于参 与抵抗外荷载作用的土体剪阻不同造成的。刚性基础 的作用范围较之柔性基础在广度和深度上都要大。同 样地，在刚性基础条件下，随内摩擦角φ的增大基础 作用范围也逐渐增大。 致谢 本文所有数值计算结果后处理均基于GPLGeneral Public License的 Gnuplot 4.0 软件完成。在此，对软件的开发人员 Thomas Williams 和 Colin Kelley 两位先生的慷慨致以敬意 参考文献 [1] 黄齐武, 黄茂松, 王贵和. 基于下限有限元法的条形浅基 础极限承载力分析[J]. 岩土工程学报, 2007, 294 572– 579. HUANG Qi-wu, HUANG Mao-Song, WANG Gui-he. Calculation of bearing capacity of strip footings using lower bound limit [J]. Journal of Chinese Geotechnical Engineering, 2007, 294 572–579. in Chinese [2] SLOAN S W. Lower bound limit analysis using finite elements and linear programming[J]. International Journal for Numerical and Analytical s in Engineering, 1988, 12 61–77. [3] LYAMIN A V, SLOAN S W. Lower bound limit analysis using nonlinear programming[J]. International Journal for Numerical s in Engineering, 2002, 555 573–611. [4] YANG H J, SHEN Z J, WANG J H. 3D lower bound bearing capacity of smooth rectangular surface footings[J]. Mechanics Research Communications, 2003, 30 481–492. [5] MAKRODIMOPOULOS A, MARTIN C M. Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[J]. International Journal for Numerical s in Engineering, 2006, 664 604–634. [6] LEU S Y. A finite-element limit analysis for materials with pressure dependent yield behavior[D]. Ann Arbor University of Michigan, 1998. [7] HEITZER M. Basis reduction technique for limit and shakedown problems[R]// Numerical s for Limit and Shakedown Analysis-Deterministic and Probabilistic Problems Part I, NIC 15, 2003. [8] KOBAYASHI S. Development of hybrid rigid plastic finite element based on primal-dual interior point [J]. Journal of Applied Mechanics, 2003, 6 96–106 JSCE, in Japanese. [9] MEHROTRA S. On the implementation of a primal-dual interior point [J]. SIAM Journal on Optimization, l992, 2 575–601. [10] LESHCHINSKY D, MARCOZZI G F. Bearing capacity of shallow foundations rigid versus flexible models[J]. Journal of Geotechnical Engineering, 1990, 11611 1750–1756. [11] GARNIER D, MAGHOUS S. Feasibility of a numerical for computing the ultimate loads of three dimensional structures[J]. Mechanics Research Communications, 1995, 223 279–288. [12] NARITA K, YAMAGUCHI H. Three-dimensional bearing capacity analysis of foundations by use of a of slices[J]. Soils and Foundations, 1992, 324 143–155. [13] SALENON J, MATAR M. Bearing capacity of circular shallow foundations[M]// PILOT G. Foundation Engineering, Paris Presses de l’ENPC France, 1982 159–168.