0.   引言

随着计算机技术的发展，数值模拟在工程中的应用越来越广泛，有限元法作为其中最具代表性的方法已经相对成熟，在固体力学、流体力学等领域发挥着重要的作用。传统的有限元方法基于小变形理论，采用有限的单元对结构进行离散，根据能量原理提出基本方程进行求解，具有较高的精度和计算效率。但对于滑坡、基础贯入、触探等岩土工程问题，其破坏过程常常伴随着较大变形的发生，采用传统基于小变形理论的有限元法会出现网格扭曲、畸变，导致计算误差较大甚至不收敛。

为了克服传统有限元方法在模拟岩土大变形问题时的局限性，学者们提出了许多大变形数值方法。Hirt等提出了任意拉格朗日-欧拉法（Arbitrary Lagrangian Eulerian, ALE），这种方法使用网格来表征结构，网格可以采用拉格朗日格式，也可以采用欧拉格式，或是任意可以减小网格畸变的形式^[1]。采用拉格朗日格式时，网格与物质固连，在大变形分析过程中需要通过网格更新、变量映射技术实现计算结果在不同分析步间的传递，这种技术又称为RITSS（remeshing and interpolation technique with small strain）技术^[2]，已被用于基础贯入^[3]等岩土大变形问题中；采用欧拉格式的网格时，计算过程中网格不产生变形，通过网格的充填状态来反映物质流动，多用于流体问题^[4]及贯入问题^[5]。不同格式的ALE方法的共性问题在于该方法中网格节点的拓扑关系固定，无法模拟自由面剧烈变化的问题。另一种常用的大变形方法是无网格法^[6]，与ALE方法不同，无网格法用粒子来表征结构，这种处理优化了粒子间的拓扑关系，适用于大变形、复杂接触等问题。无网格法中一点的状态是通过支持域内其他点构造逼近函数计算得到的，不同的点支持域内包含的粒子数量不同，逼近函数及其导数也不相同，导致计算效率较低^[7]，同时无网格法采用质点积分而非高斯积分，计算精度也较低。为了同时利用网格描述和粒子描述的优点，一类新的数值模拟方法被提出，其特点是采用粒子来表征物质，计算过程在网格上进行，典型的方法有物质点法（material point method, MPM）和粒子有限元法（particle finite element method, PFEM）。MPM法中物质点携带所有的信息，计算在背景欧拉网格上执行，计算开始时信息首先由物质点映射到网格节点上，计算完成后再由网格节点映射回物质点^[8]。在岩土工程领域，物质点法已被用于滑坡等特大变形问题的模拟^[9-10]，该方法采用质点积分，不可避免的会产生积分误差，在应用时宜考虑其适用性。PFEM法中每一增量分析计算都通过标准的有限元求解器进行，继承了有限元方法在计算效率和精度上的优势；同时由于采用粒子表征物质，在自由面的刻画、复杂边界问题时也具备类似无网格方法的优势。PFEM法首先由Oñate等提出，运用于流体-结构相互作用问题^[11]。张雪等^[12]将粒子有限元方法应用于固体力学领域，模拟了颗粒柱坍塌、土质滑坡^[13]、海底滑坡^[14]等问题。Li等^[15]引入约束Delaunay算法在ABAQUS中实现了PFEM流程，对颗粒柱塌落及滑坡问题进行了模拟。Yuan等编写python脚本在ABAQUS中实现了PFEM计算过程，模拟了T-bar贯入、地基承载力、海底管道等问题^[16]。

另一方面，岩土体在大变形过程中往往具有应变软化特性和应变局部化现象，会导致有限元数值求解的不适定及网格依赖性等问题。为保持问题求解的适定性，需要在本构方程中引入正则化机制，而采用Cosserat连续体理论是引入正则化机制的有效方法之一^[17]。为此，本文将PFEM计算方法与Cosserat连续体理论相结合，发展了Cosserat-PFEM方法。此方法既具有PFEM法在处理大变形问题、刻画自由面变化方面的优势，又能保持问题求解的适定性，适用于岩土体渐进破坏大变形的模拟。

1.   Cosserat-PFEM关键技术实现

1.1   Cosserat连续体

以二维Cosserat连续体为例简要介绍其基本理论。与经典连续体相比，在Cosserat连续体中引入了旋转自由度，使得平面内每个点具有3个自由度（即两个平动自由度，一个转动自由度），相应的位移、应力和应变定义如下:

式中：${w_z}$为旋转自由度；${k}_{zx}，{k}_{zy}$为微曲率；${m}_{zx}，{m}_{zy}$为对偶应力；${l_{\text{c}}}$为内部长度参数。

对于均质各向同性弹性体，Cosserat连续体应力、应变满足胡克定律：

式中：$\lambda$为拉梅常数，且有$\lambda = 2G\nu /(1 - 2\nu )$，$G，\nu$分别为经典连续体中的剪切模量和泊松比；${G_{\text{c}}}$为Cosserat剪切模量。根据文献[17]对${G_{\text{c}}}$和${l_{\text{c}}}$的取值研究，取${G_{\text{c}}}$/G≥0.5；从细观本质上讲，内部长度参数${l_{\text{c}}}$的大小反映了土体颗粒尺寸及有关参数的影响，在宏观有限元计算时，取值原则是在保证正则化效果和数值计算收敛性的前提下尽量取小值，经验取值为0.01~0.1倍的模型尺寸。

平衡方程和应变-位移关系为

式中：L为微分算子矩阵，

由上述理论描述可知，在Cosserat连续体本构描述中引入了内部长度参数，由此保持了应变局部化问题求解的适定性，可以模拟剪切带破坏有关的岩土工程问题，在此不再详述。

1.2   Cosserat-PFEM方法实现

PFEM法采用拉格朗日描述，通过边界识别、网格重新划分、变量映射实现信息在新旧网格间的传递。本文将Cosserat连续体理论与PFEM结合，对ABAQUS软件进行二次开发实现相关算法，其计算流程如图 1所示。首先建立Cosserat连续体有限元模型，进行常规的小变形分析，分析完成后导出变形后节点（即粒子）坐标；然后利用边界识别程序对导出的粒子集进行边界识别获得新的计算域，在新计算域内重新进行网格划分；再通过变量映射技术将信息从旧网格映射到新网格上，施加荷载和边界条件进行下一次小变形计算直到分析完成。下面对关键技术实现方法进行介绍。

图  1  Cosserat-PFEM分析流程图

Figure  1.  Flow chart of Cosserat-PFEM analysis

下载: 全尺寸图片 幻灯片

1.3   边界识别技术

边界识别是PFEM方法中最重要的技术，识别的边界可以很好的反映分析过程中自由面的变化，也给荷载、边界条件的施加带来了便利。本文采用α-shape程序进行边界识别^[18]，α-shape方法最早被应用于计算机图形学中，本文采用的是一种较为直接的实现形式。以二维情况为例，α-shape法的基本法则如下：考虑二维空间中一系列具有特征间隔h的点集合，如果能够找到一个半径大于αh的空心圆使得某些点均位于圆上（其中α是给定的参量，其值一般介于1.3到1.6之间），那么这些点就认为是边界点。针对同一组粒子集，不同α值识别得到的边界也不相同，会导致边界拾取的精确度不同，如图 2所示的粒子集（后面3.1节中的地基承载力计算模型），当选取不同的α值时，边界拾取的效果也不相同。当α取0.4时，识别的边界包含了不属于计算域内的部分；随着α值的增大，识别的边界逐渐趋近于真实的边界，当α取1.3时边界最为理想；α继续增大，开始逐渐有计算域内的点被舍弃，当α取1.5时，已有很大部分区域遗失。因此，在应用过程中应注意α的合理取值，以避免发生较大的物质损失。

图  2  α取不同值时的边界识别结果

Figure  2.  Boundary recognition results with different α

下载: 全尺寸图片 幻灯片

1.4   用户自定义单元及网格划分技术

网格划分是大变形分析中重要的一环，可通过ABAQUS内置网格生成方法来实现。需要注意的是前述边界识别技术是基于粒子间距离长短进行的，粒子分布的疏密对边界识别质量有很大影响，因此在网格划分时应尽量保证节点分布均匀。

由上可见，本文将边界识别和网格划分作为两个完全独立的过程进行，这摈弃了前人工作中PFEM法仅可采用三角形单元的限制，为采用多种单元类型打下了基础。

由于ABAQUS中没有内置的Cosserat单元，需要基于前述Cosserat连续体基本理论进行二次开发，即编写用户自定义单元子程序UEL。考虑到三角形单元对应变局部化问题模拟具有一定的倾向性，本文使用了开发的四边形八节点Cosserat单元。

1.5   变量映射技术

变量映射过程是新旧网格间信息传递的桥梁，在大变形问题求解过程中，经常需要进行多次变量映射，应选用恰当的映射方法保证映射精度。

由于采用了用户自定义单元子程序UEL，引入了新的自定义变量，变量映射无法直接通过ABAQUS自动进行。本文采用与Map Solution技术类似的映射原理，首先将旧单元积分点处的场变量计算结果外插到单元节点上，对于多个单元共用节点的情况进行平均化处理；随后进行网格重新划分，判断新网格高斯点在旧网格中的位置，再在旧网格中通过形函数插值得到新网格高斯点处的变量值^[19]。映射过程中所需要的数据文件有：新旧网格的节点信息、新旧网格的单元信息以及旧网格上的计算结果。

2.   案例分析

2.1   塌落问题分析

为验证本文开发程序的正确性，首先模拟了与文献[12]类似的塌落问题。本文的工作将Cosserat连续体理论与PFEM结合，是基于对ABAQUS软件进行二次开发实现的。为了分析一般的岩土材料，而不仅仅是无黏性的颗粒材料，选择使用了Drucker-Prager本构模型，为了使计算分析可行，依据ABAQUS软件的特点，在参数的选择上与参考文献[12]有所区别，其中主要参数（如内摩擦角、密度）一样，只是增加了很小的一个黏聚力值0.1 kPa，这对结果有点影响，但不大。

如图 3所示，计算模型长L=8 m，高H=4 m，材料参数如表 1所示。初始时刻，约束模型左右两侧边界的水平位移和底部的竖向位移，随后移除右侧的水平约束，模型在自重的作用下开始发生变形。本节考虑移除约束后0~1 s内模型的变形状态，全过程共进行18次小变形分析，每次小变形分析完成后进行一次边界识别和网格重划分。

图  3  计算模型

Figure  3.  Schematic diagram of computational model

下载: 全尺寸图片 幻灯片

表  1  材料参数表

Table  1.  Parameters of materials

弹性模量 E/kPa	泊松比 $\nu$	黏聚力 c/kPa	内摩擦角 $\varphi$/(°)	密度 $\rho$/(kg$\cdot$m-3)
30000	0.35	0.1	19.0	1490

下载: 导出CSV 

| 显示表格

文献[12]显示，移除右端约束后，颗粒柱右端随即开始坍塌，刚开始短时间内变形较小，上表面存在一部分未扰动区；随着坍塌的进行，上表面的未扰动区逐渐减小；在整个坍塌过程中，颗粒柱左下角存在一个区域始终保持静止状态；当颗粒材料达到稳定状态时，最终构型的前段倾斜角非常接近内摩擦角。

本文不同时刻的变形和位移云图如图 4所示。由图可见，右侧约束移除后的0.1 s内，变形并不明显，这段时间主要是由起始的稳定状态向不稳定状态进行转变；0.25 s时，模型的右下方有了较明显的位移；0.5 s时，模型右下角有明显位移，达1.54 m，右上表面有相应的沉陷；0.75 s时右下角已产生2.98 m的位移，此时已没有上、右边界的区分，两者已连成一条几近平滑的曲线；1.0 s时，右下角峰值位移已达到4.66 m，原来的上表面和右侧面变为一条平滑的自由面。这个过程反映了塌落破坏的全过程，这个过程与文献[12]的模拟结果高度相似，最后自由面的倾角约21°，略大于内摩擦角，这与非零的黏聚力有关。因此，总体上看变形过程、滑移线状态与文献[12]结果高度相似，说明了本文所提出的方法和开发的程序的可靠性。

图  4  不同时刻变形和位移s云图

Figure  4.  Deformation and displacement contours at different moments

下载: 全尺寸图片 幻灯片

2.2   地基承载力问题模拟

为验证本文实现的粒子有限元方法及不同单元的有效性，对地基承载力问题进行了模拟，并与前人的大变形研究结果^[20-21]进行了对比。

计算模型如图 5所示，根据对称性取半结构进行模拟，简化为二维平面应变问题，模型左右两侧施加水平向约束，底部施加水平和竖向约束。基础宽0.5 m，地基长6 m，深4 m。土体采用理想弹塑性模型，不考虑自重，材料参数如表 2所示，为了与已有的工作进行比较说明，计算参数选为与文献[20，21]一致。采用位移加载模式，加载过程中对基础施加竖直向下的位移，基础每下沉0.1 m，采用上述实现的PFEM方法进行边界识别和网格划分，直至基础下沉深度达到0.5 m。

图  5  计算模型示意图

Figure  5.  Schematic diagram of computational model

下载: 全尺寸图片 幻灯片

表  2  地基土材料参数

Table  2.  Parameters of soil

参数	弹性模量E/kPa	泊松比$\nu$	抗剪强度Su/kPa
值	100	0.49	1

下载: 导出CSV 

| 显示表格

对于地基承载力问题，Nazem等^[20]和Tian等^[21]分别采用ALE方法和RITSS技术进行过研究，将其结果绘制于图 6中。这里V代表竖向荷载，B是半基础宽，w是基础的沉降，s_u是土的抗剪强度，其中B和s_u取值都为1，为了与其他文献中变化的B和s_u的归一化处理作对比，本文也采用了类似归一化处理的形式。Nazem等^[20]的结果中，基础下沉约0.22 m时，发生计算不收敛现象；Tian等^[21]使用RITSS技术计算到基础下沉深度为0.3 m。以三角形单元为例，本文的计算下沉深度达到了预设的位移0.5 m，体现了PFEM法模拟大变形问题的优势。

图  6  承载力-贯入深度曲线图

Figure  6.  Bearing capacity-penetration depth curves

下载: 全尺寸图片 幻灯片

现有的PFEM成果大多采用常规的三角形单元，这是由于PFEM最初被提出时，边界识别与网格划分存在于同一个过程中，边界识别前首先要生成Delaunay三角剖分，再对Delaunay三角进行处理得到边界，两者有着密不可分的关系，因此限制了单元类型的选择。张雪等^[12]指出，对于固体力学大变形问题，在边界识别过程中生成的网格质量较差，难以满足计算要求。

本文将边界识别和网格划分视为两个独立的过程，边界识别仅完成对离散粒子集边界的识别，得到边界后再像常规有限元方法一样对新计算域进行网格划分。分别使用ABAQUS内置的线性三角形单元、二阶三角形单元、线性四边形单元、二阶四边形单元对前述地基承载力问题进行模拟。不同单元类型下承载力-位移曲线如图 7所示，结果表明各种单元计算结果差不多。

图  7  不同单元类型对应承载力曲线图

Figure  7.  Corresponding bearing capacity curves of different element types

下载: 全尺寸图片 幻灯片

比较而言，线性三角形和四边形单元采用了选择减缩积分，总体精度较差，映射误差也更大一些，导致每次重新划分网格后承载力跳跃回落较大一些；而采用高阶单元时精度较高，曲线更光滑，同时可以克服网格自锁对结果的影响，更有利于岩土大变形和应变局部化问题的模拟。从承载力曲线看，与三角形单元相比，四边形单元的结果更接近高阶单元的结果，计算精度更高，两者比较采用四边形单元比三角形单元精度提高约6.5%。另外，从应变局部化的模拟上，采用三角形单元计算时倾向于发展与三角形单元本身倾斜方向一致的剪切带，也就是三角形单元网格排列的倾向会影响剪切带的发展，影响剪切带模拟的客观性。因而，本文后续算例采用了所发展的四边形八节点单元进行模拟。

另外，图 6，7中显示承载力位移曲线在重新划分网格后有台阶跳跃现象，这与重划网格后变量映射产生的误差有关，三角形单元更大一些。可以通过采用高阶单元提高精度、改进映射手段、减小重划步长进一步减小误差，当然，这样可能会增加计算工作量。

2.3   边坡渐进破坏分析

（1）边坡分析模型

边坡分析模型如图 8所示，在坡顶的基础上作用一竖向加载，边坡底部固定，右侧边界水平约束，竖向自由，材料参数如表 3所示。此算例考虑黏聚力软化效应如下：

图  8  计算模型示意图

Figure  8.  Schematic diagram of computatioanl model

下载: 全尺寸图片 幻灯片

表  3  材料参数表

Table  3.  Parameters of soil

材料参数	值		材料参数	值
弹性模量E/kPa	50000		内摩擦角$\varphi$/(°)	25
泊松比$\nu$	0.30		剪胀角$\psi$/(°)	0
初始黏聚力c0/ kPa	50		Cosserat剪切模量Gc/ kPa	10000
黏聚力软化模量hp/ kPa	-30		内部长度参数lc/ m	0.06

下载: 导出CSV 

| 显示表格

式中：${c_0}$为初始黏聚力，${\varepsilon _{\text{p}}}$为等效塑性应变，${h_{\text{p}}}$为黏聚力随塑性变形而软化的软化模量。

分别采用3种尺寸的网格，即0.2，0.3，0.4 m进行计算，以考虑网格尺寸对计算结果的影响。

（2）计算结果分析

计算预设指定向下的位移为0.5 m，首先采用Cosserat-PFEM方法进行分析，图 9给出了3种网格下的荷载位移曲线，图 10给出了3种网格下坡顶指定位移为0.5 m时的等效塑性应变图。可见，3种网格下Cosserat-PFEM法分析得到的结果较为一致，显示了与网格无关的结果。文献[17]中做了类似的边坡渐进破坏算例分析，只不过采用的是小变形的分析方法，这里的Cosserat-PFEM大变形模拟得到剪切带渐进破坏过程及克服了网格依赖性的特性与文献[17]的结果是一致的。为了显示Cosserat-PFEM方法的优越性，图 11，12给出了常规FEM和PFEM方法分析得到的0.2 m和0.4 m两种网格尺寸下的荷载位移曲线，可见具有应变软化和非关联塑性的岩土材料，会表现出过早的分叉，数值计算较早的遇到困难，在本例条件下难以越过承载力峰值完成完整破坏过程的分析，而且数值解具有明显的网格依赖性。比较而言，Cosserat连续体理论由于在本构方程中引入了正则化机制，可以计算达到0.5 m的收敛位移，这充分显示在岩土体大变形分析中引入正则化机制的重要性。当然，图 11，12也表明，与FEM相比，总体上PFEM能计算到更大的变形。

图  9  Cosserat-PFEM法荷载位移曲线

Figure  9.  Load-displacement curves of Cosserat-PFEM method

下载: 全尺寸图片 幻灯片

图  10  Cosserat-PFEM计算终止时等效塑性应变图

Figure  10.  Diagram of equivalent plastic strain at the end of calculation by Cosserat-PFEM

下载: 全尺寸图片 幻灯片

图  11  FEM法荷载位移曲线

Figure  11.  Load-displacement curves of FEM method

下载: 全尺寸图片 幻灯片

图  12  PFEM法荷载位移曲线

Figure  12.  Load-displacement curves of PFEM method

下载: 全尺寸图片 幻灯片

3.   结论

本文将PFEM计算方法与Cosserat连续体理论相结合，发展了Cosserat-PFEM方法。同时，将边界识别和网格划分视作是两个独立的过程，采用了多种单元形态网格进行了大变形分析模拟，得到以下两点结论。

（1）提出并实现的Cosserat-PFEM方法既具有PFEM法解决大变形问题的优势又具有保持应变软化和应变局部化问题适定性的能力，能较好的适用于大变形渐进破坏问题的模拟。

（2）实现的算法中可以采用四边形等多形态单元，与传统的三角形单元相比，采用高阶四边形单元可以提高计算精度，克服三角形单元对应变局部化问题模拟的倾向性。

应该指出，由于映射误差的存在，承载力位移曲线在重划网格和变量映射之后会出现跳跃，除了采用高阶单元提高精度、减小重划步长进一步减小误差外，应采取相应措施改进映射手段进一步减小误差，这也是今后努力的方向。