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

1、精選優(yōu)質文檔-傾情為你奉上 2012.01-2013.06工作總結2013.06.27-基于相場理論的固態(tài)相變模擬的有限元計算目錄1、相變機理相變是材料從高能量相轉變?yōu)榈湍芰肯?，總能量減小的過程。在材料中，各相能量的高低是由其內部狀態(tài)（如應力、溫度等）所決定的。例如對于NiTi形狀記憶合金，在低溫時，馬氏體為低能量相，當溫度從高溫降為低溫時，奧氏體轉變?yōu)轳R氏體；在高溫時，奧氏體為低能量相，當溫度從低溫升為高溫時，馬氏體轉變?yōu)閵W氏體。相的變體是指晶體結構相同，取向不同的同一種相1。例如對于NiTi形狀記憶合金，奧氏體的晶體結構為高對稱的立方體（ high symmetry&#

2、160;cubic），馬氏體的晶體結構為四方體，在不同的朝向下，共有三種變體，如圖1-1。圖1-1 奧氏體與馬氏體的晶體結構示意圖在一定的內部狀態(tài)下，體系生成各變體后，總能量的大小是不同的。體系會自發(fā)地選擇生成使得體系能量降低最多的變體。如圖1-2所示，各變體首先在使得體系能量降低最多的區(qū)域形成，后生長形成特定的位向關系。圖1-2 變體產生及位向關系形成示意圖2、相場理論2.1相場理論概述相場理論是一種用來模擬相變微觀組織的演化過程的理論，它的主要特點是以Ginzburg-Laudau相變理論為基礎，基于擴散界面模型，引入連續(xù)變量來描述新舊兩相，用微分方程描述系統(tǒng)自由能的變化。 相場理論的特點

3、具體如下：（1）離散變量的連續(xù)化：在相變中，各相或者相的各變體是離散的物理量，各相或者相的各變體之間不存在連續(xù)變化的過程。相場模型將這些離散變量連續(xù)化，用序參量表示，便于在數值上獲取導數等信息。例如，對于新相含有n個變體的相變，采用序參量1=2=n=0表示體系此處為母相所占有；序參量p=1或-1且k=0（k=1,2.p-1,p+1,n）表示此處對應為新相的第p個變體。由于連續(xù)化的處理，各相或者各變體間的界面是彌散界面（diffuse interface），而不是尖銳界面（Sharp interface），如圖2-1。（2）微分方程：建立相場內的平衡微分方程，描述各相或各變體的演化。（3）界面能

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

5、關，同一相中的各變體具有相同的體積化學能密度。兩相間的體積化學能密度之差促使相變從能量高的相向能量低的相轉變。在相變中，體積化學自由能具有以下兩點性質（1）在母相處，體積化學自由能密度為0，在新相處，體積化學自由能密度為-fchem；（2）當某點為母相或新相時，該點的體積化學能驅動力為0，處于穩(wěn)定態(tài)或亞穩(wěn)定態(tài)。這兩點性質要求體積化學自由能密度函數滿足： （2-4）在不同的文獻中，體積化學能密度函數chem主要有以下兩種形式：（1）2-4-6次多項式2 （2-5）A1，A2，A3為體系膨脹系數，一般取正數；n為變體總數；序參量1=2=n=0表示體系此處為母相所占有；序參量p=±1且k=

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

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

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

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

10、變三場耦合示意圖在實際模擬中，一般忽略相變過程中的發(fā)熱、吸熱，假定在相變過程中，溫度不發(fā)生改變，原三場耦合的問題簡化為兩場耦合，如圖2-7所示。圖2-7 簡化后的馬氏體相變兩場耦合示意圖對于相場與位移場的耦合，本研究計算采用順序耦合方法求解。具體求解策略是按照順序進行位移場和相場的分析，采用上一時刻的相場分布，先求解位移場，得到上一時刻位移場分布，將位移場所得結果應用到相場求解中，再求解這一時刻的相場分布。3、有限元計算3.1有限元的離散化3.1.1 相場控制微分方程的化簡對體積化學能、界面能和彈性應變能三種能量之和的總能量泛函求變分，可得： （3-1）將式（3-1）代入式（2-1）的Ginz

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

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

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

14、流程對于位移場和相場的耦合問題，采用的求解策略是在根據t時刻序參量分布值，先在t時刻先求解位移場，再根據得到t時刻的應力分布值，求解t+t時刻相場，得到t+t時刻的序參量分布值，具體流程如圖2-8所示： 圖3-1 相變模擬計算流程圖4、已完成工作總結4.1相場模型的研究總結馬氏體相變的相場模型包含三種能量：體積化學能、界面能、彈性應變能，三種能量相互作用，改變其中任何一種，都將影響到最終結果。4.1.1 彈性應變能的影響根據彈性應變能的表達式（2-12）與彈性應變能驅動力的表達式（2-14），彈性應變能的組成分為三部分：（1）應力：應力的大小影響著相變的演化速率和演化方向。當應力的方向與相變應

15、變一致時，應力促進相變演化；反之，抑制相變演化。從應力產生的原因來分，應力可分為外荷載作用下引起的應力和相變應變引起的應力。由于相變應變引起的應力一定與相變應變相反，因此，相變應變引起的應力總是抑制相變進一步演化。（2）彈性模量：在線彈性下，彈性模量與相變應變引起的應力大小成線性關系，因此，彈性模量越大，由相變引起的應力大小越大，相變的演化越慢，反之，相變演化越快。（算例驗證）（3）泊松比：泊松比的大小與相變應變引起的應力大小相關，其影響與彈性模量相似，只是對不同的應變分量的影響程度不同。（算例驗證）（4）相變應變插值函數：相變應變的插值函數影響彈性應變能在序參量值介于0-1之間的大小，進而影

16、響彈性應變能驅動力的大小，影響相變演化進程與最終結果。在選取相變應變插值函數時，需注意以下準則：I.必須滿足：為保證當序參量為0時，相變應變?yōu)?；當序參量為1時，相變應變?yōu)?，彈性應變能插值函?)必須滿足： （4-1） II.非必須滿足：為了使序參量位于0或者1時，材料處于穩(wěn)定或者亞穩(wěn)定態(tài)，其能量驅動力應為0，由于體積化學能驅動力在序參量位于0或者1時為0，彈性應變能驅動力也應為0。因此，相變應變插值函數()須滿足： （4-2）按照以上選取準則，已進行過對比的相變插值函數如下： 表達式形式一次導數表達式線性二次三次2-3次2-4次2-3-4次其中，在材料內部無缺陷且無外力作用下，對比線性、二次

17、、三次插值函數，線性的結果好于二次和三次，原因是隨著插值函數次數增大，在序參量靠近1的部分彈性應變能驅動力增大，最終穩(wěn)定時，序參量離1的平均距離變遠。與此同時 ，體積化學能驅動力在1附近的增大幅度大于彈性應變能驅動力，因此，應變能驅動力在總驅動力的比例降低，最終形成的位向關系變差。（算例驗證）2-3次、2-4次、2-3-4次是在材料內部無缺陷且無外力作用時，結果無明顯差別。對于材料受到較大外力作用，或因為自身缺陷引起應力集中，造成相變區(qū)域為大應力區(qū)域的情況，還未進行算例分析。4.1.2 體積化學能的影響在溫度不變的情況下，體積化學能只與自身的插值函數相關，與應力大小無關。體積化學能是否促進相變

18、演化取決于在選取的溫度下，母相的體積化學能是否較大。4.1.3 界面能的影響由于相場模型本身是將離散的相通過連續(xù)的序參量值來表示，界面能采用梯度能的描述方式，與實際情況并不相符。實際各相或者各變體間是否存在界面能，界面能的大小與什么參數相關，還需進一步閱讀文獻。單就相場模型中的界面能，可以從以下角度分析：（1）數學表達式：界面能表示由于相場不均勻分布所產生的能量。在相變中，界面能的釋放將將促使相場分布的均勻化。（2）模型中的作用：通過梯度描述的界面能，在相場模型中，建立相變區(qū)域中相鄰各點的聯系，使得新相或者變體能夠通過“擴散”的方式生長。（3）界面能系數對模擬結果影響：界面能系數的大小影響相變

19、模擬結果中各變體的平均尺寸和彌散界面的寬度，界面能系數越大，各變體的平均尺寸越大，彌散界面的寬度越大。（算例驗證）4.1.4 動力學系數在式（2-1）的Ginzburg-Landau方程中，動力學系數控制著相變演化的速率。由于相變演化本身還與位移場相互耦合，相變演化的速率影響慣性力的大小，決定位移場分析是采用動力分析還是靜力分析。4.1.5 隨機噪聲相場模型是一個描述相變演化過程的模型，無法描述相變的形核機理，通常要求其化學能驅動力與彈性應變能驅動力為0，在新相形核前，體系處于亞穩(wěn)定狀態(tài)。引入隨機的噪聲，主要目的是引起序參量值的微小改變，從而使得相變能夠演化下去。因此，在隨機噪聲很小的情況下，

20、與空間坐標獨立的隨機噪聲抽樣方式和持續(xù)時間不影響相變結果的位向關系和變體的平均尺寸。（算例驗證）4.1.6 穩(wěn)定結果的非均勻性根據相場的控制方程，當結果穩(wěn)定時，原方程可化為： （4-11） （4-12）假設變體中間區(qū)域，序參量為0或者1的均勻分布，界面能為0，方程進一步簡化為： （4-13）在變體邊界，序參量不均勻分布，存在界面能，彈性應變能驅動力與體積化學能驅動力之和不為0： （4-14）由于式（4-13）與式（4-14）互斥，滿足式（4-13）的序參量值必然不滿足式（4-14）。變體中間區(qū)域的序參量值必然不等于變體邊界上的序參量值，因此，中間區(qū)域的序參量值與邊界上的序參量值存在一定的梯度，

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

22、部區(qū)域結果的不同，但所有結果整體上都符合位向關系。同時在合理步長下，步長較大，雖然計算結果收斂時對應的相變實際時間較長，但在模擬中所需要的計算步數較少，需要的數值計算時間較短。4.2.2 單元劃分在位移場中，提高網格精度，可以更準確地描述區(qū)域內應力分布；在相場中，提高網格精度，可以更準確地描述區(qū)域內序參量的梯度分布，更好地描述界面能的大小，計算結果更精確合理。同時，提高網格精度，在圖形顯示上，可以更好地描述變體間的邊界，避免出現鋸齒型區(qū)域。4.2.3 單元階次在二維分析中，對于位移場，相比Q4單元，Q8單元可以更準確地描述單元內應力分布。對于相場，Q4單元采用的雙線性插值可以保證單元內部任一點的值不超出節(jié)點值的范圍，而Q8單元采用的二次插值則不能。當節(jié)點序參量值落在0-1之間時，Q8單元插值給出的