一种基于晶体滑移机制的各向异性线弹性本构的建立方法

摘要

本发明公开了本发明涉及一种基于晶体滑移机制的各向异性线弹性本构的建立方法，属于材料力学性能有限元数值计算领域。该方法先建立有限元网格模型，并确定边界条件和加载条件；采用投影运算方法，求解有限元网格模型的坐标系与晶粒坐标系之间的转化矩阵；通过张量运算，得到有限元网格模型坐标系下的单晶体弹性刚度矩阵和施密特因子；根据设置的边界条件和加载条件进行加载；根据Hill‑Hutchinson多晶体弹塑性变形理论，求解晶粒在塑性阶段的瞬时刚度矩阵，将瞬时刚度矩阵替换单晶体弹性刚度矩阵，重复加载直至设定的最大加载步数，得到瞬时刚度矩阵即为晶粒的各向异性线弹性本构。

说明书

技术领域

本发明涉及一种基于晶体滑移机制的各向异性线弹性本构的建立方法，属于材料力学性能有限元数值计算领域。

背景技术

在利用基于真实微观组织的建模方法模拟多晶材料宏、微观力学性能方面，材料本构模型的选用和参数获取对模拟结果具有决定性的影响，同时也是一个亟待解决的薄弱环节。国内外在研究材料本构的问题上，通常是采用各向同性本构模型，即用唯象的描述方法考虑介观层次的相组织模型。获得相本构的方法主要是认为材料微观变形特征与宏观近似相同，进而将具有均一组织结构材料的整体性能赋予该相；而对于多相材料中的各个组成相，则是采用具有近似成分的单相材料的宏观性能进行替代。大型商用有限元软件LS-DYNA是一种通用显式动力分析软件，在工程应用领域具有较高的计算可靠性。其具有显/隐式求解功能，求解稳定；同时具有在复杂边界加载条件下的非线性动力分析和静力分析功能。后处理功能更是可以实现内部多种求解参数的可视化，但目前其应用主要是基于唯象学，鲜有考虑材料变形的内在机制。

然而为了考虑材料变形过程中存在的内部物理机制，各向同性的相组织模型已无法满足要求，因此所建立的模型已逐渐发展为考虑滑移机制的晶粒组织模型，即由唯象描述方发展为晶体塑性描述方法。其中用于研究材料加载变形内部机制的弹塑性自洽模拟方法已得到了广泛的应用，该方法可以用于模拟材料整体的宏观应力-应变响应；研究材料弹塑性转变过程及过程中表征参量的拟合；拟合晶格演变过程；变形过程中的织构演化等。自洽模拟方法的主要理论基础和假设有：(1)采用Hill-Hutchinson多晶体弹塑性变形机制、Vocé硬化定律及硬化系数矩阵；(2)晶粒的各向异性的实现方法是通过将晶粒设置成具有各向异性弹性常数矩阵的椭圆球体，并且给每个晶粒赋予对应其晶体结构类型的滑移系及特征参数；(3)晶粒间的相互作用是通过弹塑性Eshelby张量建立起每一个晶体与假定的均质基体之间的相互作用函数；(4)采用增量加载模式，每一加载步为小变形加载，即认为应力-应变关系为瞬时线弹性。

该方法不仅考虑了多晶材料中晶体的实际取向，还可以考虑多相材料不同的晶体结构类型和塑性变形滑移机制、硬化效应和取向变化等。但其存在一些缺陷，如：需要配合原位拉伸同步辐射X射线衍射、中子衍射等实验来获得与之相映证的数据，而这些实验不仅周期长而且费用昂贵；其相互作用是通过“限制张量”在算法内部实现的，无法进行可视化，并且无法考虑真实微观组织中相邻晶粒间的相互作用关系；可设置的加载边界条件比较简单，通常是由应力或应变控制的以一定速率的单轴拉伸或压缩等。

发明内容

本发明的目的是提供一种基于晶体滑移机制的各向异性线弹性本构的建立方法，所述方法延续了自洽模型中采用的描述晶体塑性变形的理论机制，同时与有限元方法相结合，不仅考虑了材料的微观组织形貌，同时用有限元网格的单元和节点化来实现晶粒间的相互作用，弥补了自洽模型中的不足，并可充分发挥有限元模拟计算的优势，实现复杂的边界加载条件和多样的可视化后处理操作。

本发明的目的有以下技术方案实现：

一种基于晶体滑移机制的各向异性线弹性本构的建立方法，所述方法包括以下步骤：

步骤一：

基于待求各向异性线弹性本构材料的晶粒组织图像，采用有限元网格剖分方法，建立所述晶粒组织的有限元网格模型，并确定所述有限元网格模型的边界条件和加载条件；

其中，所述加载条件包括载荷和最大加载步数，其中，初始加载步数为0；

步骤二：

采用投影运算方法，根据晶粒组织的本征参数、初始应力应变状态分量和晶粒取向的初始欧拉角求解晶粒取向欧拉角增量并将所述欧拉角增量与初始欧拉角相加，得到由于应变引起晶粒取向变化后的欧拉角根据变化后的欧拉角获得有限元网格模型的坐标系与晶粒坐标系之间的转化矩阵；

其中，所述本征参数包括晶体结构类型参数、单晶体弹性刚度矩阵、晶粒内的滑移系参数；

所述滑移系参数包括滑移系的名称、等效滑移系的数量和滑移系的Vocé硬化指数；所述Vocé硬化指数包括滑移系开动的临界分切应力值τ₀、滑移系从开动到累积切应变不再增加时的应力增量τ₁、滑移系的初始硬化率θ₀以及滑移系的饱和硬化率θ₁；

所述晶粒取向欧拉角采用Bunge方法定义，即将晶粒坐标系绕其z轴，在xy平面内转过角；再绕其x轴，在yz平面内转过φ角；最后再一次绕z轴，在xy平面内转过角，最终与有限元网格模型坐标系完全重合，由此得到每个晶粒的三个空间取向欧拉角晶粒取向的初始欧拉角利用电子背散射衍射(EBSD)实验测量得到；

所述应力应变状态分量分别按照σ_xx，σ_yy，σ_zz，σ_yz，σ_xz，σ_xy和ε_xx，ε_yy，ε_zz，ε_yz，ε_xz，ε_xy的顺序给出，初始值均设为0；

其中，所述σ_xx，σ_yy和σ_zz分别表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内沿X轴、Y轴和Z轴方向的正应力；

所述ε_xx，ε_yy和ε_zz分别表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内沿X轴、Y轴和Z轴方向的正应变；

所述σ_yz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在Y轴法向平面内平行于X轴的切应力；

所述σ_xz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴法向平面内平行于Z轴的切应力；

所述σ_xy表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴法向平面内平行于Y轴的切应力；

所述ε_yz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在Y轴和Z轴两方向上微小线段夹角的切应变；

所述ε_xz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴和Z轴两方向上微小线段夹角的切应变；

所述ε_xy表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴和Y轴两方向上微小线段夹角的切应变；

步骤三：

基于步骤二所述转化矩阵，通过张量运算，得到有限元网格模型坐标系下的单晶体弹性刚度矩阵和施密特因子α，其中，所述α＝cosλcosΦ，Φ为滑移面与外力F中心轴的夹角，所述λ为滑移方向与外力F的夹角；

步骤四：

判断所述欧拉角增量是否为(0，0，0)，若进入步骤五；若进入步骤六；

步骤五：

根据步骤一所述有限元网格模型的边界条件和加载载荷对当前有限元网格模型进行加载，然后针对当前有限元网格模型坐标系下的单晶体弹性刚度矩阵，采用有限元计算方法计算得到加载后晶粒的应力应变状态分量；将所述加载后晶粒的应力应变状态分量替换步骤二中的初始应力应变状态分量，重复步骤二至步骤四；其中，每次加载后加载步数加一；步骤六：

判断晶粒内的潜在要开动和已经开动滑移系的总个数n是否大于零；若n＝0，重复步骤五；若n＞0，进入步骤七；

潜在要开动和已经开动滑移系的判断方法：

若则该滑移系被标记为潜在要开动，若则该滑移系被标记为已经开动；其中，σ_c表示任意滑移系所在晶粒当前的应力值，表示第i个滑移系的临界分切应力值；所述αⁱ表示第i个滑移系的施密特因子；

步骤七：

根据Hill-Hutchinson多晶体弹塑性变形理论，先求解n个滑移系应变线性方程组的基向量矩阵；再根据所述基向量矩阵，求解分配给每个滑移系的切变量，进而得到晶粒在塑性阶段的瞬时刚度矩阵；

步骤八：

判断当前加载步数是否等于设定的最大加载步数，如果等于最大加载步数，则步骤七得到的瞬时刚度矩阵即为晶粒的各向异性线弹性本构；如果小于最大加载步数，则将步骤七的瞬时刚度矩阵替换步骤二中的单晶体弹性刚度矩阵，重复步骤二至步骤七，直至加载步数等于步骤一中设定的最大加载步数。

有益效果

(1)本发明所述方法延续了自洽模型中采用的描述晶体塑性变形的理论机制，同时与有限元方法相结合；采用基于晶粒组织的网格模型，不仅考虑了材料的微观组织形貌，同时用有限元网格的单元和节点化来实现晶粒间的相互作用，弥补了自洽模型中的不足。

(2)本发明所述方法利用有限元求解器中的重启动技术，基于晶粒组织建立的有限元网格模型以一定边界条件和加载条件进行增量式分步重启动加载，每一加载步都继承了上一加载步结束时各单元和节点的应力应变状态，从而实现分步增量加载；通过采用有限元求解软件中提供的各向异性线弹性本构模型来体现多晶体材料中各个晶粒的取向性，同时结合晶体塑性变形理论，可以根据每一步加载后晶粒的应力应变状态，从而达到在塑性变形机制层面对本构参数(即瞬时模量)实现实时的更新改进；充分发挥有限元模拟计算的优势，实现复杂的边界加载条件和多样的可视化后处理操作。

(3)本发明所述方法可以更进一步地考虑材料的真实微观组织形貌，从相组织模型发展为晶粒组织模型，从而可以更深层次的考虑材料变形过程中存在的内部物理机制；采用本发明提供的方法可以计算得到晶粒组织有限元网格模型中每个晶粒的随着加载过程不断变化的各向异性线弹性本构，对再现材料组织演变过程、分析材料局部受力状态，以及预测材料失效形式等提供了更为准确的研究方法。

附图说明

图1为本发明所述方法的流程示意图；

图2为所述晶粒组织的有限元网格模型的坐标系(左)和晶粒的坐标系(右)示意图；

图3为实施例1中所述Vocé硬化指数含义示意图；

图4为实施例1中所述施密特因子定义示意图。

具体实施方式

下面结合附图和具体实施例来详述本发明，但不限于此。

实施例1

如图1所示，一种基于晶体滑移机制的各向异性线弹性本构的建立方法，所述方法包括以下步骤：

步骤一：

基于待求各向异性线弹性本构材料的晶粒组织图像，基于有限元软件LS-DYNA，采用有限元网格剖分方法，建立所述晶粒组织的有限元网格模型，并确定所述有限元网格模型的边界条件和加载条件；

所述加载条件包括载荷和加载步数，其中，初始加载步数为0；

所述晶粒组织的有限元网格模型的坐标系和晶粒的坐标系示意图如图2所示；

步骤二：

采用Fortran程序读取晶粒组织的本征参数、初始应力应变状态分量和晶粒取向的初始欧拉角并根据应变增量的投影关系，所述欧拉角增量由晶粒的应变增量按以下公式计算得出：

并将所述欧拉角增量与初始欧拉角相加，得到由于应变引起晶粒取向变化后的欧拉角根据变化后的欧拉角获得有限元网格模型的坐标系OXYZ与晶粒坐标系oxyz之间的转化矩阵R；

其中，所述本征参数包括晶体结构类型参数、单晶体弹性刚度矩阵、晶粒内的滑移系参数；

所述晶体结构类型参数包括晶体结构类型和晶胞的轴比c/a；所述晶体结构类型用ihcp表示，ihcp＝1表示为密排六方晶体结构(hexagonal close-packed,HCP)，其晶向指数和晶面指数均用四坐标表示；ihcp＝0表示为立方晶体结构(CUBIC)，其晶向指数和晶面指数用三坐标；所述晶胞的轴比c/a用rca表示，以TC6钛合金中的α和β两相为例，对于其中的α相，rca＝1.597；对于其中的β相，rca＝1.0。

所述单晶体弹性刚度矩阵C_ij，为6*6的矩阵，单位为GPa，对于密排六方晶体结构具有对称性C₁₁＝C₂₂，C₁₃＝C₂₃，C₄₄＝C₅₅＝C₆₆；而对于立方晶体结构则为C₁₁＝C₂₂＝C₃₃，C₁₂＝C₁₃＝C₂₃，C₄₄＝C₅₅＝C₆₆；其中，TC6钛合金中的α相，其晶体结构类型参数为：hexagonal>

对于TC6钛合金中的β相，其晶体结构类型参数为：CUBIC，晶体结构类型ihcp＝1，晶胞的轴比rca＝1.0，单晶体弹性刚度矩阵为：

所述滑移系参数包括滑移系的名称、等效滑移系的数量和滑移系的Vocé硬化指数；晶粒内所有需要考虑的滑移系，包括文件中滑移系的总数量、需要读取的滑移系的总数量(Num_slip)；

滑移系的名称：对于密排六方晶体结构来说常见的滑移系有四个分别是：{0002}<11-20>基面滑移系、{10-10}<11-20>柱面滑移系、{10-11}<11-20>一阶锥面滑移系和{10-11}<11-23>二阶锥面滑移系，对于体心立方晶体结构来说常见的滑移系有三个：{1-10}<111>，{11-2}<111>和{12-3}<111>；

等效滑移系的数量用(m_slip)表示；

滑移系的Vocé硬化指数：其包括滑移系开动的临界分切应力值τ₀、滑移系从开动到累积切应变不再增加时的应力增量τ₁、滑移系的初始硬化率θ₀、滑移系的饱和硬化率θ₁、需要读取的滑移系间的相互作用系数h^ij(i＝j表示自硬化系数，i≠j表示互硬化)；具体含义由应力应变曲线表示，如图3所示；在Fortran程序中将所有需要读取的滑移系的所有等效滑移系顺次标号；以滑移系{0001}<11-20>为例，有6个等效滑移系(包含负向滑移)，如表1所示；

表1晶体内滑移系参数设置

所述晶粒取向欧拉角采用Bunge方法定义，即将晶粒坐标系绕z轴，在xy平面内转过角；再绕x轴，在yz平面内转过φ角；最后再一次绕z轴，在xy平面内转过角，最终与有限元网格模型坐标系完全重合，由此得到每个晶粒的三个空间取向欧拉角晶粒取向的初始欧拉角利用EBSD实验测量得到；

所述应力应变状态分量分别按照σ_xx，σ_yy，σ_zz，σ_yz，σ_xz，σ_xy和ε_xx，ε_yy，ε_zz，ε_yz，ε_xz，ε_xy的顺序给出，初始值均设为0；该文件中包含晶粒个数用变量n_grain存储；

其中，所述σ_xx，σ_yy和σ_zz分别表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内沿X轴、Y轴和Z轴方向的正应力；

所述ε_xx，ε_yy和ε_zz分别表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内沿X轴、Y轴和Z轴方向的正应变；

所述σ_yz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在Y轴法向平面内平行于X轴的切应力；

所述σ_xz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴法向平面内平行于Z轴的切应力；

所述σ_xy表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴法向平面内平行于Y轴的切应力；

所述ε_yz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在Y轴和Z轴两方向上微小线段夹角的切应变；

所述ε_xz表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴和Z轴两方向上微小线段夹角的切应变；

所述ε_xy表示在步骤一所述有限元网格模型的坐标系中，晶粒所受到的空间内在X轴和Y轴两方向上微小线段夹角的切应变；

步骤三：

基于步骤二所述转化矩阵，通过张量运算，得到有限元网格模型坐标系下的单晶体弹性刚度矩阵Cij′和施密特因子α，其中，所述α＝cosλcosΦ，Φ为滑移面与外力F中心轴的夹角，所述λ为滑移方向与外力F的夹角，如图4所示；

所述张量运算公式为：Cij′＝R^-1·Cij；

步骤四：

判断欧拉角增量是否为(0，0，0)，若进入步骤五；若进入步骤六；

步骤五：

根据步骤一所述有限元网格模型的边界条件和加载载荷对当前有限元网格模型进行加载，然后针对当前有限元网格模型坐标系下的单晶体弹性刚度矩阵，采用有限元计算方法计算得到加载后晶粒的应力应变状态分量；将所述加载后晶粒的应力应变状态分量替换步骤二中的初始应力应变状态分量，重复步骤二至步骤四；其中，每次加载后加载步数加一；步骤六：

判断晶粒内的潜在要开动和已经开动滑移系的总个数n是否大于零；若n＝0，表示该晶粒仍处于弹性变形阶段，则重复步骤五；若n＞0，表示该晶粒进入塑性变形阶段，则进入步骤七；

潜在要开动和已经开动滑移系的判断方法：

若则该滑移系被标记为潜在要开动，若则该滑移系被标记为已经开动；其中，σ_c表示滑移系所在晶粒当前的应力值，表示第i个滑移系的临界分切应力值；所述αⁱ表示第i个滑移系的施密特因子；

步骤七：

根据Hill-Hutchinson多晶体弹塑性变形理论，认为每个晶粒的塑性应变是所有滑移系应变量(γⁱ)的累加值：

根据胡克定律，每一个加载状态都对应着一个联系晶粒应力与应变的瞬时模量(L_c)，即晶粒的瞬时应力(σ_c)与应变(ε_c)关系满足：

其中，晶粒的总应变(ε_c)为弹性应变和塑性应变的加和，代入(1)和(2)得：

其中M_c为L_c的逆矩阵。

再根据Hill(1996)的理论滑移系的临界开动应力值(τⁱ)与滑移系应变量(γⁱ)存在如下关系：

将(1)～(4)代入分切应力的关系式：

τⁱ＝σ_cαⁱ(5)

中得：

其中，最终得到的n个滑移系应变线性方程组的基向量矩阵X^ij为：

X^ij＝h^ij+αⁱL_cα^j(9)

将(9)带入(8)，求解分配给每个潜在滑移系的切变量γⁱ：

γⁱ＝fⁱε_c(10)

其中，

Y^ij为X^ij的逆。

最终在原有的晶粒应力应变瞬时模量(L_c)的基础上求解出考虑了塑性变形阶段的瞬时刚度矩阵

步骤八：

判断当前加载步数是否等于设定的最大加载步数，如果等于最大加载步数，则步骤七得到的瞬时刚度矩阵即为晶粒的各向异性线弹性本构；如果小于最大加载步数，则将步骤七的瞬时刚度矩阵替换步骤二中的单晶体弹性刚度矩阵，重复步骤二至步骤七，直至加载步数等于步骤一中设定的最大加载步数。

(1)假设取向的8晶粒模型，其中4个为α相晶粒，4个为β相晶粒。有限元以10^-3应变速率每步加载0.1％的应变，经过10步加载后可以判断出α相晶粒中的4个晶粒全部进入塑性变形阶段，β相中有3个晶粒进入塑性变形阶段，从输出结果可首先对每个晶粒在当前加载步下是否进入塑性阶段做出判断：“进入塑性阶段”或者“没有进入塑性阶段”；其次可以得知晶粒中哪些滑移系已开动以及每个滑移系对应的标号及开动量等信息，所有晶粒内滑移系的开动信息如表2和表3所示：

表2假设取向的8晶粒模型中钛合金α相晶粒内滑移系的开动结果信息

表3假设取向的8晶粒模型中钛合金β相晶粒内滑移系的开动结果信息

其中以α相的晶粒1为例，其中共有6个滑移系开动，分别是第3个、第5个、第7个、第29个、第39个和第40个，每个滑移系对应的开动量分别为：0.1155E-01、0.4345E-02、0.1872E-01、0.2221E-01、0.3231E-02和0.1375E-01；

同时还给出了每个晶粒的各向异性线弹性本构，即6*6的Cij矩阵(单位为GPa，见表4和表5)，经对比后可以发现，滑移系开动数量越多、开动量越大的晶粒，各向异性线弹性本构变化越大。

表4假设取向的8晶粒模型相中钛合金α相每个晶粒的C_ij矩阵

表5假设取向的8晶粒模型相中钛合金α相每个晶粒的C_ij矩阵

(2)假设取向的108晶粒模型，其中70个为α相晶粒，38个为β相晶粒。有限元以10^-3应变速率每步加载0.1％的应变，经过10步加载后可以判断出α相晶粒中的61个晶粒进入塑性变形阶段，β相中有28个晶粒进入塑性变形阶段。与(1)类似，可首先对每个晶粒在当前加载步下是否进入塑性阶段做出判断，模型中的晶粒数量越多，配分于每个晶粒上的应力越接近于真实受力情况，此时晶粒除了有“进入塑性阶段”和“没有进入塑性阶段”两种情况外，还有可能出现某个晶粒在上一步加载时进入了塑性阶段，但该晶粒中的滑移系硬化效应致使当前加载步增加的应力不足以使其中的任何一个滑移系再次开动，在计算过程中体现为基向量方程组无解，即在当前加载步下该晶粒的状态为“由塑性回到弹性”；其次，可以得知晶粒中哪些滑移系已开动及每个滑移系对应的标号及开动量等所有晶粒内滑移系的开动信息；最终还可以给出每个晶粒的各向异性线弹性本构，即6*6的C_ij矩阵(单位为GPa)，列举α相第27个晶粒的滑移系开动结果信息和各向异性线弹性本构C_ij矩阵信息如表6所示。

表6假设取向的108晶粒模型中α相第27个晶粒的滑移系开动结果信息和各向异性线弹性本构C_ij矩阵

本发明包括但不限于以上实施例，凡是在本发明精神的原则之下进行的任何等同替换或局部改进，都将视为在本发明的保护范围之内。