基于线性规划模型的极限分析上限有限元的实现

极限分析上限有限元法常利用三角形常应变率单元和摩尔-库仑屈服函数线性化的方法,以形成较易求解的线性规划模型。然而为满足计算精度,需引入大量的优化变量,增加了计算和存储的难度。为此,基于线性规划模型,利用MATLAB编制极限分析上限有限元程序,针对线性规划模型等式约束矩阵的高度稀疏性的特点,以稀疏矩阵的方式存储,从而在一定程度上解决了上述问题,使得上限有限元法能处理较大规模的岩土工程稳定性问题。以条形基础地基承载力课题为例进行算例分析,验证该方法的有效性,同时讨论了模型网格单元和塑性乘子数目对计算结果精度的影响。     

边坡稳定分析刚体有限元上限法的锥规划模型

极限分析方法是土边坡稳定性分析的重要方法之一。刚体有限元上限法是其中的一类,此类方法仍旧存在一些关键问题需要完善。由于单元的刚性假设,系统的塑性变形内能耗散仅发生在单元间的界面上,故此类方法的性能主要取决于界面的布局,即采用非结构化三角形单元计算往往精度较差。为此,提出了基于滑动面摄动的刚体有限元上限法及临界滑动面的搜索方法。首先,在考虑刚体转动的基础上构造刚体有限元上限法的二阶锥规划模型,用于确定在给定试滑动面条件下的运动许可速度场。其次,将试滑动面的控制参数视为决策变量,建立搜索临界滑动面的非线性非凸优化问题模型,并采用非线性单纯形方法和粒子群方法求解此优化问题找出临界滑动面。通过经典边坡稳定问题的分析求解,验证了所提出的新方法,进一步证实了网格类型（即界面的布局）是影响刚体有限元上限法计算精度的主要因素。经过计算结果的对比发现,在刚体有限元上限法中考虑刚体转动是非常必要的,不仅可以提高刚体有限元上限法的计算精度,还可以克服此方法对界面布局的依赖性。     

基于二阶锥规划有限元增量加载法的条形浅基础极限承载力分析

地基极限承载力分析是土力学研究中的一个经典课题。基于Hellinger-Reissner混合变分原理和有限元方法,将岩土体弹塑性问题构造成基于有限元框架的二阶锥规划（second-order cone programming,SOCP）问题,进而提出一种基于二阶锥规划理论的增量有限元法,即FEM-SOCP法。将岩土体弹塑性问题构造成二阶锥规划的数学优化问题,可以避免采用传统弹塑性计算中复杂的应力点积分等算法和屈服面棱角的平滑处理。此外,对于二阶锥规划问题,可以采用具有原始?对偶内点求解法的标准数学规划求解器MOSEK进行求解。将增量加载FEM-SOCP法应用于经典的基底粗糙的条形浅基础地基极限承载力分析中,分别考虑了关联和非关联塑性条件下的Mohr-Coulomb屈服准则。数值结果表明：所提出的增量加载FEM-SOCP法获得的地基承载力系数及地基承载力与传统FEM计算结果基本一致,而与常规有限元计算结果相比,基于增量加载的FEM-SOCP法所获得的屈服区更加平滑。     

断裂问题的扩展有限元法研究

扩展有限元（extended finite element method,XFEM）是近年来发展起来的、在常规有限元框架内求解不连续问题的有效数值计算方法,其基于单位分解的思想,在常规有限元位移模式中加入能够反映裂纹面不连续性的跳跃函数及裂尖渐进位移场函数,避免了采用常规有限元计算断裂问题时需要对裂纹尖端重新加密网格造成的不便。在推导扩展有限元算法的基础上,分析了应力强度因子的J积分计算方法及积分区域的选取。采用XFEM对I型裂纹进行了计算,有限元网格独立于裂纹面,无需在裂纹尖端加密网格;分析了积分区域、网格密度对应力强度因子计算精度的影响,指出了计算应力强度因子的合适参数,验证了此方法的可靠性和准确性。     

浅埋隧道掌子面稳定性二维自适应上限有限元分析

为研究浅埋隧道掌子面稳定性及获取精细化的破坏模式,提出了一种上限有限元非结构化网格自适应加密策略。以单元耗散能权重指标作为网格自适应加密评判准则,该策略同时兼顾了单元尺度与塑性应变。应用高阶的6节点三角形单元并建立上限有限元线性规划模型,以多次反复计算和网格加密的方式实现了二维自适应上限有限元分析并编制了计算程序。利用条形基础地基极限承载力课题,从上限解精度和网格加密形态方面验证了该程序的有效性。针对浅埋隧道掌子面稳定性问题,展开多参数条件下的自适应上限有限元计算,分析了网格加密过程中单元总数与上限解精度的关系,列出不同隧道埋深和内摩擦角对应的隧道掌子面稳定性临界值的上限解,揭示出掌子面稳定性变化规律及精细化的破坏模式。     

基于塑性极限分析上限法理论的土质边坡可靠度分析

将极限分析的上限定理、有限元离散思想、随机规划理论和蒙特卡洛方法这四者结合起来,提出了一种土质边坡可靠度分析的上限数值方法。首先采用三节点有限单元离散土质边坡,然后将土体的抗剪参数设为随机变量,根据上限定理构建同时满足三角形单元的塑性流动约束条件、单元公共边的塑性流动约束条件和单元速度边界条件的机动许可速度场,并根据内功功率等于外功功率条件建立目标函数,构建土质边坡可靠度分析的上限法随机规划模型。采用蒙特卡洛方法求解上限法随机规划模型,同时提出了一种基于上限法速度场的边坡失效风险系数估算方法,该方法特别适用于具有多种失效模式的边坡风险分析。对2个经典算例进行了深入分析,验证了方法的正确性。     

基于多重网格法的边坡稳定分析

边坡稳定分析的刚体极限平衡法不能很好反映岩体中实际的应力分布,而基于有限元的边坡稳定分析方法很难得到一个明确的滑动安全系数.为了解决这些问题,采用多重网格法,分别建立用于有限元计算的结构网格和用于计算滑面稳定安全系数的滑面网格,可以方便地获得任意滑面或滑块的稳定安全系数,从而将非线性有限元方法和极限平衡分析结合起来,充分考虑了变形和岩体弹塑性应力调整对边坡稳定的影响.将该算法集成到三维非线性有限元程序TFINE中,在分析网格密度对结果精度的影响的基础上,应用于锦屏高边坡的稳定分析中.结果表明,该方法是一种有效的边坡稳定分析方法,更符合实际情况.     

水平受荷桩的弹性有限元虚拟加载上限分析

基于将塑性上限分析等效为弹性迭代计算的总量虚拟加载上限分析理论,在商业化有限元软件ABAQUS中实现了弹性有限元虚拟加载上限方法（弹性有限元T-EMSD）。应用弹性有限元T-EMSD法分析了不排水黏土中的二维水平受荷桩,其获得的荷载-位移曲线与弹塑性有限元分析结果一致,其极限承载力与塑性解相近。在极限位移加载量下弹性有限元T-EMSD法对应的上限机构从弹性始速度场开始随迭代逐渐演化,迭代收敛后的速度场和解析塑性破坏机构相似。与其他基于可变强度概念（MSD）的方法相比,弹性有限元T-EMSD法对水平受荷桩桩身的分析具有更高的精度。弹性有限元T-EMSD法最大的优势在于可在计算中自然地获得塑性机构,因而可被用于研究一些塑性机构难以构造的复杂问题,并对弹塑性数值方法进行验证。     

大应变静力触探数值模拟及锥形因子影响因素分析

以模拟圆锥在不排水黏土中的静力贯入和分析其影响因素为目标,假设土体为均质弹性—完全塑性材料且服从Mises屈服准则,采用任意拉格朗日—欧拉(ALE)网格划分技术确保锥尖土体在大应变条件下的网格质量,进行大应变有限元数值模拟,并分析了稳定状态下土体刚度指数、原位应力状态和锥尖粗糙程度对塑性区半径与锥形因子的影响,获得了锥形因子表达式。模拟结果表明:塑性区随着刚度指数的增大而增大,锥尖周围塑性区的径向扩张处于柱形孔扩张和球形孔扩张之间,更接近于球形孔扩张;锥形因子随土体刚度指数、锥尖粗糙程度的增大而增大,随土体原位应力状态参数的增大而减小;得到的锥形因子表达式可以量化土体刚度指数、原位应力状态和锥尖粗糙程度的影响,具有较高的精确度。     

带岛屿约束的三角形空间剖分算法的改进方法

本文针对数值模拟中三角形有限元网格的具体要求，对经典的Bowyer-watson算法作了局部改进，利用具有岛屿约束的Delaunay剖分方法来构建三角形网格，使其适用于有岛屿、关键点和特征线约束的区域剖分的要求。此方法在（苏州，无锡，常州）以下简称苏锡常地区地下水流数值模拟中得到了应用，结果表明此方法自动生成带岛屿约束的三角形网格是可行的。     

Upper bound limit analysis using finite elements and linear programming

This paper describes a technique for computing rigorous upper bounds on limit loads under conditions of plane strain. The method assumes a perfectly plastic soil model, which is either purely cohesive or cohesive-frictional, and employs finite elements in conjunction with the upper bound theorem of classical plasticity theory. The computational procedure uses three-noded triangular elements with the unknown velocities as the nodal variables. An additional set of unknowns, the plastic multiplier rates, is associated with each element. Kinematically admissible velocity discontinuities are permitted along specified planes within the grid. The finite element formulation of the upper bound theorem leads to a classical linear programming problem where the objective function, which is to be minimized, corresponds to the dissipated power and is expressed in terms of the velocities and plastic multiplier rates. The unknowns are subject to a set of linear constraints arising from the imposition of the flow rulé and velocity boundary conditions. It is shown that the upper bound optimization problem may be solved efficiently by applying an active set algorithm to the dual linear programming problem. Since the computed velocity field satisfies all the conditions of the upper bound theorem, the corresponding limit load is a strict upper bound on the true limit load. Other advantages include the ability to deal with complicated loading, complex geometry and a variety of boundary conditions. Several examples are given to illustrate the effectiveness of the procedure.     

An element for modeling potential channel flows between impervious underground structures in an urban regional groundwater flow field

In an urban regional groundwater flow field, the presence of thousands of channels between impervious structures makes the field difficult to simulate using the finite-element method (FEM), because the scale of the field is usually several orders of magnitude larger than that of the channels. To overcome this problem, a simple element for the simulation of potential channel flows has been developed. This element works with linear triangular elements and can be easily implemented in a finite-element code to simulate the channel flows with a sparse mesh without a loss of global accuracy. The transmissivity matrix of this element is deduced from the analytical solution for channel flow. The application of the element is discussed, and the accuracy of the element is assessed. The element makes it easy to merge small structures that are close to each other into a larger one that can be modeled with a sparse mesh.     

复杂区域四边形网格生成的一种改进方法

提出了一种改进的四边形网格生成算法。该法在行波法生成平面三角形网格的基础上，将一个三角形单元分解为三个四边形单元。经过优化处理，获得优良的计算域四边形网格。针对岩土工程结构问题，还提出了网格剖分中尖灭不连续面和锚杆及锚索的处理方案。地下洞室和边坡算例表明了这种方法的可行性和可靠性。     

边坡稳定分析的自适应性局部网格加密解法

常规的有限元强度折减法在确定边坡塑性区时,不能直观地反映出塑性区的形状和滑裂面位置,计算结果精度低。在塑性区范围内进行网格局部加密,不仅使安全系数的计算精度大幅度提高,而且可求出滑裂面的形状和位置。针对求解塑性区搜索滑裂面这一具体问题,提出塑性应变界限值网格加密准则,采用自适应性网格局部加密技术求解边坡的滑裂面和稳定安全系数。数值计算结果表明,所提出的加密准则能正确地确定滑裂面的形状和位置,计算所得稳定安全系数和理论稳定安全系数仅相差0.2 %。     

Kinematic modelling of shear band localization using discrete finite elements

Modelling shear band is an important problem in analysing failure of earth structures in soil mechanics. Shear banding is the result of localization of deformation in soil masses. Most finite element schemes are unable to model discrete shear band formation and propagation due to the difficulties in modelling strain and displacement discontinuities. In this paper, a framework to generate shear band elements automatically and continuously is developed. The propagating shear band is modelled using discrete shear band elements by splitting the original finite element mesh. The location or orientation of the shear band is not predetermined in the original finite element mesh. Based on the elasto‐perfect plasticity with an associated flow rule, empirical bifurcation and location criteria are proposed which make band propagation as realistic as possible. Using the Mohr–Coulomb material model, various results from numerical simulations of biaxial tests and passive earth pressure problems have shown that the proposed framework is able to display actual patterns of shear banding in geomaterials. In the numerical examples, the occurrence of multiple shear bands in biaxial test and in the passive earth pressure problem is confirmed by field and laboratory observations. The effects of mesh density and mesh alignment on the shear band patterns and limit loads are also investigated. Copyright © 2003 John Wiley & Sons, Ltd.     

Limit analysis of 2-D and 3-D structures based on an ellipsoid yield criterion

In this paper, a nonlinear numerical technique is developed to calculate the limit load and failure mode of structures obeying an ellipsoid yield criterion by means of the kinematic limit theorem, nonlinear programming theory and displacement-based finite element method. Using an associated flow rule, a general yield criterion expressed by an ellipsoid equation can be directly introduced into the kinematic theorem of limit analysis. The yield surface is not linearized and instead a nonlinear purely kinematic formulation is obtained. The nonlinear formulation has a smaller number of constraints and requires less computational effort than a linear formulation. By applying the finite element method, the kinematic limit analysis with an ellipsoid yield criterion is formulated as a nonlinear mathematical programming problem subject to only a small number of equality constraints. The objective function corresponds to the dissipation power which is to be minimized and an upper bound to the plastic limit load of a structure can then be calculated by solving the minimum optimization problem. An effective, direct iterative algorithm has been developed to solve the resulting nonlinear programming formulation. The calculation is based purely on kinematically admissible velocities. The stress field does not need to be calculated and the failure mode of structures can be obtained. The proposed method can be used to calculate the bearing capacity of clay soils in a direct way. Some examples are given to illustrate the validity and effectiveness of the proposed method.     

Upper bound limit analysis of cohesive soils using mesh-free method

A mesh-free based limit analysis approach is proposed, to determine the upper bound solutions for the collapse loads associated with cohesive soils, under plane strain conditions. In the presented technique, the geometry of problem is just simulated by nodes and there is no need of mesh in the traditional sense. The process of finding an upper bound solution consists of combining limit analysis theory and a mesh-free numerical technique as a discretisation tool. To satisfy the required conditions for the admissibility of the discretised velocity field at the entire problem domain, a strain rate smoothing technique has been adopted. The outcome of proposed combination is a nonlinear optimisation problem which is solved by a direct iterative technique. The solution found by an iterative algorithm is an upper bound for limit load of the stability problem. The efficiency of the proposed method is demonstrated by solving different example problems in the soil mechanics engineering field, at the end of the paper.