黃永剛的案例
黃永剛晶體塑性UMAT及VUMAT理論及程序詳細解讀
黃永剛院士編寫的單晶晶體塑性UMAT，主要用于在Abaqus有限元軟件中進行單晶及多晶晶體塑性變形的計算。許多科研人員通過改寫其晶體塑性UMATs，建立自己的本構(gòu)子程序并發(fā)表論文。利用該本構(gòu)子程序進行材料模擬的研究方向很廣泛，包括但不限于修改材料滑移、孿晶系，修改硬化方程，加入損傷，將隱式分析的UMATs修改為顯式分析的VUMATs等，理解其基本理論并將公式與代碼對應(yīng)是改寫的關(guān)鍵。
晶體塑性的理論和UMAT自學(xué)難度較高，對于初學(xué)者自己讀懂代碼基本需要半年以上的時間，本課程旨在講解黃永剛晶體塑性UMAT的理論、公式及代碼，有助于初學(xué)者在兩周之內(nèi)熟悉和掌握晶體塑性的基本理論和子程序，加快代碼改寫進度。
本課程課件PPT長達90頁+，課程形式為一對一線上講解，時長約為3-4小時。主要面向高校研究生，需要具備張量分析的基本知識。可以根據(jù)學(xué)生的基礎(chǔ)適當(dāng)增加UMAT和VUMAT編程的講解，F(xiàn)ortran基本語法等，亦可根據(jù)需求針對性的分析UMATs修改方法。課程的內(nèi)容主要包括如下方面：
本構(gòu)模型推導(dǎo)
主要包含了應(yīng)變分解、本構(gòu)方程、硬化方程、本構(gòu)時間積分方法及雅克比矩陣等推導(dǎo)。
雅克比矩陣推導(dǎo)過程
本構(gòu)子程序逐行解讀
主要包括了UMATs基本功能、UMATs結(jié)構(gòu)、函數(shù)解讀、主程序逐行解讀等；
主程序代碼逐行展開解讀
3.UMATs改寫VUMATs方法
UMAT主要應(yīng)用于隱式分析，而對于大變形接觸問題，隱式分析往往計算效率較低。對于接觸、碰撞、沖擊等問題采用VUMAT往往具有更高的計算效率和收斂速度。該部分主要對UMATs和VUMATs的區(qū)別進行講解，介紹UMATs改寫VUMATs的要點。
展開 擴展黃永剛原始晶體塑性程序加入AF背應(yīng)力模擬金屬疲勞問題 ￥800
擴展黃永剛原始晶體塑性程序加入AF背應(yīng)力模擬金屬疲勞問題
參考文獻：《Low-cycle fatigue life prediction of a polycrystalline nickel-base superalloy using crystal plasticity modelling approach》
在原始程序中修改流動方程，加入背應(yīng)力項，引入運動硬化項，從而可以描述多晶金屬循環(huán)加載中的包辛格效應(yīng)
背應(yīng)力的演化遵循AF模型
并使用原始的PAN模型描述滑移系統(tǒng)的硬化行為
為了表征多晶的疲勞壽命，引入兩類疲勞指示因子分別為
一：累計塑性滑移
二：累計能量耗散
基于黃永剛umat的FCC多晶在簡單剪切變形下的有限元分析------案例二 ￥99
基于黃永剛umat的FCC多晶在簡單剪切變形下的有限元分析
1，使用neper軟件生成包含500個晶粒的幾何模型
2，賦予相應(yīng)的材料參數(shù)（基于腳本完成材料的批量賦值）
3，施加相應(yīng)的邊界條件（X0方向完全固定，X1方向施加Z1方向0.15mm的位移）
4，結(jié)果與后處理
5，模型文件見附錄
6，提供材料屬性賦予腳本
關(guān)于黃永剛晶體塑性的一些初學(xué)者資料整理
適用于新手入門介紹
材料參數(shù)：
狀態(tài)變量：
相關(guān)文檔：
（1）A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program, Mech Report 178
講解子程序和各個子函數(shù)作用以及整體的程序框架并包含原始的1991年代碼和inp文件（pdf版本）
（2）Addendum to 'A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program, Mech Report 178
1997年針對原始的程序中存在的個別問題進行了修改，修改后的文件以CFIXA開頭，主要是修改了硬化相關(guān)的部分
（3）umat_documentation 翻譯
主要對第一個pdf內(nèi)容進行了翻譯
（4）黃永剛程序迭代部分解釋
該pdf文檔主要對umat子程序NR迭代部分進行解釋
（5）晶體塑性硬化公式
該pdf主要介紹常見的幾類唯象的硬化模型
相關(guān)程序：
（1）原始程序可參考文檔一
（2）huang_umat_97對應(yīng)于97版本修改后的子程序
（3）考慮孿晶效應(yīng)的本構(gòu)模型參考HCP材料多晶滑移系統(tǒng)與Miller指數(shù) (qq.com)，晶體塑性每日文章推薦（六） (qq.com)
（4）位錯密度模型參考晶體塑性每日文章推薦（三） (qq.com)
（5）疲勞模型參考晶體塑性每日文章推薦（五） (qq.com)，闞前華老師的《非線性
本構(gòu)關(guān)系在ABAQUS中的實現(xiàn)》
（6）顯式模型參考晶體塑性每日文章推薦
展開 
擴展黃永剛原始晶體塑性程序加入AF背應(yīng)力模擬金屬疲勞問題
參考文獻：《Low-cycle fatigue life prediction of a polycrystalline nickel-base superalloy using crystal plasticity modelling approach》
在原始程序中修改流動方程，加入背應(yīng)力項，引入運動硬化項，從而可以描述多晶金屬循環(huán)加載中的包辛格效應(yīng)
背應(yīng)力的演化遵循AF模型
并使用原始的PAN模型描述滑移系統(tǒng)的硬化行為
為了表征多晶的疲勞壽命，引入兩類疲勞指示因子分別為
一：累計塑性滑移
二：累計能量耗散
以文獻的例，驗證修改模型的準(zhǔn)確性，其中文獻作者的幾何模型和材料參數(shù)如下
依據(jù)該模型，作者模擬得到單調(diào)拉伸以及循環(huán)加載下材料的宏觀應(yīng)力應(yīng)變響應(yīng)為
微觀響應(yīng)結(jié)果為
基于兩類疲勞指示因子，作者通過線性外推得到了基于模擬的壽命預(yù)測結(jié)果:
基于作者提供的思路和參數(shù)，對黃永剛原始程序進行修改，考慮背應(yīng)力效應(yīng)，并進行簡單的數(shù)值驗證
1，建立包含200晶粒的二維多晶模型（0.1*0.03mm），并使用四節(jié)點平面應(yīng)變單元進行網(wǎng)格劃分，如下圖
2，施加正弦形式的循環(huán)拉壓的位移載荷（1%），引力比為-1
3，模擬結(jié)果如下：
第一個滑移系統(tǒng)的背應(yīng)力：
累計塑性剪切：
累計能量耗散：
宏觀應(yīng)力應(yīng)變響應(yīng)：
展開 黃永剛晶體塑性模型耦合相場方法模擬多晶斷裂
通過和黃永剛晶體塑性模型進行耦合可以實現(xiàn)介觀尺度下，多晶材料的完整彈-塑-損傷力學(xué)行為分析，并且相比與其他損傷模型耦合方式而言，耦合相場法物理含義更加清晰，數(shù)值實現(xiàn)格式簡介，處理雅可比矩陣方便且易于收斂。因此逐漸受到介觀尺度分析材料損傷分析學(xué)者的青睞。
這里通過耦合常用的晶體塑性模型（黃-umat（修改取向到狀態(tài)變量））和斷裂相場方法，剛度和應(yīng)力退化使用二次退化函數(shù)形式。模擬包含200個晶粒的多晶模型，使用平面應(yīng)變簡化，拉伸變形30%，模擬拉伸過程中裂紋的產(chǎn)生和發(fā)展，其中斷裂總能量包含彈性變形能和塑性耗散功兩部分，模擬結(jié)果如下
初始多晶模型：
網(wǎng)格劃分（CPE4網(wǎng)格）：
相場分布（0：材料完好，1：材料完全失效）：
退化程度分布：
歐拉角（phi）分布：
可以看到耦合相場的晶體塑性模型具有潛在的預(yù)測裂紋萌生和發(fā)展的能力，其準(zhǔn)確程度取決于斷裂能參數(shù)的選擇，與更精細的實驗對比，如原位的ebsd拉伸將成為良好的校核手段。這可能成為介觀尺度下斷裂力學(xué)的應(yīng)用提供良好的參考
展開 黃永剛原始程序?qū)?yīng)的滑移系統(tǒng)順序
只考慮一組的FCC單晶滑移系統(tǒng)序號，法線方向和滑移方向分別為：
（一）(111) [110] 滑移系統(tǒng)：
1 ( 1 1 1 ) [ 0 -1 1 ]
2 ( 1 1 1 ) [ 1 0 -1 ]
3 ( 1 1 1 ) [ -1 1 0 ]
4 ( -1 1 1 ) [ 1 0 1 ]
5 ( -1 1 1 ) [ 1 1 0 ]
6 ( -1 1 1 ) [ 0 -1 1 ]
7 ( 1 -1 1 ) [ 0 1 1 ]
8 ( 1 -1 1 ) [ 1 1 0 ]
9 ( 1 -1 1 ) [ 1 0 -1 ]
10 ( 1 1 -1 ) [ 0 1 1 ]
11 ( 1 1 -1 ) [ 1 0 1 ]
12 ( 1 1 -1 ) [ -1 1 0 ]
只考慮一組的BCC單晶滑移系統(tǒng)序號，法線方向和滑移方向分別為：
（一）(110) [111] 滑移系統(tǒng)：
1 ( 0 1 1 ) [ 1 -1 1 ]
2 ( 0 1 1 ) [ 1 1 -1
展開 基于遺傳算法的晶體塑性參數(shù)自動標(biāo)定
在使用晶體塑性理論進行分析時，材料參數(shù)的標(biāo)定往往是一個枯燥繁瑣卻十分重要的工作，但由于模型考慮了滑移孿晶相變等眾多的微觀因素，造成了本構(gòu)模型包含了大量的待確定參數(shù)，目前主流的方案依然以試錯法為主，但該方案往往效率十分低下，且需要對每個參數(shù)的影響趨勢去做出準(zhǔn)確判斷，才能給出相對合理的參數(shù)更改，一些研究人員使用特定的優(yōu)化算法可以做到參數(shù)的高效標(biāo)定工作，如：蟻群算法，遺傳算法，機器學(xué)習(xí)，神經(jīng)網(wǎng)絡(luò)等，這里以黃永剛唯象的本構(gòu)模型為例，通過遺傳算法的引入，實現(xiàn)參數(shù)的自動標(biāo)定，在遺傳算法中每個設(shè)計點都被視為一個具有特定適應(yīng)度值的個體，該適應(yīng)度值基于目標(biāo)函數(shù)和約束懲罰的值。目標(biāo)函數(shù)值和懲罰值越大的個體，其適應(yīng)度值就越高。假設(shè)在模擬中待確定的材料參數(shù)為Tau_0，Tau_s，H_0，并通過黃永剛初始的材料參數(shù)Tau_0=60.9，Tau_s=109.5，H_0=540.5得到初始的拉伸曲線作為目標(biāo)函數(shù)，并給定參數(shù)對應(yīng)的區(qū)間，Tau_0【30,80】,Tau_s【100,150】,H_0【200,1000】作為待定函數(shù)的區(qū)間，給定初始測試值為Tau_0=50，Tau_s=125，H_0=350，作為初始試探值提供給遺傳算法作為初始值，將遺傳算法得到的不同參數(shù)值對應(yīng)的力-位移曲線和原始黃永剛參數(shù)的力-位移曲線的標(biāo)準(zhǔn)差作為目標(biāo)函數(shù)對參數(shù)進行優(yōu)化。優(yōu)化效果如下圖示：
在使用遺傳算法進行22次的嘗試過程中，遺傳算法給出的參數(shù)以及對應(yīng)目標(biāo)函數(shù)的值為
可以看到參數(shù)均落在了給定的初始區(qū)間中，隨機迭代次數(shù)的增加，對應(yīng)的目標(biāo)函數(shù)逐漸下降。
展開 適合HCP mg合金的黃永剛Vumat的晶體塑性（之前推文對應(yīng)的vumat版本，同時附上對應(yīng)的材料參數(shù)）
subroutine vumat(nblock, ndir, nshr, nstatev, nfieldv, nprops,
* lanneal, steptime, totaltime, dt, cmname, coordmp, charlength,
* props, density, straininc, relspininc, tempold, stretchold,
* defgradold, fieldold, stressold, stateold, enerinternold,
* enerinelasold, tempnew, stretchnew, defgradnew, fieldnew,
* stressnew, statenew, enerinternnew, enerinelasnew)
c
include 'vaba_param.inc'
c modified from harewood vumat- jg:20/07/12
dimension stressold(nblock,ndir+nshr),
1 stressnew(nblock,ndir+nshr),
1 stateold(nblock,nstatev),statenew(nblock,nstatev),
1 straininc(nblock, ndir+nshr)
dimension slpdir(3,18),slpnor(3,18),slpdef(6,18),slpspn(3,18),
1 dspdir(3,18),dspnor(3,18),d(6,6),fslip(18),dfdxsp(18),
1 ddemsd
展開 huang隱式程序修改為顯式及計算案例
黃永剛原始晶體塑性具有良好的收斂性，以及高效的計算效率，在一般變形下無需修改，即可直接使用。然而一些特殊的工況，如切削，軋制，沖壓等隱式存在收斂性問題。因此通常使用顯示程序進行計算。但從頭完成顯式晶體塑性構(gòu)造對于一般學(xué)者顯然難度過高，一個簡單的想法就是直接將現(xiàn)成的黃永剛隱式程序改成顯式。abaqus里這是可以實現(xiàn)的。其基本的步驟是：
1，加入vumat接口程序（見附錄abaqus官網(wǎng)有）
2，對nblock進行循環(huán)，計算應(yīng)力和狀態(tài)變量
3，更新應(yīng)力與狀態(tài)變量，重復(fù)計算直到增量結(jié)束。
值得注意的是，umat與vumat程序里面剪應(yīng)力分量定義順序與應(yīng)力不同
umat：12，13，23（工程剪應(yīng)變）
vumat：12，23，13（2*工程剪應(yīng)變）
同時采用該方法計算時計算效率顯著高于完全顯式，并允許較大的時間增量。為評估模型計算效率，采用1000個晶粒80000個單元的二維模型進行20%的壓縮模擬。耗時3小時，計算結(jié)果與隱式結(jié)果類似。
展開 晶體塑性每日文章推薦（七）
然而，晶粒取向?qū)暧^應(yīng)力-應(yīng)變曲線的影響很小
基于作者提供的思路和參數(shù)，對黃永剛原始程序進行修改，考慮背應(yīng)力效應(yīng)，并進行簡單的數(shù)值驗證
1，建立包含200晶粒的二維多晶模型（0.1*0.03mm），并使用四節(jié)點平面應(yīng)變單元進行網(wǎng)格劃分，如下圖
2，施加正弦形式的循環(huán)拉壓的位移載荷（1%），引力比為-1
3，模擬結(jié)果如下：
第一個滑移系統(tǒng)的背應(yīng)力：
累計塑性剪切：
累計能量耗散：
宏觀應(yīng)力應(yīng)變響應(yīng)：
整體上于作者的模擬情況比較接近，做疲勞分析的同學(xué)可以進行類似嘗試
精品abaqus學(xué)習(xí)教材分享
后續(xù)計劃上傳黃永剛程序的一般使用案例。
基于密西西比州立大學(xué)晶體塑性模型預(yù)測不同變形下織構(gòu)演化案例教學(xué)------案例八 ￥99
基于密西西比州立大學(xué)晶體塑性模型預(yù)測不同變形下織構(gòu)演化
官方使用原始案例
案例一，單向壓縮75%（FCC）
加載條件
織構(gòu)演化結(jié)果
1， 案例二，單向拉伸75%（FCC）
加載條件
織構(gòu)演化結(jié)果
密西西比州立大學(xué)晶體塑性有限元代碼和黃永剛院士的程序一樣，均是開源代碼，可免費獲得，并且同時可以考慮FCC，BCC，HCP的滑移和孿晶變形，有著廣泛的應(yīng)用，目前該代碼已經(jīng)集成到FEPX計算軟件中，支持并行運算，計算效率很高。本案例采用該代碼，研究FCC，BCC兩種結(jié)構(gòu)在單向拉伸，壓縮，平面應(yīng)變壓縮等75%的變形量下織構(gòu)的演變（需要注意的是，這個代碼的輸出使用的Kocks輸出表示取向，為了使用方便，已經(jīng)在程序中修改bunge標(biāo)號，可用MTEX直接繪制極圖）
使用包含500個隨機取向的單元預(yù)測取向演化
初始隨機取向
一：FCC織構(gòu)演化
單向壓縮75%的取向分布
單向拉伸75%的取向分布
平面應(yīng)變壓縮75%取向分布
一：BCC織構(gòu)演化
單向壓縮75%取向分布
單向拉伸75%取向分布
平面應(yīng)變壓縮75%取向分布
其中FCC和壓縮和拉伸與官網(wǎng)所提供的案例保持一致，F(xiàn)CC，BCC的平面應(yīng)變壓縮與已有文獻的典型織構(gòu)一致，完全正確。
展開 在huang子程序中進行晶體取向更新
由于網(wǎng)上更新取向的通常只是公式的簡單介紹，這里嘗試結(jié)合lingzhi的鎂合金更新取向的方案對黃永剛的原始程序進行修改，并以狀態(tài)變量的形式保存。將三個歐拉角分別存儲為SDV123，SDV124，SDV125。
模擬FCC沿著ND方向壓縮50%之后的取向分布，并與之前介紹的基于matlab腳本生成歐拉角的方法對比。
輸入的初始極圖
基于Matlab腳本實現(xiàn)的極圖
寫入fortran狀態(tài)變量對應(yīng)的極圖
可以看到兩者幾乎一致，因此更推薦在Fortran程序內(nèi)部進行取向更新
這里展示部分修改的Fortran程序，大家可以嘗試參考lingzhi的程序在黃程序的基礎(chǔ)上進行取向更新的加入
優(yōu)先建議大家自行嘗試修改加入，操作出錯無法解決可以取技術(shù)鄰找到我修改后的Fortran程序，鏈接如下：
https://www.yqgqt.org.cn/post/1936806
展開 基于abaqus的Huang晶體塑性UMAT改VUMAT
黃永剛院士編寫的單晶晶體塑性UMAT，主要用于在Abaqus有限元仿真中進行單晶及多晶晶體塑性變形的計算，是許多科研工作者學(xué)習(xí)晶體塑性模擬的教學(xué)資源。可以在其基礎(chǔ)上對硬化模型進行修改，甚至引入損傷。
UMAT主要應(yīng)用于隱式分析，而對于大變形接觸問題，隱式分析往往計算效率較低。對于接觸、碰撞、沖擊等問題采用VUMAT往往具有更高的計算效率和收斂速度。本文旨在將Huang編寫的UMAT改寫為VUMAT，并進行對比驗證。
將UMAT改寫為VUMAT需要從以下方面考慮：（1）UMAT是在積分點上調(diào)用的，而VUMAT一次調(diào)用會計算很多個積分點上的變量，需要對子程序接口形式進行修改；（2）Huang本構(gòu)中的轉(zhuǎn)動張量DROT是用于對滑移面和滑移方向進行旋轉(zhuǎn)的，在UMAT中，Abaqus會提供轉(zhuǎn)動張量DROT，在VUMAT中，子程序接口沒有提供DROT，需要通過VUMAT傳入的變形梯度更新滑移面和滑移方向；（3）顯示分析采用了Green-Naghdi率，而隱式分析采用Jaumann率，需要對應(yīng)變率進行修改。
采用了兩個多晶模型進行一致性的驗證，第一個模型是125個網(wǎng)格的單位長度代表體積單元，每25個網(wǎng)格設(shè)置1個取向。第二個模型是采用Voronoi方法獲得的15個不同取向晶粒的多晶模型。
（1）15個不同取向晶粒的多晶模型
15個不同取向晶粒的多晶模型，采用狗骨單軸拉伸試件進行數(shù)值試驗，有限元模型如下圖所示。開展單軸拉伸，UMAT采用隱式分析，VUMAT采用顯式分析。
展開 