Development of an implicit material point method for geotechnical applications

An implicit material point method (MPM), a variant of the finite element method (FEM), is presented in this paper. The key feature of MPM is that the spatial discretisation uses a set of material points, which are allowed to move freely through the background mesh. All history-dependent variables are tracked on the material points and these material points are used as integration points similar to the Gaussian points. A mapping and re-mapping algorithm is employed, to allow the state variables and other information to be mapped back and forth between the material points and background mesh nodes during an analysis. In contrast to an explicit time integration scheme utilised by most researchers, an implicit time integration scheme has been utilised here. The advantages of such an approach are twofold: firstly, it addresses the limitation of the time step size inherent in explicit integration schemes, thereby potentially saving significant computational costs for certain types of problems; secondly, it enables an improved algorithm accuracy, which is important for some constitutive behaviours, such as elasto-plasticity. The main purpose of this paper is to provide a unified MPM framework, in which both quasi-static and dynamic analyses can be solved, and to demonstrate the model behaviour. The implementation closely follows standard FEM approaches, where possible, to allow easy conversion of other FEM codes. Newton’s method is used to solve the equation of motion for both cases, while the formation of the mass matrix and the required updating of the kinematic variables are unique to the dynamic analysis. Comparisons with an Updated Lagrangian FEM and an explicit MPM code are made with respect to the algorithmic accuracy and time step size in a couple of representative examples, which helps to illustrate the relative performance and advantages of the implicit MPM. A geotechnical application is then considered, illustrating the capabilities of the proposed method when applied in the geotechnical field.     

A GPU parallel computing strategy for the material point method

The material point method (MPM), which is a combination of the finite element and meshfree methods, suffers from significant computational workload due to the fine mesh that is required in spite of its advantages in simulating large deformations. This paper presents a parallel computing strategy for the MPM on the graphics processing unit (GPU) to boost the method’s computational efficiency. The interaction between a structural element and soil is investigated to validate the applicability of the parallelisation strategy. Two techniques are developed to parallelise the interpolation from soil particles to nodes to avoid a data race; the technique that is based on workload parallelisation across threads over the nodes has a higher computational efficiency. Benchmark problems of surface footing penetration and a submarine landslide are analysed to quantify the speedup of GPU parallel computing over sequential simulations on the central processing unit. The maximum speedup with the GPU used is ∼30 for single-precision calculations and decreases to ∼20 for double-precision calculations.     

A finite element approach to solve contact problems in geotechnical engineering

This paper presents a finite element approach to solve geotechnical problems with interfaces. The behaviours of interfaces obey the Mohr–Coulomb law. The FEM formulae are constructed by means of the principle of virtual displacement with contact boundary. To meet displacement compatibility conditions on contact boundary, independent degrees of freedom are taken as unknowns in FEM equations, instead of conventional nodal displacements. Examples on pressure distribution beneath a rigid strip footing, lateral earth pressure on retaining walls, behaviours of axially loaded bored piles, a shield‐driven metro tunnel, and interaction of a sliding slope with the tunnels going through it are solved with this method. The results show good agreement with analytical solutions or with in situ test results. Copyright © 2005 John Wiley & Sons, Ltd.     

非完全粗糙平板锚旋转安装过程中的丢失埋深研究

吸力贯入平板锚被用于系泊深海浮式结构,实际应用中必须预估平板锚旋转安装过程中的丢失埋深。采用大变形有限元方法探索非完全粗糙平板锚在正常固结黏土中的旋转过程。大变形有限元法通过网格重剖避免平板锚大幅值平动和转动引起的土体单元扭曲。根据平板锚旋转到达的埋深,实时更新“锚-土”界面上的摩擦剪切强度。将数值模拟结果与离心机模型试验进行对比验证,表明高岭土中平板锚表面的粗糙系数约为0.3。平板锚丢失埋深随粗糙系数的减小而增大,但粗糙系数对丢失埋深的影响受平板锚厚度比和拉力偏心比的耦合作用。厚度比和偏心比越小,粗糙系数对丢失埋深的影响越显著。     

Numerical investigation of pile installation effects in sand using material point method

The installation of displacement piles in sand leads to severe changes in the stress state, density and soil properties around the pile tip and shaft, and therefore has a significant influence on the pile bearing capacity. Most current numerical methods predicting pile capacity do not take installation effects into account, as large deformations can lead to mesh distortion and non-converging solutions. In this study, the material point method (MPM) is applied to simulate the pile installation process and subsequent static pile loading tests. MPM is an extension of the finite element method (FEM), which is capable of modelling large deformations and soil-structure interactions. This study utilizes the moving mesh algorithm where a redefined computational mesh is applied in the convective phase. This allows a fine mesh to be maintained around the pile tip during the installation process and improves the accuracy of the numerical scheme, especially for contact formulation. For the analyses a hypoplastic constitutive model for sand is used, which takes into account density and stress dependent behaviour. The model performs well in situations with significant stress level changes because it accounts for very high stresses at the pile tip. Numerical results agree with centrifuge experiments at a gravitational level of 40 g. This analysis confirms the importance of pile installation effects in numerical simulations, as well as the proposed numerical approach’s ability to simulate installation and static load tests of jacked displacement piles.     

An edge-based strain smoothing particle finite element method for large deformation problems in geotechnical engineering

To solve large deformation geotechnical problems, a novel strain-smoothed particle finite element method (SPFEM) is proposed that incorporates a simple and effective edge-based strain smoothing method within the framework of original PFEM. Compared with the original PFEM, the proposed novel SPFEM can solve the volumetric locking problem like previously developed node-based smoothed PFEM when lower-order triangular element is used. Compared with the node-based smoothed PFEM known as “overly soft” or underestimation property, the proposed SPFEM offers super-convergent and very accurate solutions due to the implementation of edge-based strain smoothing method. To guarantee the computational stability, the proposed SPFEM uses an explicit time integration scheme and adopts an adaptive updating time step. Performance of the proposed SPFEM for geotechnical problems is first examined by four benchmark numerical examples: (a) bar vibrations, (b) large settlement of strip footing, (c) collapse of aluminium bars column, and (d) failure of a homogeneous soil slope. Finally, the progressive failure of slope of sensitive clay is simulated using the proposed SPFEM to show its outstanding performance in solving large deformation geotechnical problems. All results demonstrate that the novel SPFEM is a powerful and easily extensible numerical method for analysing large deformation problems in geotechnical engineering.     

吸力式筒形基础沉贯过程的大变形有限元模拟

吸力式筒形基础在海洋工程中已获得越来越广泛地应用,其安装过程的数值模拟对指导工程实践具有重要意义。在大型通用有限元软件ABAQUS平台上建立了二维轴对称模型,基于ALE（任意拉格朗日-欧拉法）技术模拟了黏土中吸力筒的大变形沉贯过程。模拟过程利用了子程序VUFIELD控制土体的不排水抗剪强度和弹性模量随土体深度变化。参照离心机试验及理论计算,对模型进行验证。利用已验证模型分析不同吸力下沉贯阻力、土塞高度,并讨论了筒壁摩擦特性。数值计算结果表明,ALE技术能有效地模拟吸力筒贯入过程,避免网格畸变。贯入方式对贯入阻力影响很大,吸力式贯入阻力明显低于压力式贯入阻力。进一步研究发现,随着最终吸力值的增大,沉贯阻力会显著降低,土塞高度会显著提高。对内壁摩擦特性的研究表明,内壁摩擦阻力是导致沉贯阻力改变的主要因素,并且相比吸力式贯入方式,筒壁摩擦特性会对压力式贯入造成更大的影响。     

边坡稳定有限元分析中的三个问题

对大型有限元软件ABAQUS进行二次开发,提出了自动搜索安全系数的边坡稳定有限元分析模型。在此基础上进行大量变动参数研究,探讨了迭代不收敛、塑性区贯通及等效塑性应变贯通等3种边坡失稳判据的内在联系与适用性,其中迭代不收敛判据易于自动搜索的编程实现,且较少依赖研究者的经验。大多数情况下单元阶次不影响安全系数的确定,但一阶单元有时可能高估安全系数,建议采用二阶单元。基于更新拉格朗日格式实施边坡稳定大变形有限元分析,结果表明迭代不收敛准则不适用于大变形分析。     

模拟摩擦接触问题的新型无网格数值方法

针对摩擦接触问题,给出一种新型的无网格数值方法。该方法基于单位分解思想,在标准无网格Galerkin法的位移模式中嵌入不连续函数和裂尖奇异函数,分别反映接触面的不连续性和接触面端点的奇异性;结合接触摩擦定律,继而构造出求解摩擦接触问题的无网格线性互补模型。在该方法中,接触面方程使用接触点对离散,全局离散系统方程则转化为标准的线性互补问题,可以很方便地使用Lemke算法求解。算例分析证明,本文方法的正确性和有效性。     

Deformation analysis of shallow penetration in clay

A new method of analysis is described for estimating the deformations and strains caused by shallow undrained penetration of piles and caissons in clay. The formulation combines previous analyses for steady, deep penetration, with methods used to compute soil deformations due to near-surface ground loss, and is referred to as the Shallow Strain Path Method (SSPM). Complete analytical solutions for the velocity and strain rates are given for a planar wall, an axisymmetric, closed-ended pile and unplugged, open-ended pile geometries. In these examples, the analyses consider a single source penetrating through the soil at a constant rate, generating a family of penetrometers with rounded tips, referred to as simple wall, pile and tube geometries. Soil deformations and strains are obtained by integrating the velocity and strain rates along the particle paths. The transition from shallow to deep penetration is analysed in detail. Shallow penetration causes heave at the ground surface, while settlements occur only in a thin veneer of material adjacent to the shaft and in a bulb-shaped region around the tip. The size of this region increases with the embedment depth. Deformations inside an open-ended pile/caisson are affected significantly by details of the simple tube wall geometry. © 1997 by John Wiley & Sons, Ltd.     

岩土参数随机场离散的三角形单元局部平均法

将不确定性岩土参数建模为随机场而非传统意义上的随机变量,基于随机场的局部平均理论,提出了用于二维随机场离散的三角形单元局部平均法。通过面积坐标变换和高斯数值积分,给出了三角形单元局部平均随机场协方差矩阵的解析计算方法和数值计算方法。采用算例再现了所提方法的分析过程和有效性,并与传统二维随机场四边形单元离散法进行了对比。结果表明:提出的二维随机场三角形单元离散法能与有限元三角形单元离散法完美结合,随机场单元与有限元单元的对应关系清晰,易于随机有限元程序的编制;对于随机场单元的均值,传统四边形单元离散法与所提方法的计算结果相同;对于随机场单元的方差,传统四边形单元离散法计算结果偏小,所提方法显得更加科学、合理。     

二维接触问题的互补算法及工程应用

接触问题是一种常见的非线性问题,如何能够很好地模拟接触面的变形及受力特性,以及实现对变形体间的接触问题的真实模拟是该领域研究的难点问题。基于二维接触问题的实际物理意义,分别在法向和切向建立等价的互补模型。用非线性互补函数（NCP）中的Fischer-Burmeister（FB）函数将互补函数模型转化为非光滑方程组表达,用常规的Newton法求解。同时,基于高斯积分法可以用较少的积分点达到较高的精度,为了进一步提高求解精度,改善不连续的通病,对面-面接触模型在高斯点上对接触面进行处理,可通过调节积分点数目对求解精度进行控制,方法易于理解,实现方便。在此基础上建立二维接触有限元模型,通过一系列工程算例验证该方法的可行性与有效性。结果表明,与ABAQUS有限元的计算结果相比,该方法有着较高的精度,更真实地反映问题的实际。     

Application of the third medium method for frictionless contact problems in geomechanics

Traditional approaches in contact mechanics demand complicated search algorithms at the interface between the contacting bodies. Recently, a new contact method based on the concept of a third medium has been developed, which overcomes the drawbacks of conventional contact mechanics techniques. This new scheme is based on a space filling mesh, in which the contacting bodies can move and interact. The ability and accuracy of this method in predicting displacements, as well as the contact forces, is validated by solving selected numerical examples. The potential merits of this method for analysing geotechnical problems by the finite element method are addressed.     

A new model of edge-to-edge contact for three dimensional discontinuous deformation analysis

A robust contact theory can be regarded as key to three dimensional discontinuous deformation analyses (3D-DDA). Not only must this theory provide an efficient algorithm to judge the type and location of contacts but also it must be able to present comprehensive formulations for every kind of contact (open, sliding and locked contact). There are six types of contact in three dimensional discontinuous deformation analyses (vertex-to-vertex, vertex-to-edge, vertex-to-face, edge-to-edge, edge-to-face and face-to-face) that can be converted to vertex-to-face and edge-to-edge contacts. This paper presents a new model of edge-to-edge contact to three dimensional discontinuous deformation analyses (3D-DDA). This new model considers both kinds of edge-to-edge contact (cross-over and parallel edge-to-edge contact) and presents a criterion for inter-penetration. Sub matrices of normal and shear spring and friction force are derived by geometrical analysis and penalty method. This new model is implemented in a 3D-DDA computer programme, and the example results demonstrate the validity of the model.     

New point-to-face contact algorithm for 3-D contact problems using the augmented Lagrangian method in 3-D DDA

This paper presents a new point-to-face contact algorithm for contacts between two polyhedrons with planar boundaries. A new discrete numerical method called three-dimensional discontinuous deformation analysis (3-D DDA) is used and formulations of normal contact submatrices based on the proposed algorithm are derived. The presented algorithm is a simple and efficient method and it can be easily coded into a computer program. This approach does not need to use an iterative algorithm in each time step to obtain the contact plane, unlike the ‘Common-Plane’ method applied in the existing 3-D DDA. In the present 3-D DDA method, block contact constraints are enforced using the penalty method. This approach is quite simple, but may lead to inaccuracies that may be large for small values of the penalty number. The penalty method also creates block contact overlap, which violates the physical constraints of the problem. These limitations are overcome by using the augmented Lagrangian method that is used for normal contacts in this research. This point-to-face contact model has been programmed and some illustrative examples are provided to demonstrate the new contact rule between two blocks. A comparison between results obtained by using the augmented Lagrangian method and the penalty method is presented as well.     

接触处理Munjiza方法在数值流形法中的应用

数值流形法（NMM）在接触处理过程包括：接触方式判断,包括点-点和点-线接触;开-闭迭代,确定块体系统的约束状态;以及接触传递,将有效的接触对传递到下一个时步。过程稍显繁琐和耗时,且与真实的物理接触状态有异。由Munjiza所提出的NBS（no binary search）接触检测算法,将单元映射到规则格子中,以链表结构将其有效地连接在一起,只在单元所在格子以及周围格子内部进行接触判断,接触检测效率大为提高,计算量仅随单元数线性增长,内存需求也很低。在计算接触力时以所定义的势为媒介,用重叠面积来衡量接触力的大小,属于分布式接触力,更接近于实际,避免了原NMM接触处理过程的繁琐。因此,将其作为一种平行的接触处理方法引入到NMM中,并以算例验证,证实了算法的可行性。     

Solution of quasi‐static large‐strain problems by the material point method

A procedure for solving quasi‐static large‐strain problems by the material point method is presented. Owing to the Lagrangian–Eulerian features of the method, problems associated with excessive mesh distortions that develop in the Lagrangian formulations of the finite element method are avoided. Three‐dimensional problems are solved utilizing 15‐noded prismatic and 10‐noded tetrahedral elements with quadratic interpolation functions as well as an implicit integration scheme. An algorithm for exploiting the numerical integration procedure on the computational mesh is proposed. Several numerical examples are shown. Copyright © 2010 John Wiley & Sons, Ltd.     

Simulation of geomembrane response to settlement in landfills by using the material point method

Because of the multiple layers of dissimilar materials and large deformations involved in the subsidence of a landfill system, large-scale computer simulation of the geomechanical response to subsidence with the use of conventional numerical methods are problematic. The Material Point (MPM), which was recently developed for dynamic problems such as penetration and perforation, is a newly emerging numerical method. The MPM is modified in this paper to simulate the geomechanical response of a landfill cover system that includes a geomembrane under quasi-static loading conditions. Sample problems, for which an analytical solution is available with certain assumptions, are considered to demonstrate the proposed solution procedure. Future work is discussed based on current research results. Copyright © 1999 John Wiley & Sons, Ltd.     

均质黏土中圆形平板锚的抗拉承载力分析

基于网格重新生成和场变量映射的大变形有限元模型,探索了立即脱离和无脱离两种典型条件下均质黏土中圆形平板锚的抗拉承载力。与小变形有限元比较,大变形分析克服了锚周围土体初始网格畸变的不利影响,能够追踪平板锚整个拔出过程中抗拉力的变化。通过具体算例,考察平板锚表面摩擦性质和上覆土重等因素对立即脱离工况承载力的影响程度,指出有重土中深锚的承载力小于无重土中对应的承载力与上覆土重之和,其上限是无脱离条件下的承载力。计算结果表明：土重对无脱离条件下的承载力影响很小,进而给出了无脱离承载力系数与初始埋深的关系曲线。