基于相場理論的固態(tài)相變模擬的有限元計(jì)算

1、 工作總結(jié)-基于相場理論的固態(tài)相變模擬的有限元計(jì)算目錄1、相變機(jī)理12、相場理論32.1相場理論概述32.2相場的控制微分方程32.3體積化學(xué)自由能42.4界面能72.5彈性應(yīng)變能72.6多物理場耦合83、有限元計(jì)算93.1有限元的離散化93.1.1 相場控制微分方程的化簡93.1.2 相場與位移場的離散化103.2相場與位移場的動力求解113.2.1 相場的動力求解113.2.1 位移場的動力求解123.3相變模擬計(jì)算流程134、已完成工作總結(jié)144.1相場模型的研究總結(jié)144.1.1 彈性應(yīng)變能的影響144.1.2 體積化學(xué)能的影響164.1.3 界面能的影響164.1.4 動力學(xué)系數(shù)17

2、4.1.5 隨機(jī)噪聲174.1.6 穩(wěn)定結(jié)果的非均勻性174.2有限元算法的研究總結(jié)184.2.1 時域積分184.2.2 單元劃分184.2.3 單元階次194.2.4 高斯點(diǎn)選取194.2.5 計(jì)算效率191、相變機(jī)理相變是材料從高能量相轉(zhuǎn)變?yōu)榈湍芰肯?，總能量減小的過程。在材料中，各相能量的高低是由其內(nèi)部狀態(tài)（如應(yīng)力、溫度等）所決定的。例如對于NiTi形狀記憶合金，在低溫時，馬氏體為低能量相，當(dāng)溫度從高溫降為低溫時，奧氏體轉(zhuǎn)變?yōu)轳R氏體；在高溫時，奧氏體為低能量相，當(dāng)溫度從低溫升為高溫時，馬氏體轉(zhuǎn)變?yōu)閵W氏體。相的變體是指晶體結(jié)構(gòu)相同，取向不同的同一種相1。例如對于NiTi形狀記憶合金，奧氏體

3、的晶體結(jié)構(gòu)為高對稱的立方體（ high symmetry cubic），馬氏體的晶體結(jié)構(gòu)為四方體，在不同的朝向下，共有三種變體，如圖1-1。圖1-1 奧氏體與馬氏體的晶體結(jié)構(gòu)示意圖在一定的內(nèi)部狀態(tài)下，體系生成各變體后，總能量的大小是不同的。體系會自發(fā)地選擇生成使得體系能量降低最多的變體。如圖1-2所示，各變體首先在使得體系能量降低最多的區(qū)域形成，后生長形成特定的位向關(guān)系。圖1-2 變體產(chǎn)生及位向關(guān)系形成示意圖2、相場理論2.1相場理論概述相場理論是一種用來模擬相變微觀組織的演化過程的理論，它的主要特點(diǎn)是以Ginzburg-Laudau相變理論為基礎(chǔ)，基于擴(kuò)散界面模

4、型，引入連續(xù)變量來描述新舊兩相，用微分方程描述系統(tǒng)自由能的變化。 相場理論的特點(diǎn)具體如下：（1）離散變量的連續(xù)化：在相變中，各相或者相的各變體是離散的物理量，各相或者相的各變體之間不存在連續(xù)變化的過程。相場模型將這些離散變量連續(xù)化，用序參量表示，便于在數(shù)值上獲取導(dǎo)數(shù)等信息。例如，對于新相含有n個變體的相變，采用序參量1=2=n=0表示體系此處為母相所占有；序參量p=1或-1且k=0（k=1,2.p-1,p+1,n）表示此處對應(yīng)為新相的第p個變體。由于連續(xù)化的處理，各相或者各變體間的界面是彌散界面（diffuse interface），而不是尖銳界面（Sharp interface），如圖2-1

5、。（2）微分方程：建立相場內(nèi)的平衡微分方程，描述各相或各變體的演化。（3）界面能：基于擴(kuò)散界面模型建立界面能，描述相場內(nèi)部序參量的變化。圖2-1 尖銳界面與彌散界面示意圖2.2相場的控制微分方程在相場理論中，相變微觀組織演化通過求解如下的Ginzburg-Landau方程： （2-1）式中，L為相場動力學(xué)系數(shù)，控制著相場的演化速率；G為系統(tǒng)總的自由能，在馬氏體相變中，一般來說，包括體積化學(xué)能、界面能和彈性應(yīng)變能；為高斯隨機(jī)噪聲，滿足以下條件： （2-2） （2-3）其中，為相場熱起伏對空間和時間的平均；kB為Boltzmann常數(shù)，T為溫度，為Kronecker Delta函數(shù)。2.3體積化學(xué)

6、自由能體積化學(xué)能是指材料本身具有的一種化學(xué)勢能，各相之間的體積化學(xué)能密度與溫度相關(guān)，同一相中的各變體具有相同的體積化學(xué)能密度。兩相間的體積化學(xué)能密度之差促使相變從能量高的相向能量低的相轉(zhuǎn)變。在相變中，體積化學(xué)自由能具有以下兩點(diǎn)性質(zhì)（1）在母相處，體積化學(xué)自由能密度為0，在新相處，體積化學(xué)自由能密度為-fchem；（2）當(dāng)某點(diǎn)為母相或新相時，該點(diǎn)的體積化學(xué)能驅(qū)動力為0，處于穩(wěn)定態(tài)或亞穩(wěn)定態(tài)。這兩點(diǎn)性質(zhì)要求體積化學(xué)自由能密度函數(shù)滿足： （2-4）在不同的文獻(xiàn)中，體積化學(xué)能密度函數(shù)chem主要有以下兩種形式：（1）2-4-6次多項(xiàng)式2 （2-5）A1，A2，A3為體系膨脹系數(shù)，一般取正數(shù)；n為變體總

7、數(shù)；序參量1=2=n=0表示體系此處為母相所占有；序參量p=±1且k=0（k=1,2.p-1,p+1,n）表示此處對應(yīng)為新相的第p個變體。當(dāng)k=0（k=1,2.p-1,p+1,n）時，體積化學(xué)能密度chem與序參量p的關(guān)系曲線見圖2-2與圖2-3。由圖所示，該插值函數(shù)有三個穩(wěn)定點(diǎn)：p=0,±1。在該插值函數(shù)下，相變從母相到新相轉(zhuǎn)化時，需先從外界吸收一部分能量，達(dá)到躍遷所需的能量后（如圖2-3中p=±0.1），才能向外界放出能量，完成相變轉(zhuǎn)化。圖2-2 2-4-6多項(xiàng)式體積化學(xué)能密度插值函數(shù)與序參量的關(guān)系曲線圖2-3 2-4-6多項(xiàng)式體積化學(xué)能密度插值函數(shù)與序參量的

8、關(guān)系曲線（靠近0部分）（2）2-3-4次多項(xiàng)式3 （2-6） 與2-4-6次多項(xiàng)式類似，A1，A2，A3為體系膨脹系數(shù)，一般取正數(shù)；n為變體總數(shù)；序參量1=2=n=0表示體系此處為母相所占有；序參量p=1且k=0（k=1,2.p-1,p+1,n）表示此處對應(yīng)為新相的第p個變體。當(dāng)k=0（k=1,2.p-1,p+1,n）時，體積化學(xué)能密度chem與序參量p的關(guān)系曲線見圖1-3與圖1-4。由圖所示，該插值函數(shù)有兩個穩(wěn)定點(diǎn)：p=0,1。與2-4-6次多項(xiàng)式類似，在該插值函數(shù)下，相變從母相到新相轉(zhuǎn)化時，需先從外界吸收一部分能量，達(dá)到躍遷所需的能量后（如圖1-4中p=0.011），才能向外界放出能量，完

9、成相變轉(zhuǎn)化。 圖2-4 2-3-4多項(xiàng)式體積化學(xué)能密度插值函數(shù)與序參量的關(guān)系曲線圖2-5 2-3-4多項(xiàng)式體積化學(xué)能密度插值函數(shù)與序參量的關(guān)系曲線（靠近0部分）2.4界面能在相場模型中，界面能表示的是由于相場不均勻分布所產(chǎn)生的能量。在相變中，界面能的釋放將將促使相場分布的均勻化。界面能與序參量的梯度相關(guān)，其表達(dá)式為： （2-7）其中，為界面梯度能系數(shù)；一般為了簡化問題，將材料考慮為各向同性，且各變體的界面梯度能系數(shù)相等，界面能表達(dá)式可化簡為： （2-8）2.5彈性應(yīng)變能在相變中，由于晶體結(jié)構(gòu)的改變，從母相到新相，材料將發(fā)生相變變形。在變形協(xié)調(diào)下，相變變形使得材料的變形重分布，改變材料的彈性應(yīng)變

10、能。對于線彈性材料，彈性應(yīng)變能密度的表達(dá)式為： （2-9）在相變中，總應(yīng)變包含彈性應(yīng)變與相變應(yīng)變： （2-10）其中，相變應(yīng)變在相場模型中采用序參量構(gòu)造插值函數(shù)()來表示： （2-11）式中，為第p個變體由于晶體結(jié)構(gòu)改變產(chǎn)生的相變應(yīng)變，是材料常數(shù)。彈性應(yīng)變能密度的表達(dá)式可展開為： （2-12）因此，彈性應(yīng)變能驅(qū)動力為： （2-13）2.6多物理場耦合相變是從高能量相轉(zhuǎn)變?yōu)榈湍芰肯嗟倪^程，由于對于材料，相的能量高低是由內(nèi)部狀態(tài)所決定的，相變中微觀組織演化將受到其它物理場的影響。同時，相變中微觀組織的演化也會改變其它物理場的分布。因此，相變本身是一個多物理場耦合的問題。在馬氏體相變中，相互耦合的三

11、個物理場分別是相場、位移場與溫度場，三者的關(guān)系如圖2-6所示。圖2-6 馬氏體相變?nèi)龍鲴詈鲜疽鈭D在實(shí)際模擬中，一般忽略相變過程中的發(fā)熱、吸熱，假定在相變過程中，溫度不發(fā)生改變，原三場耦合的問題簡化為兩場耦合，如圖2-7所示。圖2-7 簡化后的馬氏體相變兩場耦合示意圖對于相場與位移場的耦合，本研究計(jì)算采用順序耦合方法求解。具體求解策略是按照順序進(jìn)行位移場和相場的分析，采用上一時刻的相場分布，先求解位移場，得到上一時刻位移場分布，將位移場所得結(jié)果應(yīng)用到相場求解中，再求解這一時刻的相場分布。3、有限元計(jì)算3.1有限元的離散化3.1.1 相場控制微分方程的化簡對體積化學(xué)能、界面能和彈性應(yīng)變能三種能量之

12、和的總能量泛函求變分，可得： （3-1）將式（3-1）代入式（2-1）的Ginzburg-Landau方程，可得： （3-2）式（3-2）與非穩(wěn)態(tài)的溫度場控制方程類似，在有限元離散化中，只需保證C0連續(xù)，只需采用常規(guī)的有限元插值。3.1.2 相場與位移場的離散化將相場中的變量序參量用有限元離散化，如圖3-1，采用節(jié)點(diǎn)變量表示序參量，將式（3-2）離散化，可得： （3-3）（1）C由單元矩陣c組裝后所得 （3-4）（2）K由單元矩陣k組裝所得 （3-5）（3）Qk包含相變驅(qū)動力和隨機(jī)項(xiàng)k （3-6）圖3-1 有限元離散化示意圖同理，位移場的動力方程或者靜力方程可離散化為： （3-7）或：（1）M

13、為單元質(zhì)量矩陣m組裝后所得 （3-8）（2）K為單元剛度矩陣k組裝所得 （3-9）（3）R包含外荷載fext和相變引起的等效荷載fequi3.2相場與位移場的動力求解3.2.1 相場的動力求解相場中ti+1時刻的序參量與ti時刻的序參量的關(guān)系通過下式表示： （3-10）將上述遞推公式代入離散后的相場方程，整理可得 （3-11）其中等效剛度陣和等效右端項(xiàng)為： （3-12） （3-13）該求解方法數(shù)值穩(wěn)定的最大時間步為： （3-14）為方程中的最大特征值。當(dāng)=0，該方法為時域顯式積分法，式（3-8）和式（3-9）可化為 （3-15） （3-16）3.2.1 位移場的動力求解（1）Newmark-法

14、：在Newmark-法中ti+1時刻的位移、速度、加速度與ti時刻的位移、速度、加速度的關(guān)系通過下式表示： （3-17） （3-18） 將式（3-17）和式（3-18）遞推公式位移場的動力方程，整理后得逐步積分公式為： （3-19）其中，等效剛度陣和等效荷載向量分別為： （3-20） （3-22） Newmark-法的穩(wěn)定性條件為： （3-23）當(dāng)=1/2，=1/4時，算法為無條件穩(wěn)定，同時具有二階精度。（2）中心差分法：在中心差分法中ti時刻的速度、加速度與ti+1時刻和ti-1時刻的位移的關(guān)系通過下式表示： （3-24） （3-25）將遞推公式代入力平衡方程，得： （3-26）其中： （3

15、-27） （3-28）中心差分法的穩(wěn)定性條件為： （3-29）3.3相變模擬計(jì)算流程對于位移場和相場的耦合問題，采用的求解策略是在根據(jù)t時刻序參量分布值，先在t時刻先求解位移場，再根據(jù)得到t時刻的應(yīng)力分布值，求解t+t時刻相場，得到t+t時刻的序參量分布值，具體流程如圖2-8所示： 圖3-1 相變模擬計(jì)算流程圖4、已完成工作總結(jié)4.1相場模型的研究總結(jié)馬氏體相變的相場模型包含三種能量：體積化學(xué)能、界面能、彈性應(yīng)變能，三種能量相互作用，改變其中任何一種，都將影響到最終結(jié)果。4.1.1 彈性應(yīng)變能的影響根據(jù)彈性應(yīng)變能的表達(dá)式（2-12）與彈性應(yīng)變能驅(qū)動力的表達(dá)式（2-14），彈性應(yīng)變能的組成分為三

16、部分：（1）應(yīng)力：應(yīng)力的大小影響著相變的演化速率和演化方向。當(dāng)應(yīng)力的方向與相變應(yīng)變一致時，應(yīng)力促進(jìn)相變演化；反之，抑制相變演化。從應(yīng)力產(chǎn)生的原因來分，應(yīng)力可分為外荷載作用下引起的應(yīng)力和相變應(yīng)變引起的應(yīng)力。由于相變應(yīng)變引起的應(yīng)力一定與相變應(yīng)變相反，因此，相變應(yīng)變引起的應(yīng)力總是抑制相變進(jìn)一步演化。（2）彈性模量：在線彈性下，彈性模量與相變應(yīng)變引起的應(yīng)力大小成線性關(guān)系，因此，彈性模量越大，由相變引起的應(yīng)力大小越大，相變的演化越慢，反之，相變演化越快。（算例驗(yàn)證）（3）泊松比：泊松比的大小與相變應(yīng)變引起的應(yīng)力大小相關(guān)，其影響與彈性模量相似，只是對不同的應(yīng)變分量的影響程度不同。（算例驗(yàn)證）（4）相變應(yīng)變

17、插值函數(shù)：相變應(yīng)變的插值函數(shù)影響彈性應(yīng)變能在序參量值介于0-1之間的大小，進(jìn)而影響彈性應(yīng)變能驅(qū)動力的大小，影響相變演化進(jìn)程與最終結(jié)果。在選取相變應(yīng)變插值函數(shù)時，需注意以下準(zhǔn)則：I.必須滿足：為保證當(dāng)序參量為0時，相變應(yīng)變?yōu)?；當(dāng)序參量為1時，相變應(yīng)變?yōu)?，彈性?yīng)變能插值函數(shù)()必須滿足： （4-1） II.非必須滿足：為了使序參量位于0或者1時，材料處于穩(wěn)定或者亞穩(wěn)定態(tài)，其能量驅(qū)動力應(yīng)為0，由于體積化學(xué)能驅(qū)動力在序參量位于0或者1時為0，彈性應(yīng)變能驅(qū)動力也應(yīng)為0。因此，相變應(yīng)變插值函數(shù)()須滿足： （4-2）按照以上選取準(zhǔn)則，已進(jìn)行過對比的相變插值函數(shù)如下： 表達(dá)式形式一次導(dǎo)數(shù)表達(dá)式線性二次三次

18、2-3次2-4次2-3-4次其中，在材料內(nèi)部無缺陷且無外力作用下，對比線性、二次、三次插值函數(shù)，線性的結(jié)果好于二次和三次，原因是隨著插值函數(shù)次數(shù)增大，在序參量靠近1的部分彈性應(yīng)變能驅(qū)動力增大，最終穩(wěn)定時，序參量離1的平均距離變遠(yuǎn)。與此同時 ，體積化學(xué)能驅(qū)動力在1附近的增大幅度大于彈性應(yīng)變能驅(qū)動力，因此，應(yīng)變能驅(qū)動力在總驅(qū)動力的比例降低，最終形成的位向關(guān)系變差。（算例驗(yàn)證）2-3次、2-4次、2-3-4次是在材料內(nèi)部無缺陷且無外力作用時，結(jié)果無明顯差別。對于材料受到較大外力作用，或因?yàn)樽陨砣毕菀饝?yīng)力集中，造成相變區(qū)域?yàn)榇髴?yīng)力區(qū)域的情況，還未進(jìn)行算例分析。4.1.2 體積化學(xué)能的影響在溫度不變的

19、情況下，體積化學(xué)能只與自身的插值函數(shù)相關(guān)，與應(yīng)力大小無關(guān)。體積化學(xué)能是否促進(jìn)相變演化取決于在選取的溫度下，母相的體積化學(xué)能是否較大。4.1.3 界面能的影響由于相場模型本身是將離散的相通過連續(xù)的序參量值來表示，界面能采用梯度能的描述方式，與實(shí)際情況并不相符。實(shí)際各相或者各變體間是否存在界面能，界面能的大小與什么參數(shù)相關(guān)，還需進(jìn)一步閱讀文獻(xiàn)。單就相場模型中的界面能，可以從以下角度分析：（1）數(shù)學(xué)表達(dá)式：界面能表示由于相場不均勻分布所產(chǎn)生的能量。在相變中，界面能的釋放將將促使相場分布的均勻化。（2）模型中的作用：通過梯度描述的界面能，在相場模型中，建立相變區(qū)域中相鄰各點(diǎn)的聯(lián)系，使得新相或者變體能夠

20、通過“擴(kuò)散”的方式生長。（3）界面能系數(shù)對模擬結(jié)果影響：界面能系數(shù)的大小影響相變模擬結(jié)果中各變體的平均尺寸和彌散界面的寬度，界面能系數(shù)越大，各變體的平均尺寸越大，彌散界面的寬度越大。（算例驗(yàn)證）4.1.4 動力學(xué)系數(shù)在式（2-1）的Ginzburg-Landau方程中，動力學(xué)系數(shù)控制著相變演化的速率。由于相變演化本身還與位移場相互耦合，相變演化的速率影響慣性力的大小，決定位移場分析是采用動力分析還是靜力分析。4.1.5 隨機(jī)噪聲相場模型是一個描述相變演化過程的模型，無法描述相變的形核機(jī)理，通常要求其化學(xué)能驅(qū)動力與彈性應(yīng)變能驅(qū)動力為0，在新相形核前，體系處于亞穩(wěn)定狀態(tài)。引入隨機(jī)的噪聲，主要目的是

21、引起序參量值的微小改變，從而使得相變能夠演化下去。因此，在隨機(jī)噪聲很小的情況下，與空間坐標(biāo)獨(dú)立的隨機(jī)噪聲抽樣方式和持續(xù)時間不影響相變結(jié)果的位向關(guān)系和變體的平均尺寸。（算例驗(yàn)證）4.1.6 穩(wěn)定結(jié)果的非均勻性根據(jù)相場的控制方程，當(dāng)結(jié)果穩(wěn)定時，原方程可化為： （4-11） （4-12）假設(shè)變體中間區(qū)域，序參量為0或者1的均勻分布，界面能為0，方程進(jìn)一步簡化為： （4-13）在變體邊界，序參量不均勻分布，存在界面能，彈性應(yīng)變能驅(qū)動力與體積化學(xué)能驅(qū)動力之和不為0： （4-14）由于式（4-13）與式（4-14）互斥，滿足式（4-13）的序參量值必然不滿足式（4-14）。變體中間區(qū)域的序參量值必然不等于

22、變體邊界上的序參量值，因此，中間區(qū)域的序參量值與邊界上的序參量值存在一定的梯度，變體中間區(qū)域的序參量值不滿足式（4-13），與假設(shè)矛盾。4.2有限元算法的研究總結(jié)4.2.1 時域積分（1）顯式積分與隱式積分：在完全積分的情況下，由于求解相場驅(qū)動力與相變引起的等效荷載所需的計(jì)算時間較長，同線性回代時間差不多，因此，當(dāng)每步求解過程中，剛度陣不發(fā)生改變時，無論是在相場計(jì)算還是位移場計(jì)算（動力分析）中采用顯式積分不能明顯地提高效率。由于顯式積分的等效剛度陣為對角陣，相比隱式積分未經(jīng)優(yōu)化的變帶寬等效剛度陣，占用內(nèi)存較少。應(yīng)用顯式積分，可以大幅降低內(nèi)存使用。（2）積分步長：步長是否合理的判斷：I是否滿足相

23、場數(shù)值計(jì)算的穩(wěn)定性要求；II序參量是否大比例超限。在合理步長下，會造成相變結(jié)果局部區(qū)域結(jié)果的不同，但所有結(jié)果整體上都符合位向關(guān)系。同時在合理步長下，步長較大，雖然計(jì)算結(jié)果收斂時對應(yīng)的相變實(shí)際時間較長，但在模擬中所需要的計(jì)算步數(shù)較少，需要的數(shù)值計(jì)算時間較短。4.2.2 單元劃分在位移場中，提高網(wǎng)格精度，可以更準(zhǔn)確地描述區(qū)域內(nèi)應(yīng)力分布；在相場中，提高網(wǎng)格精度，可以更準(zhǔn)確地描述區(qū)域內(nèi)序參量的梯度分布，更好地描述界面能的大小，計(jì)算結(jié)果更精確合理。同時，提高網(wǎng)格精度，在圖形顯示上，可以更好地描述變體間的邊界，避免出現(xiàn)鋸齒型區(qū)域。4.2.3 單元階次在二維分析中，對于位移場，相比Q4單元，Q8單元可以更準(zhǔn)確地描述單元內(nèi)應(yīng)力分布。對于相場，Q4單元采用的雙線性插值可以保證單元內(nèi)部任一點(diǎn)的值不超出節(jié)點(diǎn)值的范圍，