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

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

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

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

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

1 數(shù)學模型

1.1 相場方程和溶質(zhì)擴散方程

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

(1)

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

(2)

內(nèi)部的實際成分C為

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

(3)

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

(4)

關(guān)于ψ的表達式q為

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

(5)

1.2 數(shù)值模擬

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

1.3 自適應(yīng)有限元法

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

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

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

3) 計算后驗誤差估計，根據(jù)后驗誤差估計對網(wǎng)格進行自適應(yīng)調(diào)整；

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

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

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

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

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

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

2 結(jié)果和討論

2.1 相場模型的自適應(yīng)有限元法求解分析

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

圖2 采用自適應(yīng)有限元法網(wǎng)格圖

圖3 采用自適應(yīng)有限元法和均勻網(wǎng)格法CPU使用時間

圖4 均勻網(wǎng)格法和自適應(yīng)有限元法的加速比

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

(6)

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

在定向凝固的研究過程中，首先模擬了熱流方向與晶體擇優(yōu)方向平行時，固液界面按豎直方向生長的海藻狀枝晶，如圖5所示，上半部分為模擬定向凝固的相場圖，下半部分為溶質(zhì)場圖，這是枝晶定向凝固中最為普遍的例子.當熱流方向與晶體擇優(yōu)方向不平行時，設(shè)定枝晶優(yōu)選生長的錯向角φ0=π/6，冷卻速率分別設(shè)定為R=0.032 K/s，R=0.12 K/s，R=0.24 K/s，其他參數(shù)均見表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，且出現(xiàn)了第三枝晶臂.當冷卻速率R=0.12 K/s時，枝晶的實際生長方向偏離了晶體的優(yōu)選生長方向，向溫度梯度方向靠攏，偏向角度為π/9.當冷卻速率達到R=0.24 K/s，枝晶的實際生長的角度更加靠近溫度梯度的方向，基本與溫度梯度角度相吻合，偏向角度為π/36，枝晶形貌呈現(xiàn)出較為規(guī)則的樹枝形貌.總體而言，枝晶的原始生長角度與冷卻速率有很大的關(guān)系.界面過冷度是控制枝晶生長的一個關(guān)鍵因素，隨著冷卻速率不斷增大，界面過冷度也隨之增大，從而控制了枝晶生長角度的偏差，枝晶的生長角度從晶體的優(yōu)選生長方向π/6逐漸向π/36偏移.

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

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

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

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

對傾斜枝晶實際生長方向與熱流方向的夾角α和晶體擇優(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有很大的關(guān)系，并且隨著將晶體擇優(yōu)取角的數(shù)值調(diào)整大，枝晶尖端的分裂更加劇烈.以上現(xiàn)象與Deschamps和Pocheau等[13-14]實驗給出的結(jié)果相吻合.

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

(7)

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

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

3 結(jié)論

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

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

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

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