晶體塑性有限元仿真入門(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
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
