晶体塑性有限元仿真入门(3)--开源代码平台EVOCD
晶体塑性有限元仿真入门(3)--开源代码平台EVOCD
晶体塑性有限元初学者较为熟知的两个工具Huang's UMAT以及DAMASK平台,这篇文章介绍另外一个晶体塑性有限元方法(CPFEM)的开源代码平台EVOCD,讲解如何使用这些开源代码进行材料的塑性变形模拟以及模拟变形过程中晶体取向的变化(织构)。
图1 EVOCD的CPFEM流程图
(E.B. Marin, Center for Advanced Vehicular Systems, Mississippi State University)
我们在网上搜索晶体塑性的关键字''CPFEM''时,会发现搜索引擎的网页排名第一是马普所(MPI, 大名鼎鼎的DAMASK就是他们团队的成果)的研究成果,其次是密西西比州立大学先进车辆系统中心(Center for Advanced Vehicular Systems, Mississippi State University)的开源代码平台EVOCD,第三是基于Huang的晶体塑性有限元方法,由此可见EVOCD在晶体塑性有限元方法中的重要性。
图2 CPFEM搜索结果
(从上到下分别是马普所 (dierk-raabe.com) 、密西西比州立大学 (msstate.edu) 、哈佛大学 (columbia.edu) 的相关研究成果)
国内的晶体塑性有限元初学者,最主要的还是使用Huang's UMAT以及DAMASK平台,而对密西西比州立大学的开源代码平台EVOCD不太常用。这篇文章将讲解该平台的使用方法以及如何使用该平台进行晶体塑性有限元变形模拟。
YouTube有相关教学视频:https://youtu.be/UPbUiaj8us8
B站也有搬运相关教程:https://www.bilibili.com/video/BV1SK4y1d74z
图3 EVOCD平台视频教学
本文将参考EVOCD平台的官网教学以及以上的视频教学,全文包括以下几个部分:
1) 开源代码平台EVOCD的介绍
2) 子程序与ABAQUS进行关联
3) 子程序材料文件下载
4) 子程序材料文件修改
5) 输入文件
6) 输出文件
7) 参考资料
1. 开源代码平台EVOCD的介绍
如下图所示,EVOCD平台里简介了晶体塑性有限元法,以及如何使用该平台进行简单的模拟,这些文件都能在平台里免费开源下载。
图4 EVOCD平台以及教程截图
根据教学文档的参考文献,我们可以找到许多Marin博士已发表的理论方面的文献。
图5 EVOCD平台教程相关的参考文献
https://icme.hpc.msstate.edu/mediawiki/images/1/10/Introduction_to_CPFEM_manual-1.pdf
虽然这些教程都是免费开源的,我曾经也推荐过这个平台,但对于初学者来说,难点在于如何跟着教程快速上手案例,并在此基础上进行改进,实现较为复杂模型的构建。
图6 EVOCD平台教程的模拟结果
https://icme.hpc.msstate.edu/mediawiki/images/8/8e/CPFEM_Simulation_of_Aluminum_V2.pdf
下面,我将按照教程文件的顺序,一步步对晶体塑性的变形进行教学。
2. 子程序与ABAQUS进行关联
ABAQUS的材料库里没有晶体塑性材料模型,因此需要使用用户材料子程序(User-Defined Material, UMAT),而所有编写好的UMAT能够成功使用的前提是ABAQUS与UMAT子程序进行关联。因此,第一步是安装子程序的相关平台,并与ABAQUS进行关联。
图7 EVOCD平台教程的第一步:子程序与ABAQUS进行关联
ABAQUS要是想调用子程序,需要以下软件进行关联:Microsoft Visual Studio、Intel Visual Fortran Composer。需要注意安装顺序是先VS,再IVF,具体安装过程可参考下列B站视频或百度文库。
B站视频教学:https://www.bilibili.com/video/av584237026/
百度文库教学:https://wenku.baidu.com/view/3ba92eefb9d528ea80c7794b
图8 子程序与ABAQUS进行关联步骤
3. 子程序材料文件下载
安装好子程序的相关平台,并与ABAQUS进行关联后,第二步是下载EVOCD平台的子程序材料文件,需下载EVOCD平台提供的下列6个子程序材料文件才能进行晶体塑性的模拟。
图9 EVOCD平台教程的第二步:子程序材料文件下载
注1:如果单击网站的文件后不会自动下载文件,只会打开文件,需复制粘贴所有内容到文本编辑器里,然后再另存为相应格式的文件。
图10 EVOCD平台支撑文件下载(FCC晶格材料)
https://icme.hpc.msstate.edu/mediawiki/index.php/Code:_ABAQUS_CPFEM.html
注2:网页的文字显示有左图和右图两种格式,默认是左图,fcc.sx、test.xtali、params_xtal.inc应该使用显示源代码调整为右图显示效果再进行复制粘贴。
图11 EVOCD平台支撑文件显示(fcc.sx)
4. 子程序材料文件修改
这六个文件的功能分别如下:
(1) umat_xtal.f
constitutive model - polycrystal average model
UMAT subroutine for crystal plasticity - SNL & MSU
umat_xtal.f文件下载后需根据目标目录对路径进行修改:
data filePath
& /'/cavs/cmd/data1/users/qma/abaqus_xtalplas/oneelement/'/
图12 umat_xtal.f文件修改目录位置
(2) texture.txti
initial orientation distribution
texture.txti文件包括了许多欧拉角的初始取向,以案例中的数据为例,各数据的含义是:The number 500 means that there are 500 orientations in this texture file. The last line number 555 is the seed number it should be larger than the total orientation number (500 in this case). The three columns of data are the corresponding three Euler angles in Kocks convention
(3) fcc.sx
single crystal parameters
fcc.sx is another key input file. It includes all the deformation modes and the corresponding hardening parameters. When you change these parameters you will get different output results. The parameters of these deformation modes can be obtained in the published literature or by calibration using experimental data.
(4) test.xtali
control for the time step and deformation
This input file controls which structure you will use and how many orientations you will use in your simulation.
(5) params_xtal.inc
number of slip systems
The params_xtal.inc input file defines some control parameters. Usually, users can not change this input file. But, users should change the parameter like below when you simulate 2D or 3D problems.
parameter (NTENS=6, NDI=3, NSHR=3) // this parameter for 3D structure
c parameter (NTENS=4, NDI=3, NSHR=1) // this parameter for 2D structure
(6) numbers.inc
numerical constants
This file lists a number of numerical constants used in the UMAT. This file does not need to be changed.
以上复制了教程里对各参数含义的解释,这部分的关键是umat_xtal.f文件下载后需根据目标目录对路径进行修改,否则运算会报错。
5. 输入文件
如下所示,建立一个单元的变形模型,提交Job生成inp文件。
图13 oneelement变形模型
生成inp之后,结合之前保存的umat_xtal.f文件作为材料子程序,提交仿真计算。下图红框为所有仿真计算的输出结果,其中两个重要输出文件为tension.odb和texture.txto,tension.odb保存了所有变形场量随时间变化的结果,texture.txto保存了所有取向随时间变化的结果。
图14 oneelement计算结果
6. 输出文件
(1) tension.odb
使用ABAQUS打开ODB文件,如下图所示查看应力应变结果并导出数据。
图15 oneelement计算结果 (应力和应变)
将导出的应力应变进行结合,绘制应力-应变曲线:
图16 oneelement计算的应力S33-应变L33曲线
下图左边是EVOCD平台拉伸教程的模拟结果,右边是本教程的模拟结果,可以发现二者的结果一致。
图6 EVOCD平台拉伸教程的模拟结果对比
(2) texture.txto
texture.txto文件由以下两部分构成:
一是变形前的初始取向部分:
*----- Initial Assigned Orientations -----*
二是变形中的变形取向部分
*----- Euler Angles at incr 20 -----*
... ...
*----- Euler Angles at incr 180 -----*
如下所示,每部分都记录了所有晶粒的取向(三个欧拉角)随增量步变化的变化情况。
图17 texture.txto记录了所有晶粒取向的变化
texture.txto文件记录了变形过程中所有晶粒取向变化的原始数据,如需以极图或者其他的方式呈现取向变化的规律,需将取向数据进行转换,通常是使用MTEX软件,以下是将欧拉角转换为极图的结果。
图18 变形前的极图和反极图
图19 变形后的极图和反极图
MTEX画图代码:
clc;close all;clear;
Euler = load('texture_euler.txt');
cs = crystalSymmetry('cubic'); %晶体对称性
ss = specimenSymmetry('triclinic'); %试样对称性
ori = orientation.byEuler(Euler(:,1)*degree,Euler(:,2)*degree,Euler(:,3)*degree,cs,ss);
plotPDF(ori,Miller({1,0,0},{1,1,0},{1,1,1},cs),'all') %绘制极图
plotIPDF(ori,[vector3d.X,vector3d.Y,vector3d.Z],'all')%绘制反极图
下图是EVOCD平台拉伸教程的模拟结果,与本教程的结果一致。
图6 EVOCD平台拉伸教程的模拟结果
https://icme.hpc.msstate.edu/mediawiki/images/8/8e/CPFEM_Simulation_of_Aluminum_V2.pdf
7. 参考资料
https://icme.hpc.msstate.edu/mediawiki/images/1/10/Introduction_to_CPFEM_manual-1.pdf
https://icme.hpc.msstate.edu/mediawiki/images/8/8e/CPFEM_Simulation_of_Aluminum_V2.pdf
B站视频教学:https://www.bilibili.com/video/av584237026/
百度文库教学:https://wenku.baidu.com/view/3ba92eefb9d528ea80c7794b
查看更多评论 >