CN108629147B - 一种多晶体几何建模方法 - Google Patents

Abstract

本发明公开了一种多晶体几何建模方法，包括：利用MATLAB生成多晶体预织构，求取并获得晶粒各晶面的顶点坐标值及其排列顺序索引值；根据所获得的顶点坐标数据及其排列顺序数组索引值在ANSYS中构建多晶体的几何模型；进而采用晶粒取向随机因子法生成取向随机的各晶粒，并赋予其材料参数后得到体心立方多晶微观结构的力学模型。本发明避免了在建模过程中划分单元归属种子点时考虑相邻晶粒的单元归属的协调性等问题；考虑植入有厚度的内聚力层时更加方便，扩大了模型的使用范围，为后续的BCC金属多晶体方面的研究与优化奠定基础。

Description

一种多晶体几何建模方法

技术领域

本发明涉及多晶体的结构分析，具体地说是一种基于Voronoi拓扑原理，采用MATLAB和ANSYS混合编程的多晶微观结构的几何建模方法。

背景技术

因为多晶体材料的塑性变形常常伴随材料的应变局部化以及塑性损伤，其内部的应变及应力多方向变化造成其塑性变形的不均匀，采用解析法去研究存在一定困难。1970年左右发展了连续介质滑移理论后，即出现了有限元框架下有限变形晶体塑性的本构模型。早期模拟计算金属塑性变形的模型多数以各向同性与率无关模型为主，未考虑金属的大变形等内容，要想深入研究金属材料属性就需要构建更加精确的模型来完成模拟计算，如基于内变量模型，采用屈服准则或应变硬化等理论来合理解释金属的塑性流动等方面的内容。为了更加透彻的研究金属塑性变形，目前材料研究已经逐渐从宏观的模拟计算转变为多尺度材料模拟计算。然而，与金属塑性变形相关的有限元计算模拟中，最基础的也是最重要的工作之一是建立能够表征金属多晶体细观结构的单晶或多晶体的几何模型。与此有关的很多国内外文献多是建立在单晶体面心立方晶格的塑形基础上，再去研究多晶体。对于多晶体的研究，目前国内外学者多以面心立方金属为主，而针对体心和密排六方晶格的金属塑性模拟研究相对较少。此外，由于传统的X射线检测出的宏观的织构往往不能构建出直观的微观结构和晶粒的取向，所以在研究多晶体过程中，对晶粒的几何外形缺乏一个直观的认识和描述，以致从晶体塑性理论的诞生到现在的整个多晶体模拟过程中所建晶粒的外形上发生了很大变化：

1).比较早的多晶体研究方法是通过把一个模型划分单元后每个单元作为一个单晶体的多晶体预织构,如皮华春利用Taylor模型和有限元模型建立了面心立方晶格模型并在ABAQUS中进行了轧制织构的模拟演化进行了分析对比；王琦是以纯铝为研究对象，平面四边形单元为单晶体并构建了umat程序去模拟不同加载方式下以及晶粒尺寸、不同晶粒取向组成的多晶体内部织构受载后的演化情况等。

2).先把模型划分网格，每个单元上的高斯积分点被用作单颗晶粒并最终组成多晶体，如李大永，张少睿，彭颖红等以Mindlin壳单元模型为研究对象，该模型是把每一个高斯积分点作为单颗晶粒，构建了铜织构和铝板材料属性利用率相关晶体塑性理论并结合切线系数法求解出了剪切应变率。

3).构建规则六面体或六边形作为单颗晶粒去模拟多晶体，如Borg为了研究颗粒尺寸对多晶体材料的屈服和流动应力的影响使用了应变梯度晶体塑性理论进行计算。最终的数值解主要是通过有限元方法模拟了由40个晶粒组成的多晶体，且每个晶粒有三个滑移系。晶粒的尺寸主要通过控制d-n大小来实现，其中在初始屈服阶段n的取值范围是0.82～1.25，进入对数应变阶段时n的取值变为0.77～1.09。

4).使用电子设备制备金相图如扫描电镜照片或金相图，通过电子设备扫描样品的表面信息，后期通过处理这些信息就能合成相应的晶体模型。但是这种建模方法数据采集工作量大，最终得到的数据构建出的模型与真实模型相比还是有一点区别。

5).利用元胞自动机生成多晶体模型。

6).利用Voronoi图拓扑单晶几何外形，许多学者开始借助于Voronoi虚拟拓扑的种子点相互作用生成单个Voronoi晶粒，可以较准确的模拟出与实际晶粒外形相似的晶粒。如，Lehmann，Loehnert，Wriggers通过Delaunay三角剖分算法计算出预先给定的随机种子点得到任意尺寸的单颗Voronoi晶粒并组建多晶体结构，用于研究BCC多晶体材料的板材成型过程中的塑性力学特性，并利用体积平均的计算方法对能代表该多晶体的单元RVE进行均化；Paggi， Lehmann，Weber利用Voronoi拓扑方法模拟2D/3D的BCC多晶体微观结构，在有限元软件FEAP中模拟了晶体的塑性变形以及利用内聚力单元模拟颗粒间的结合力，发现两者之间的相互作用是非线性的；Lehmann，Faβmann，Loehnert 提出结晶行为是由晶体塑性框架下的一致连续介质力学表征的，并且把晶体塑性行为转化为基于变形梯度的弹塑性变形的乘法分解，以及由于在滑移系中的滑移率和临界切应力的本构关系下得到的粘塑性正则化。该文利用EBSD实验手段检测样本DC04钢的材料属性，并把所得到织构的数据转化为欧拉角生成 Voronoi多晶体微观模型并将模拟数据用Schmid因子来表征，最终验证其模拟的多晶体的塑性演化行为同实验获得的实验数据对比且两者较为吻合。 V.V.C.Wan,J.Jiang,D.W.Maclachlan.etc利用EBSD实测一些参数并利用 VGrain构建Voronoi结合模型的前处理数据，之后利用ABAQUS/CAE完成几何建模以及分析等工作，重点研究了NI合金的实验数据与模拟结果对比，以能量作为参考值，以及Voronoi生成晶粒去模拟断裂即研究织构的外形对材料的影响等。郑战光，冯强，袁帅等认为在Voronoi模型的构造上先根据种子点分布划定晶粒边界，都要借助于种子点的对应单元形心坐标去计算晶粒归属问题需要考虑晶粒间的协调性，此外，由于不同晶粒的参考系坐标系不同以及根据种子点对应单元的节点坐标求取形心坐标需要较大的运算量，构建多晶体几何模的效率不佳，因此在此基础上提出了一种新的基于Voronoi拓扑原理，结合MATLAB和ANSYS混合编程的多晶体几何建模方法。

综上所诉，目前CAE多晶体几何模型种类比较多，为规避建模及求解过程中的困难，多晶体模型会简化为二维模型模拟，这样会释放掉模型的一些空间自由度，且往往是以失去解的精度为代价的；构建规则外形的三维模型又没有考虑到晶粒的外形分布呈现出很强的任意性，因此目前如何构建多晶体几何模型仍然是研究多晶体的关键因素。此外，方法1),2)晶体模拟方法中在进行多晶体织构演变过程中不能够演化晶粒内部的变化情况；方法3)采用规则多面体或者是多边形未能考虑到不规则的晶界对晶粒的塑性的影响；方法4)是建模过程中工作量比较大，实现过程相对繁琐；方法5)多晶体模型依赖于元胞自动机方法；方法6)目前使用的比较多且能够较真实的反应晶粒的外形，但是目前默认为晶界是零厚度的。

目前，用于模拟多晶体织构演化所构建的几何模型主要是采用Voronoi图，其外形与真实晶粒较为接近，而且便于控制晶粒不规则外形和尺寸量级。然而采用Voronoi作为晶粒去研究多晶体的几何建模方法大同小异，部分方法主要是借MATLAB初始化Voronoi图种子点信息，之后将这些数据导入ABAQUS/CAE 中进行种子点归属区划分，其划分准则是依次求取种子点间的距离，依据“就近关联”原则将空间或二维平面剖分为多个区域且每个区域只包含一个种子点，这些独立区域就变成单颗晶粒，随后利用Python语言中mdb.models.Sketch.Line把归属区域的顶点连线并写入inp文件完成建模；还有将归属区编号并建立集合set，借助inp文件完成建模。这些建模方法虽然能够建立Voronoi模型，但是学习并使用inp文件过程较慢；需要先将整个模型划分网格，依据网格与种子点的距离来定义单元归属区，需要考虑许多因素如确定单元归属区时相邻晶粒间的协调性等因而计算步骤相比较繁琐；此外，如果多晶体模型变化频繁或者是差异性大时，需要重新划分网格，重复修改inp文件，增加了许多重复性工作量，耽误大量科研工作。

发明内容

为解决现有技术中存在的上述缺陷，本发明的目的在于提供一种以Voronoi 拓扑图作为晶胞模型，首先利用MATLAB构建多晶体预织构模型，并且在MATLAB 中完成数据的整理封存；之后利用这些数据在ANSYS中进行实体化几何建模；最后用Python对ABAQUS进行二次开发完成多晶体几何模型的力学分析，从而验证该方法所建几何模型的正确性和有效性，为后续的BCC金属多晶体方面的研究与优化奠定基础。

本发明是通过下述技术方案来实现的。

一种多晶体几何建模方法，该方法包括以下步骤：

1)基于Vonoroi拓扑原理，利用MATLAB生成多晶体预织构，求取并获得晶粒各晶面的顶点坐标值及其排列顺序索引值；

2)根据所获得的顶点坐标数据及其排列顺序数组索引值在ANSYS中构建多晶体的几何模型；

将晶粒各晶面的顶点坐标值及其排列顺序数组索引值并结合所需ANSYS命令流一并写入TXT文档中，在ANSYS经典界面中读取该TXT文档，ANSYS就会构建出多晶体几何模型；

3)获得多晶体微观结构的力学模型；

当多晶微观结构总体尺寸不变但包含不同数量晶粒时，采用晶粒取向随机因子法生成取向随机的各晶粒，并赋予其材料参数后得到体心立方多晶微观结构的力学模型。

进一步，所述步骤1)中，编写程序初始化三维Voronoi种子点，生成多晶体预组织，其具体步骤如下：

1a)初始化晶胞数量的种子点，且这些种子点都是等间距的，将三维空间平均分成n＝n_x×n_y×n_z个大小相等的子空间，且各个子空间形心坐标等间距；

1b)赋予每一形心坐标一个微小增量(UI_x,UI_y,UI_z)，使形心坐标偏离原形心坐标一个偏移量来控制晶粒颗粒的外形；

1c)将不同偏移量的形心坐标输入到MATLAB的mpt工具箱中，按照就近关联原则，生成一个种子点控制一个Voronoi区域，即得到不同外形三维的 Voronoi拓扑图。

进一步，所述步骤1)中，求取并获得晶粒各晶面的顶点坐标值及其排列顺序索引值，其具体步骤如下：

1d)编写MATLAB程序初始化三维的Voronoi种子点，调用MATLAB的mpt 工具箱生成多晶体预织构，获取各个晶粒上多面体顶点的坐标值并逐个封装在 MATLAB预先定义好的元胞数组中；

1e)逐个调用封存晶粒数据的元胞数组，并利用空间点法式平面方程找到存放晶粒上位于同一晶面的顶点坐标数组索引值，并将位于同一晶面上的顶点在数组中的索引值按单个晶面、顶点无序存储；

1f)把所有共晶面的顶点依照逆或顺时针排序，依次求解每个平面法向量同X、Y、Z轴的夹角，若法向量同X轴夹角最小，就会将该晶面顶点投影在 YOZ平面上；若同Y或Z轴的夹角最小，则投影在XOZ或XOY平面上；

1g)求解投影在YOZ平面上多边形的形心坐标值，并将该形心平移到坐标原点处，将多边形的其余顶点等距平移，求出各顶点与原点连线，以及该连线与Y或Z轴夹角余弦值，再依据投影面上该顶点所在象限求出该顶点与坐标轴的夹角，并将这些夹角排序，依据投影面上各顶点排序信息逆推出投影前晶面各顶点的空间排列顺序索引值。

进一步，所述步骤3)中，采用随机因子法生成取向随机的各晶粒，是为了使晶体模型更接近于实际，考虑多晶体内各晶粒取向的随机性对多晶材料织构演变和力学性质的影响，以及晶粒内部切应力沿着滑移方向作用从而造成的滑移位错，使晶粒内部乃至多晶材料塑性变形的不均匀性的影响。

进一步，所述随机因子法的具体步骤如下：

3a)针对多晶体内各晶粒取向的随机性对多晶体材料织构演变和力学性质的影响，定义晶粒取向随机因子方法模型如下：

其中，(α,β,γ)和(α,β,γ)_i分别是多晶集合体整体取向的欧拉角和第i个晶粒取向的欧拉角；随机因子r∈[0,1]；r_i∈[0,1]为随机生成的一组服从标准正态分布的随机数；N为多晶集合体中的晶粒个数；

3b)给定(α,β,γ)和随机因子r后，从随机生成的随机数中任取一个r_i，利用晶粒取向随机因子法公式即可确定出(α,β,γ)_i，从而逐一完成各晶粒随机取向的设置；

3c)对于晶粒内部乃至多晶体塑性变形的不均匀性，通过Schmid因子相对发生率进行预测；

若单晶内α滑移系上的Schmid因子大小为i，其中，i∈[0,0.5]；则将单晶中Schmid因子i出现的频率定义为Schmid因子i的相对发生率

其中：

∑S_i ^α为α滑移系上Schmid因子为i的频次；N是单晶中α滑移系包含的滑移面总数，M是单晶中滑移系的个数；

的大小表明了单晶乃至多晶集合体内各滑移系开动的相对难易程度。

本发明同现有技术相比具有以下几个优点：

本发明属于材料力学和结构力学领域的研究，结合了现代计算机科学的应用，具有明显的前瞻性和现实的应用前景。就本发明来说，它具有以下特点：

本发明避免了在建模过程中划分单元归属种子点时考虑相邻晶粒的单元归属的协调性等问题；

考虑植入有厚度的内聚力层时更加方便，扩大了模型的使用范围。

附图说明

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

图1是多晶体建模过程的流程图；

图2(a)-(d)分别是单个空间晶面顶点顺序求解示意图；

图3(a)-(d)分别是不同增量下的Voronoi拓扑图；其中：

图3(a)是在增量(UI_x＝1,UI_y＝1,UI_z＝1)下，正六面体晶胞多晶体的Voronoi 拓扑图；图3(b)是在增量(UI_x＝0.2,UI_y＝0.2,UI_z＝1)下，yoz半对称晶胞多晶体的Voronoi拓扑图；图3(c)是在增量(UI_x＝0.09,UI_y＝0.06,UI_z＝0.05)下，不规则晶体的Voronoi拓扑图；图3(d)是在增量(UI_x＝0.06,UI_y＝0.03,UI_z＝0.09) 下，不规则晶体的Voronoi拓扑图；

图4是MATLAB程序处理单颗Voronoi数据点实现ANSYS自下而上建模处理流程图；

图5是MATLAB程序实现Voronoi数据点处理并得到ANSYS建模命令流的 txt文档流程图；

图6(a)-(f)分别是建模过程演化和多晶体模型对比验证；其中：

图6(a)、(b)分别是ANSYS中初始化得到的单晶胞顶点(左)，逐次寻找各晶面顶点(右)的过程；图6(c)、(d)分别是ANSYS中由各晶面围成的单晶粒壳体(左)及进一步得到的单晶粒实体化模型(右)；图6(e)、(f)分别是 MATLAB中预织构的64个Voronoi晶粒与在ANSYS中完成参数化建模得到的64 个晶粒多晶体几何模型；

图7(a)、7(b)分别是多晶体微观结构有限元模型的应力以及不同变形量下微观结构在滚轧方向上的变形图；其中：

图7(a)是XOY平面上10％变形下被拉伸后的微观结构；图7(b)是30％变形下被拉伸后的微观结构；

图8(a)-(f)分别是多晶体微观结构在不同变形量不同主平面上的极图。

具体实施方式

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

本发明主要有三个流程。首先基于拓扑Voronoi拓扑原理，利用MATLAB构建多晶体预织构模型，并且在MATLAB中完成数据的整理封存；之后利用这些数据在ANSYS中进行实体化几何建模；最后用Python对ABAQUS进行二次开发完成多晶体几何模型的力学分析，从而验证该方法所建几何模型的正确性，为后续的BCC金属多晶体方面的研究与优化奠定基础。具体技术方案如下：

步骤1，利用MATLAB生成多晶体预织构模型，求取并获得晶粒各晶面的顶点坐标值及其排列顺序索引值。

编写MATLAB程序初始化3D种子点，利用MATLAB中mpt工具箱生成Voronoi 图并利用MATLAB生成多晶体预织构，编写程序初始化三维Voronoi种子点，生成多晶体预组织，其具体方法为：

1a)初始化晶粒数量的种子点，且这些种子点都是等间距的，将三维空间平均分成n＝n_x×n_y×n_z个大小相等的子空间且各个子空间形心坐标等间距；

1b)赋予每一形心坐标一个微小增量(UI_x,UI_y,UI_z)，使形心坐标偏离原形心坐标一个距离来控制颗粒的外形；

1c)将不同偏移量的形心坐标输入到MATLAB的mpt工具箱中，按照就近关联原则，生成一个种子点控制一个Voronoi拓扑区域，即得到不同外形三维的Voronoi拓扑图如图3(a),图3(b),图3(c)和图3(d)所示。

求取并获得晶胞各晶面的顶点坐标值及其排列顺序索引值，见图4所示，其具体步骤如下：

1d)编写MATLAB程序初始化三维的Voronoi种子点，调用MATLAB的mpt 工具箱生成多晶体预织构，获取各个晶粒上多面体顶点的坐标值并逐个封装在 MATLAB预先定义好的元胞数组中；

1e)逐个调用封存晶粒数据的元胞数组并利用空间点法式平面方程找到存放晶粒上位于同一晶面的顶点坐标数组索引值，并将位于同一晶面上的顶点在数组中的索引值按单个晶面存储(顶点无序)如图2(a)所示；

1f)把所有共晶面的顶点依照逆或顺时针排序，依次求解每个平面法向量同X、Y、Z轴的夹角，若法向量同X轴夹角最小就会将该晶面顶点投影在YOZ 平面上；若同Y或Z轴的夹角最小则投影在XOZ或XOY平面上，如图2(b) 所示；

1g)求解投影在YOZ平面上多边形的形心坐标值并将该形心平移到坐标原点处，将多边形的其余顶点等距平移，求出各顶点与原点连线与Y或Z轴夹角余弦值，再依据投影面上该顶点所在象限求出该顶点与坐标轴的夹角并把这些夹角排序，依据投影面上各顶点排序信息逆推出投影前晶面各顶点的空间排列顺序如图2(c),图2(d)所示。

得到各晶粒顶点坐标值并封存于MATLAB的一个元胞数组中，并基于元胞数组求得晶粒各晶面的顶点坐标值及其排列顺序索引值。其具体实现过程如流程图图4和图5所示。

图4中应用MATLAB程序处理单颗Voronoi数据点实现ANSYS自下而上建模处理流程如下：

①编写MATLAB程序，初始化3D种子点；

②基于Voronoi图，利用MATLAB中mpt生成晶粒顶点坐标值，并将其封装于MATLAB的一个元胞数组中；

③利用空间点平面方程确定出晶粒中位于同一晶面的所有顶点，并将顶点数组的索引值按晶面存储；

④将各晶面顶点按顺序排列并求其法向量与坐标轴的夹角，将该晶面投影到夹角最小的投影面上；经晶面投影图形的形心平移到坐标原点，其余顶点等距平移，依据一定规则将晶面投影图形各顶点排序，并依据该排序信息逆推从而确定出投影面各顶点的相应顺序；

⑤将晶粒各晶面顶点坐标值及其排序索引值结合ANSYS命令VA，ALL流写入txt文本中；

⑥读取txt文本并在ANSYS界面显示构建的多晶体几何模型。

图5中通过MATLAB程序实现Voronoi数据点处理并得到ANSYS建模命令流的TXT文档流程如下：

①初始化Voronoi的参数，并将发生元的位置按照偏移量开始重新生成新的发生元坐标点；

②调用MATLAB中mpt工具箱预织构出Voronoi多晶体织构；

③依照空间平面方程A*x+B*y+C*z+D＝0找到所有晶粒的所有晶面上的共面顶点，3个顶点为一组存入数组；

④调用a_face子程序将3个为一组的共面点构建平面，找到位于这个晶面上其他的属于该晶粒的其它所有共面点；

⑤调用rank_pointsl子程序将上面找到的共面顶点按照顺时针或逆时针排序，并调用unique_face子程序将重复的面除掉；

⑥将处理后的数据结合ANSYS命令流一并写入txt文档中用于单颗Voronoi 模型建立，并将所有颗粒装配的命令流逐个写入txt文档，用于多晶体装配。

步骤2，利用ANSYS建立晶胞模型并验证

根据所获得的顶点坐标数据及其排列顺序数组索引值在ANSYS中构建多晶体的几何模型；

将晶粒各晶面的顶点坐标值及其排列顺序索引值结合ANSYS命令流写入TXT 文本后，在ANSYS中读取该文本并显示出构建好的多晶体几何模型。ANSYS建模过程如图6(a)-(f)所示。MATLAB中预织构多晶体与ANSYS环境中建立的多晶体几何模型对比为图6(e)、(f)所示。图中重叠点的部分是由于ANSYS建模过程中空间点投影到一个平面上导致的，其中，图6(a)、6(b)中上部数字“22”与“44”重叠，下部数字“29”与“32”重叠，右下方数字“19”、“20”和“21”有重叠。图6(c)、6(d)同理。

步骤3，多晶体模型的力学分析模拟

当多晶微观结构总体尺寸不变但包含不同数量晶粒时，采用晶粒取向随机因子法生成取向随机的各晶粒，并赋予其材料参数后得到体心立方多晶微观结构的力学模型。

采用随机因子法生成取向随机的各晶粒，是为了使晶体模型更接近于实际，考虑多晶体内各晶粒取向的随机性对多晶材料织构演变和力学性质的影响，以及晶粒内部切应力沿着滑移方向作用从而造成的滑移位错，使晶粒内部乃至多晶材料塑性变形的不均匀性的影响，该随机因子法的具体内容如下：

3a)针对多晶体内各晶粒取向的随机性对多晶体材料织构演变和力学性质的影响，本发明提出了一种晶粒取向随机因子法，定义晶粒取向随机因子法模型如下：

其中(α,β,γ)和(α,β,γ)_i分别是多晶集合体整体取向的欧拉角和第i个晶粒取向的欧拉角；随机因子r∈[0,1]，一般给定；r_i∈[0,1]为随机生成的一组服从标准正态分布的随机数；N为多晶集合体中的晶粒个数；

3b)给定(α,β,γ)和随机因子r后，从随机生成的随机数中任取一个r_i，利用式(1)即可确定出(α,β,γ)_i，从而逐一完成各晶粒随机取向的设置。从式(1) 可知：随机因子r越大，与整体欧拉角(α,β,γ)取向不一致的晶粒越多，即具有随机取向的晶粒越多，多晶微观结构形态的随机性越强；

3c)对于晶粒内部乃至多晶体塑性变形的不均匀性，本发明提出了一种 Schmid因子相对发生率来预测这种不均匀性。若单晶内α滑移系上的Schmid因子大小为i(i∈[0,0.5])，则将单晶中Schmid因子i出现的频率定义为Schmid因子i 的相对发生率

其中：

∑S_i ^α为α滑移系上Schmid因子为i的频次；N是单晶中α滑移系包含的滑移面总数，M是单晶中滑移系的个数。

的大小表明了单晶乃至多晶集合体内各滑移系开动的相对难易程度。

下面通过具体实施例来进一步说明本发明。

在ABAQUS中输入DC04钢的材料属性见表1

表1为DCO4材料属性

通过体积模量(K)和剪切模量(G)可以计算出杨氏模量E和泊松比μ

带入表1中的数值，求解式子(2)、(3)可得到DC04钢的杨氏模量 E＝202.7GPa，相应泊松比μ＝0.278

将用本方法在ANSYS中构建的含有279颗晶粒的多晶体模型导入ABAQUS中，采用四面体十节点单元进行网格划分后，模拟了滚轧和横向两个方向上的双轴拉伸测试，并对滚轧方向和横向方向上的应力和应变等力学性能进行了分析。施加的变形量分别为相应变形方向上整个多晶体微观结构外形尺寸的10％和 30％，多晶体微观结构有限元模型的应力以及不同变形量下微观结构在滚轧方向上的变形如图7(a)、图7(b)所示。

为了更好地表征多晶体材料在外载荷下其微观织构的演变如晶粒的滑移以及晶格的旋转等，使用Python对ABAQUS进行二次开发得到如图8(a)-(f) 所示的多晶体微观结构在不同变形量不同主平面上的极图。其中，图8(a)为横向(100)上施加10%变形，图8(b)为法向(010)上施加10 ％变形，图8(c) 为滚轧方向(001)上施加10 ％变形，图8(d)为横向(100)上施加30 ％变形，图8 (e)为法向(010)上施加30％变形，图8(f)为滚轧方向(001)上施加30 ％变形。

由图8(a)-(f)可见：横向、法线方向和滚轧方向上各晶粒的初始取向相对分散，当变形量达到30％时，各晶胞的取向越来越集中到某一或某几个区域。这是因为：多晶体在变形过程中，首先少量晶粒出现滑移位错和晶格旋转，随着变形的增加，少数晶粒的变形逐渐演变为更多晶粒发生变形，最终达到稳定状态。

表2(100)和(001)两个分向上的应力应变变数据

结合图8和表2给出了施加的变形量达到30％过程中，微观结构的应力与应变关系。可见随着变形量的不断增大，参与变形以及发生旋转的晶粒数量增加，此时对应的应力数值明显变大。由此可知，织构演变过程中参与变形的晶粒越多，多晶体的硬化过程中的应力值也明显变大。

从以上实施例可以看出，本发明针对多晶体微观结构受载发生不均匀塑性变形的情况构建多晶体模型，采用了一种新的建模方法进行建模，该方法结合了MATLAB与ANSYS的APDL参数化建模，首先利用MATLAB构建多晶体预织构模型，并且在MATLAB中完成数据的整理；之后利用这些数据在ANSYS中进行实体化几何建模；最后用Python对ABAQUS进行二次开发完成多晶体几何模型的力学分析，从而验证该方法所建几何模型的正确性，为后续的BBC金属多晶体方面的研究与优化奠定基础。

本发明并不局限于上述实施例，在本发明公开的技术方案的基础上，本领域的技术人员根据所公开的技术内容，不需要创造性的劳动就可以对其中的一些技术特征做出一些替换和变形，这些替换和变形均在本发明的保护范围内。

Claims (3)

1.一种多晶体几何建模方法，其特征在于，该方法包括以下步骤：

1)基于Vonoroi拓扑原理，利用MATLAB生成多晶体预织构，求取并获得晶粒各晶面的顶点坐标值及其排列顺序索引值；

2)根据所获得的各晶粒顶点坐标数据及其排列顺序数组索引值在ANSYS中构建多晶体的几何模型；

将晶粒各晶面的顶点坐标值及其排列顺序数组索引值并结合所需ANSYS命令流一并写入TXT文档中，在ANSYS经典界面中读取该TXT文档，ANSYS就会构建出多晶体几何模型；

3)获得多晶体微观结构的力学模型；

当多晶微观结构总体尺寸不变但包含不同数量晶粒时，采用晶粒取向随机因子法生成取向随机的各晶粒，并赋予其材料参数后得到体心立方多晶微观结构的力学模型；

所述步骤1)中，编写程序初始化三维Voronoi种子点，生成多晶体预织构，其具体步骤如下：

1a)初始化晶粒数量的种子点，且这些种子点都是等间距的，将三维空间平均分成n＝n_x×n_y×n_z个大小相等的子空间，且各个子空间形心坐标等间距；

1b)赋予每一形心坐标一个微小增量(UI_x,UI_y,UI_z)，使形心坐标偏离原形心坐标一个偏移量来控制晶胞颗粒的外形；

1c)将不同偏移量的形心坐标输入到MATLAB的mpt工具箱中，按照就近关联原则，生成一个种子点控制一个Voronoi区域，即得到不同外形三维的Voronoi拓扑图；

所述步骤1)中，求取并获得晶粒各晶面的顶点坐标值及其排列顺序索引值，其具体步骤如下：

1d)编写MATLAB程序初始化三维的Voronoi种子点，调用MATLAB的mpt工具箱生成多晶体预织构，获取各个晶粒上多面体顶点的坐标值并逐个封装在MATLAB预先定义好的元胞数组中；

1e)逐个调用封存晶粒数据的元胞数组，并利用空间点法式平面方程找到存放晶粒上位于同一晶面的顶点坐标数组索引值，并将位于同一晶面上的顶点在数组中的索引值按单个晶面、顶点无序存储；

1f)把所有共晶面的顶点依照逆或顺时针排序，依次求解每个平面法向量同X、Y、Z轴的夹角，若法向量同X轴夹角最小，就会将该晶面顶点投影在YOZ平面上；若同Y或Z轴的夹角最小，则投影在XOZ或XOY平面上；

1g)求解投影在YOZ平面上多边形的形心坐标值，并将该形心平移到坐标原点处，将多边形的其余顶点等距平移，求出各顶点与原点连线，以及该连线与Y或Z轴夹角余弦值，再依据投影面上该顶点所在象限求出该顶点与坐标轴的夹角，并将这些夹角排序，依据投影面上各顶点排序信息逆推出投影前晶面各顶点的空间排列顺序索引值；

所述采用晶粒取向随机因子法生成取向随机的各晶粒，过程如下：

针对多晶体内各晶粒取向的随机性对多晶体材料织构演变和力学性质的影响，定义晶粒取向随机因子方法模型；

利用晶粒取向随机因子法模型即可确定出多晶集合体整体取向的欧拉角(α,β,γ)_i，从而逐一完成各晶粒随机取向的设置；

通过Schmid因子相对发生率进行预测；若单晶内α滑移系上的Schmid因子大小为i，其中，i∈[0,0.5]，则将单晶中Schmid因子i出现的频率定义为Schmid因子i的相对发生率。

2.根据权利要求1所述的多晶体几何建模方法，其特征在于，所述步骤3)中，采用随机因子法生成取向随机的各晶粒，是为了使晶体模型更接近于实际，考虑多晶体内各晶粒取向的随机性对多晶材料织构演变和力学性质的影响，以及晶粒内部切应力沿着滑移方向作用从而造成的滑移位错，使晶粒内部乃至多晶材料塑性变形的不均匀性的影响。

3.根据权利要求2所述的多晶体几何建模方法，其特征在于，所述随机因子法的具体步骤如下：

3a)针对多晶体内各晶粒取向的随机性对多晶体材料织构演变和力学性质的影响，定义晶粒取向随机因子方法模型如下：

其中，(α,β,γ)和(α,β,γ)_i分别为多晶集合体整体取向的欧拉角和第i个晶粒取向的欧拉角；随机因子r∈[0,1]；r_i∈[0,1]为随机生成的一组服从标准正态分布的随机数；N为多晶集合体中的晶粒个数；

3b)给定(α,β,γ)和随机因子r后，从随机生成的随机数中任取一个r_i，利用晶粒取向随机因子法公式即可确定出(α,β,γ)_i，从而逐一完成各晶粒随机取向的设置；

3c)对于晶粒内部乃至多晶体塑性变形的不均匀性，通过Schmid因子相对发生率进行预测；

若单晶内α滑移系上的Schmid因子大小为i，其中，i∈[0,0.5]，则将单晶中Schmid因子i出现的频率定义为Schmid因子i的相对发生率

其中：

∑S_i ^α为α滑移系上Schmid因子为i的频次；N为单晶中α滑移系包含的滑移面总数，M为单晶中滑移系的个数；

的大小表明了单晶乃至多晶集合体内各滑移系开动的相对难易程度。

Images