趙達(dá)文，李秋書，郝維新

（太原科技大學(xué)材料學(xué)院，山西 太原 030024）

在凝固過程中，液-固界面上的一些物理現(xiàn)象如界面能、界面動力學(xué)具有不同程度的各向異性。理論和實(shí)驗(yàn)研究均表明，各向異性對晶體生長形態(tài)和生長行為有著決定性的影響［1-3］。通過數(shù)值方法模擬凝固組織演化是凝固學(xué)研究的熱點(diǎn)之一，其中液-固界面追蹤是數(shù)值求解的難點(diǎn)。近年來出現(xiàn)的相場模型通過引入序參量場，把凝固模型從尖銳界面模型轉(zhuǎn)化為數(shù)學(xué)上等價的相場模型，從而避免了數(shù)值求解中追蹤界面的困難，在凝固組織模擬中得到了廣泛應(yīng)用［4］。

相場模型由一個（或組）序參量場、溫度場、濃度場的偏微分方程組成，通常通過有限差分或者自適應(yīng)有限元方法進(jìn)行數(shù)值求解。數(shù)值求解過程中首先要把求解區(qū)域劃分為多個單元組成的網(wǎng)格，這在計算中往往引入額外的各向異性，一般稱之為網(wǎng)格各向異性［5］。網(wǎng)格各向異性的存在限制了計算所采用空間步長的大小，同時影響了計算結(jié)果的可靠性。因此，必須通過計算確定網(wǎng)格各向異性與計算中各個因素的關(guān)系。

由于有限元和有限差分方法微分方程的離散、網(wǎng)格剖分等方面不同，因此二者的網(wǎng)格各向異性行為也存在差異。對于有限差分法的網(wǎng)格各向異性已有研究報道［5］，而對于自適應(yīng)有限元法求解時的網(wǎng)格各向異性尚缺乏研究。因此，本文中通過計算晶體的平衡形態(tài)，對自適應(yīng)有限元求解相場模型時的網(wǎng)格各向異性進(jìn)行研究。

1 相場模型及求解

這里以純物質(zhì)凝固為研究對象，相場模型為［5］：

方程（1）是序參量方程，（2）是溫度場方程。其中，φ為序參量，u為無量綱溫度D為熱擴(kuò)散系數(shù)，λ為耦合常數(shù)，τ（θ）是界面法線方向與x軸夾角。（θ）為原子弛豫時間，W（θ）為彌散界面厚度。τ（θ）和W（θ）分別對應(yīng)界面能和界面動力學(xué)效應(yīng)，二者都是與界面法向角θ相關(guān)的函數(shù)，其具體形式取決于各向異性函數(shù)的設(shè)定。

這里采用自適應(yīng)有限元方法求解相場模型，并選擇四分樹作為動態(tài)自適應(yīng)網(wǎng)格數(shù)據(jù)結(jié)構(gòu)來建立、調(diào)整自適應(yīng)網(wǎng)格。采用Galerkin加權(quán)余量法對相場控制方程進(jìn)行空間域離散，分別采用向前差分和中心差分格式對序參量場和溫度場進(jìn)行時間域離散。采用非零元素按行壓縮方法對其進(jìn)行存儲，同時以ICCG方法迭代求解。具體求解過程見［6］。

2 晶體平衡形態(tài)數(shù)值算法

處于液-固平衡狀態(tài)的晶體，其外形由各向異性界面能函數(shù)γ（θ）決定。假設(shè)晶體中心位于坐標(biāo)系原點(diǎn)，其外形由以下參數(shù)方程給出:

從上式消去 x，y，得：

其中，R為界面與原點(diǎn)距離，γθ為界面能函數(shù)對θ的導(dǎo)數(shù)。

把晶體平衡形態(tài)模擬結(jié)果代入（3）式，即可得出實(shí)際各向異性以及網(wǎng)格各向異性。

由于晶體處于液-固平衡狀態(tài)，液相和固相溫度相同，因此計算中只需要求解相場方程（1）。具體計算過程為：

1）設(shè)定輸入各向異性系數(shù)；

(4)讓學(xué)生組隊進(jìn)行小組討論。在課程進(jìn)行到一定階段時,組織學(xué)生就某一領(lǐng)域的新進(jìn)展,某一課程重點(diǎn)難點(diǎn)內(nèi)容,或某個科研成果的案例進(jìn)行小組討論,全部過程使用英文,激發(fā)學(xué)生的學(xué)習(xí)熱情,也提高學(xué)生用英文進(jìn)行學(xué)術(shù)交流的能力。例如,對2016年諾貝爾化學(xué)獎授予的分子機(jī)器研究成果分組討論,激發(fā)了學(xué)生發(fā)散的充滿想象力和創(chuàng)新性的點(diǎn)子和思路,課堂氣氛熱烈有趣,充滿生機(jī)。

2）設(shè)定初始圓形晶核半徑；

3）設(shè)定計算區(qū)域內(nèi)所有節(jié)點(diǎn)溫度為-Δ；

4）求解相場方程（1），并計算固/液界面速度V；

5） 如果 V＞0，則減小過冷度 Δ 值；如果 V＜0，則增大過冷度Δ值；

6）判斷V是否小于給定的收斂速度Vcon。如果V＜Vcon，則輸出平衡界面形態(tài)；否則返回繼續(xù)執(zhí)行步驟（4）；

7）由平衡界面形態(tài)依據(jù)（3）式求出實(shí)際的各向異性值。

計算表明上述算法具有良好的收斂性。圖1a）是具有四重對稱界面能的晶體形態(tài)演化過程。在迭代過程中，〈110〉晶向界面向內(nèi)收縮，〈100〉晶向界面則向外生長。經(jīng)過約500次迭代后，晶體形狀收斂于平衡形態(tài)。圖1b）為計算過程中x軸方向的界面速度Vx變化情況，可見隨迭代進(jìn)行Vx迅速趨于0。

圖1 晶體平衡形態(tài)相場模擬

3 網(wǎng)格各向異性與影響因素

計算中許多參數(shù)的取值都可能影響網(wǎng)格各向異性，其中最主要的影響因素是自適應(yīng)網(wǎng)格中最細(xì)的單元尺寸Δx/W0以及晶核半徑R0取值。這里以四重對稱性的界面能為例進(jìn)行討論，其界面能函數(shù)為 γ（θ）=γ0（1+ε4cos4θ）。 計算區(qū)域大小為 409.6W0×409.6 W0。

3.1 Δx/W0對網(wǎng)格各向異性的影響

設(shè)置 R0=20，ε4=0.05， Δx/W0分別為 0.4，0.8，1.2，1.6，分別計算出實(shí)際各向異性系數(shù)ε4＇。圖2為有效各向異性系數(shù)ε4＇。與空間步長 Δx/W0關(guān)系。Δx/W0=0.4，0.8時，二者符合的很好；隨著 Δx/W0增大，ε4＇和 ε4的差值（網(wǎng)格各向異性）逐漸增大。

圖2 實(shí)際各向異性系數(shù)ε4＇與空間步長Δx/W0關(guān)系

3.2 晶體半徑的影響

設(shè)置 ε4=0.05，Δx/W0=0.8，晶體半徑 R0分別為10，15，20，25，30，35，40，45，50，由上節(jié)計算方法求得實(shí)際各向異性系數(shù)ε4＇。圖3為實(shí)際各向異性系數(shù)ε4＇與晶體半徑R0關(guān)系。在有限差分方法中，R0較小時二者偏差較大，隨著 R0的增加，ε4＇并逐漸趨近 ε4［5］；而在自適應(yīng)網(wǎng)格中對不同R0時ε4＇與ε4非常接近。可見，采用自適應(yīng)有限元方法可以降低網(wǎng)格各向異性。

圖3 有效各向異性系數(shù)ε4＇與初始晶體半徑R0關(guān)系

3.3 網(wǎng)格各向異性與各向異性系數(shù)關(guān)系

設(shè)置 R0=20， Δx/W0=0.8，ε4為 0.00～0.06，由上節(jié)計算方法求得不同各向異性值對應(yīng)的實(shí)際各向異性 ε4＇。 表 1為不同 ε4取值計算得到的實(shí)際 ε4＇和網(wǎng)格各向異性 ε4＇-ε4。 在 ε4=0.01 時， （ε4＇-ε4） /ε4大小為18%，此時網(wǎng)格各向異性的影響必須加以考慮； 隨 ε4增加 （ε4＇-ε4） /ε4＇迅速減小， 當(dāng) ε4＞0.02 時， （ε4＇-ε4） /ε4低于 5%， 網(wǎng)格各向異性的影響已經(jīng)可以忽略。

表1 計算輸入各向異性系數(shù)和網(wǎng)格各向異性

4 結(jié) 論

通過相場模型模擬晶體平衡狀態(tài)，標(biāo)定了自適應(yīng)網(wǎng)格各向異性以及對相場模擬的影響。計算結(jié)果表明，隨著空間步長增大，實(shí)際各向異性系數(shù)和輸入的各向異性差值增大；網(wǎng)格各向異性與晶體大小無關(guān)；低各向異性下網(wǎng)格各向異性對計算的影響較大，當(dāng)各向異性大于0.03時可以忽略。

［1］Kessler D A，Levine H.Velocity selection in dendritic growth ［J］.Phys Rev B，1986，33：7867-7870.

［2］Amar M B，Pomeau Y.Theory of Dendritic Growth in a Weakly Undercooled Melt［J］.Euro Phys Lett，1986（2）：307-314.

［3］Langer J S.Existence of needle crystals in local models of solidification ［J］.Phys Rev A，1986，33：435-441.

［4］Boettinger W J，Warren J A，Karma A.Phase-field simulation of solidification ［J］.Annu Rev Mater Res，2002，32：163-194.

［5］Karma A，Rappel W J.Quantitative phase-field modeling of dendritic growth in two and three dimensions ［J］.Phys Rev E，1998，57：4323-4349.

［6］趙達(dá)文，楊根倉，王錦程.過冷純物質(zhì)凝固的相場法的自適應(yīng)有限元模擬 ［J］.自然科學(xué)進(jìn)展，2006，16：1009-1015.