趙達(dá)文

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

凝固是典型的一級(jí)相變過(guò)程。在純物質(zhì)的凝固過(guò)程中，首先通過(guò)起伏形成晶核，其后沿著各個(gè)方向生長(zhǎng)從而形成各種凝固組織。理論分析表明[1，2]，這些生長(zhǎng)方向與液固界面上的各向異性密切相關(guān)。各向異性是指液固界面的一些物理性質(zhì)是晶向的函數(shù)，常見(jiàn)的有界面能各向異性和動(dòng)力學(xué)各向異性。當(dāng)各向異性效應(yīng)較顯著時(shí)，晶體沿著各向異性決定的方向以枝晶方式生長(zhǎng)。在無(wú)、弱各向異性條件下不存在優(yōu)勢(shì)生長(zhǎng)方向，因此各個(gè)晶向上都以相同速度生長(zhǎng)，當(dāng)界面曲率半徑超過(guò)臨界值后，會(huì)在擾動(dòng)的作用下發(fā)生失穩(wěn)現(xiàn)象[3]。弱各向異性條件下界面形態(tài)演化以及生長(zhǎng)方式等都是凝固學(xué)的重要內(nèi)容。對(duì)于顯著各向異性條件下生長(zhǎng)行為已有眾多研究報(bào)道，而對(duì)弱各向異性下生長(zhǎng)缺乏相應(yīng)的研究。

由于數(shù)學(xué)上的困難，凝固理論無(wú)法預(yù)測(cè)弱各向異性下界面形態(tài)演化過(guò)程；同時(shí)，由于界面各向異性值很難精確測(cè)量、各向異性值無(wú)法自由調(diào)整等原因，導(dǎo)致實(shí)驗(yàn)研究也存在較大困難，所以通常用數(shù)值模擬方法進(jìn)行研究。相場(chǎng)模型通過(guò)引入序參量場(chǎng)從而避免了顯式追蹤界面的困難，在凝固組織模擬中得到了廣泛應(yīng)用[5，6]。

本文以自適應(yīng)有限元方法求解相場(chǎng)模型，分別對(duì)界面能和動(dòng)力學(xué)弱各向異性條件下的晶體生長(zhǎng)過(guò)程進(jìn)行模擬，從而對(duì)其生長(zhǎng)規(guī)律進(jìn)行深入研究。

1 相場(chǎng)模型

這里采用Karma等提出的純物質(zhì)凝固相場(chǎng)模型[4]。該模型優(yōu)點(diǎn)是數(shù)值求解時(shí)保證精度的前提下可以采用較大的界面厚度，并且所引入的網(wǎng)格各向異性遠(yuǎn)小于傳統(tǒng)相場(chǎng)模型。其控制方程由序參量方程和溫度場(chǎng)控制方程組成：

式中：為序參量；u≡（T-Tm）/（L/cp）是無(wú)量綱溫度；T是熱力學(xué)溫度；Tm是熔點(diǎn)；L是結(jié)晶潛熱；cp是等壓熱容；D代表熱擴(kuò)散系數(shù)；λ代表耦合常數(shù)；θ為界面法線與x軸夾角。τ（θ）為原子弛豫時(shí)間，W（θ）為彌散界面厚度，二者分別對(duì)應(yīng)界面能和動(dòng)力學(xué)效應(yīng)，在不同各向異性條件下有不同形式，具體形式見(jiàn)模擬結(jié)果部分。

2 求解方法與計(jì)算參數(shù)

采用自適應(yīng)有限元方法求解相場(chǎng)模型，使用自適應(yīng)網(wǎng)格可以在計(jì)算精度不變前提下使求解計(jì)算量降低一階。在空間域上采用伽遼金加權(quán)余量法對(duì)方程（1）、（2）進(jìn)行離散，在時(shí)間域上分別采用向前差分和C-N格式對(duì)二者進(jìn)行離散，采用ICCG方法求解線性方程。

計(jì)算中所采用矩形單元會(huì)引入額外的網(wǎng)格各向異性，當(dāng)各向異性值較弱時(shí)，網(wǎng)格各向異性的影響必須予以考慮。計(jì)算中采用如下措施來(lái)降低網(wǎng)格各向異性影響[5～7]：1）采用薄界面漸近的相場(chǎng)模型；2）采用較小的空間步長(zhǎng)；3）通過(guò)計(jì)算各種參數(shù)下的網(wǎng)格各向異性，對(duì)所輸入各向異性值進(jìn)行校準(zhǔn)。

3 模擬結(jié)果與討論

如上所述，液固界面上常見(jiàn)各向異性有界面能和動(dòng)力學(xué)各向異性。在相場(chǎng)模型中，二者分別體現(xiàn)在界面厚度W（θ）和原子的弛豫時(shí)間τ（θ）上。本節(jié)中分別對(duì)弱界面能、動(dòng)力學(xué)各向異性作用下的晶體生長(zhǎng)進(jìn)行模擬。計(jì)算中用到的參數(shù)為W0=1.0，τ0=1.0，D=1.0，d0/W0=0.577，△x/W0=0.4，△t/τ0=0.008，無(wú)量綱過(guò)冷度△=0.65。計(jì)算區(qū)域大小為409.6W0×409.6W0，初始晶核半徑R0=10，計(jì)算區(qū)域邊界為絕熱邊界條件。

3.1 弱動(dòng)力學(xué)各向異性控制下晶體生長(zhǎng)

界面動(dòng)力學(xué)是指液相原子穿過(guò)界面最終吸附在固相的過(guò)程，在此過(guò)程中消耗的驅(qū)動(dòng)力即為動(dòng)力學(xué)過(guò)冷。這里設(shè)定界面厚度W（θ）各向同性、弛豫時(shí)間τ（θ）具有四重對(duì)稱性

來(lái)模擬界面動(dòng)力學(xué)各向異性作用下的生長(zhǎng)過(guò)程，其中耦合系數(shù)λ=5.77.

模擬結(jié)果表明，當(dāng)動(dòng)力學(xué)各向異性εk≤0.02時(shí)，固相不再以枝晶方式生長(zhǎng)。圖1a）為εk=0.02時(shí)的液固界面形態(tài)演變過(guò)程，圖中輪廓線為間隔1萬(wàn)時(shí)間步長(zhǎng)的液固界面；圖1b）為x軸方向界面前沿的生長(zhǎng)速度曲線。在7萬(wàn)步左右時(shí)，界面前沿的曲率半徑達(dá)到極大值，在M-S不穩(wěn)定性的調(diào)制下，界面開(kāi)始發(fā)生分叉（最內(nèi)側(cè)粗實(shí)線）；此時(shí)界面生長(zhǎng)速度位于第一個(gè)波谷位置。緊接著產(chǎn)生分叉。在12萬(wàn)步左右，界面分叉過(guò)程已經(jīng)完成，此時(shí)界面前沿的曲率半徑達(dá)到極小值（圖1a）次內(nèi)側(cè)粗實(shí)線）；而V達(dá)到第一個(gè)波峰位置的極大值。隨著生長(zhǎng)過(guò)程的進(jìn)行，在16萬(wàn)步左右界面曲率半徑再次達(dá)到極大值（圖1a）中次外側(cè)粗實(shí)線），開(kāi)始第二次分叉過(guò)程；此時(shí)生長(zhǎng)速度V為極小值（圖1b）中的第二個(gè)波谷位置）。在21萬(wàn)步左右，第二次分叉過(guò)程完成，界面曲率半徑再次取極小值（圖1a）外側(cè)粗實(shí)線），生長(zhǎng)速度到達(dá)圖1b）中的第二個(gè)波峰位置。可見(jiàn)，隨著時(shí)間的推移生長(zhǎng)前沿重復(fù)性的發(fā)生分叉現(xiàn)象，生長(zhǎng)速度隨時(shí)間的推移而波動(dòng)。

3.2 弱各向異性界面能控制下晶體生長(zhǎng)

為了單獨(dú)考察界面能各向異性的影響，這里設(shè)定W（θ）=W0（1+εccos4θ），τ（θ）=τ0（1+εccos4θ）2，并通過(guò)選擇耦合常數(shù)λ來(lái)消除界面動(dòng)力學(xué)效應(yīng)。當(dāng)界面能各向異性系數(shù)εc取≥0.02時(shí)，固相以枝晶方式生長(zhǎng)；而當(dāng)值εc取0.00001時(shí)，晶體不再以枝晶方式生長(zhǎng)。圖2a）是εc=0.01液固界面演化過(guò)程，圖中相鄰界面時(shí)間間隔為1萬(wàn)步；圖2b）為x軸附近界面前沿的生長(zhǎng)速度。在生長(zhǎng)過(guò)程中界面前沿周期性的發(fā)生分叉，相應(yīng)的生長(zhǎng)速度持續(xù)性的波動(dòng)，同樣不存在穩(wěn)定生長(zhǎng)狀態(tài)。

綜合以上模擬結(jié)果，可以發(fā)現(xiàn)在弱各向異性下具有相同的生長(zhǎng)方式：界面前沿總是持續(xù)性分叉、生長(zhǎng)速度周期性波動(dòng)、不存在與枝晶生長(zhǎng)類似的穩(wěn)定生長(zhǎng)狀態(tài)。究其原因，是由于在生長(zhǎng)過(guò)程中由于不存在優(yōu)勢(shì)生長(zhǎng)方向，使得界面每個(gè)位置以相近或相同的速度向各個(gè)方向生長(zhǎng)；當(dāng)界面前沿曲率半徑大于臨界半徑值時(shí)，在M-S不穩(wěn)定性的調(diào)制下界面失穩(wěn)發(fā)生分叉[3]；當(dāng)前沿的曲率半徑再度超過(guò)臨界半徑時(shí)，生長(zhǎng)前沿再次發(fā)生分叉。通常將這種生長(zhǎng)方式稱之為分形生長(zhǎng)[8]。注意分形生長(zhǎng)與有限空間內(nèi)的doublon生長(zhǎng)模式不同[9]，雖然二者的生長(zhǎng)形態(tài)類似，但是分形生長(zhǎng)不存在穩(wěn)定生長(zhǎng)形態(tài)和生長(zhǎng)速度，而doublon存在穩(wěn)定生長(zhǎng)狀態(tài)。

4 結(jié) 論

使用相場(chǎng)模型對(duì)弱各向異性條件下晶體生長(zhǎng)進(jìn)行了模擬。模擬結(jié)果表明，在各向異性較弱時(shí)固相均以分形方式生長(zhǎng)。由于不存在優(yōu)勢(shì)生長(zhǎng)方向，生長(zhǎng)過(guò)程中界面前沿總是持續(xù)的發(fā)生分叉、生長(zhǎng)速度周期性的波動(dòng)，不存在類似枝晶穩(wěn)態(tài)生長(zhǎng)的生長(zhǎng)狀態(tài)。

［1］Langer J S，Hong D C.Solvability Conditions for Dendritic Growth in the Boundary-Layer Model with Capillary Anisotropy［J］.Phys Rev A，1986，34（2）：1462-1471.

［2］Langer J S.Recent Developments in the Theory of Pattern-Formation［J］.Physica A，1986，140（1-2）：44-50.

［3］Mullins WW，Sekerka R F.Morphological Stability of a Particle Growing By Diffusion of Heat Flow ［J］.J Appl Phys，1963，34：323-329.

［4］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.

［5］趙達(dá)文.過(guò)冷熔體凝固的相場(chǎng)法自適應(yīng)有限元模擬［D］.西安：西北工業(yè)大學(xué)，2006.

［6］趙達(dá)文，李金富.相場(chǎng)法模擬動(dòng)力學(xué)各向異性對(duì)過(guò)冷熔體中晶體生長(zhǎng)的影響［J］.金屬學(xué)報(bào)，2009，45（10）：1237-1241.

［7］趙達(dá)文.自適應(yīng)有限元方法凝固問(wèn)題模擬的網(wǎng)格各向異性研究［J］.鑄造設(shè)備與工藝，2011（2）：20-22.

［8］Johnson B K，SekerkaR F，F(xiàn)oley M P.Scaling of Fractal Aggregates［J］.Phys Rev E.1995，52（1）：796-800.

［9］Brener E，Müller-Krumbhaar H，Saito Y，et al.Crystal-Growth in a Channel-Numerical Study of the One-Sided Model［J］.Phys Rev E，1993，47（2）：1151-1155.