文章名稱：《Comparison of the implicit and explicit finite element methods using crystal plasticity》

doi：10.1016/j.commatsci.2006.08.002

推薦理由：作為顯式晶體塑性與隱式晶體塑性模型比較，以及適用性討論最經典的文章之一，詳細介紹了顯式與隱式求解器的區別，其研究結果表明：對于更簡單的加載條件，隱式方法的求解時間更短。在涉及接觸的載荷條件下，顯式方法被證明是優選的選擇。與使用多個處理器解決分析的隱式方法相比，顯式方法顯示出持續高水平的并行化效率。

眾所周知：在隱式方法中通常使用牛頓迭代方法，解涉及迭代，直到每個增量都滿足收斂標準。因此計算收斂于精確解（隱式：基于t+Δt時刻確定t+Δt時刻的狀態變量），而顯式方法中的有限元方程使用的通常是歐拉向前方法，并且在這種形式下，它們可以直接求解，以在增量結束時確定解，當增量步小于穩定時間增量時，總是可以保證計算的進行，但結果不一定收斂于精確解，（顯式：基于t時刻確定t+Δt時刻的狀態變量）因此準靜態問題使用隱式往往更加高效，而涉及到接觸則使用顯式可以保證計算穩定

由于晶體塑性模型考慮了介觀尺度的塑性變形的真實物理過程，因此被廣泛用于已在模擬金屬和金屬基材料中的介觀尺度下大變形和應變局部化，然而其高度非線性的積分過程，相較于傳統模型，其數值成本往往很高，因此基于復雜的晶體塑性模型更能體現顯式于隱式積分的差異，這里作者的討論顯式程序和隱式程序使用huang的亞彈性框架進行

由于隱式程序在大量的博士論文可以找到詳細內容，這里不做贅述，這里主要提到文獻關于顯式的一些內容。

顯式計算存在臨界穩定時間的概念，當最大增量步大于該時刻時會導致計算結果發散，從而使得計算結果失去意義。

其中Le是單元特征長度，Cd是膨脹波速，公式中可以看出Le越大，（單元尺寸越大）穩定時間越大，密度愈大，穩定時間越大（這意味著可以通過質量縮放加快計算時間，需要根據能量評估質量縮放的影響），需要注意的是顯式時間為真實變形時間，因此使用晶體塑性這類粘塑性本構模型顯式求解器并不適用于準靜態分析，計算時間太長，并且人為修正時間會影響本構的響應情況，當利用顯式求解器進行準靜態使用質量縮放是合理的備選方案（保證動能于總能量之比小于5%）

作者開發的顯式晶體塑性框圖

作者使用的是晶體塑性里面經典的PAN模型，材料為316L，晶體塑性參數為：

C11 = 205 GPa, C12 = 138 GPa, and C44 = 126 GPa，g0 = 50 MPa; g1 = 330 MPa; h0 = 225 MPa;gamma_dot - 0:001，n=20

注意：顯式求解過程涉及大量增量，每個增量都很容易求解。因此，顯式有限元計算程序非常適合由多個處理器分解和求解。考慮到這一點，VUMAT是用矢量化接口構建的。在每次增量開始時，應力和狀態變量數據以二維陣列的形式傳入。陣列中的每一列都包含與材料的積分點相關的信息。當使用多個處理器執行模擬時，可以將分析數據拆分為塊并獨立求解。因此，在子程序的編寫中保留了矢量化，以便可以實現最佳的處理器并行化。使用ABAQUS/explicit不需要相同的過程，因為時間增量足夠小。為了確定ABAQUS/explicit中初始時間增量的大小，在分析開始時將一組偽的微小應變增量傳遞給VUMAT。根據材料的應力響應，計算出穩定時間增量的保守值。有限元求解器要求材料在初始檢查時具有彈性。由于晶體塑性子程序計算密集，必須定義彈性應力-應變響應，以確保相對較小的時間增量。在這個初始時間增量計算階段，材料響應被定義為具有與用于描述晶體塑性子程序主體中的彈性相同的彈性特性。這確保了采用相對有效的時間增量。如果需要較小的增量，則可以使用更硬的彈性模量，盡管解決時間會更長。這不會影響材料在分析過程中的響應；它純粹是為了計算時間增量。

同時需要注意ABAQUS/標準采用Jaumann應力率，而ABAQUS/explicit采用Green–Naghdi應力率

Ω表示變形梯度極分解對應的旋轉部分，

在ABAQUS/標準中，材料被視為基于固定的全局坐標系。增量旋轉在每次增量開始時傳遞給UMAT。這個數組dR是應力和應變數組在上一次增量結束和當前增量開始之間旋轉的量。這樣處理：

其中σt+Δt是當前柯西應力，dσ是在先前增量中計算的應力增量。每次增量開始時傳遞的值是模型坐標系中的柯西應力。如圖所示，它在每次增量結束時都會“旋轉”。dR變量使用Hughes–Winget算法計算

其中I是單位矩陣，ΔW是增量速度梯度的反對稱部分，即W的增量形式

在VUMAT的情況下，材料位于同旋坐標系上，即坐標系隨材料旋轉。由于每次增量都很小，所以假設所有材質旋轉都是剛體。因此，對于每個增量，問題都可以看作是一個小應變問題。ABAQUS/explicit中使用的旋轉是Green–Naghdi自旋速率。Hughes–Winget算法也被使用，但形式略有不同

注意;由于晶體塑性子程序不包含運動硬化，兩種方法產生的結果相同。

作者比較了二維和三維典型的變形情況以及接觸問題，如下圖所示

二維拉伸：(比較了流動應力局部累計剪切和計算時間)

三維：

以及顯式和隱式處理器并行問題：

其研究結論為：

由于接觸和材料響應引起的非線性的嚴重性使得所需的時間步長遠小于實際可用于用隱式方法生成解的時間步長。顯式求解器更適合處理復雜的接觸和滑動條件，特別是在大單元變形的情況下。在包括剛體的接觸分析中，使用隱式代碼的運行時間短于顯式代碼，同時顯式可以處理近似率無關情況，隱式則難以計算，顯式對于多核的處理效率的提升顯著高于隱式

通過文獻的思路可以實現將黃隱式程序通過接口文件轉化為顯式的vumat子程序，單晶與多晶對應的案例如下（參數均為原始的Cu參數）：

案例一：顯式與隱式單晶含圓孔薄板的簡單拉伸對比（顯式開啟了質量縮放）

使用位移邊界條件，沿著X方向施加10%工程應變的拉伸變形，其中等效應力，累計剪切應變為對比指標，結果如下：

等效應力分布：

累計剪切應變分布：

案例二：顯式與隱式多晶的簡單壓縮對比（包含30個晶粒）（顯式開啟了質量縮放）

使用位移邊界條件，沿著X方向施加10%工程應變的拉伸變形，其中等效應力，累計剪切應變以及滑移系統的當前強度為對比指標，結果如下：

等效應力分布：

累計剪切應變分布：

滑移系統的當前強度：

可以看到這種基于接口的子程序可以顯式與隱式的一致性良好（部分差異在于顯式使用了過大的質量縮放）