多尺度物理本构模型嵌入ABAQUS的方法与流程

本发明涉及一种本构模型二次开发方法，特别是涉及一种多尺度物理本构模型嵌入abaqus的方法。

背景技术：

文献“基于vumat的固体推进剂材料本构模型二次开发(《机械强度》2013年第35期391—394页”公开了一种本构模型二次开发方法，该方法只针对唯相模型，本构模型不涉及微观组织变化，内变量之间无互相耦合，可直接进行求导处理后通过fortran编写子程序即可实现嵌入。现阶段，材料微观组织演变的研究已成为热门，多尺度物理本构模型的内变量较多，相互之间耦合，由于本构方程内变量之间的耦合关系复杂，传统的插值方法求解的运算量较大，成为嵌入的难点，因此，非常有必要研究建立一种方法，并解决多尺度物理本构模型的嵌入问题。

技术实现要素：

为了克服现有本构模型二次开发方法复杂的不足，本发明提供一种多尺度物理本构模型嵌入abaqus的方法。该方法对于多尺度物理本构模型，计算子程序所需的屈服应力yield(i)以及多尺度物理本构模型中的晶粒尺寸d、动态再结晶分数x和位错密度ρ；对已知应变率下的屈服应力对塑性应变进行拟合，针对选取的本构方程，选取三种应变率下的应力应变离散数据进行多项式拟合；在不同应变率下的a0、a1、a2、a3、a4、a5，对进行线性拟合；通过两次拟合，得到拟合方程。本发明通过建立屈服应力与应变与应变率的联系，直接求得vuhard接口所需的各个变量。使后续的计算更准确，减少了计算的复杂性。由于采用对目标应变与应变率下的屈服应力进行拟合的方法，解决了多尺度本构模型内变量耦合过多，无法直接求得接口变量的技术问题。

本发明解决其技术问题所采用的技术方案是：一种多尺度物理本构模型嵌入abaqus的方法，其特点是包括以下步骤：

步骤一、对于多尺度物理本构模型，子程序所需的屈服应力yield(i)以及多尺度物理本构模型中的晶粒尺寸d、动态再结晶分数x、位错密度ρ通过式(1)计算：

式中，表示晶粒尺寸变化率，为位错密度变化率，为再结晶率，为等效塑性应变增量。

步骤二、对已知应变率下的屈服应力对塑性应变进行拟合，针对选取的本构方程，选取三种应变率下的应力应变离散数据进行多项式拟合，拟合方程如下：

式中，yield^a(i)，yield^b(i),yield^c(i)为三种应变率下的屈服应力。

步骤三、在不同应变率下的a0、a1、a2、a3、a4、a5，对进行线性拟合，拟合方程为：

式中，为等效塑性应变率。

步骤四、通过两次拟合，得到的拟合方程如下：

本发明的有益效果是：该方法对于多尺度物理本构模型，计算子程序所需的屈服应力yield(i)以及多尺度物理本构模型中的晶粒尺寸d、动态再结晶分数x和位错密度ρ；对已知应变率下的屈服应力对塑性应变进行拟合，针对选取的本构方程，选取三种应变率下的应力应变离散数据进行多项式拟合；在不同应变率下的a0、a1、a2、a3、a4、a5，对进行线性拟合；通过两次拟合，得到拟合方程。本发明通过建立屈服应力与应变与应变率的联系，直接求得vuhard接口所需的各个变量。使后续的计算更准确，减少了计算的复杂性。由于采用对目标应变与应变率下的屈服应力进行拟合的方法，解决了多尺度本构模型内变量耦合过多，无法直接求得接口变量的技术问题。

下面结合具体实施方式对本发明作详细说明。

具体实施方式

本发明多尺度物理本构模型嵌入abaqus的方法具体步骤如下：

1、子程序嵌入模式及其所需接口变量的求解

子程序中只需要定义4个接口变量就可以实现与abaqus的仿真交互过程，包括：材料变形的屈服应力yield(i)以及超塑性屈服应力对应变eqps(i)、应变率eqpsrate(i)的导数dyielddeqps(i,1)、dyielddeqps(i,2)。子程序接口如下：

在abaqus嵌入的过程中，需要在usercoding部分求得接口所需的屈服应力yield(i)，屈服应力对应变的导数dyielddeqps(i,1)、屈服应力对应变率的导数dyielddeqps(i,2)。

选取“神秘”材料采用多尺度物理本构模型进行嵌入，选取某多尺度物理本构模型，其通用形式如下：

其中yield(i)为屈服应力，为等效塑性应变，为应变率，d为晶粒尺寸，x为动态再结晶分数、ρ为位错密度，为晶粒尺寸变化率，为动态再结晶分数变化率、为位错密度变化率。可以看出，本构模型内变量之间互相耦合，形式复杂，在本构模型的运算过程中，无法直接通过求导求得子程序接口所需的屈服应力yield(i)以及对应变的导数dyielddeqps(i,1)、屈服应力对应变率的导数dyielddeqps(i,2)。本文提出一种于仿真过程中计算中间变量的新方法如下：

式中：表示晶粒尺寸变化率，为位错密度变化率，为再结晶率，为等效塑性应变增量，通过子程序接口提供的状态变量state(i,j)记录上一增量步的等效塑性应变与当前增量步的等效塑性应变做差求得等效塑性应变增量。

通过实验获得“神秘”材料在10^-2、10^-3、10^-4三种应变速率下的材料性能曲线如下表：

表一10^-2、10^-3、10^-4三种应变速率下的材料性能

直接求得偏导的在分析过程中工作量较大，影响运行速度，为了提高计算效率，求取，选取10^-2、10^-3、10^-4三种应变率下的应力应变离散数据进行拟合，拟合方程如下：

表二多项式方程参数表

在不同应变率下的a0、a1、a2、a3、a4、a5，对进行线性拟合，拟合方程为：

线性方程参数表三所示：

表三线性方程拟合参数表

此时即可通过直接对多项式方程求导即可得出屈服应力yield(i)对应变的导数dyielddeqps(i,1)如下：

屈服应力对应变率的导数dyielddeqps(i,2)如下：

此时即可通过frotran编写子程序求得子程序所需的变量，实现子程序的嵌入。