自適應有限元法對合金定向傾斜枝晶生長的相場法模擬

朱昶勝, 高宏偉, 馬芳蘭, 馮 力, 雷 鵬

(1. 蘭州理工大學 計算機與通信學院, 甘肅 蘭州 730050; 2. 蘭州理工大學 省部共建有色金屬先進加工與再利用國家重點實驗室, 甘肅 蘭州 730050; 3. 蘭州理工大學 網絡與信息中心, 甘肅 蘭州 730050)

近年來，對于凝固過程的界面形貌演化，不管是在實驗、理論還是數值模擬，人們已經做了大量的研究并獲得了重要的研究成果[1].一直以來，在傳統的研究方法中，合金的定向凝固是研究非平衡體系中由于枝晶形態(tài)不穩(wěn)定而引起競爭效應的標準范式[2].在現代先進的材料制備工藝中，定向凝固技術已被廣泛應用，其組織的演化規(guī)律對于實際生產具有重要的指導作用.與自由枝晶生長不同的是，定向凝固過程中晶體生長的方向性一方面取決于傳熱方向，另一方面取決于晶體學上的擇優(yōu)生長方向.例如在均勻過冷熔體中生長的晶體，由于向四周散熱速度是一樣的，因此其生長方向完全是由晶體界面能的各向異性來決定，晶體沿著界面能最低的方向能夠快速生長，最終形成等軸形狀的晶體.當熱流方向與晶體擇優(yōu)方向不平行時，枝晶的整體結構將會改變其豎直生長的特性，固液界面將會演化成非對稱的傾斜枝晶結構[3].Pettersen等[4]首先研究了定向凝固過程中冷卻速率對其枝晶生長的影響，其主要目的是研究合金生長的優(yōu)選方向.Xing等[5]研究了樹枝狀結構的不對稱形態(tài)以及錯向角對主晶距選擇機制的影響,將實驗結果與Ghmadh等[6]的相場模擬結果進行對比，表明三維模擬有助于獲得更準確的定量結果.Tourret和Karma[7]通過針對錯向角進行大量的二維相場模擬，得出通過主晶距的變化，從而導致熱梯度對枝晶的生長方向產生影響，且表面張力的各向異性強度也會影響合金陣列生長方向的選擇和枝晶形態(tài)的變化.

研究了退化海藻向傾斜枝晶轉變的機理及其生長動力學，特別研究了溫度梯度、抽拉速度和各向異性對凝固方式選擇的影響.通過將葉尖分裂不穩(wěn)定性與側枝之間發(fā)生的關系，討論了退化海藻傾斜枝晶轉變的機理以及其生長動力學.Chen等[9]將射線攝影的實時觀察結果與數值模擬結果進行比較，為Al-Cu合金凝固過程中樹狀枝晶或細胞結構向海藻結構的轉變提供了直接證據.其數值結果表明，當噪聲足夠大時，即使立方原鋁固體的晶體首選生長方向與熱梯度方向一致，枝晶到海藻的轉變也會發(fā)生.雖然上述的研究中海藻的結構是處于瞬態(tài)的，并且最終實現了穩(wěn)態(tài)樹狀枝晶生長的微觀結構，但足夠證明在實際的特殊情況下可以獲得對海藻結構有利的條件.Zhu等[10]重點研究了生長方向規(guī)定在豎直方向的枝晶生長動力學和定向凝固過程中的形態(tài)轉變，特別是分析了擾動和各向異性強度對凝固組織的影響，以及定向凝固過程中界面形態(tài)的生長機理和海藻組織的轉化機理.

由于對所建立相場模型的研究越來越復雜，在計算區(qū)域、計算量以及計算精度上的要求也逐漸增加，導致CPU和計算機內存空間的嚴重負荷，限制了模擬計算的精度，所以采用了自適應有限元方法進行模擬，根據有限元分析結果和事后誤差估計來自動地實現網格的加密或粗化，進行網格的自適應調整，以實現用盡可能少的網格數目獲得較高精度的解，從而進一步提高了計算效率和求解精度[19].本文采用自適應有限元的方法對Al-Cu(w(Cu)=4%)合金定向凝固過程中傾斜枝晶生長形貌進行深入研究，通過控制固液界面的溫度梯度、枝晶尖端的冷卻速率、主晶間距以及抽拉速度等參數的變化，得到不同的結晶形態(tài).

1 數學模型

1.1 相場方程和溶質擴散方程

對于Al-Cu(w(Cu)=4%)合金傾斜生長的定量模擬，由于界面寬度選擇的影響，相場法模型模擬的尺度都過小，實現起來會有諸多困難，所以將采用薄界面理論分析[11]，使界面寬度的選擇大于實際界面的寬度.但是當界面寬度增加后，會產生溶質陷落，這種情況會帶來一定程度上的計算誤差，為了減小誤差，將在溶質場方程中引入Karma等提出的反溶質陷落流[20].以下為所使用的相場模型，其中ψ為1時為固相，ψ為-1時為液相.相場方程為

(1)

通過求解溶質守恒方程的方式得到溶質過飽度U,溶質場方程為

(2)

內部的實際成分C為

C/C0=[1+1-k)U][(1+k)-(1-k)ψ]/2

(3)

方程(2)中的反溶質陷落流jat為

(4)

關于ψ的表達式q為

q(ψ)=(1-ψ)+k(1+ψ)Ds/Dl

(5)

1.2 數值模擬

表1 定向凝固過程模擬使用的計算參數

1.3 自適應有限元法

自適應有限元法是一種以誤差估計與自適應網格改進技術為核心，通過后驗誤差估計進行自動調整算法以改進求解過程的數值方法，具有高效率和高可靠性的特點.在自適應有限元方法中，采用各單元基函數的線性組合來逼近單元的實際解.這種方法的優(yōu)點是它可以將復雜的幾何圖形離散成單元的網格形狀，可以有多種類型，如三角形、矩形和四邊形.本文所研究的傾斜枝晶形貌較為復雜，將采用H型自適應有限元法完成網格的加密與粗化的自適應調整.H型自適應有限元法首先在一個較為粗糙的網格上進行求解，根據計算得到的數值結果判斷解是否具有要求的精度，然后結合事后的誤差估計結果對誤差較大的區(qū)域進行局部網格加密和粗化，接著計算機重新求解，直至滿足精度要求.自適應有限元算法求解相場模型的大致步驟如下：

1) 創(chuàng)建相場模型，在當前網格上對偏微分方程進行求解；

2) 設置相場和溶質場的初始值以及邊界條件；

3) 計算后驗誤差估計，根據后驗誤差估計對網格進行自適應調整；

4) 若計算結果不滿足要求精度，則重復以上步驟.

為了建立高效的自適應網格算法，使用幾何等級遺傳樹法來統一管理網格，完成網格的加密與粗化，具體細化模型如圖1，主要有以下操作.如圖1a將計算域中每個網格Γ0細化，得到一個均勻的三角形網格，連續(xù)細化n次后得到網格{Γn}，最后得到一個以Γ0為根節(jié)點的網格四叉樹.將每個網格單元記作τ，每個網格單元的指示子記作Eτ，如圖1b具體自適應算法為[12]：

1) 局部加密.一個三角形單元τ,如果Eτ>2N+αΞ，其中N和α為控制常數，Ξ為加密粗化偏差，那么該網格τ將被細化，其所有子網格的指示子將被設置為Eτ/2N+α；

4) 粗化.在幾何等級樹中由底部向上遍歷，對于任意結點，如果它的指示子滿足Eτ<2N+αΞ，那么它的所有子孫都將被刪除.如圖1c所示，由于計算結果的精度要求，需要在求解過程中對初始網格半正則化處理，這樣在加密的過程中會出現懸掛點A，B.接著對調整后的網格進行正則化處理，使得整個網格中的單元的每條邊都沒有懸掛點, M1和M2為雙生三角形.在經過正則化處理的網格上建立有限元空間，最終會實現偏微分方程的高效求解.具體的半正則化和正則化算法參考文獻[12].

圖1 幾何等級遺傳樹管理網格的細化模型Fig.1 Refinement model of geometric hierarchical genetic tree management grid

采用自適應有限元方法求解相場方程,這種方式既可以根據方程需求對區(qū)域進行合理形狀的網格劃分，又可以根據指定場的函數需求將節(jié)點進行排列，所以對區(qū)域形狀有良好的適應性.

2 結果和討論

2.1 相場模型的自適應有限元法求解分析

數值模擬方法的不同會對方程的處理、計算模型結果以及計算效率有很大影響.采用基于自適應有限元法求解相場模型，如圖2a為初始網格，圖2b為生長時間t=0.3 s時的枝晶生長形貌，圖2c為框處放大后的結果.可以發(fā)現，在固液界面區(qū)域A處網格密度相對稠密，而遠離界面的B、C、D區(qū)域網格相對稀疏.通過這種稀疏分布的網格方式，導致自適應有限元網格節(jié)點相對較少.隨著求解時間的推移，節(jié)點的數量隨即增加，在t=0.3 s時，節(jié)點數為2.79×105，遠小于采用均勻網格法計算時1.6×106的節(jié)點數量[9]，二者相比少了一個數量級，所以采用自適應有限元法大大降低了計算量，有效提高了計算效率.此外，采用非零元素存儲法存儲系數矩陣時，自適應有限元法對應系數矩陣所占用的存儲空間為M1～L，而均勻網格占用的存儲空間M1～L2，相比較，自適應有限元法所需要的存儲空間低了一個數量級，這對于較大規(guī)模的數值模擬計算效率會有提高.將自適應有限元法的CPU使用時間Ta和均勻網格法的CPU使用時間Tu進行對比，如圖3所示，在相同計算域LB大小的情況下，自適應有限元法的CPU使用時間遠小于均勻網格法，可知自適應有限元法的計算效率要遠高于均勻網格法.圖4為兩種數值模擬所得加速比，可以看出，隨著計算域的增大，自適應有限元法的計算效率越高.

圖2 采用自適應有限元法網格圖

圖3 采用自適應有限元法和均勻網格法CPU使用時間

圖4 均勻網格法和自適應有限元法的加速比

自適應有限元法的CPU使用時間Ta，均勻網格法的CPU使用時間Tu以及加速比Tu/Ta公式為

(6)

2.2 冷卻速率對界面形貌的影響

在定向凝固的研究過程中，首先模擬了熱流方向與晶體擇優(yōu)方向平行時，固液界面按豎直方向生長的海藻狀枝晶，如圖5所示，上半部分為模擬定向凝固的相場圖，下半部分為溶質場圖，這是枝晶定向凝固中最為普遍的例子.當熱流方向與晶體擇優(yōu)方向不平行時，設定枝晶優(yōu)選生長的錯向角φ0=π/6，冷卻速率分別設定為R=0.032 K/s，R=0.12 K/s，R=0.24 K/s，其他參數均見表1，圖6a～6c分別為3種不同冷卻速率下樹枝狀枝晶的瞬態(tài)演化，可以看出，在不同的冷卻速率下，枝晶的形貌有很大不同.

圖5 熱流方向與晶體擇優(yōu)方向平行時固液界面變?yōu)榘簇Q直方向生長的海藻狀枝晶

圖6 Al-Cu(w(Cu)=4%)合金在3種不同冷卻速率下的樹狀枝晶形態(tài)Fig.6 Dendritic morphology of Al-Cu(w(Cu)=4%) alloys at three different cooling rates

當R=0.032 K/s時，枝晶的生長方向更加傾向于晶體的優(yōu)選生長方向π/6，且出現了第三枝晶臂.當冷卻速率R=0.12 K/s時，枝晶的實際生長方向偏離了晶體的優(yōu)選生長方向，向溫度梯度方向靠攏，偏向角度為π/9.當冷卻速率達到R=0.24 K/s，枝晶的實際生長的角度更加靠近溫度梯度的方向，基本與溫度梯度角度相吻合，偏向角度為π/36，枝晶形貌呈現出較為規(guī)則的樹枝形貌.總體而言，枝晶的原始生長角度與冷卻速率有很大的關系.界面過冷度是控制枝晶生長的一個關鍵因素，隨著冷卻速率不斷增大，界面過冷度也隨之增大，從而控制了枝晶生長角度的偏差，枝晶的生長角度從晶體的優(yōu)選生長方向π/6逐漸向π/36偏移.

枝晶由細胞狀向樹狀演變的過程中，圖7a～7c為從胞狀向樹狀枝晶過渡和尖端分裂的進化序列，觀察圖7所表示的枝晶的初始狀態(tài)，當冷卻速率R=0.032 K/s時，此時界面過冷度很小，圖7a中黑色虛線標注上方，枝晶尖端處出現了明顯的分裂現象.當冷卻速率達到R=0.12 K/s時，從圖7b可以看出，枝晶尖端的分裂現象減少，只有個別尖端出現分裂.隨著冷卻速率的不斷增加，界面過冷度也隨之增大，枝晶中靠向溫度梯度方向的尖端受到青睞.如圖7c冷卻速率達到R=0.24 K/s時，枝晶的尖端沒有了明顯的分裂現象.Amoorezaei等[18]曾報道了界面過冷度和溶質之間的相互作用將會大幅度的影響尖端分裂，以上現象與他的理論是一致的.

圖7 定向凝固組織形態(tài)在冷卻速度分別為R=0.032 K/s，R=0.12 K/s，R=0.24 K/s的變化Fig.7 The changes of the directional solidification structure at the cooling rate of R=0.032 K/s, R=0.12 K/s, R=0.24 K/s

通過觀察圖7d溶質場中枝晶的生長形貌，可以看出，枝晶內部的溶質濃度要明顯低于枝晶間隙的溶質濃度，這是由于溶質偏析的影響(合金中各組成元素在結晶時分布不均勻導致).當冷卻速率較低時，固液界面的移動速率非常緩慢，溶質的擴散在凝固的過程中發(fā)揮著主要的作用，因此枝晶出現了不穩(wěn)定性，并且沿著高對稱的晶軸發(fā)生觸發(fā)，所以枝晶的生長不易發(fā)生改變，更傾向于晶體的原始方向.但是隨著冷卻速率的不斷增大，界面的過冷度也隨即快速增大，此時熱釋放主導了界面演化，從而影響枝晶生長方向和界面速度的變化，枝晶的角度由原始的擇優(yōu)晶體方向逐漸向溫度梯度方向轉變.

如圖6a～6c所示，冷卻速率從R=0.032 K/s轉變到R=0.24 K/s的同時，枝晶的生長方向也由原始的擇優(yōu)晶體方向φ0=π/6轉變到φ0=π/36.由此可見，冷卻速率很小時，具有控制枝晶生長有關鍵作用的過冷度也很小，此時溶質的擴散是完全的，所以熱流方向的枝晶沿著擇優(yōu)取向不斷生長.隨著冷卻速率的增加，尖端過冷度隨即不斷增大，熱萃取控制了取向角的偏差，使枝晶的生長角度向溫度梯度方向發(fā)生偏轉.

2.3 間距(Péclect數)對界面形貌的影響

對傾斜枝晶實際生長方向與熱流方向的夾角α和晶體擇優(yōu)取角φ0之間的變化規(guī)律進行模擬.由圖8可以看出，在Al-Cu(w(Cu)=4%)合金的定向凝固過程中，當晶體的擇優(yōu)生長方向與溫度梯度方向不一致時，枝晶的實際生長方向位于晶體的擇優(yōu)方向與溫度梯度方向之間.當晶體的擇優(yōu)生長方向分別為φ0=10°，20°，30°時，枝晶的實際生長方向為α=7.85°，13.92°，22.08°.這與枝晶的間距λ和枝晶的生長速度V有很大的關系，并且隨著將晶體擇優(yōu)取角的數值調整大，枝晶尖端的分裂更加劇烈.以上現象與Deschamps和Pocheau等[13-14]實驗給出的結果相吻合.

圖8 V=110 μm，λ=140 μm時枝晶的實際生長方向α隨晶體擇優(yōu)方向φ0的變化 Fig.8 When V=110 μm and λ=140 μm, the actual growth direction α of the dendrite varies with the crystal preferred direction φ0

枝晶之間的間距Péclect數，在理論上定義為Pe=λV/Dl.設定三個晶體擇優(yōu)方向，φ0分別為10°、20°和30°，并計算枝晶實際生長角度α與晶體擇優(yōu)取向角φ0之間的比值，通過記錄三個角度隨著Pe變化的α/φ0值，觀察其變化的情況.Pocheau等[15]通過大批實驗數據得出關于主晶間距Pe、晶體擇優(yōu)取向角φ0和枝晶實際生長角α之間的經典公式：

(7)

其中：a和b均為擬合參數，通過以上經典公式并擬合10°、20°和30°三個擇優(yōu)角度的計算數值繪制曲線，計算得出當φ0=10°時擬合參數為a=0.63,b=1.32；φ0=20°時擬合參數為a=0.49,b=1.66；φ0=30°時擬合參數為a=0.31，b=1.98，可知模擬結果基本吻合公式(7).通過圖9可以看出，隨著Pe的逐漸增大，10°、20°和30°對應的α/φ0值也逐漸增大，慢慢趨近于1.以上模擬結果表明，隨著枝晶間距Pe值的逐漸增大，α值逐漸趨近于φ0值，即隨著主晶間距和抽拉速度的增大，枝晶的實際生長角度也逐漸從溫度梯度方向向晶體擇優(yōu)方向偏轉，這與Pocheau等[14-15]給出的實驗結論一致.

圖9 φ0分別為30°，20°，10°時α/φ0隨Péclect數變化規(guī)律Fig.9 When φ0 is 30°, 20°, 10°, α/φ0varies with Pe number

3 結論

使用了自適應有限元法求解相場模型，解決了傳統方法解決相場模型計算量大、計算效率低等缺點.研究了Al-Cu(w(Cu)=4%)合金的定向傾斜枝晶在溫度梯度固定的情況下，冷卻速率、主晶間距對形貌演變的影響，結論如下：

1) 隨著冷卻速率的增大，熱萃取控制了取向角的偏差，枝晶錯向角向溫度梯度方向偏轉.不同的冷卻速率會影響樹狀枝晶尖端分裂的不穩(wěn)定性，導致枝晶形貌較為復雜，向海藻狀過渡.

2) 傾斜枝晶的實際生長方向位于溫度梯度與晶體擇優(yōu)取向角之間，并且隨著一次主晶間距Pe值和抽拉速度的逐漸增大，枝晶的實際生長角度也逐漸從溫度梯度方向向晶體擇優(yōu)方向偏轉.

3) 自適應有限元法相較于均勻網格法在網格節(jié)點數、網格占用的存儲空間和CPU的運行時間上均低了一個數量級，并且隨著計算域的增大，自適應有限元法的計算效率越高.