CN109446731A - 一种基于abaqus的岩土工程数值模拟方法 - Google Patents

Abstract

本发明涉及一种基于ABAQUS的岩土工程数值模拟方法，该方法基于Euler积分算法，利用Fortran语言编制HS本构模型的子程序，实现ABAQUS软件的二次开发；将考虑材料硬化效应的HS模型编入ABAQUS软件的UMAT子程序中，实现土体的数值模拟；克服了ABAQUS自带土体相关模型的缺陷，在一定程度上弥补了ABAQUS在岩土工程计算分析方面的不足。

Description

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

【技术领域】

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

【背景技术】

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

目前常用于岩土工程的土体本构模型有：邓肯-张(DC)模型、莫尔-库仑(MC)模型、修正剑桥(MCC)模型、硬化土(hardening soil)模型(简称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为材料参数，λ为压缩参数，κ为回弹参数，e₀为参照应力对应的孔隙比；

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

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

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

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

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

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

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

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

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

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

其中，E₅₀为加载模量，即：

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

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

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

【附图说明】

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

图1是HS模型的屈服面。

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

【具体实施方式】

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

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

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

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

其中，E₅₀为加载模量，即：

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

当q≥q_f时，土体进入塑性阶段，产生塑性变形，随着硬化参数的变化，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为材料参数，决定ρ的发展速度。λ为压缩参数，κ为回弹参数，e₀为参照应力对应的孔隙比。

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

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

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

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

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

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

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

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

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

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

Claims (5)

1.一种基于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为材料参数，λ为压缩参数，κ为回弹参数，e₀为参照应力对应的孔隙比；

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

2.根据权利要求1所述的方法，其特征在于，所述步骤二包括以下实现步骤：

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

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

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

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

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

3.根据权利要求1-2任意一项所述的方法，其特征在于，设q为偏应力，q_f为土体强度，则

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

4.根据权利要求1-3任意一项所述的方法，其特征在于，当q＜q_f时，土体处于弹性阶段，竖向应变ε₁和偏应力q之间满足双曲线关系：

其中，E₅₀为加载模量，即：

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

5.根据权利要求1-4任意一项所述的方法，其特征在于，当q≥q_f时，土体进入塑性阶段，产生塑性变形。