CN109063275B - 基于feap的三维多晶微观结构材料模型的构建方法 - Google Patents

基于feap的三维多晶微观结构材料模型的构建方法 Download PDF

Info

Publication number
CN109063275B
CN109063275B CN201810765570.8A CN201810765570A CN109063275B CN 109063275 B CN109063275 B CN 109063275B CN 201810765570 A CN201810765570 A CN 201810765570A CN 109063275 B CN109063275 B CN 109063275B
Authority
CN
China
Prior art keywords
model
feap
umat
crystal
dimensional
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810765570.8A
Other languages
English (en)
Other versions
CN109063275A (zh
Inventor
王云飞
马娟
高静丽
孙逸
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201810765570.8A priority Critical patent/CN109063275B/zh
Publication of CN109063275A publication Critical patent/CN109063275A/zh
Application granted granted Critical
Publication of CN109063275B publication Critical patent/CN109063275B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明公开了一种基于FEAP的三维多晶微观结构材料模型的构建方法，包括：预织构出三维多晶Voronoi几何模型；对几何模型实体化、网格划分，将得到的inp文件转译成FEAP软件下的input文件，完成三维多晶微观结构几何建模；基于Fortran程序运行平台，将编写的St.Venant Kirchhoff晶体模型和晶体弹塑性模型的用户材料子程序UMAT添加到FEAP软件的模型库中，与FEAP进行集成安装，得到三维多晶微观结构的材料模型；给不同晶粒指定不同的晶格方向，得到三维多晶微观结构的力学模型；经有限元力学模拟分析验证所编写用户材料子程序UMAT的有效性。本发明完成了两种材料模型用户材料子程序UMAT的编写，并进行有效性验证，为多晶材料的均化研究提供了理论基础。

Description

基于FEAP的三维多晶微观结构材料模型的构建方法
技术领域
本发明属于材料力学和结构力学领域的研究，具体是基于FEAP的三维多晶微观结构材料模型的创建方法。本发明是对BCC结构的多晶体建立的可描述其微观结构和变形机制的数值模型，可从微观尺度对BCC结构金属材料进行力学分析，可用于对BCC结构金属材料的成型进行预测，具有一定的工程实际意义。
背景技术
金属材料通过结晶形成，其在诸多工程领域中发挥着不可替代的作用。由于复杂结晶方式导致不同晶核控制生成的晶体有一定差异，故又称为多晶材料。多晶材料的复杂结晶方式也使其在微观尺度上呈现出较强的几何不规则性和复杂的变形机制，进而使得其宏观力学响应受到诸多因素如结晶几何形状、晶格方向和滑移系统等的影响。因此，想要通过计算机方法研究金属材料的力学性质，就必须从其微观结构出发，构建真实有效的模型。正是由于非均质材料组成成分和结构的多样性，使得对其力学性质的研究较均质材料更为复杂和困难。
目前，对于非均质材料多用均化方法对其材料性质和变形机制进行相关研究，该方法是通过创建一个介观结构即表征体积单元(RVE：Representative Volume Element)去联系微观结构和宏观响应之间的关系，主题思想是先创建微观几何模型，进而给微观几何模型赋予应有的材料属性，得到其力学模型，并对多组不同尺寸下的非均质材料微观模型进行力学测试，其力学响应收敛于容差范围内的非均质材料微观模型尺寸就是目标材料RVE模型的尺寸，其力学响应就是宏观下的该非均质材料的力学响应。对于有限变形下非均质材料均化问题的研究，由于其存在强烈的非线性，无法构建显示方程去描述，故Kouznetsova等人提出了数值方法，对于随机参数非均质材料的均化问题的研究，Vel和Goupee通过随机形态描述方程给出了其随机微观结构的建模方法，Cottereau提出了一种数值耦合法去降低数值模拟的计算量，Feng和Li提出了一种高效的基于高斯随机场的非线性变换算法，将随机问题与统计学相结合。对于随机微观结构在力学和热学上的多尺度均化分析，马娟等不仅考虑了随机微观几何形状的影响，还涉及了随机参数的对非均质材料宏观响应的影响。对于热弹性方面的均化研究，Rosen和Hashin做出了较早的尝试，Francfort则首先提出了线性热弹性周期性介质的均化方法。对充分考虑微观结构不确定性的非均质材料的均化分析，已发表的文献相对较少，仅有少部分学者做过相应的研究。
为了更加准确的描述晶体结构的材料响应，除了构建合理的材料模型之外，其几何模型的构建对最终的响应结果也有着至关重要的影响。迄今为止，晶体几何建模发展了几十年，不外乎以下几种方式：(1)将每个单元作为一个单晶体；这种建模方式最为简单，但单元划分往往为六面体的规则单元，且这些单元构成的多晶体都是规则形状，无法从几何外形上突出晶体真实的微观结构；例如皮华春王琦就是通过这种方式来模拟晶体几何形状的。(2)取每个单元上相同位置的高斯积分点来表示一组晶粒；其主要是通过对高斯积分点的积分去求解问题，因此它无法反映晶粒内部发生的诸如晶格旋转等的变化，具有很强的局限性。(3)通过规则的几何形状去描述单晶；这种方法通过正六边形的排列去描述晶粒，除了与(1)具有相同的局限性之外，它还无法反应晶粒之间的不规则晶界，例如Borg就是通过这种方式去构建40颗粒的多晶体模型从而研究晶体的屈服和塑性流动。(4)通过元胞自动生成器去构建晶粒的几何外形，这种方法可以生成近似于晶体结晶的几何外形，但是由于其自动生成，使得科研工作者无法对其进行自如的控制，因此也具有很大的局限性。(5)通过扫描电镜照片和金相图片得到试样表面的模型，然后仿照得到的模型去构建相似的晶粒模型；这种方式获得的晶粒几何形状与真实最为接近，但是其不具有一般性，而且建模方式复杂、工作量大，无法得到推广。(6)通过Voronoi图拓扑的方法去获得能近似描述晶粒几何外形的微观结构；Voronoi图是计算机图形学中重要的数据结构之一，自从其问世以来就得到了快速的发展，直到上世纪80年代达到了辉煌阶段，它不仅应用于结晶学，在工程优化和气象学领域也有着极其广泛的应用。Voronoi图拓扑出的晶体颗粒，不仅能从几何外形上接近真实晶粒，也使得使用者更容易对其形状进行控制和干涉，因此在现在的晶体结构材料力学响应研究中，其建模方式最为准确，例如司良英、Lehmann都是通过这种方法去构建晶体的几何模型的。
发明内容
基于上述问题，本发明采用具备完善均化模块的FEAP软件进行三维多晶微观结构的构建，首先基于Voronoi拓扑原理，结合Matlab预织构出三维多晶Voronoi几何模型，在ANSYS下对其实体化，然后结合ABAQUS软件完成FEAP下的三维多晶微观结构的几何建模；然后进行FEAP有限元分析平台的集成与创建及St.Venant Kirchhoff晶体模型和晶体弹塑性模型用户材料子程序UMAT的编写，进一步完成基于FEAP的三维多晶微观结构的材料模型的建模，最后对所得几何模型指定不同材料属性，得到三维多晶微观结构的力学模型，经有限元力学模拟分析验证所编写用户材料子程序UMAT的有效性，为后续将要开展的多晶金属材料的均化研究提供了参考依据。
本发明是通过下述技术方案来实现的。
本发明一种基于FEAP的三维多晶微观结构材料模型的构建方法，包括以下步骤：
(1)基于Voronoi拓扑原理，使用MATLAB控制各种子点空间位置的方法预织构出三维多晶Voronoi几何模型；
(2)在ANSYS下对生成的三维多晶Voronoi几何模型实体化，然后通过Python语言脚本导入到ABAQUS中并进行网格划分，得到相应的inp文件；
(3)将得到的inp文件通过MATLAB转译成FEAP软件下的input文件，进而完成基于FEAP软件的三维多晶微观结构的几何建模；
(4)创建有限元分析软件FEAP的Fortran程序运行平台；
(5)建立St.Venant Kirchhoff晶体模型和晶体弹塑性模型两种材料模型，并基于FEAP软件完成了这两种材料模型的用户材料子程序UMAT的编写；
(6)将编写的两种材料模型的用户材料子程序UMAT添加到FEAP软件的模型库中，与FEAP进行集成安装，得到三维多晶微观结构的材料模型；
(7)在FEAP软件下对所得几何的三维多晶微观结构的几何模型指定不同材料属性，得到三维多晶微观结构的力学模型；
(8)对三维多晶微观结构的力学模型进行有限元力学模拟分析，验证所编写材料模型的用户材料子程序UMAT的有效性。
进一步，所述步骤(1)是通过程序控制生成种子点的数量和位置，然后调用MATLAB中的多参数工具箱MPT集成的程序包完成三维Voronoi多晶几何模型的预织构。
进一步，所述步骤(2)的过程为：
(2a)通过MATLAB对三维多晶Voronoi几何模型所有单晶各个顶点进行排序和读写操作，结合ANSYS参数化设计语言APDL，对MATLAB写出的包含模型顶点信息的命令流语言文本进行读取完成单颗晶粒的建模，依次循环完成所有组成多晶微观结构的单颗晶粒的建模，从而完成实体化建模，并以iges格式写出模型信息；
(2b)通过MATLAB下写出的Python脚本命令将iges模型导入ABAQUS中，在ABAQUS下生成该三维多晶体的几何模型，并对包容盒边界上不完整的Voronoi多面体进行去边处理；
(2c)对去边处理后的Voronoi多晶模型进行网格划分，得到相应的inp文件。
进一步，所述步骤(3)中，将步骤(2)中得到的inp文件通过MATLAB转译成FEAP软件下的input文件，从而完成FEAP下的三维多晶微观结构的几何建模。
进一步，所述步骤(4)有限元分析软件FEAP的Fortran程序运行平台的创建的过程为：
(4a)通过批处理方式，将所有FEAP的Fortran文件以及C语言头文件复制到编译文件夹中，为Visual Studio的编译工作做好准备；
(4b)使用Visual Studio 2010创建一个新的工程，工程模板选择Inter(R)VisualFortran-Library-Static Library，指定工程名称和路径；选择发布版本的解决方案配置和对应的Windows解决方案平台；通过添加已存在的项目，将步骤(4a)中编译文件夹所有的Fortran程序文件添加到工程项目中；选择工程配置属性，最后，选择生成解决方案，便在工程发布文件夹下生成了所有程序的源代码二进制目标文件；
(4c)使用Visual Studio 2010再创建一个新的工程，工程模板选择Inter(R)Visual Fortran-QuickWin Application-Static Library，指定工程名称和路径，路径要与步骤(4b)在同一路径下；
其他步骤与步骤(4b)相同，只是在添加存在的项目时选择的是步骤(4b)中生成的工程名二进制目标文件库文件和Fortran主程序文件。
进一步，所述步骤(5)中，St.Venant Kirchhoff晶体模型和晶体弹塑性模型两种材料模型用户子程序UMAT的编写，是通过向FEAP软件中添加UMATLn和UMATIn两个接口程序完成，其中，n为0～9之间的整数，UMATLn是用户自己编写的，材料本构关系包含在其中，当其被单元调用时，便会接收传递进来的变形信息，完成每个单元积分点的计算并返回应力和切线刚度模量；UMATLIn完成宏命令的指派，并定义模型所需参数，从而使得可以通过输入文件实现对UMAT的调用。
进一步，对于UMATIn接口程序的架构包括：
type，用于指定宏命令名，当输入文件中使用该宏命令时，FEAP便调用对应的UMAT程序进行计算；
vv，为实型数组，当UMAT被调用时，FEAP从输入文件UMAT宏命令名后读取材料参数并存储到vv实型数组中；
d，为FEAP标准程序模块参数；
n1和n3，用于指定用户自定义材料本构程序中的历史变量数；
ud，为用户材料参数数组，通过vv数组对其赋值从而参与程序的计算，通过UMATIn结构程序，完成用户自定义UMAT的参数读取和赋值。
进一步，对于UMATLn接口子程序，其架构不同于UMATIn只做参数的定义、赋值和传递，UMATLn则是通过接收传入参数，在单元调用时进行有限元计算，然后返回应力值和切线刚度模量，也即是真正参与模拟运算的子程序。
进一步，所述步骤(5)中，对St.Venant Kirchhoff晶体模型和晶体弹塑性模型两种材料模型用户子程序UMAT的编写，通过UMATLn和UMATIn两个接口程序的编译；其中，对晶体弹塑性模型用户子程序UMAT的编译通过ABAQUS用户子程序进行改写，其主要过程如下：
(5a)完成变形梯度的极分解，将n+1时刻和n时刻的变形梯度分别通过极分解得到相应时刻的旋转张量和拉伸张量；
(5b)通过两个时刻的旋转张量得到n到n+1时刻的旋转增量，再将两个时刻的拉伸张量通过计算得出相应的对数应变，然后得到n到n+1时刻的应变增量；
(5c)处理ABAQUS下历史变量数组STATEV到FEAP下历史变量数组hn和h1的转换问题，使得不同时刻下FEAP能成功完成历史变量的更新。
而对于St.Venant Kirchhoff晶体模型则主要是处理立方晶体材料晶格方向的构建以及通过左、右Cauchy-Green变形张量得到Green-Lagrange应变张量，从而得到模型整体的切线刚度矩阵和Cauchy应力的计算输出。
进一步，所述步骤(6)中，主要是通过Visual Studio 2010以及Parallel StudioXE 2011使自编写的St.Venant Kirchhoff晶体模型和晶体弹塑性模型与FEAP进行集成。
进一步，所述步骤(7)中，通过在FEAP下对不同Voronoi单晶进行组集，指定不同的诸如晶格方向的材料属性，进而得到三维多晶微观结构的力学模型。
进一步，所述步骤(8)中，通过运用FEAP自身的Hill材料模型的力学响应与自编写的St.Venant Kirchhoff晶体材料模型和晶体弹塑性材料模型的力学响应进行对比，以验证模型的有效性；以体心立方BCC结构金属DC04钢的材料属性为基础，结合上述几何建模和材料建模，对三维多晶微观结构进行Hill材料模型、St.Venant Kirchhoff晶体材料模型和晶体弹塑性材料模型三种材料模型下的有限元模拟分析，对比三种材料模型的相通性和差异性。
本发明与现有技术相比，具有以下几个特点：
1.本发明属于材料力学和结构力学领域的研究，结合了现代计算机科学的应用，具有明显的前瞻性和现实的应用前景。
2.本发明可从微观尺度对BCC结构金属材料进行力学分析，可以对BCC结构金属材料的成型进行预测，具有一定的工程实际意义。
3.本发明以BCC结构的金属材料为研究对象，充分考虑了其几何和材料两个层次上的客观因素，具有较强的理论和实际意义。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的不当限定，在附图中：
图1为本发明方法流程图；
图2为Voronoi图的算法实现流程图；
图3(a)-(d)分别为Voronoi多晶微观结构几何模型建模过程图；其中：
图3(a)为MATLAB下预织构出Voronoi多晶模型；
图3(b)为ANSYS下生成的Voronoi多晶模型；
图3(c)为ABAQUS下生成的Voronoi多晶模型；
图3(d)为ABAQUS下去边处理后的Voronoi多晶模型；
图4为ABAQUS inp文件到FEAP input文件转译算法流程图；
图5为FEAP PostScript下的几何模型图；
图6为FEAP运行窗口界面；
图7为UMATIn接口子程序架构示意图；
图8为UMATLn接口子程序架构示意图；
图9为44颗粒三维Voronoi多晶微观结构几何模型及边界条件示意图；
图10(a)-(f)为44颗粒三维Voronoi多晶微观结构三种材料模型下S33应力云图；其中：
图10(a)为随机晶格方向下晶体弹塑性模型S33应力云图；
图10(b)为相同晶格方向下晶体弹塑性模型S33应力云图；
图10(c)为随机晶格方向下St.Venant Kirchhoff晶体模型S33应力云图；
图10(d)为相同晶格方向下St.Venant Kirchhoff晶体模型S33应力云图；
图10(e)为随机晶格方向下Hill模型S33应力云图；
图10(f)为相同晶格方向下Hill模型S33应力云图。
具体实施方式
下面将结合附图以及具体实施例来详细说明本发明，在此本发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。
参照图1，本发明基于FEAP的三维多晶微观结构材料模型的构建方法，具体步骤如下：
步骤1，基于Voronoi拓扑原理，使用MATLAB控制各种子点空间位置的方法预织构出三维多晶Voronoi几何模型
MATLAB中集成了生成二维和三维Voronoi图的算法，可以通过调用多参数工具箱MPT实现，因此只需通过程序控制生成种子点的数量和位置，直接调用MPT工具箱集成的程序包便可以完成三维Voronoi多晶几何模型的预织构。
Voronoi算法的实现流程图如图2所示，通过该算法便可以在一个指定大小的六面体空间中预织构出三维Voronoi模型。如图3(a)所示，为MATLAB预织构出的含有121个颗粒的Voronoi多晶模型。
步骤2，在ANSYS下对生成的三维多晶Voronoi几何模型实体化，然后通过Python语言脚本导入到ABAQUS中并进行网格划分，得到相应的inp文件
(2a)通过MATLAB对三维多晶Voronoi几何模型所有单晶各个顶点进行排序和读写操作，结合ANSYS参数化设计语言APDL，对MATLAB写出的包含模型顶点信息的命令流语言文本进行读取完成单颗晶粒的建模，依次循环完成所有组成多晶微观结构的单颗晶粒的建模，从而完成实体化建模，并以iges格式写出模型信息，如图3(b)所示；
(2b)通过MATLAB下写出的Python脚本命令将iges模型导入ABAQUS中，在ABAQUS下生成该三维多晶体的几何模型，如图3(c)所示，由于预织构时包容盒边界的限制，使得边界上的Voronoi多面体不完整，故须进行去边处理，如图3(d)所示为121个多晶颗粒去边处理后在ABAQUS下生成的模型；
(2c)对去边处理后的Voronoi多晶模型进行网格划分，得到相应的inp文件。
步骤3，将得到的inp文件通过MATLAB转译成FEAP软件下的input文件，进而完成基于FEAP的三维多晶微观结构的几何建模
MATLAB转译算法的实现流程如图4所示，在FEAP中生成的去边处理后三维多晶几何模型如图5所示，其通过FEAP的PostScript宏命令创建。将步骤2中得到的inp文件通过MATLAB转译成FEAP软件下的input文件，通过ABAQUS的inp文件转译得到FEAP下的数字化几何模型后，还需要在FEAP输入文件中编译受载和约束等边界条件。
步骤4，有限元分析软件FEAP的Fortran程序运行平台的创建
有限元分析软件FEAP的Fortran程序运行平台的创建过程如下：
(4a)通过批处理方式，将所有FEAP的Fortran文件以及C语言头文件复制到编译文件夹中，为Visual Studio的编译工作做好准备。
(4b)使用Visual Studio 2010创建一个新的工程，工程模板选择Inter(R)VisualFortran-Library-Static Library，指定工程名称和路径；选择发布版本的解决方案配置和对应的Windows解决方案平台；通过添加已存在的项目，将步骤(4a)中编译文件夹所有的Fortran程序文件添加到工程项目中；选择工程配置属性，最后，选择生成解决方案，便在工程发布文件夹下生成了所有程序的源代码二进制目标文件。
(4c)使用Visual Studio 2010再创建一个新的工程，工程模板选择Inter(R)Visual Fortran-QuickWin Application-Static Library，指定工程名称和路径，路径要与步骤(4b)在同一路径下。其他步骤与步骤(4b)相同，只是在添加存在的项目时选择的是步骤(4b)中生成的工程名二进制目标文件库文件和Fortran主程序文件。
经过上述过程后便生成了FEAP运行窗口主界面，如图6所示。
步骤5，建立St.Venant Kirchhoff晶体模型和晶体弹塑性模型两种材料模型，并基于FEAP完成了这两种材料模型的用户材料子程序UMAT的编写
St.Venant Kirchhoff晶体模型和晶体弹塑性模型两种材料模型用户材料子程序UMAT的编写，是通过向FEAP软件中添加UMATLn和UMATIn两个接口程序完成，其中，n为0～9之间的整数，UMATLn是用户自己编写的，材料本构关系包含在其中，当其被单元调用时，便会接收传递进来的变形信息，完成每个单元积分点的计算并返回应力和切线刚度模；UMATLIn完成宏命令的指派，并定义模型所需参数，从而使得可以通过输入文件实现对UMAT的调用。对于UMATIn接口程序的架构如图7所示，其包含五部分内容：
(a)type用于指定宏命令名，当输入文件中使用该宏命令时，FEAP便调用对应的UMAT程序进行计算。
(b)vv为实型数组，当UMAT被调用时，FEAP从输入文件UMAT宏命令名后读取材料参数并存储到vv数组中。
(c)d为FEAP标准程序模块参数。
(d)n1和n3用于指定用户自定义材料本构程序中的历史变量数。
(e)ud为用户材料参数数组，通过vv数组对其赋值从而参与程序的计算。通过UMATIn结构程序，完成用户自定义UMAT的参数读取和赋值。
如图8所示为有限变形模式下UMATLn接口子程序架构示意图，不同于UMATIn只做参数的定义、赋值和传递，UMATLn则是通过接收传入参数，在单元调用时进行有限元计算，然后返回应力值和切线刚度模量，也即是真正参与模拟运算的子程序。
进一步，步骤5中，通过UMATLn和UMATIn两个接口程序的编译，完成对St.VenantKirchhoff晶体模型和晶体弹塑性模型两种材料模型用户子程序UMAT的编写；其中，对晶体弹塑性模型用户子程序UMAT的编译主要是通过ABAQUS用户子程序进行改写的，其主要过程如下：
(5a)完成变形梯度的极分解，将n+1时刻和n时刻的变形梯度分别通过极分解得到相应时刻的旋转张量和拉伸张量；
(5b)通过两个时刻的旋转张量得到n到n+1时刻的旋转增量，再将两个时刻的拉伸张量通过计算得出相应的对数应变，然后得到n到n+1时刻的应变增量；
(5c)处理ABAQUS下历史变量数组STATEV到FEAP下历史变量数组hn和h1的转换问题，使得不同时刻下FEAP能成功完成历史变量的更新。
而对于St.Venant Kirchhoff晶体模型则主要是处理立方晶体材料晶格方向的构建以及通过左、右Cauchy-Green变形张量得到Green-Lagrange应变张量，从而得到模型整体的切线刚度矩阵和Cauchy应力的计算输出。
步骤6，将编写的两种材料模型的用户材料子程序UMAT添加到FEAP软件的模型库中，与FEAP进行集成安装，得到三维多晶微观结构的材料模型
通过Visual Studio 2010以及Parallel Studio XE 2011将编写的St.VenantKirchhoff晶体模型和晶体弹塑性模型与FEAP进行集成。
步骤7，在FEAP下对所得几何模型指定不同材料属性，得到三维多晶微观结构的力学模型
通过在FEAP下对不同Voronoi单晶进行组集，指定不同的诸如晶格方向的材料属性，进而得到三维多晶微观结构的力学模型。
步骤8，对三维多晶微观结构的力学模型进行有限元力学模拟分析，验证所编写材料模型UMAT的有效性
通过运用FEAP自身的Hill材料模型的力学响应与编写的St.Venant Kirchhoff晶体材料模型和晶体弹塑性材料模型的力学响应进行对比，以验证编写的用户材料子程序UMAT的有效性；本发明以体心立方BCC结构金属DC04钢的材料属性为基础，结合上述几何建模和材料建模，对三维多晶微观结构的力学模型进行Hill材料模型、St.Venant Kirchhoff晶体材料模型和晶体弹塑性材料模型三种材料模型下的有限元模拟分析，对比了三种材料模型的相通性和差异性。
下面通过具体实例来进一步说明本发明。
如表1所示为BCC结构的DC04钢的力学属性表，其中K为体积模量，G为剪切模量，τ0为初始剪切屈服应力，h0为初始硬化模量，σs和σb为宏观上的屈服极限与强度极限。
表1 DC04钢的力学属性表
Figure GDA0002412994150000141
对于Hill模型，假设其弹性阶段是根据材料的弹性模量E和泊松比μ去表征的，这两个物理参数可根据表1由下式求得：
Figure GDA0002412994150000142
Figure GDA0002412994150000143
此外，在塑性阶段，Hill模型在FEAP中算法的实现是通过Hill参数去描述各个方向的屈服状态的，对于Hill参数的求解，可通过式(3)和(4)计算得出：
Figure GDA0002412994150000144
Figure GDA0002412994150000145
其中，r0、r45和r90分别为单轴拉伸下切取试样与滚轧方向成0°、45°和90°时的Hill塑性应变比，根据文献资料得r0＝2.057、r45＝1.852和r90＝2.612，代入上式可得R22＝1.044，R33＝1.047，R12＝0.618。
而对于FEAP下自编写的St.Venant Kirchhoff模型和晶体弹塑性模型，其弹性阶段的三个弹性常数可由材料的弹性模量E和泊松比μ通过解析的方法获得，如式(5)至(7)所示：
Figure GDA0002412994150000146
Figure GDA0002412994150000147
Figure GDA0002412994150000148
对于晶体弹塑性模型塑性阶段的初始剪切应力τ0和极限剪切应力τs，采用由FCC结构晶体给出的解析公式近似替代：
τs＝1.8τ0 (8)
则τs＝144MPa，对于剪切应变率和率敏感系数，Voronoi多晶模拟取为
Figure GDA0002412994150000151
m＝1。
故可得三种材料模型下的材料参数分别如表2所示：
表2算例中三种材料模型下的材料参数
Figure GDA0002412994150000152
基于以上三种材料模型及其参数，用本发明所述的建模方法创建一个包含44个颗粒的三维Voronoi多晶微观结构模型，如图9所示，边界条件采用一端三向固定，一端两向固定并在第三向三向拉伸，图红色标识即为所有约束节点。
在上述边界条件下，随机出一组立方晶格方向向量(材料方向)创建多晶模型，并基于St.Venant Kirchhoff晶体模型和晶体弹塑性模型进行纵向拉伸；其中晶体弹塑性模型下为10％左右的变形量，St.Venant Kirchhoff晶体模型采用其收敛情况下的最大变形量，作出三种模型下的S33应力云图如图10(a)-(f)所示。
由图10(a)-(f)可知，对于同一个Voronoi多晶微观结构模型，其在三种材料模型下，当每个单晶晶格方向(材料主轴方向)随机且不相同时，多晶材料的力学响应与在相同晶格方向下(材料主轴方向)的力学响应存在显著区别。但从模型表面应力分布云图图10(e)和图10(f)来看，在Hill模型下图10(f)表面应力所在的两个区间均处图10(e)表面应力分布的区间中，因此Hill模型下材料主轴变化对应力响应的影响较低。此外，晶体弹塑性模型下其表面应力分布变化大于St.Venant Kirchhoff晶体模型的应力变化，这是因为晶格方向不同导致了滑移系统的两个方向向量在全局坐标系中的不同。
从上述实例可以看出，本发明以体心立方结构为基础的DC04钢为研究对象进行力学响应的研究，充分考虑了其几何和材料两个层次上的客观因素，完成了基于FEAP的三维多晶微观结构的力学模型的构建，并对不规则Voronoi多晶微观结构的力学模型进行有限元力学模拟分析，验证了所编写的用户材料子程序UMAT的有效性，为后续将要开展的BCC结构的多晶金属材料的均化研究奠定了基础。
本发明并不局限于上述实例，在本发明公开的技术方案的基础上，本领域的技术人员根据所公开的技术内容，不需要创造性的劳动就可以对其中的一些技术特征作出一些替换和变形，这些替换和变形均在本发明的保护范围内。

Claims (9)

1.基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，该方法包括以下步骤：
(1)基于Voronoi拓扑原理，使用MATLAB控制各种子点空间位置的方法预织构出三维多晶Voronoi几何模型；
(2)在ANSYS下对生成的三维多晶Voronoi几何模型实体化，然后通过Python语言脚本导入到ABAQUS中并进行网格划分，得到相应的inp文件；
(3)将得到的inp文件通过MATLAB转译成FEAP软件下识别的input文件，进而完成基于FEAP软件的三维多晶微观结构的几何建模；
(4)创建有限元分析软件FEAP的Fortran程序运行平台；
(5)建立St.Venant Kirchhoff晶体模型和晶体弹塑性模型两种材料模型，并基于FEAP软件完成了这两种材料模型的用户材料子程序UMAT的编写；
所述步骤(5)中，通过UMATLn和UMATIn两个接口程序的编译，完成对St.VenantKirchhoff晶体模型和晶体弹塑性模型两种材料模型用户材料子程序UMAT的编写；其中，对晶体弹塑性模型用户材料子程序UMAT的编译通过ABAQUS用户子程序进行改写，其主要过程如下：
(5a)完成变形梯度的极分解，将n+1时刻和n时刻的变形梯度分别通过极分解得到相应时刻的旋转张量和拉伸张量；
(5b)通过两个时刻的旋转张量得到n到n+1时刻的旋转增量，再将两个时刻的拉伸张量通过计算得出相应的对数应变，然后得到n到n+1时刻的应变增量；
(5c)处理ABAQUS下历史变量数组STATEV到FEAP下历史变量数组hn和h1的转换问题，使得不同时刻下FEAP能成功完成历史变量的更新；
而对于St.Venant Kirchhoff晶体模型则主要是处理立方晶体材料晶格方向的构建以及通过左、右Cauchy-Green变形张量得到Green-Lagrange应变张量，从而得到模型整体的切线刚度矩阵和Cauchy应力的计算输出；
(6)将编写的两种材料模型的用户材料子程序UMAT添加到FEAP软件的模型库中，与FEAP进行集成安装，得到三维多晶微观结构的材料模型；
(7)在FEAP软件下对所得的三维多晶微观结构的几何模型指定不同材料属性，得到三维多晶微观结构的力学模型；
(8)对三维多晶微观结构的力学模型进行有限元力学模拟分析，验证所编写两种材料模型的用户材料子程序UMAT的有效性。
2.根据权利要求1所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，所述步骤(1)是通过程序控制生成种子点的数量和位置，然后调用MATLAB中的多参数工具箱MPT集成的程序包完成三维Voronoi多晶几何模型的预织构。
3.根据权利要求1所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，所述步骤(2)的过程为：
(2a)通过MATLAB对三维多晶Voronoi几何模型所有单晶各个顶点进行排序和读写操作，结合ANSYS参数化设计语言APDL，对MATLAB写出的包含模型顶点信息的命令流语言文本进行读取完成单颗晶粒的建模，依次循环完成所有组成多晶微观结构的单颗晶粒的建模，从而完成实体化建模，并以iges格式写出模型信息；
(2b)通过MATLAB下写出的Python脚本命令将iges模型导入ABAQUS中，在ABAQUS下生成该三维多晶体的几何模型，并对包容盒边界上不完整的Voronoi多晶体颗粒去除处理；
(2c)对去毛边处理后的Voronoi多晶模型进行网格划分，得到相应的inp文件。
4.根据权利要求1所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，所述步骤(4)有限元分析软件FEAP的Fortran程序运行平台的创建的过程为：
(4a)通过批处理方式，将所有FEAP的Fortran文件以及C语言头文件复制到编译文件夹中，为Visual Studio的编译工作做好准备；
(4b)使用Visual Studio 2010创建一个新的工程，工程模板选择Inter(R)VisualFortran-Library-Static Library，指定工程名称和路径；选择发布版本的解决方案配置和对应的Windows解决方案平台；通过添加已存在的项目，将步骤(4a)中编译文件夹所有的Fortran程序文件添加到工程项目中；选择工程配置属性，最后，选择生成解决方案，便在工程发布文件夹下生成了所有程序的源代码二进制目标文件；
(4c)使用Visual Studio 2010再创建一个新的工程，工程模板选择Inter(R)VisualFortran-QuickWin Application-Static Library，指定工程名称和路径，路径要与步骤(4b)在同一路径下；
其他步骤与步骤(4b)相同，只是在添加存在的项目时选择的是步骤(4b)中生成的工程名二进制目标文件库文件和Fortran主程序文件。
5.根据权利要求1所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，所述步骤(5)中，编写St.Venant Kirchhoff晶体模型和晶体弹塑性模型两种材料模型的UMAT子程序，是通过向FEAP软件中添加UMATLn和UMATIn两个接口程序完成，其中，n为0～9之间的整数，UMATLn是用户自己编写的，材料本构关系包含在其中，当其被单元调用时，便会接收传递进来的变形信息，完成每个单元积分点的计算并返回应力和切线刚度模量；UMATLIn完成宏命令的指派，并定义模型所需参数，从而使得可以通过输入文件实现对UMAT的调用。
6.根据权利要求5所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，对于UMATIn接口程序的架构包括：
type，用于指定宏命令名，当输入文件中使用该宏命令时，FEAP便调用对应的UMAT程序进行计算；
vv，为实型数组，当UMAT被调用时，FEAP从输入文件UMAT宏命令名后读取材料参数并存储到vv实型数组中；
d，为FEAP标准程序模块参数；
n1和n3，用于指定用户自定义材料本构程序中的历史变量数；
ud，为用户材料参数数组，通过vv数组对其赋值从而参与程序的计算，通过UMATIn结构程序，完成用户自定义UMAT的参数读取和赋值。
7.根据权利要求5所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，对于UMATLn接口子程序，其架构不同于UMATIn只做参数的定义、赋值和传递，UMATLn则是通过接收传入参数，在单元调用时进行有限元计算，然后返回应力值和切线刚度模量，也即是真正参与模拟运算的子程序。
8.根据权利要求1所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，所述步骤(6)中，通过Visual Studio 2010以及Parallel Studio XE 2011使自编写的St.Venant Kirchhoff晶体模型和晶体弹塑性模型与FEAP进行集成。
9.根据权利要求1所述的基于FEAP的三维多晶微观结构材料模型的构建方法，其特征在于，所述步骤(8)中，通过运用FEAP自身的Hill材料模型的力学响应与自编写的St.Venant Kirchhoff晶体材料模型和晶体弹塑性材料模型的力学响应进行对比，以验证所编写材料模型的用户材料子程序UMAT的有效性；以体心立方BCC结构金属DC04钢的材料属性为基础，结合上述几何建模和材料建模，对三维多晶微观结构的力学模型进行Hill材料模型、St.Venant Kirchhoff晶体材料模型和晶体弹塑性材料模型三种材料模型下的有限元模拟分析，对比三种材料模型的相通性和差异性。
CN201810765570.8A 2018-07-12 2018-07-12 基于feap的三维多晶微观结构材料模型的构建方法 Active CN109063275B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810765570.8A CN109063275B (zh) 2018-07-12 2018-07-12 基于feap的三维多晶微观结构材料模型的构建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810765570.8A CN109063275B (zh) 2018-07-12 2018-07-12 基于feap的三维多晶微观结构材料模型的构建方法

Publications (2)

Publication Number Publication Date
CN109063275A CN109063275A (zh) 2018-12-21
CN109063275B true CN109063275B (zh) 2020-07-17

Family

ID=64816166

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810765570.8A Active CN109063275B (zh) 2018-07-12 2018-07-12 基于feap的三维多晶微观结构材料模型的构建方法

Country Status (1)

Country Link
CN (1) CN109063275B (zh)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107403037B (zh) * 2017-07-04 2021-05-07 清华大学 一种开源有限元求解及优化分析方法
CN110092646B (zh) * 2019-05-30 2021-11-12 陕西理工大学 一种类晶格点阵增强相陶瓷复合材料的制备方法
CN110222442B (zh) * 2019-06-12 2020-04-21 四川大学 面心立方材料疲劳过程晶体塑性本构模型建立方法
CN111159927B (zh) * 2019-11-24 2021-09-17 浙江大学 基于体素矩阵的三维不规则形状颗粒投放的数值建模方法
CN111027244B (zh) * 2019-12-03 2024-03-12 天津大学 一种百亿级颗粒模型的构建方法
CN111028899B (zh) * 2020-01-03 2023-03-24 河南理工大学 一种建立多晶体几何模型的方法
CN111539071B (zh) * 2020-04-27 2023-06-02 武汉工程大学 一种差厚板晶体塑性本构模型建立方法、系统及电子设备
CN112001106B (zh) * 2020-08-26 2024-09-10 上海大学 一种基于有限元法的压敏电阻微观结构电热特性仿真方法
CN113358678B (zh) * 2021-05-11 2022-08-09 哈尔滨工业大学（深圳） α钛变形过程介观应力和织构的半定量预测及可视化方法
CN113221416B (zh) * 2021-05-14 2022-06-24 上海工程技术大学 一种构建颗粒增强复合材料二维微观构型的方法
CN113536623B (zh) * 2021-06-24 2022-03-25 河海大学 一种材料不确定性结构稳健性拓扑优化设计方法
CN115831271B (zh) * 2022-09-23 2023-08-08 哈尔滨工业大学 一种CuZrAl多晶模型的分子动力学模拟构建方法
CN116644646B (zh) * 2023-05-17 2024-03-22 中山大学 岩石的细观重构方法、装置、设备及存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104361167A (zh) * 2014-11-04 2015-02-18 南京航空航天大学 一种基于相场法分析的含电极铁电单晶有限元预测方法
JP2017211887A (ja) * 2016-05-26 2017-11-30 ファイフィット株式会社 有限要素法解析方法、有限要素法解析装置、解析サービスシステムおよび有限要素法解析プログラムを記録した記録媒体
CN108182335A (zh) * 2018-01-26 2018-06-19 山东科技大学 一种基于abaqus的岩石力学试验数值仿真方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104361167A (zh) * 2014-11-04 2015-02-18 南京航空航天大学 一种基于相场法分析的含电极铁电单晶有限元预测方法
JP2017211887A (ja) * 2016-05-26 2017-11-30 ファイフィット株式会社 有限要素法解析方法、有限要素法解析装置、解析サービスシステムおよび有限要素法解析プログラムを記録した記録媒体
CN108182335A (zh) * 2018-01-26 2018-06-19 山东科技大学 一种基于abaqus的岩石力学试验数值仿真方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
一种基于 Voronoi 图的多晶体有限元建模方法;郑战光等;《广西大学学报( 自然科学版)》;20160430;第41卷(第2期);第460-469页 *
多晶体材料加工的细观塑性有限元模拟;汪凯;《中国优秀硕士学位论文全文数据库 工程科技Ⅰ辑》;20140215(第2期);B022-5 *

Also Published As

Publication number Publication date
CN109063275A (zh) 2018-12-21

Similar Documents

Publication Publication Date Title
CN109063275B (zh) 基于feap的三维多晶微观结构材料模型的构建方法
Vantyghem et al. VoxelPrint: A Grasshopper plug-in for voxel-based numerical simulation of concrete printing
CN109313670B (zh) 在计算机辅助设计应用中生成晶格建议的方法和系统
CN108629147B (zh) 一种多晶体几何建模方法
US7275023B2 (en) System and method of interactively generating a family of mesh models
Nguyen et al. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems
Park et al. Structural optimization based on CAD–CAE integration and metamodeling techniques
US8831913B2 (en) Method of design optimisation
San-Vicente et al. Cubical mass-spring model design based on a tensile deformation test and nonlinear material model
CN108875138A (zh) 考虑了制造引发状态的增材制造零件的结构优化
US6947879B2 (en) Mesh generation system, design support system, analysis system, analysis method, mesh generation method, and storage medium and program transmission apparatus therefor
KR20080103504A (ko) 가상 테스트에 근거한 파라미터화 재료 및 성능 특성
Chen et al. Triangulated manifold meshing method preserving molecular surface topology
US10318675B2 (en) Post-processing system for finite element analysis
Eason et al. A parallelized multi-degrees-of-freedom cell mapping method
CN114996858A (zh) 飞行器仿真方法、装置、终端设备和存储介质
Delgado et al. Automated generation of structural solutions based on spatial designs
WO2021044891A1 (ja) 曲面のフィッティング処理方法、フィッティング処理装置およびフィッティング処理プログラム、並びに、該フィッティング処理プログラムを記憶したコンピュータ読取可能な記憶媒体
Henrich et al. DRAGen–A deep learning supported RVE generator framework for complex microstructure models
US20110191068A1 (en) Multiscale substructures in finite element analysis
Hestroffer et al. Xtalmesh toolkit: High-fidelity mesh generation of polycrystals
Panagiotopoulos et al. A group-based space-filling design of experiments algorithm
Bai et al. VCMM: A visual tool for continuum molecular modeling
Areias et al. Implicit solutions with consistent additive and multiplicative components
JP2009064164A (ja) 曲面形状作成装置、曲面形状作成方法、及び曲面形状作成プログラム

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant