CPFEM的案例
基于晶體塑性有限元（CPFEM）的鈦合金圓棒拉伸過程模擬
作者：辭殤
關鍵詞：CPFEM；鈦合金；單軸拉伸；織構極圖；孿晶
晶體塑性有限元是一種結合了晶體塑性理論和有限元方法的數值模擬技術?。這種方法考慮了晶體材料的各向異性、滑移系統的開動和相互作用、以及變形過程中的硬化效應。它主要用于分析和預測晶體材料的塑性變形行為，特別是在微觀尺度上的變形機制。
晶體塑性有限元在材料科學和工程領域有著廣泛的應用，特別是在金屬加工、航空航天、汽車制造和生物醫學等領域。通過這種技術，研究人員和工程師可以更好地理解材料的力學行為，從而開發出更輕、更強、更耐用的材料和產品。此外，晶體塑性有限元仿真還能夠考慮材料的微觀結構特征，如晶粒取向、晶界、相分布以及滑移系統的活動，從而能夠預測材料在細觀尺度上的織構演化。
利用CPFEM方法對鈦合金圓棒拉伸過程進行模擬，使用UMAT子程序以及Abaqus有限元軟件作為晶體塑性有限元分析的實現方式。并且，在一些復雜工藝條件下如切削、軋制、沖壓等，CPFEM方法同樣適用，能夠模擬材料變形過程中的非線性行為和動態響應。
在晶體塑性有限元中，首先在Abaqus中建立了單軸拉伸有限元模型如圖1所示，材料被建模為包含大量晶粒的集合體如圖2所示，每個晶粒都有其特定的晶體取向，并且每個晶粒的變形過程均考慮了滑移和孿晶的變形機制。
圖1 單軸拉伸有限元模型示意圖
圖2 單軸拉伸晶體塑性模型示意圖
通過有限元方法，可以計算出在給定拉伸載荷下，這些晶粒如何相互作用，以及它們如何隨時間變形。這種方法能夠提供關于晶體材料內部應力、應變和變形機制的詳細信息，有助于理解材料在受力時的響應，并優化材料的設計和加工過程。圖3所示為單軸拉伸過程應力云圖，圖4所示為單軸拉伸過程孿晶云圖。
展開 晶體塑性每日文章推薦（十四）
同時取向變化速率關系為：
（5）Cube和Goss取向的晶粒相比于RCube和RZ取向，對相鄰取向的晶粒更加敏感，取向變化對周圍晶粒取向的敏感性排序為：
（6）Cube和RZ取向的晶粒顯示相比于Goss更高的晶內旋轉梯度，而RC晶粒則分裂為兩個相等的coppor組分
作者的研究思路：
一：得到和實驗等效取向分布分取向數據和材料參數，這里作者研究一退火良好的鋁為例子，如下圖所示
二：研究taylor模型CPFEM模型以及CPFEM模型預測結果于網格的相關性，以平面應變壓縮為典型變形條件，模擬結果如下，這顯示：全場模型的預測能力高于taylor模型，同時CPFFT模型在預測織構演化，應變分布時幾乎分辨率無關
三，研究CPFEM和CPFFT在預測雙晶變形時預測晶內應變梯度的差異，結果顯示CPFFT預測的晶粒旋轉幾乎于晶內的單元個數無關，良好的收斂。并且預測的晶內取向旋轉于CPFEM結果存在較大差異
四，研究了典型取向的單晶演化速率受周圍晶粒取向影響的情況。
展開 晶體塑性有限元仿真入門(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不太常用。這篇文章將講解該平臺的使用方法以及如何使用該平臺進行晶體塑性有限元變形模擬。
展開 晶體塑性有限元仿真入門(2)--BCC、FCC、HCP晶格材料以及多相材料的有限元模擬
晶體塑性有限元仿真入門(2)--BCC、FCC、HCP晶格材料以及多相材料的有限元模擬
這篇文章講解如何使用晶體塑性有限元方法(CPFEM)進行不同晶格材料以及多相材料的變形模擬，CPFEM是基于商業有限元軟件ABAQUS完成的建模，晶體塑性本構模型是使用的開源的UMAT用戶子程序(源碼和inp文件見附件)。采用CPFEM模擬了面心立方結構(FCC)、體心立方結構(BCC)和密排六方結構(HCP)的單晶、多晶及多相材料受到外部載荷時的力學響應。基于滑移原理的晶體變形理論，隨著變形的進行各滑移系統的臨界剪應力都會增大，CPFEM將捕捉到材料的力學響應(應力-應變曲線)。這些應力-應變數據有助于從根本上理解晶粒尺度下金屬變形的性質。
首先我們從一個簡單的FCC晶格材料的例子入手，講解如何進行有限元模型的創建，從完全新手的角度出發，一步步講解如何建模，賦予材料和處理仿真結果。
本文章包括以下八個部分：
1) FCC晶格材料的變形模擬-單晶體
2) FCC晶格材料的變形模擬-多晶體
3) BCC晶格材料的變形模擬-單晶體
4) BCC晶格材料的變形模擬-多晶體
5) HCP晶格材料的變形模擬-單晶體
6) HCP晶格材料的變形模擬-多晶體
7) 多相材料的變形模擬
8) 參考資料
1. FCC晶格材料的變形模擬-單晶體
幾何模型
如圖1.1a在草圖里繪制R0.015mm的圓形，拉伸0.05mm，最后得到如圖1.1b所示的圓柱體(R0.015mm&H0.05mm)。
展開 
晶體塑性耦合連續損傷本構框架
結果表明，具有損傷模型的CPFEM可以合理地預測頸縮行為和頸縮方向，而沒有任何初始缺陷。
晶體塑性每日文章推薦（六）
《Simulation of texture evolution and macroscopic properties in Mg alloys using the crystal plasticity finite element method》
給出了結合PTR方法處理的TWIN-CPFEM方法的詳細方案，同時與宏觀實驗進行了對比，為HCP結構多晶模擬提供了相對簡單的解決思路，以及合理的參數用于AZ31B鎂合金，可以較為容易的在黃永剛程序上進行修改易于數值實現。
這里簡要介紹一下作者理論框架和數值模擬結果：
理論框架基于經典的亞彈性本構框架，不做贅述
流動方程使用經典的冪律流動方程
硬化方程：
其中作者考慮四組滑移＋一組孿晶
對應的材料參數為
另外所有滑移和孿晶系統使用統一的：參考剪切應變率0.001，率相關系數20
孿晶處理方案基于PTR方法
即孿晶體積分數達到臨界體積分數，晶粒整體發生旋轉，孿晶體積分數定義為：
即孿晶的剪切應變/孿晶特征剪切應變（例：AZ31B：拉伸孿晶特征剪切0.129），PTR模型假設允許每個方向相對于最活躍的孿晶系統中的鏡面的法線方向重新定向，如由CPFEM確定的。矩陣中晶格取向和孿晶區中晶格取向之間的變換矩陣T可以定義為
數值案例及實驗對比：
樣品初始取向：
沿著RD方向壓縮實驗結果：
作者的有限元邊界條件設置以及模擬結果
模擬的變形后取向
變形過程中各個系統的相對活躍情況
孿晶演化情況
模擬結果表明，單軸加載下的織構演變、孿晶體積分數和宏觀性能與加載f方向密切相關。變形模式的相對活動性可以從理論上解釋AZ31鎂合金板的拉壓屈服的不對稱行為和明顯的R值。
展開 晶體塑性有限元仿真入門(5)—歐拉角與晶體取向
sqrt(sum(cross(SDV_Input_hkl,SDV_Input_uvw).^2));
g_1 = [transpose(SDV_Input_uvw),transpose(SDV_Input_rst),transpose(SDV_Input_hkl)];
g_ij_new = G_1 * inv(g_1);
Euler = Rotation_to_Euler(g_ij_new);
SDV數據批處理
% SDV數據批處理，本程序適用于將輸出SDV數據轉換為歐拉角(角度制0~180)
clc
clear
close all
% Data = load('D:\CPFEM_Abaqus_V2\Course3_3\Out_end_rolling.txt');
Data = load('D:\CPFEM_Abaqus_V2\Course3_3\Out_end_compression.txt');
HKLandUVW = reshape(Data,length(Data)/6,6);
data = HKLandUVW(1:1:end,:);
Unique_data = data(1:20:end,:);
Euler = zeros(length(Unique_data(:,1)),3);
for i = 1:length(Unique_data(:,1))
Euler(i,:) = SDV_to_Euler(Unique_data(i,1:3),Unique_data(i,4:6));
end
Unique_Euler = Euler;
plot(Unique_Euler(:,1),'r.')
展開 abaqus晶體塑性有限元分析后處理晶界顯示插件分享
abaqus后處理中顯示晶界可以是多晶塑性分析更加直觀，但abaqus未內置此功能，需要通過二次開發實現，這里分享一個插件用于實現該功能，插件源于一位法國讀博士老哥的分享，將該插件放入到abaqus plug-in中即可輕松的實現后處理晶界的顯示問題
軟件用戶界面：
得到的效果圖如下：
如果您在文章中使用了該插件，請引用該作者對應的兩篇文獻：
1，A physically-based mixed hardening model for the prediction of the ductility limits of thin metal sheets using a CPFE approach
2，Investigation of the effect of morphological and crystallographic textures on the ductility limits of thin metal sheets using a CPFEM-based approach
插件壓縮包：
鏈接：https://pan.baidu.com/s/1tlAH2su4DiY5zGnxBFMfUA?pwd=cxbi
提取碼：cxbi
此外如果使用該插件存在任何疑問，可以加入我的知識星球，并提出使用過程的疑問。會根據你遇到的問題進行統一解答。
展開 超彈性與亞彈性，顯式與隱式的HCP多晶滑移+孿晶（主導孿晶重定向（PTR））計算效率比較
參考文獻：《Influence of texture distribution in magnesium welds on their non-uniform mechanical behavior: A CPFEM study》
主導孿晶重定向（PTR）方案作為目前處理HCP晶格結構的多晶材料孿晶模擬中最常使用的方案被廣泛討論，然而晶體取向旋轉過程可能會造成模擬的收斂性問題，選擇一個相對穩定的本構框以及迭代變量對模擬計算效率的提升是有意義的。文章中公式以及其對應的參數總結如下：
這里使用文章的模型和參數對超彈性和亞彈性PTR方案進行比較。
二維（200個晶粒X方向壓縮20%）
以下各個圖中左圖為超彈性結果，右圖為亞彈性結果：
應力分布云圖
應變分布云圖：
孿晶分布云圖：
這里使用文章的模型和參數對顯示和隱式PTR方案進行比較
二維（200個晶粒Y方向剪切變形20%）
以下各個圖中左圖為顯示結果，右圖為隱式結果：
應力分布云圖
應變分布云圖：
孿晶分布云圖
拉壓非對稱與織構演化方面超彈性與亞彈性保持一致：
初始極圖：
RD拉伸20%：
RD壓縮20%：
應力應變曲線
模擬的結果建議，使用PTR方案，超彈性建議使用PK2應力和當前強度為迭代變量，并使用雙重迭代方案，亞彈性建議使用柯西應力為迭代變量，兩者在模擬過程中，計算效率相差較小，無論是局部晶粒的應力應變響應，整體的流動應力，以及變形后的織構結果幾乎保持一致。同時涉及到接觸，碰撞問題，修改為顯式對于收斂性的提升是必要的。
展開 晶體塑性耦合元胞自動機模擬熱壓縮過程中的再結晶行為
參考文獻：《A straightforward 3D polycrystal plasticity finite element method for dynamic/static recrystallization simulation》
文章doi：10.1016/j.jmst.2024.09.005
在這個文章中，作者提出了一種直接在 CPFEM 中實現 DRX/SRX 的方法，以位錯密度為核心變量，利用 UMAT 進行應力積分、在 UEXTERNALDB 中執行形核與晶界遷移，在 FE 網格上同步更新取向、位錯密度、等狀態變量，同一仿真中先熱壓縮（出現 DDRX 多循環與流動應力峰值震蕩），隨后“原地”退火，追蹤 SRX 繼續演化，避免傳統方法里變形場到顯微組織模型的跨尺度映射難題。
需要注意的是，作者模擬的兩類再結晶（動態和靜態）使用的機制是相同的，只是材料參數不同。
本構理論分成晶體塑性和再結晶兩部分，其中晶體塑性部分公式如下：
流動方程（經典的唯象流動）：
硬化方程使用的taylor位錯模型
位錯密度的演化使用經典的KM方程：
再結晶部分公式包含形核和晶界遷移兩部分，其中形核的理論公式是
晶界遷移速度為：
整體數值實現框架示意圖如下：
作者以OFHC銅為研究對象，對775K和875K的熱壓縮進行了研究，分析了溫度對再結晶的影響，以及定向形核和生長選擇兩類機制的差異，同時模擬了順序耦合的 DRX→SRX（退火）過程及異常晶粒長大（AGG），模擬效果如下：
根據作者提供的思路（相對簡單清晰），可以編寫對應的子程序，完整晶體塑性和元胞自動機的完全耦合，同樣使用隱式umat實現。
展開 基于有限元-元胞自動機法（CAFE）的增材制造過程組織模擬
通過CA模型捕獲的晶粒結構，特別是其??晶體取向??信息，可以直接為后續的??晶體塑性有限元（CPFEM）?? 模擬提供微觀結構輸入。CPFEM利用每個晶粒的取向和滑移系來預測其在外部載荷下的塑性變形行為和各向異性力學響應，從而建立從制造過程→凝固組織→宏觀力學性能的完整預測鏈條。這使得CAFE方法成為理解和優化增材制造過程-結構-性能關系不可或缺的核心工具。</p><div contenteditable="false" width="100%" class="ql-align-justify">
<img src="https://p3-sign.toutiaoimg.com/tos-cn-i-6w9my0ksvp/e6ada7111e6d44bb843fd19ba475eac8~tplv-tt-origin-web:gif.jpeg?_iz=58558&from=article.pc_detail&lk3s=953192f4&x-expires=1758680646&x-signature=Oc4dHm8CFXOfll2FVsVJ9kwEkPI%3D" img_width="561" img_height="323" image_type="1" mime_type="image/png" web_uri="tos-cn-i-6w9my0ksvp/e6ada7111e6d44bb843fd19ba475eac8" class="syl-page-img" data-src="https://p3-sign.toutiaoimg.com/tos-cn-i-6w9my0ksvp/e6ada7111e6d44bb843fd19ba475eac8~tplv-tt-origin-web:gif.jpeg?
展開 
考慮GND的大變形冷軋模擬
參考文獻：《Rolling deformation mechanism of dual-phase NiTiNb shape memory alloy thin strip based on crystal plasticity finite element method》
該文章聚焦雙相 NiTiNb 形狀記憶合金薄帶在 20% 軋制壓下量下的微觀變形機制，作者用晶體塑性有限元（CPFEM）在三維多晶尺度上解析 NiTi（硬相）與 β-Nb（軟相）在軋制過程中的“分工”，并用實驗應力–應變曲線對模型進行了驗證。
其中流動方程使用經典的唯象流動方程：
硬化模型則使用了同時考慮SSD和GND的位錯密度硬化模型：
作者構建了包含 520 個晶粒的三維 RVE（NiTi 基體晶粒 442 個、β-Nb 晶粒 78 個），并在 ABAQUS 中進行單道次軋制變形20%的模擬。
軋制模型如下所示：
變形后的SSD的GND分布如下：
作者的研究表明：
1）多晶變形不均勻性來源：軋制過程中應變分布的非均勻性主要由晶粒取向差異及相/結構差異共同導致；滑移優先在晶界與自由表面萌生，并沿晶內逐步擴展形成明顯的滑移帶。
2）織構演化特征：NiTi 相在軋制區形成以 Cube 織構 {001}<100> 為主的織構組分（約 9 MRD），同時伴隨較弱的 γ-fiber、α-fiber 等；β-Nb 相的變形以 {110}<111> 與 {112}<111> 滑移體系主導，并呈現 λ-fiber 等織構特征。
3）“硬/軟相協同”力學響應：NiTi 相表現出較高的承載能力與較低的剪切應變率，其屈服強度約 1027 MPa，且 SSD 累積更高；β-Nb 相更易發生滑移，屈服強度約 364 MPa，剪切應變率更高而 SSD 相對更低。
展開 晶體塑性模擬中的大變形網格重劃分
相關做法完美的集中到damask3.0版本里面，然而需要指出的是：DAMASK/譜方法更偏向規則網格與RVE范式，而工程里經常需要：任意幾何與復雜邊界（非周期、接觸、局部細化等），以及不同工藝路徑（多道次、換向、局部約束），Abaqus CPFEM（UMAT/VUMAT）在這些方面更“通用”，所以把“remesh + 狀態變量映射”做成一套工作流，就能把大變形晶體塑性更穩地推進到更高壓縮/更大應變階段。
因此結合作者提供的思路，嘗試把相關方法遷移到abaqus，并初步實現了理想的效果。
這里展示模擬的案例的效果，初始模型尺寸0.1*0.03mm的二維模型，并沿著RD方向壓縮40%.vs.20%（remesh）+20%使用簡單的唯象模型測試
初始模型如下圖所示：
壓縮20%后應力分布如下：
累計剪切滑移如下：
晶粒旋轉角度：
在20%變形后進入網格重劃分，重劃分后的變量傳遞：
累計剪切滑移分布如下：
晶粒旋轉角度如下圖：
可以看到所有相關變量良好的映射到規則網格上面。
接下來對比單次壓縮40%（左側）和20%remesh（右側）之后再壓縮20%的結果對比：
應力分布結果：
累計剪切滑移分布：
晶粒旋轉角度分布：
累計剪切滑移------應力曲線分布
重劃分后應力略低于不劃分單次壓縮的結果，其余結果網格重劃分和原始模型基本一致，驗證了作者提出方案的準確性。做成型和大變形相關內同可以參考文章進行對應的嘗試
展開 任意多晶微觀結構生成，GUI操作，模型直接下載
多格式導出： 生成的模型支持導出為坐標數據、拓撲連接信息等，方便后續導入 ABAQUS、ANSYS 或自編的有限元/晶體塑性（CPFEM）程序中。
【操作流程：三步搞定】
第一步：設定全局參數。 在左側面板選擇晶粒總數及 RVE 尺寸。
第二步：精修幾何特征。 調整權重系數（Weights）和偏度，生成不規則或特定分布的晶粒形狀。
第三步：導出與應用。 預覽滿意后，點擊下載按鈕獲取幾何模型文件。
五大基礎案測試例如下：
一：2D單相對數正態分布1000個晶粒模型（尺寸1*1）（0.1s）
二：2D雙相對數正態分布500個晶粒模型（0.1s）（體積分數0.2：0.8，晶粒個數450：50）
自動生成對應的統計信息：
點擊下載可以直接獲取所有相關的晶粒信息
三：隨機的500個晶粒的3D單相模型（0.1s）
四：隨機的1000個晶粒的3D雙相模型（0.5s）二次相分布于晶界
五：隨機的2000個晶粒的周期性單相模型（2s）
同時內置豐富的可視化，比如調整透明度，切片顯示，修改顏色等等
功能非常豐富直接登錄網站即可使用：
https://david-bourne.shinyapps.io/synthetmic-gui/
展開 一種新穎的多尺度晶體塑性實現方案-------Direct FE2
參考文獻:《A novel concurrent multiscale method based on the coupling of Direct FE2 and CPFEM》
文章doi：10.1016/j.tws.2024.112610
Direct FE2（Finite Element Square）多尺度模擬方案是一種常用于材料科學中的數值模擬方法，特別適用于研究復雜的材料行為（如晶體塑性、微觀結構演化等）。它結合了微觀和宏觀兩個尺度的有限元分析，能夠在多尺度上捕獲材料的細觀力學響應。
Direct FE2 的原理：
Direct FE2 是一種強耦合多尺度方法，其核心思想是在每個宏觀積分點處嵌套一個細觀有限元模型，以實現宏觀響應與微觀行為之間的直接相互作用。
宏觀有限元（FE）模型：
描述宏觀結構的整體變形和力學行為。
宏觀模型的單元中積分點的材料行為由微觀模型提供。
宏觀模型中的應力、應變是由微觀模型的細觀響應直接計算得到。
微觀有限元（RVE）模型：
微觀模型（通常為代表性體積單元，Representative Volume Element, RVE）用于描述材料微觀結構，如晶粒、裂紋、孔洞等。
通過微觀有限元模擬，計算微觀應力和應變場，以及其等效宏觀響應。
微觀模型會受到來自宏觀模型的邊界條件（如位移或應力）的驅動。
宏微耦合：
宏觀積分點傳遞的變形梯度作為邊界條件輸入到微觀模型。
微觀模型計算得到的等效應力或切線剛度返回到宏觀模型，用于更新宏觀求解。
這種雙向耦合的計算流程在每個加載步中反復迭代。
Direct FE2 的計算流程
在宏觀有限元模型中，計算每個單元的變形梯度 FFF。
展開 