买专利卖专利找龙图腾，真高效！ 查专利查商标用IPTOP,全免费！专利年费监控用IP管家,真方便！

申请/专利权人：合肥哈工热气球数字科技有限公司

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

主权项：1.一种基于ABAQUS的岩土工程数值模拟方法，其特征在于，包括以下步骤：步骤一：基于Euler积分算法，利用Fortran语言编制HS本构模型的子程序，实现ABAQUS软件的二次开发；步骤二：将考虑材料硬化效应的HS模型编入ABAQUS软件的UMAT子程序中，实现土体的数值模拟；所述HS模型包括剪切屈服和帽盖屈服函数；其中，步骤一中的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将应变分量、等效塑性应变等更新之后放在状态变量STATEV1：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-σ3ult。

全文数据：一种基于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将应变分量、等效塑性应变等更新之后放在状态变量STATEV1：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-σ3ult。进一步地，当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-σ3ult。当q≥qf时，土体进入塑性阶段，产生塑性变形，随着硬化参数的变化，HS模型屈服面也在不断的变化。HS模型为双屈服函数，包括剪切屈服和帽盖屈服函数，其屈服面如图1所示。本发明基于上述HS模型，提供了基于ABAQUS的岩土工程数值模拟方法，具体包括以下步骤。步骤一：基于Euler积分算法，利用Fortran语言编制HS本构模型的子程序，实现ABAQUS软件的二次开发。所述Euler积分算法是UMAT子程序编写的关键，它的基本思路是已知第n增量步的所有变量值，给定时间步长增量和总应变增量，通过数学算法求得第n+1增量步的满足本构方程的应力解。应力积分算法非常重要，可以影响计算的精度、效率和稳定性。具体的，参见附图2，假设第n增量步得到的应力为σn，塑性应变为总应变为εn，则求取第n+1增量步的相关参数的具体实现步骤如下：1通过输入的应力状态参数计算应力状态不变量、屈服函数对应力的偏导数和弹性矩阵De；2进行弹性试算，假设应变增量dε都是弹性应变，那么试探应力σtrial可表示为：σtrial＝σn+Dedε；3计算试探应力对应的应力状态不变量，然后根据应力状态不变量计算塑性乘子Λ，根据塑性乘子的正负号判定加卸载状态。若塑性乘子大于0，则处于加载状态，若塑性乘子小于0，则处于卸载状态。4根据步骤3中的判定结果进行选择，若处于卸载状态，则不进行塑性修正，且将雅可比矩阵D设为弹性矩阵De，然后跳转到步骤7；若处于加载状态，则执行步骤5进行塑性修正。5塑性修正。加载状态时产生塑性变形，需要进行塑性修正，修正后的第n+1增量步的应力σn+1、塑性应变和总应变εn+1的更新式如下：6更新并返回雅克比矩阵。获得第n+1步的应力值后，要进行雅克比矩阵的更新，以供ABAQUS主程序进行下一个增量步的计算。雅可比矩阵的更新公式如下：其中，Dep是更新后的雅克比矩阵，σii为主应力，σm为平均有效应力，ρ为新增的状态变量，与超固结比相关，a为材料参数，决定ρ的发展速度。λ为压缩参数，κ为回弹参数，e0为参照应力对应的孔隙比。7将应变分量、等效塑性应变等更新之后放在状态变量STATEV1：NSTATEV中存储，以便在后处理中输出查看。步骤二：将考虑材料硬化效应的HS模型编入ABAQUS软件的UMAT子程序中，实现土体的数值模拟。具体的，步骤二包括以下实现步骤：1提取岩土工程结构的几何参数，在ABAQUS软件中建立几何模型；2在ABAQUS软件的材料输入界面，使用关键词“USERMATERIAL”，表示定义用户材料属性，根据试验结果，将材料参数输入到数值模型中；3在ABAQUS的LOAD模块中，对模型的边界条件进行设置，并施加初始荷载；4在ABAQUS中选取合适的单元类型划分网格，网格划分方式应同时满足计算精度和计算效率；5将编译完成的UMAT子程序通过ABAQUS主程序接口嵌入到ABAQUS有限元模型中，提交模型并使其运算完成。基于上述方法，将土体的HS模型通过编译用户材料子程序UMAT的方式，嵌入到ABAQUS软件中，最终实现了岩土工程中土体本构模拟的数值计算。本发明克服了ABAQUS自带土体相关模型的缺陷，在一定程度上弥补了ABAQUS在岩土工程计算分析方面的不足。通过三轴试验模拟，验证了所编制的UMAT子程序思路正确，能合理反映土体应力应变的非线性关系，对指导工程实际问题具有重要的科学价值。以上所述仅是本发明的较佳实施方式，故凡依本发明专利申请范围所述的构造、特征及原理所做的等效变化或修饰，均包括于本发明专利申请范围内。

权利要求：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是更新后的雅克比矩阵，σii为主应力，σm为平均有效应力，ρ为新增的状态变量，与超固结比相关，a为材料参数，λ为压缩参数，κ为回弹参数，e0为参照应力对应的孔隙比；1.7将应变分量、等效塑性应变等更新之后放在状态变量STATEV1：NSTATEV中存储，以便在后处理中输出查看。2.根据权利要求1所述的方法，其特征在于，所述步骤二包括以下实现步骤：2.1提取岩土工程结构的几何参数，在ABAQUS软件中建立几何模型；2.2在ABAQUS软件的材料输入界面，使用关键词“USERMATERIAL”，表示定义用户材料属性，根据试验结果，将材料参数输入到数值模型中；2.3在ABAQUS的LOAD模块中，对模型的边界条件进行设置，并施加初始荷载；2.4在ABAQUS中选取合适的单元类型划分网格，网格划分方式应同时满足计算精度和计算效率；2.5将编译完成的UMAT子程序通过ABAQUS主程序接口嵌入到ABAQUS有限元模型中，提交模型并使其运算完成。3.根据权利要求1-2任意一项所述的方法，其特征在于，设q为偏应力，qf为土体强度，则其中，c为土体粘聚力，为土体内摩擦角，σ3为小主应力。4.根据权利要求1-3任意一项所述的方法，其特征在于，当q＜qf时，土体处于弹性阶段，竖向应变ε1和偏应力q之间满足双曲线关系：其中，E50为加载模量，即：其中，σ1为大主应力；σref为相关应力；为相关应力σref时的加载模量；m为幂指数；qa为极限偏应力，即：qa＝σ1-σ3ult。5.根据权利要求1-4任意一项所述的方法，其特征在于，当q≥qf时，土体进入塑性阶段，产生塑性变形。

百度查询： 合肥哈工热气球数字科技有限公司 一种基于ABAQUS的岩土工程数值模拟方法