静力触探过程中的大变形弹塑性有限元分析

吴鹏祁敢

中铁六局集团北京铁路建设有限公司北京100036

摘要：本文利用有限元软件ABAQUS为工具，采用任意拉格朗日-欧拉（ALE）运动学方法描述，考虑土体的大变形特性，模拟静力触探过程中贯入阻力及沿探杆长度的应力场和位移场。结果表明，贯入阻力随着贯入深度的增加线增大而后基本保持恒定，随土层深度的增加而增大，与土的内摩擦角和内聚力正相关。通过对比大变形和小变形计算结果，表明对于静力触探，采用考虑大变形的ALE网格重划分技术可以较好的模拟锥尖处的应力集中现象。

关键词：静力触探；拉格朗日-欧拉（ALE）方法；大变形

引言：本文以怀柔火工区铁路专用线一期工程、二期工程、张家口大唐电力铁路专用线工程、河北宣化钢铁集团铁路专用线工程、天津静海东方红路下穿京沪铁路立交桥工程的路基检测工程实践为课题基础，联合北京交大、首钢地质勘察院、北京金地通检测技术有限公司等多家单位进行研究，利用有限元软件ABAQUS为工具，采用任意拉格朗日-欧拉（ALE）方法描述的有限元方法，考虑土体的大变形特性，模拟静力触探过程中贯入阻力及沿探杆长度的应力场和位移场，并分析了各土性参数对应力场和位移场的影响。

1有限元大变形问题的ALE描述简介

在大变形分析中，位移与应变一般呈现非线性关系，使得以往的的小位移小应变分析方法不能适用于大变形分析计算。目前，大变形分析的方法一般采用增量分析法，即在一般弹塑性分析方程中加入增量部分，随着加载变化而进行构型变化。所以其求解方法以增量法为基础主要分为三种：Lagrange法、Euler法和Lagrange-Euler混合方法。

在拉格朗日描述（即物质描述或L描述）中，计算网格固定在物体上随物体一起运动，即网格点与物质点在物体的变形过程中始终保持重合，因此物质点与网格点之间不存在相对运动（即迁移运动，也称对流运动）。对于物质的变形，一般采用应力和应变增量表示。这大大地简化了控制方程的求解过程，而且能准确描述物体的移动界面，并可跟踪质点的运动轨迹。但在涉及到特大变形的问题中，物质的扭曲将导致计算网格的畸形而使得计算失败。Lagrange法又可分为全拉格朗日列式法（TL法，TotallagrangianFormulation）和更新的拉格朗日列式法（UL法，Updated-LangrangianFormulation）。全拉格朗日列式法（TL法）以初始（0时刻）构型作为应力与应变描述的参考构型，因而，采用柯西应力（增量）和格林应变（增量）。更新的拉格朗日列式法（L法）以t时刻（即现时构型）为参考构型，也就是说参考构型是变化的，因而，采用现时柯西应力（增量）和现时格林应变（增量）。

Euler法是着眼于空间点，设在空间中的每一个点上描述出流体运动随时间的变化状况。在Euler法有限元计算中，单元在空间上是固定的，材料在空间流动，单元不会随着材料的运动而扭曲，形成了一个应力应变场。网格固定在空间中，即计算网格在物体的变形过程中保持不变，易于处理物质的扭曲。但对运动界面需要引入非常复杂的数学映射，将可能导致较大的误差计算。但是，由于材料的流动，本构方程的处理和更新复杂。同时，Euler单元在处理移动边界和相互作用的问题上有较大限制。

纯拉格朗日和纯欧拉方法都存在严重的缺陷，但也具有各自的优势。如果能将二者有机地结合，充分吸收各自的优势，克服各自的缺点，则可解决一大批只用纯拉格朗日和纯欧拉描述所解决不了的问题。任意拉格朗日-欧拉（ALE）方法就是基于此目的最早由Noh（1964）[5]以耦合欧拉-拉格朗日的术语提出的，并用有限差分法求解带有移动边界的二维流体动力学问题。在Noh的研究工作中，网格点可以随物质点一起运动，但也可以在空间中固定不动，甚至网格点可以在一个方向上固定，而在另一个方向上随物体一起运动，因此ALE描述也被称为耦合欧拉-拉格朗日描述。例如在液体表面波的传播问题中，网格点在垂向随物质点一起运动，而在水平方向上固定不动。这样可很容易描述液体表面的运动，而且网格不会发生扭曲。在ALE描述中计算网格可以在空间中以任意的形式运动，即可以独立于物质坐标系和空间坐标系运动。这样通过规定合适的网格运动形式可以准确地描述物体的移动界面，并维持单元的合理形状。类似于E描述，在ALE描述下的控制方程中也将出现对流项，因此也可能得到振荡解，需要进行相应的数值处理。纯拉格朗日和纯欧拉描述实际上是ALE描述的两个特例，即当网格的运动速度等于物体的运动速度时就退化为拉格朗日描述，而当网格固定于空间不动时就退化为欧拉描述。

ALE最早是为了解决流体动力学问题而引入的，并且使用有限差分法。由于核反应堆结构安全分析的需要，Donea[6]、Belytschko[7]等人分别将ALE法引入有限元法中，用以求解流体与结构相互作用问题。Hughes等人[8]建立了ALE描述的运动学理论，并使用有限元法解决了粘性不可压缩流体流动和自由表面流动问题。

ALE描述中，参考构形是已知的，而初始构形和现时构形都是待求解的。因此ALE法尤其适合于在初始构形和现时构形都未知的问题中使用，如接触问题等。目前，ALE法已被应用于固体力学领域中求解大变形问题，如碰撞、接触和弹性断裂力学、路径相关材料（如弹塑性材料等）、加工成型等。另外对网格的运动以及控制方程的求解策略等问题也进行了多方面的研究。

2ABAQUS软件中Mohr-Coulomb模型简介

ABAQUS中莫尔一库仑模型，认为材料屈服服从莫尔一库仑屈服准则，采用不关联流动法则，流动势面在子午应力平面上成双曲线函数，而在偏应力平面(即兀平面)上，为分段的椭圆曲线，如图1~图3所示。该模型反映拉压不同屈服强度，而且可以考虑料剪胀和软化的特性。正是由于莫尔一库仑模型能较好地描述岩土类材料的材料特性，且使用方便，在岩土工程界获得了广泛的应用，基于此，本章选用该模型。本章应力以压为负，拉为正。

2.1屈服准则

在简单应力状态下，莫尔一库仑屈服准则表示为

图1莫尔-库伦屈服模型

按照上式计算的可保证塑性势面在π平面受拉和受压的角点上与屈服面相切。用户也可以指定的大小，但其范围必须为：。图5给出了不同的大小对应的塑性势面。当材料内摩擦角和剪胀角相同，且。取值较小时，在子午应力面上屈服面与流动势面几乎重合，此时流动法则即为相关联的流动法则，当材料内摩擦角和剪胀角不同时即为不相关联的流动法则，通常情况下，剪胀角与内摩擦角不相等，流动势函数与屈服面函数不是相关联的。

图4Mohr-Coulomb模型在子午平面上的塑性流动势

图5几种模型在π平面上的偏心率e

3静力触探过程有限元分析

3.1基本假定

（1）不考虑桩土体固结，采用总应力法分析；

（2）结合实际静力触探过程，为了简化计算，采用二维轴对称问题来模拟；

（3）标准锥采用解析刚体模拟；

（4）土体为连续的弹塑性体，符合Mohr-Coulomb模型，塑性流动遵从非相关联的流动法则考虑几何和材料非线性，采用四结点二维轴对称单元。

（4）考虑土体初始应力场的影响。

（5）桩–土界面采用面–面摩擦接触单元，接触面的摩擦类型为库仑摩擦。

3.2有限元计算模型

（1）计算区域

静力触探的工作过程属于空间轴对称问题，即荷载和变形都成空间轴对称，以探杆中心线为对称轴。

本文探头取Ⅱ-1型，探头角度60o，直径36mm，底面积10cm2。为了减小边界条件的影响，取土体的分析区域高度为1.5m，宽度为1.5m。由于本模型只分析土体的指标，不考虑椎体的应力应变指标。所以椎体采用了刚体，在锥顶设置以参考点，在计算过程中通过参考点的位移来模拟椎体的贯入过程。

（2）荷载及边界条件

在静力触探贯入过程中，处于锥尖处的土体首先劈裂，然后被挤压到四周，随着探头的不断贯入，探杆周围的土体还受到侧摩擦力的牵引，因此土体的受力及变形过程相当复杂，完全真实地再现静力触探过程是很困难的。本章模拟静力触探贯入土体后在垂直贯入力作用下，贯入器周围土体内的响应。考虑土体的重力场，同时在孔壁处考虑向下的摩擦力的作用。贯入器与土体在刚度与强度上存在的较大差异，在外力作用下，其接触面有可能会产生相对滑动或脱离。为此在有限元数值分析中，在两者之间设置接触面单元，并采用Coulomb摩擦接触模型描述界面的接触特性地表面边界为自由边界，地面和外侧边界为固定约束边界。

(3)计算参数

土的计算参数见表1。但在参数影响分析中，相应改变分析因子。

图6贯入过程中椎体阻力的变化

3.3.2土体应力状态

图7和图8为贯入36cm计算求得的径向应力、轴向应力和切向应力云图。由图可知，对于径向应力主要影响区域为径向为3~5倍椎体直径，轴向锥尖以下2~3倍椎体直径范围内。在探头附近存在一个应力集中区域，最大应力达到2.0MPa；对于轴向应力主要影响区域为径向为2~3倍椎体直径，轴向锥尖以下3~5倍椎体直径范围内。在探头附近存在一个应力集中区域，最大应力达到2.0MPa（本算例中没有考虑破坏应力）。

图9土体径向位移u

图10土体轴向位移v

4大变形理论和小变形理论计算结果对比

土体的基本计算参数同上，分别利用大变形理论和小变形理论针对圆锥在贯入深度为36cm时引起的土体内的应力、位移以及贯入阻力计算。

大变形给出的贯入阻力（5.44KN）明显大于小变形（1.20KN），对于大变形理论，土体最大应力出现在锥尖处（径向2MPa），大于小变形理论（径向1.45MPa），说明大变形理论更能模拟锥尖处的应力集中现象。从图11可以看出，基于小变形理论的普通网格技术，较之大变形理论的ALE网格技术单元畸变很大，小变形蓝色部分单元已经由四边形扭曲成8字形，相应的计算精度也下降。

图12土层深度对贯入阻力的影响

6土体物理参数对贯入阻力的影响

土体其他参数不变，分别改变粘土内摩擦角和内聚力C，研究其对锥体阻力的影响。

图14粘土内聚力对贯入阻力的影响

从图13和图14可以看出，随着粘土内摩擦角和内聚力增加，贯入阻力增加；但内聚力大于40kpa后对贯入阻力的影响较小。

结束语：

通过利用ABAQUS有限元分析软件对静力触探过程的数值模拟，表明贯入阻力随着贯入深度的增加线增大而后基本保持恒定，随土层深度的增加而增大，与土的内摩擦角和内聚力正相关。通过对比大变形和小变形计算结果，表明对于静力触探，采用考虑大变形的ALE网格重划分技术可以较好的模拟锥尖处的应力集中现象。

参考文献：

[1]蒋明镜,沈珠江.饱和黏土的塑性大变形有限元平面固结分析[J].河海大学学报，1998，26(1):73-77。

[2]李波,栾茂田,静力触探工作特性的大变形弹塑性数值分析.武汉理工大学学报2009，31（2），

[3]马海鹏,静力触探比贯入阻力与土体抗剪强度相关关系研究[D].中国水利水电科学研究院，2013

[4]刘斌,静力触探的机理研究[D]，天津大学，2011.12