自適應(yīng)有限元方法凝固問題模擬的網(wǎng)格各向異性研究

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

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

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

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

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

1 相場(chǎng)模型及求解

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

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

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

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

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

從上式消去 x，y，得：

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

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

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

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

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

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

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

4）求解相場(chǎng)方程（1），并計(jì)算固/液界面速度V；

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

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

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

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

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

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

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

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

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

圖2 實(shí)際各向異性系數(shù)ε4＇與空間步長(zhǎng)Δ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é)計(jì)算方法求得實(shí)際各向異性系數(shù)ε4＇。圖3為實(shí)際各向異性系數(shù)ε4＇與晶體半徑R0關(guān)系。在有限差分方法中，R0較小時(shí)二者偏差較大，隨著 R0的增加，ε4＇并逐漸趨近 ε4［5］；而在自適應(yīng)網(wǎng)格中對(duì)不同R0時(shí)ε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é)計(jì)算方法求得不同各向異性值對(duì)應(yīng)的實(shí)際各向異性 ε4＇。 表 1為不同 ε4取值計(jì)算得到的實(shí)際 ε4＇和網(wǎng)格各向異性 ε4＇-ε4。 在 ε4=0.01 時(shí)， （ε4＇-ε4） /ε4大小為18%，此時(shí)網(wǎng)格各向異性的影響必須加以考慮； 隨 ε4增加 （ε4＇-ε4） /ε4＇迅速減小， 當(dāng) ε4＞0.02 時(shí)， （ε4＇-ε4） /ε4低于 5%， 網(wǎng)格各向異性的影響已經(jīng)可以忽略。

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

4 結(jié) 論

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

［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á)文，楊根倉(cāng)，王錦程.過(guò)冷純物質(zhì)凝固的相場(chǎng)法的自適應(yīng)有限元模擬 ［J］.自然科學(xué)進(jìn)展，2006，16：1009-1015.