一种基于ABAQUS的岩土工程数值模拟方法与流程

本发明属于岩土工程计算机辅助设计技术领域，尤其涉一种基于abaqus的岩土工程数值模拟方法。

背景技术：

对岩土工程进行分析时，不大可能用解析方法来完成，只能采用实验和数值模拟计算的方法。实验研究虽然能提供大量宝贵的研究资料，但是需花费大量的人力、物力，实验周期往往也相当长，所得到实验成果往往相当有限，需进行处理才能得到可用于分析工程岩土介质的宏观力学参数。岩土工程数值模拟的基本方法是有限元法，尤其是对于大范围的工程施工效应的动态分析，有限元法是十分有效的。

目前常用于岩土工程的土体本构模型有：邓肯-张(dc)模型、莫尔-库仑(mc)模型、修正剑桥(mcc)模型、硬化土(hardeningsoil)模型(简称hs模型)等。

dc模型为非线性弹性模型，可以反映土体应力、应变的非线性特性，但却不能反映土体的塑性应变，也不能不同的应力路径；mc模型是目前岩土力学中应用最广和应用时间最长的土体模型之一，但mc模型过高地估计了岩土的抗拉强度；mcc模型参数较多且较难确定。

以上本构模型均未考虑材料硬化，即把材料当做理想塑性材料。hs模型考虑了材料的硬化，不仅可以反映土体应力、应变的非线性特性和复杂的应力路径，而且模型参数可以从常规三轴试验获得，模型参数简单，因此成为岩土工程分析应用最广泛，也是最准确的土体本构模型之一。

abaqus是国际上最先进的有限元软件之一，它具有丰富的适合岩土工程分析的本构模型，多达数百种的各种单位类型，非凡的非线性分析和耦合场分析能力，使它成为岩土工程分析领域的有力工具。

但是岩土工程中有些常用的本构模型abaqus并未包含，这就限制了abauqs在岩土工程中的应用，但abaqus提供的umat子程序可以灵活地创建自定义材料模型，采用这种方式可以更合理有效地模拟土体的本构特性。

技术实现要素：

为了解决上述问题，本发明提出了一种基于abaqus的岩土工程数值模拟方法。

本发明采用的技术方案如下：

一种基于abaqus的岩土工程数值模拟方法，包括以下步骤：

步骤一：基于euler积分算法，利用fortran语言编制hs本构模型的子程序，实现abaqus软件的二次开发；

步骤二：将考虑材料硬化效应的hs模型编入abaqus软件的umat子程序中，实现土体的数值模拟；

其中，步骤一中的euler积分算法根据已知第n增量步的所有变量值，给定时间步长增量和总应变增量，求得第n+1增量步的满足本构方程的应力解，假设第n增量步得到的应力为σn，塑性应变为总应变为εn，则求取第n+1增量步的相关参数的具体实现步骤如下：

(1.1)通过输入的应力状态参数计算应力状态不变量、屈服函数对应力的偏导数和弹性矩阵d^e；

(1.2)进行弹性试算，假设应变增量dε都是弹性应变，那么试探应力σ^trial为：

σ^trial＝σn+d^edε；

(1.3)计算试探应力对应的应力状态不变量，然后根据应力状态不变量计算塑性乘子λ，根据塑性乘子的正负号判定加卸载状态，若塑性乘子大于0，则处于加载状态，若塑性乘子小于0，则处于卸载状态；

(1.4)根据步骤1.3中的判定结果进行选择，若处于卸载状态，则不进行塑性修正，且将雅可比矩阵d设为弹性矩阵d^e，然后跳转到步骤1.7；若处于加载状态，则执行步骤1.5进行塑性修正；

(1.5)塑性修正，修正后的第n+1增量步的应力σn+1、塑性应变和总应变εn+1的更新式如下：

(1.6)更新并返回雅克比矩阵。获得第n+1步的应力值后，要进行雅克比矩阵的更新，以供abaqus主程序进行下一个增量步的计算。雅可比矩阵的更新公式如下：

其中，d^ep是更新后的雅克比矩阵，σii为主应力，σm为平均有效应力，ρ为新增的状态变量，与超固结比相关，a为材料参数，λ为压缩参数，κ为回弹参数，e0为参照应力对应的孔隙比；

(1.7)将应变分量、等效塑性应变等更新之后放在状态变量statev(1：nstatev)中存储，以便在后处理中输出查看。

所述步骤二包括以下实现步骤：

(2.1)提取岩土工程结构的几何参数，在abaqus软件中建立几何模型；

(2.2)在abaqus软件的材料输入界面，使用关键词“usermaterial”，表示定义用户材料属性，根据试验结果，将材料参数输入到数值模型中；

(2.3)在abaqus的load模块中，对模型的边界条件进行设置，并施加初始荷载；

(2.4)在abaqus中选取合适的单元类型划分网格，网格划分方式应同时满足计算精度和计算效率；

(2.5)将编译完成的umat子程序通过abaqus主程序接口嵌入到abaqus有限元模型中，提交模型并使其运算完成。

进一步地，设q为偏应力，qf为土体强度，则

其中，c为土体粘聚力，为土体内摩擦角，σ3为小主应力。

进一步地，当q＜qf时，土体处于弹性阶段，竖向应变ε1和偏应力q之间满足双曲线关系：

其中，e50为加载模量，即：

其中，σ1为大主应力；σ^ref为相关应力；为相关应力σ^ref时的加载模量；m为幂指数；qa为极限偏应力，即：qa＝(σ1-σ3)ult。

进一步地，当q≥qf时，土体进入塑性阶段，产生塑性变形。

本发明的有益效果为：本发明克服了abaqus自带土体相关模型的缺陷，如dc模型不能反映土体的塑性应变和不同的应力路径，mc模型对材料抗拉强度的过高估计等问题，通过fortran语言编译umat子程序将土体hs模型嵌入到abaqus软件中，实现了土体的数值模拟计算。

【附图说明】

此处所说明的附图是用来提供对本发明的进一步理解，构成本申请的一部分，但并不构成对本发明的不当限定，在附图中：

图1是hs模型的屈服面。

图2是本发明是求取第n+1增量步参数的流程图。

【具体实施方式】

下面将结合附图以及具体实施例来详细说明本发明，其中的示意性实施例以及说明仅用来解释本发明，但并不作为对本发明的限定。

岩石工程中的hs模型考虑了材料的硬化，设q为偏应力，qf为土体强度，由莫尔-库仑强度理论有：

其中，c为土体粘聚力，为土体内摩擦角，σ3为小主应力。

当q＜qf时，土体处于弹性阶段，竖向应变ε1和偏应力q之间满足双曲线关系：

其中，e50为加载模量，即：

上面两式中，σ1为大主应力；σ^ref为相关应力，一般取100kpa；为相关应力σ^ref时的加载模量；m为幂指数；qa为极限偏应力，即：qa＝(σ1-σ3)ult。

当q≥qf时，土体进入塑性阶段，产生塑性变形，随着硬化参数的变化，hs模型屈服面也在不断的变化。hs模型为双屈服函数，包括剪切屈服和帽盖屈服函数，其屈服面如图1所示。

本发明基于上述hs模型，提供了基于abaqus的岩土工程数值模拟方法，具体包括以下步骤。

步骤一：基于euler积分算法，利用fortran语言编制hs本构模型的子程序，实现abaqus软件的二次开发。

所述euler积分算法是umat子程序编写的关键，它的基本思路是已知第n增量步的所有变量值，给定时间步长增量和总应变增量，通过数学算法求得第n+1增量步的满足本构方程的应力解。应力积分算法非常重要，可以影响计算的精度、效率和稳定性。

具体的，参见附图2，假设第n增量步得到的应力为σn，塑性应变为总应变为εn，则求取第n+1增量步的相关参数的具体实现步骤如下：

(1)通过输入的应力状态参数计算应力状态不变量、屈服函数对应力的偏导数和弹性矩阵d^e；

(2)进行弹性试算，假设应变增量dε都是弹性应变，那么试探应力σ^trial可表示为：

σ^trial＝σn+d^edε；

(3)计算试探应力对应的应力状态不变量，然后根据应力状态不变量计算塑性乘子λ，根据塑性乘子的正负号判定加卸载状态。若塑性乘子大于0，则处于加载状态，若塑性乘子小于0，则处于卸载状态。

(4)根据步骤3中的判定结果进行选择，若处于卸载状态，则不进行塑性修正，且将雅可比矩阵d设为弹性矩阵d^e，然后跳转到步骤7；若处于加载状态，则执行步骤5进行塑性修正。

(5)塑性修正。加载状态时产生塑性变形，需要进行塑性修正，修正后的第n+1增量步的应力σn+1、塑性应变和总应变εn+1的更新式如下：

(6)更新并返回雅克比矩阵。获得第n+1步的应力值后，要进行雅克比矩阵的更新，以供abaqus主程序进行下一个增量步的计算。雅可比矩阵的更新公式如下：

其中，d^ep是更新后的雅克比矩阵，σii为主应力，σm为平均有效应力，ρ为新增的状态变量，与超固结比相关，a为材料参数，决定ρ的发展速度。λ为压缩参数，κ为回弹参数，e0为参照应力对应的孔隙比。

(7)将应变分量、等效塑性应变等更新之后放在状态变量statev(1：nstatev)中存储，以便在后处理中输出查看。

步骤二：将考虑材料硬化效应的hs模型编入abaqus软件的umat子程序中，实现土体的数值模拟。

具体的，步骤二包括以下实现步骤：

(1)提取岩土工程结构的几何参数，在abaqus软件中建立几何模型；

(2)在abaqus软件的材料输入界面，使用关键词“usermaterial”，表示定义用户材料属性，根据试验结果，将材料参数输入到数值模型中；

(3)在abaqus的load模块中，对模型的边界条件进行设置，并施加初始荷载；

(4)在abaqus中选取合适的单元类型划分网格，网格划分方式应同时满足计算精度和计算效率；

(5)将编译完成的umat子程序通过abaqus主程序接口嵌入到abaqus有限元模型中，提交模型并使其运算完成。

基于上述方法，将土体的hs模型通过编译用户材料子程序umat的方式，嵌入到abaqus软件中，最终实现了岩土工程中土体本构模拟的数值计算。本发明克服了abaqus自带土体相关模型的缺陷，在一定程度上弥补了abaqus在岩土工程计算分析方面的不足。通过三轴试验模拟，验证了所编制的umat子程序思路正确，能合理反映土体应力应变的非线性关系，对指导工程实际问题具有重要的科学价值。

以上所述仅是本发明的较佳实施方式，故凡依本发明专利申请范围所述的构造、特征及原理所做的等效变化或修饰，均包括于本发明专利申请范围内。