自適應有限元方法凝固問題模擬的網格各向異性研究

趙達文，李秋書，郝維新

（太原科技大學材料學院，山西 太原 030024）

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

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

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

1 相場模型及求解

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

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

這里采用自適應有限元方法求解相場模型，并選擇四分樹作為動態自適應網格數據結構來建立、調整自適應網格。采用Galerkin加權余量法對相場控制方程進行空間域離散，分別采用向前差分和中心差分格式對序參量場和溫度場進行時間域離散。采用非零元素按行壓縮方法對其進行存儲，同時以ICCG方法迭代求解。具體求解過程見［6］。

2 晶體平衡形態數值算法

處于液-固平衡狀態的晶體，其外形由各向異性界面能函數γ（θ）決定。假設晶體中心位于坐標系原點，其外形由以下參數方程給出:

從上式消去 x，y，得：

其中，R為界面與原點距離，γθ為界面能函數對θ的導數。

把晶體平衡形態模擬結果代入（3）式，即可得出實際各向異性以及網格各向異性。

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

1）設定輸入各向異性系數；

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

2）設定初始圓形晶核半徑；

3）設定計算區域內所有節點溫度為-Δ；

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

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

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

7）由平衡界面形態依據（3）式求出實際的各向異性值。

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

圖1 晶體平衡形態相場模擬

3 網格各向異性與影響因素

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

3.1 Δx/W0對網格各向異性的影響

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

圖2 實際各向異性系數ε4＇與空間步長Δx/W0關系

3.2 晶體半徑的影響

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

圖3 有效各向異性系數ε4＇與初始晶體半徑R0關系

3.3 網格各向異性與各向異性系數關系

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

表1 計算輸入各向異性系數和網格各向異性

4 結 論

通過相場模型模擬晶體平衡狀態，標定了自適應網格各向異性以及對相場模擬的影響。計算結果表明，隨著空間步長增大，實際各向異性系數和輸入的各向異性差值增大；網格各向異性與晶體大小無關；低各向異性下網格各向異性對計算的影響較大，當各向異性大于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］趙達文，楊根倉，王錦程.過冷純物質凝固的相場法的自適應有限元模擬 ［J］.自然科學進展，2006，16：1009-1015.