晶體塑性理論的案例
黃永剛晶體塑性UMAT及VUMAT理論及程序詳細解讀
黃永剛院士編寫的單晶晶體塑性UMAT，主要用于在Abaqus有限元軟件中進行單晶及多晶晶體塑性變形的計算。許多科研人員通過改寫其晶體塑性UMATs，建立自己的本構子程序并發表論文。利用該本構子程序進行材料模擬的研究方向很廣泛，包括但不限于修改材料滑移、孿晶系，修改硬化方程，加入損傷，將隱式分析的UMATs修改為顯式分析的VUMATs等，理解其基本理論并將公式與代碼對應是改寫的關鍵。
晶體塑性的理論和UMAT自學難度較高，對于初學者自己讀懂代碼基本需要半年以上的時間，本課程旨在講解黃永剛晶體塑性UMAT的理論、公式及代碼，有助于初學者在兩周之內熟悉和掌握晶體塑性的基本理論和子程序，加快代碼改寫進度。
本課程課件PPT長達90頁+，課程形式為一對一線上講解，時長約為3-4小時。主要面向高校研究生，需要具備張量分析的基本知識。可以根據學生的基礎適當增加UMAT和VUMAT編程的講解，Fortran基本語法等，亦可根據需求針對性的分析UMATs修改方法。課程的內容主要包括如下方面：
本構模型推導
主要包含了應變分解、本構方程、硬化方程、本構時間積分方法及雅克比矩陣等推導。
雅克比矩陣推導過程
本構子程序逐行解讀
主要包括了UMATs基本功能、UMATs結構、函數解讀、主程序逐行解讀等；
主程序代碼逐行展開解讀
3.UMATs改寫VUMATs方法
UMAT主要應用于隱式分析，而對于大變形接觸問題，隱式分析往往計算效率較低。對于接觸、碰撞、沖擊等問題采用VUMAT往往具有更高的計算效率和收斂速度。該部分主要對UMATs和VUMATs的區別進行講解，介紹UMATs改寫VUMATs的要點。
展開 基于遺傳算法的晶體塑性參數自動標定
在使用晶體塑性理論進行分析時，材料參數的標定往往是一個枯燥繁瑣卻十分重要的工作，但由于模型考慮了滑移孿晶相變等眾多的微觀因素，造成了本構模型包含了大量的待確定參數，目前主流的方案依然以試錯法為主，但該方案往往效率十分低下，且需要對每個參數的影響趨勢去做出準確判斷，才能給出相對合理的參數更改，一些研究人員使用特定的優化算法可以做到參數的高效標定工作，如：蟻群算法，遺傳算法，機器學習，神經網絡等，這里以黃永剛唯象的本構模型為例，通過遺傳算法的引入，實現參數的自動標定，在遺傳算法中每個設計點都被視為一個具有特定適應度值的個體，該適應度值基于目標函數和約束懲罰的值。目標函數值和懲罰值越大的個體，其適應度值就越高。假設在模擬中待確定的材料參數為Tau_0，Tau_s，H_0，并通過黃永剛初始的材料參數Tau_0=60.9，Tau_s=109.5，H_0=540.5得到初始的拉伸曲線作為目標函數，并給定參數對應的區間，Tau_0【30,80】,Tau_s【100,150】,H_0【200,1000】作為待定函數的區間，給定初始測試值為Tau_0=50，Tau_s=125，H_0=350，作為初始試探值提供給遺傳算法作為初始值，將遺傳算法得到的不同參數值對應的力-位移曲線和原始黃永剛參數的力-位移曲線的標準差作為目標函數對參數進行優化。優化效果如下圖示：
在使用遺傳算法進行22次的嘗試過程中，遺傳算法給出的參數以及對應目標函數的值為
可以看到參數均落在了給定的初始區間中，隨機迭代次數的增加，對應的目標函數逐漸下降。
展開 塑性工程學報：Custom450鋼拉伸的晶體塑性有限元分析
而研究微細觀尺度的變形不均勻性是新材料開發及優選的重要準則，晶體塑性有限元方法將晶體塑性理論和有限元軟件進行了恰當的融合，成為研究細觀層次塑性變形行為的強有力工具。來自華東理工大學機械與動力工程學院的艾鑫團隊，基于Voronoi方法建立了Custom450 鋼拉伸的二維晶體塑性模型，分析了初始硬化模量、參考剪切應變率、應變率敏感系數、初始屈服應力以及飽和流動應力對材料應力——應變曲線的影響，并對晶體塑性參數進行了標定。
在文獻中，作者所建立的單晶本構模型參考了HUANGY的單晶體模型的子程序UMAT，此率相關硬化晶體塑性模型需要確定的參數包括初始硬化模量h0、初始屈服應力τ0、參考剪切應變率γ，應變率敏感系數n和飽和流動應力τs，其他參數通過計算和查找文獻獲得。基于Voronoi方法，作者在有限元軟件Abaqus中建立了Custom450材料的多晶體二維幾何模型并將本構關系嵌入軟件中，進行拉伸過程的模擬。
圖1所示是微結構模型及其網格劃分，幾何模型尺寸長度為0. 2 mm，寬度為0. 5mm，共包含100個晶粒，大小和形狀隨機，且晶粒取向隨機分布。
圖1包含100個晶粒的微結構模型及其網格劃分
圖2是邊界條件的約束情況，模型的上端面和下端面的所有節點在y方向上具有均勻的位移，左側所有節點在x方向上設置約束，使其不能橫向移動，y方向自由，在右邊界施加載荷，右側的所有節點x方向上經受同等應變載荷，而在y方向上是自由的。
圖2邊界條件示意圖
對于體心立方晶體來說，3個滑移系包括1個主滑移系和2個次滑移系。分別對包含1、2、3組滑移系開動的情形進行模擬，結果如圖3所示，只有主滑移系 { 110} < 111 >啟動時，應力——應變曲線在彈塑性區間過渡的位置存在明顯拐點，并與試驗曲線吻合良好。
展開 基于晶體塑性有限元（CPFEM）的鈦合金圓棒拉伸過程模擬
作者：辭殤
關鍵詞：CPFEM；鈦合金；單軸拉伸；織構極圖；孿晶
晶體塑性有限元是一種結合了晶體塑性理論和有限元方法的數值模擬技術?。這種方法考慮了晶體材料的各向異性、滑移系統的開動和相互作用、以及變形過程中的硬化效應。它主要用于分析和預測晶體材料的塑性變形行為，特別是在微觀尺度上的變形機制。
晶體塑性有限元在材料科學和工程領域有著廣泛的應用，特別是在金屬加工、航空航天、汽車制造和生物醫學等領域。通過這種技術，研究人員和工程師可以更好地理解材料的力學行為，從而開發出更輕、更強、更耐用的材料和產品。此外，晶體塑性有限元仿真還能夠考慮材料的微觀結構特征，如晶粒取向、晶界、相分布以及滑移系統的活動，從而能夠預測材料在細觀尺度上的織構演化。
利用CPFEM方法對鈦合金圓棒拉伸過程進行模擬，使用UMAT子程序以及Abaqus有限元軟件作為晶體塑性有限元分析的實現方式。并且，在一些復雜工藝條件下如切削、軋制、沖壓等，CPFEM方法同樣適用，能夠模擬材料變形過程中的非線性行為和動態響應。
在晶體塑性有限元中，首先在Abaqus中建立了單軸拉伸有限元模型如圖1所示，材料被建模為包含大量晶粒的集合體如圖2所示，每個晶粒都有其特定的晶體取向，并且每個晶粒的變形過程均考慮了滑移和孿晶的變形機制。
圖1 單軸拉伸有限元模型示意圖
圖2 單軸拉伸晶體塑性模型示意圖
通過有限元方法，可以計算出在給定拉伸載荷下，這些晶粒如何相互作用，以及它們如何隨時間變形。這種方法能夠提供關于晶體材料內部應力、應變和變形機制的詳細信息，有助于理解材料在受力時的響應，并優化材料的設計和加工過程。圖3所示為單軸拉伸過程應力云圖，圖4所示為單軸拉伸過程孿晶云圖。
展開 
晶體塑性耦合連續損傷本構框架
經典文章推薦
《Necking behavior of AA 6022-T4 based on the crystal
plasticity and damage models
是最經典的耦合晶體塑性理論和連續損傷的文章之一，損傷力學有兩種主要方法。第一種是Gurson提出的基于微觀力學的損傷模型。在基于微觀力學的方法中，損傷演化通過孔隙成核、生長和聚結來描述。對空穴成核和生長進行了建模，必須使用實驗數據確定相關系數。另一種方法是連續損傷力學（CDM）。在CDM框架中，使用應力、壓力、溫度和應力三軸性確定斷裂應變。在這些研究之后，提出了許多改進的模型，以包括洛德角和各向異性損傷的影響，
作者在研究中使用的損傷模型基于連續損傷力學（CDM）。然而，通過結合CPFEM可以預測孔隙的萌生、生長和聚結行為。此外，材料因損傷而弱化用于描述頸縮后承載能力的突然下降，通過顯式時間積分方案進行了分析，這為通過CPFEM預測頸縮行為提供了可能性。然而，沒有預測頸縮形狀和載荷位移曲線。為了準確預測頸縮和載荷位移曲線，使用隱式時間積分方案進行了分析，可以獲得更合理的載荷位移曲線。此外，還進行了實驗，并與分析結果進行了比較。最后，新提出了四種不同的帶系數校準的損傷模型，并提出了一種最能描述頸縮行為的模型。
作者使用的四類連續損傷模型理論如下
（1）最大塑性應變損傷模型：該模型將損傷定義為當主塑性應變大于某一臨界值時開始和累積的損傷。此模型寫為：
ε1f.ini是損傷萌生塑性應變值，ε1f.ini是最大塑性應變值，D是損傷因子，M是損傷指數（通常取值大于1.0有利于流動應力平滑過渡）
（2）等效塑性應變損傷模型：該模型將損傷定義為當等效塑性應變大于某一臨界值時開始和累積的損傷。
展開 金屬學報：孿生誘發軟化與強化效應的Cu晶體塑性行為模擬
眾所周知，位錯滑移和孿生是主導多晶體材料塑性行為的主要變形機制。一方面，在孿生主導塑性條件下，孿晶激活演化過程中應力-應變曲線存在明顯的應力突降現象，即孿生軟化效應；另一方面，孿晶阻礙位錯運動使得晶體材料在塑性變形過程中表現出強化現象。為了能夠精確描述孿晶激活演化及其與位錯交互作用對宏觀塑性行為的影響，來自于天津理工大學的郭祥如和申俊杰兩人基于晶體塑性理論建立描述孿晶形核、增殖和長大的位錯密度基晶體塑性本構模型，揭示了不同晶體取向Cu單晶拉伸變形過程中位錯滑移、孿生激活及其交互作用下的宏觀塑性行為演化規律，進一步分析了Cu多晶拉伸變形過程中晶粒間交互作用對孿生軟化、應變硬化等宏觀塑性行為的影響。
為了應用該模型準確模擬材料的宏觀力學響應，必須確定該模型相關材料參數。作者結合fcc晶體材料滑移系和孿生系的晶體學特征，根據前人對Cu的研究結果，最終得出晶體塑性模型Cu單晶材料參數。建立如圖1所示的Cu單晶CPFE模型。
圖1 Cu單晶拉伸過程的晶體塑性有限元模型示意圖
為了驗證上述CPFE模型的可靠性，圖2給出了Cu單晶沿不同取向拉伸過程的力學響應模擬及實驗結果的對比情況。可以看出，模擬結果與實驗結果吻合良好。在Cu單晶沿[541]取向加載條件下，其應力-應變曲線分成明顯的3個階段，即滑移階段A、孿生階段B及位錯與孿晶交互作用階段C。為了深入揭示Cu單晶塑性變形過程中各滑移系和孿生系激活演化行為及孿晶對位錯滑移的影響，圖3給出了[541]取向下Cu單晶拉伸變形過程中各滑移系和孿生系激活演化結果。
展開 一份FEPX晶體塑性軟件使用介紹的文檔
no.ntnu_inspera_140385465_3592470.pdf
文檔介紹了
理論部分
1，極圖相關理論
2，晶體塑性有限元的理論框架
neper軟件的四個模塊
1，Generation-Module
2， Meshing-Module
3，Simulation-Module
4， Visualization-Module
fepx軟件的四個模塊
1， Code Structure
2， Running FEpX
3，Input Files
4，Output Files
并給出了linux上的fepx和neper的安裝和驗證教程
比較了fepx模擬和abaqus模擬的差異性，介紹了fepx的優缺點
并附上了邊界條件設置的腳本。對于相了解fepx使用的同學是一份難得的教程資源，可以結合官網幫助文檔一起使用。
如果本文提供的資源對您有所幫助，歡迎點贊訂閱打賞！
部分文檔內容截圖
neper安裝
Fepx安裝
文獻提供的軋制變形的案例
邊界腳本：
文檔放在技術鄰上，有興趣去搜即可
展開 晶體塑性有限元仿真入門(5)—歐拉角與晶體取向
晶體塑性有限元仿真入門(5)—歐拉角與晶體取向.pdf
晶體塑性有限元仿真入門(5)—歐拉角與晶體取向
備注：網頁排版有亂碼，建議下載附件pdf查看
晶體取向是材料學科中的重要分支，當晶粒發生擇優取向時，則導致材料性能（力學，物理和化學性能）的各向異性。各向異性會造成材料實際應用中的各種問題，如鋁合金典型的制耳現象，再如取向硅鋼中存在Goss織構時，有利于其磁學性能。在基礎研究領域，織構的形成與演變是基本的科學問題。在工業應用領域，通過織構的設計和控制可以提高材料的性能。隨著近年來EBSD和XRD等表征技術的發展，各種SCI期刊的發文都已離不開對材料晶體學取向的分析。這篇文章介紹晶體塑性有限元仿真過程中的歐拉角與晶體取向。
圖1 塑性變形過程導致的材料各向異性
全文包括以下幾個部分：
1) 材料晶體結構
2) EBSD工作原理
3) 晶體取向分析
4) 晶體塑性材料模型
5) 織構演變結果
6) 參考資料
7) 附錄
材料晶體結構
在晶體學中，晶體結構是對晶體材料中原子、離子或分子有序排列的描述。有序結構由組成粒子的內在性質產生，形成沿物質三維空間的主要方向重復的對稱模式，如圖2所示。
圖2 高分辨率透射電子顯微鏡圖片的鐵晶體，完美單晶的二維示意圖
構成這種重復圖案的材料中最小的一組粒子是結構的晶胞。晶胞完全反映了整個晶體的對稱性和結構，這是通過晶胞沿其主軸重復平移而建立的。平移向量定義布拉維點陣的節點，不同的晶體內部原子排列稱為具有不同的晶格結構。各種晶格結構可以歸納為七大晶系，各種晶系分別與十四種空間格（稱為Bravais晶格）相對應，如圖3所示。
展開 熱-彈-黏塑性晶體塑性模型文章推薦
傳統室溫本構模型通常需要依賴大量不同溫度、不同加載路徑下的實驗數據進行擬合，很難真正解釋“溫度如何影響晶體滑移和多晶塑性響應”。
Cyr 等人針對這一問題提出了一個三維熱-彈-黏塑性晶體塑性模型，即 TEV 模型，用于描述 FCC 多晶材料，特別是 AA5754 鋁合金在升溫條件下的力學行為。該模型的核心思想是：材料變形不僅包含彈性變形和晶體塑性滑移，還需要顯式考慮熱膨脹變形。因此，總變形梯度被分解為彈性/剛體轉動部分、熱變形部分和塑性變形部分。
在本構層面，作者保留了 FCC 晶體的 12 個 {111}<110> 滑移系，并采用冪律型滑移率方程描述率相關塑性流動。與常規晶體塑性模型不同的是，該模型把溫度效應系統地引入到多個關鍵物理量中：首先，單晶彈性常數 C11、C12、C44 隨溫度變化；其次，滑移阻力引入熱軟化函數，用來描述溫度升高后滑移更容易發生的現象；再次，單滑移硬化參數也被寫成溫度函數，包括參考臨界分切應力、初始硬化率和硬化指數。
這個模型的優勢在于，它不是簡單地給宏觀應力-應變曲線加一個溫度修正系數，而是從晶體滑移層面描述溫度對材料響應的影響。換句話說，它可以同時分析宏觀應力變化、微觀滑移活動、織構演化、局部應變集中和熱軟化機制。因此，它比普通經驗型熱塑性模型更適合用于多晶材料溫成形模擬。
作者首先利用 AA5754 鋁合金在 25 ℃、148 ℃、204 ℃ 和 232 ℃ 下的單軸拉伸實驗數據標定溫度相關硬化參數。隨后，又預測了 177 ℃ 和 260 ℃ 下的拉伸響應。
展開 初識彈塑性與晶體塑性——面向新手 ￥4.9
文檔以ppt形式展示，從簡單走向深入，講解了彈塑性與晶體塑性在實現時的基本邏輯，內容展示：
材料的理論斷裂強度 附晶體材料強度與斷裂微觀理論下載
但這只是實驗觀察到的現象，它與材料的理論斷裂值還有很大的區別。
假設材料的斷裂是由于原子間距被拉的太遠，超過了極限從而發生的斷裂。我們知道，原子之間的力與原子間的距離存在一定的關系，當原子靠的特別近的時候，原子間存在排斥力，當原子離的比較遠的時候，原子間存在相互吸引力，在某一距離下，原子間的作用力為0，即平衡位置。
現在我們來考慮原子間的力與應力的關系，根據應力的定義
顯然，曲線上的最大值σm即代表原子間的最大結合力——理論斷裂強度，即在理論上認為材料應力超過σm時將被拉斷。作為一級近似，該曲線可用正弦曲線表示。
而實際上，對于純鐵的抗拉強度是只有170～270MPa左右，我們熟知的Q235鋼，其抗拉極限為375～460MPa，Q345鋼的抗拉強度約是490-620MPa，遠遠低于材料的理論斷裂強度。主要原因在于公式(11)表示的是理想材料的斷裂強度，也就是說材料中沒有任何的缺陷。但這是不可能的，材料在冶金、鑄造、加工等過程中難免會產生一些初始缺陷，造成應力集中從而大大降低了材料的強度缺陷。
下載地址：晶體材料強度與斷裂微觀理論
展開 
JMPS：多主元合金塑性和應變硬化的分層多尺度晶體塑性框架
在細觀尺度上，晶體塑性有限元(CPFE)方法可以考慮相變、位錯滑移和變形孿生等多種細觀變形機制，在描述基于微觀結構演化的材料塑性行為方面具有明顯的優勢。而晶體塑性本構模型的參數通常是通過擬合宏觀實驗結果得到的，但是其缺乏亞微米變形機理，所以擬合參數可能不是唯一的，從而降低了CPFE模擬的預測精度。由于MPEAs的微觀結構是多尺度的，如原子空位和晶格畸變、微尺度位錯和中尺度晶粒等，所以需要考慮微尺度的變形機理來獲得精確的晶體塑性本構模型參數，然后開發一種從納米-微-中尺度微觀結構集成的新的模擬方法。湖南大學的Qihong Fang等人將原子模擬、離散位錯動力學和晶體塑性有限元方法結合起來，建立了一個新的框架，研究MPEAs的應變硬化行為，實現了包括納米尺度晶格畸變和微尺度位錯硬化在內的復雜跨尺度因素對塑性變形的影響，作者結合MD、DDD、CPFE模擬方法和隨機場理論(圖1)，提出了一種可捕捉MPEAs中嚴重晶格畸變的分層多尺度方法來建模MPEAs，該方法連接了三個長度尺度(納米尺度、微觀尺度和中尺度)，為深入理解納米-微米-中尺度結構相關的微尺度變形機制提供了新的思路，并為研究先進MPEAs的多尺度微結構調控相關的優越力學性能提供了可能和途徑。
圖1：用分層多尺度建模方法估計晶體塑性本構模型中的硬化參數。用MD、DDD和CPFE耦合模型預測了多晶材料在不同長度尺度下的力學響應。
圖2是通過MD模擬得到的Al0.1FeCoCrNi MPEA中邊緣位錯速度隨不同剪切應力/溫度比的變化規律。在作者測試的外加應力范圍內，位錯速度幾乎隨σ/T線性增加，這符合聲子阻尼理論。采用DDD模擬研究了邊緣和螺桿段遷移率對Al0.1FeCoCrNi單晶[001]取向應力應變曲線的影響。從圖3中可以看出，不同位錯遷移率下的應力應變曲線與相同遷移率下的應力應變曲線變化不大。
展開 ABAQUS三維Voronoi晶體幾何建模
材料晶體塑性理論與細觀尺度上晶體幾何模型相融合的模擬方法為探究材料在塑性變形過程中的行為機制以及晶體材料優化開辟了新途徑。本案例演示在CAD軟件內通過Voronoi建立晶體三維模型，并將模型導入到Abaqus CAE內，完成晶體材料的有限元建模。
在AutoCAD軟件內，采用CAD_Voronoi V1.0.1插件建立晶體結構三維模型，并將整個模型導出為.iges格式文件備用。
CAD_Voronoi V1.0.1插件
將導出的Voronoi模型文件以部件的形式導入到ABAQUS內。
插件可建立包含晶界的模型，在Abaqus內將晶格及晶界分別賦值不同材料。
也可建立無晶界模型，對不同晶格分別指派材料。
可將Voronoi晶體部件進行裝配。
及完成網格劃分等操作。
展開 晶體塑性VUMAT結合VUSDFLD實現晶體變形過程中的臨界狀態單元的刪除------案例二十八
并將返回值給予對應的變量，如下：
常用于獲取的變量如：
S：所有的應力分量
MISES:mises等效應力
TRIAX:應力三軸度
LODE:洛德角參數
PEEQ:等效塑性應變
SDEG:折減剛度系數
TEMP:積分點溫度
需要注意的是：
（1）該子函數不能應用于獲取用戶定義的狀態變量。同時注意2D與3D獲取的分量順序：
（2）需要在材料界面打開用戶自定義場
（3）做單元刪除時也要指明利用哪個狀態變量表示材料狀態，哪個狀態變量用于判斷是否符合應力刪除的特征。從而實現單元的刪除。
本案例介紹如下：
1，模型幾何尺寸20*20mm包含中心區域直徑為2.5mm的缺口
2，使用包含1514個CPE4R單元，每個單元包含代表一個單獨的晶粒
3，分別固定X0和Y0在x，y方向的自由度，并施加X1方向10%的工程應變
4，設置單元的最大Mises等效應力為100Mpa,超過100Mpa后單元失效刪除。
結果如下：
初始時刻的應力狀態
單元刪除效果展示
展開 huang晶體塑性umat耦合Johnson-cook 損傷模型，實現晶體材料彈-塑-損傷模擬分析
參考應變率：ε0
當滿足下列條件時，損傷初始化準則得以滿足：
等效塑性應變認為與應力三軸度和應變率相關聯。
θ^是無量綱溫度，表示為：
其中，θ是當前溫度，θ-melt是熔化溫度，θ_transition是指轉變溫度，在該溫度或低于該溫度時，損傷應變εpl_D的表達式不存在溫度依賴性。材料參數必須在轉變溫度或低于轉變溫度時測量。
損傷的發展可以公式化為:
公式中分母表示單元失效對應的Johnson-cook等效塑性應變，公式為：
分子表示為等效塑性應變增量，公式為：
公式中可以看到，損傷隨著塑性應變的增大不斷累積，直至材料的失效，通過損傷變量進一步與晶體材料的屈服面或者彈性性能的退化可以實現材料彈-塑-損傷的耦合模擬，當不對其進行耦合時，可以用來判斷材料的失效狀態與相關參數的關系。
參考文獻：《Crystal plasticity finite element modeling and simulation of diamond cutting of polycrystalline copper》編寫對應的材料子程序。在huang晶體塑性程序的基礎上，調用johnson-cookd損傷函數，編寫過程中，需要自定義響應的狀態變量，如等效塑性應變，等效塑性應變率，損傷變量，以及是否進行損傷單元的刪除分析。其中等效塑性應變增量的計算，通過滑移系統的分切應力與對應滑移系統剪切應變的乘積絕對值之后與等效應力的比值獲得。并最終實現損傷的表征，采用umat子程序進行編寫。
展開 