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

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

A numerical simulation method of geotechnical engineering based on ABAQUS

The present invention relates to a numerical simulation method of geotechnical engineering based on ABAQUS, which is based on Euler integral algorithm, uses Fortran language to compile a subprogram of HS constitutive model, realizes the secondary development of ABAQUS software, incorporates HS model considering material hardening effect into UMAT subroutine of ABAQUS software, realizes the numerical simulation of soil, overcomes the deficiency of ABAQUS self-contained soil related model. Subsidence, to a certain extent, makes up for the deficiency of ABAQUS in geotechnical engineering calculation and analysis.

【技术实现步骤摘要】
一种基于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)通过输入的应力状态参数计算应力状态不变量、屈服函数对应力的偏导数和弹性矩阵De；(1.2)进行弹性试算，假设应变增量dε都是弹性应变，那么试探应力σtrial为：σtrial＝σn+Dedε；(1.3)计算试探应力对应的应力状态不变量，然后根据应力状态不变量计算塑性乘子Λ，根据塑性乘子的正负号判定加卸载状态，若塑性乘子大于0，则处于加载状态，若塑性乘子小于0，则处于卸载状态；(1.4)根据步骤1.3中的判定结果进行选择，若处于卸载状态，则不进行塑性修正，且将雅可比矩阵D设为弹性矩阵De，然后跳转到步骤1.7；若处于加载状态，则执行步骤1.5进行塑性修正；(1.5)塑性修正，修正后的第n+1增量步的应力σn+1、塑性应变和总应变εn+1的更新式如下：(1.6)更新并返回雅克比矩阵。获得第n+1步的应力值后，要进行雅克比矩阵的更新，以供ABAQUS主程序进行下一个增量步的计算。雅可比矩阵的更新公式如下：其中，Dep是更新后的雅克比矩阵，σ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)通过输入的应力状态参数计算本文档来自技高网...

【技术保护点】
1.一种基于ABAQUS的岩土工程数值模拟方法，其特征在于，包括以下步骤：步骤一：基于Euler积分算法，利用Fortran语言编制HS本构模型的子程序，实现ABAQUS软件的二次开发；步骤二：将考虑材料硬化效应的HS模型编入ABAQUS软件的UMAT子程序中，实现土体的数值模拟；其中，步骤一中的Euler积分算法根据已知第n增量步的所有变量值，给定时间步长增量和总应变增量，求得第n+1增量步的满足本构方程的应力解，假设第n增量步得到的应力为σn，塑性应变为

【技术特征摘要】
1.一种基于ABAQUS的岩土工程数值模拟方法，其特征在于，包括以下步骤：步骤一：基于Euler积分算法，利用Fortran语言编制HS本构模型的子程序，实现ABAQUS软件的二次开发；步骤二：将考虑材料硬化效应的HS模型编入ABAQUS软件的UMAT子程序中，实现土体的数值模拟；其中，步骤一中的Euler积分算法根据已知第n增量步的所有变量值，给定时间步长增量和总应变增量，求得第n+1增量步的满足本构方程的应力解，假设第n增量步得到的应力为σn，塑性应变为总应变为εn，则求取第n+1增量步的相关参数的具体实现步骤如下：(1.1)通过输入的应力状态参数计算应力状态不变量、屈服函数对应力的偏导数和弹性矩阵De；(1.2)进行弹性试算，假设应变增量dε都是弹性应变，那么试探应力σtrial为：σtrial＝σn+Dedε；(1.3)计算试探应力对应的应力状态不变量，然后根据应力状态不变量计算塑性乘子Λ，根据塑性乘子的正负号判定加卸载状态，若塑性乘子大于0，则处于加载状态，若塑性乘子小于0，则处于卸载状态；(1.4)根据步骤1.3中的判定结果进行选择，若处于卸载状态，则不进行塑性修正，且将雅可比矩阵D设为弹性矩阵De，然后跳转到步骤1.7；若处于加载状态，则执行步骤1.5进行塑性修正；(1.5)塑性修正，修正后的第n+1增量步的应力σn+1、塑性应变和总应变εn+1的更新式如下：(1.6)更新并返回雅克比矩阵：获得第n+1步的应力值后，要进行雅克比矩阵的更新，以供ABAQUS主程序进行下一个增量步的计算，雅可比矩阵的更新公式如下：其中，Dep是更新后的雅克比矩阵，σ...

【专利技术属性】
技术研发人员：张洁洁，丁亮，金马，董事，李传国，甘建安，夏科睿，
申请(专利权)人：哈工大机器人合肥国际创新研究院，
类型：发明
国别省市：安徽,34