海水結(jié)冰過程中冰晶生長的相場模擬

韓端鋒， 王永魁， 鞠磊， 王慶

(哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001)

在極地復(fù)雜環(huán)境下，進行科學(xué)考察的船舶和極地保障平臺會面臨結(jié)冰問題。過冷水滴凍結(jié)累積形成的冰層將對外部設(shè)備造成危害。現(xiàn)有的船舶結(jié)冰研究多集中于水滴凍結(jié)或船舶上層建筑結(jié)冰量計算[1-5]，多忽略了冰的微觀生長過程，而冰晶生長行為的研究有助于深入揭示結(jié)冰微觀機理及過程。六角冰晶也稱枝晶，是凝固過程中最常見的一種晶體生長形式，六角冰晶粒相互凍結(jié)在一起，最終形成表面光滑的冰殼[6]。在枝晶生長問題中，一大難題是對固液相變界面的處理，界面的移動往往很難追蹤[7]，而相場法通過引入一個寬度非常小的過渡層來代替跳躍間斷或者尖銳界面，避免了復(fù)雜的界面追蹤，不僅可以通過耦合相場與外場方程將微觀與宏觀結(jié)合起來，還能模擬分析晶體生長過程中的物性參數(shù)對冰晶生長的影響。相場法最先應(yīng)用于金屬模擬領(lǐng)域，而后經(jīng)過發(fā)展，已廣泛應(yīng)用于材料科學(xué)、固體物理、電化學(xué)等領(lǐng)域。Kobayashi[8]通過數(shù)值求解了考慮各向異性的二維相場模型，并得到了過冷純金屬凝固的二維復(fù)雜枝晶生長形貌。Wheeler等[9]提出了模擬二元合金等溫凝固的相場模型(WBM模型)。Wheeler等[10]模擬計算了過冷鎳的凝固過程，結(jié)果與Ivantsov理論相一致。1999年，Kim等[11]建立了可用于合金模擬的相場模型(KKS模型)，消除了WBM模型中界面厚度限制。2005年，Qin等[12]建立了多元合金凝固模型，與實驗結(jié)果基本相符。

近年來，國內(nèi)學(xué)者開始將相場法運用到水的凝固。陳梅英等在食品存儲領(lǐng)域的冰晶生長有較多的研究[13-16]；李方方等[17]實現(xiàn)了細胞尺度下的冰晶數(shù)值模擬；鄒陽[18]模擬了制冷劑中冰晶生長過程，結(jié)果與實驗觀察相一致。但這些都沒有與真實物理場建立較好的聯(lián)系。徐立等[19]討論了無量綱過冷度對冰晶生長的影響，但沒有對Wheeler模型的準(zhǔn)確性及各主要參數(shù)意義進行討論分析。

本文將基于Wheeler模型，采用海水物理參數(shù)來實現(xiàn)海水冰晶的數(shù)值模擬，為過冷水滴凍結(jié)的再輝階段(即枝晶生長過程)提供理論基礎(chǔ)，進一步揭示結(jié)冰機理，并為模擬船舶結(jié)冰過程及評估極地船舶安全提供技術(shù)參考。

1 Wheeler 相場模型

1.1 自由能函數(shù)形式及控制方程

Wheeler相場模型能夠?qū)⑾鄨鰠?shù)與實際物理參數(shù)聯(lián)系起來，使模擬結(jié)果更加真實[20]，這也是本文選取該模型的原因。該模型以金茲堡-朗道理論為基礎(chǔ)，引入相場變量φ，φ=0表示固相，φ=1表示液相，其采用的自由能函數(shù)形式為：

(1)

式中：W為常數(shù)，表示單位體積的能量；β(T)是關(guān)于溫度的函數(shù)，表示熱力學(xué)驅(qū)動力，β(T)<1/2；p(φ)為固相分?jǐn)?shù)，1-p(φ)為液相分?jǐn)?shù)，p(φ)=φ3(10-15φ+6φ2)。f(φ,T)為雙穩(wěn)勢函數(shù)。

相場及溫度場的最終控制方程為：

(2)

(3)

(4)

表1 模型參數(shù)取值Table 1 Model parameter value

1.2 溶質(zhì)場控制方程

海水冰晶可仿照二元合金，視為純水與鹽的混合物，因此其與純冰晶相比，需要考慮濃度場影響，在模型中引入濃度主要有2種方式：1)在自由能密度函數(shù)中包含濃度參數(shù)；2)考慮濃度對過冷度的影響[21]。本文是在純冰晶Wheeler相場模型的基礎(chǔ)上，考慮濃度對過冷度的影響，相場及溫度場方程保持不變，基于Fick定律引入溶質(zhì)場方程為：

(5)

(6)

式中：C為海水中鹽的濃度(摩爾分?jǐn)?shù))；Ds為溶質(zhì)固相擴散率；Dl為溶質(zhì)液相擴散率；k0為溶質(zhì)平衡分配分?jǐn)?shù)；溶質(zhì)場方程中濃度C采用海水中鹽的摩爾分?jǐn)?shù)，已知海水中鹽的質(zhì)量分?jǐn)?shù)約為3.5%，可求得鹽的摩爾分?jǐn)?shù)為0.011。

采用顯式差分格式即可求解每個時間步內(nèi)控制區(qū)域的濃度場，并通過式(7)影響區(qū)域內(nèi)過冷度分布：

(7)

式中：Δ′為濃度影響下的實時無量綱過冷度，其值在每個時間步迭代完成后更新，進而影響相場控制方程；mL為液相線斜率，K/(mol%)；ΔT為實際過冷溫度。

1.3 噪聲強度的引入

在實際的海水凝固結(jié)冰的過程中，受到外界的干擾是不可避免的，因此在相場中引入噪聲項來模擬固液界面的干擾，其干擾會一直存在，造成界面處發(fā)生失穩(wěn)現(xiàn)象，從而導(dǎo)致二次分枝、三次以及四次分枝的出現(xiàn)[22]。本文在相場控制方程中加入噪聲項，記為An(Ar,Rn)：

(8)

式中：Ar為噪聲強度，描述引入擾動的大小；Rn為(-1,1)的隨機數(shù)，符合高斯分布。

2 數(shù)值求解

2.1 控制方程離散形式及穩(wěn)定性條件

有限差分法在均勻網(wǎng)格中進行離散求解，具有原理簡單、方程便于建立和易于編程等優(yōu)點，本文將采用有限差分法求解控制方程，其中時間離散采用一階向前差分，空間離散采用中心差分，通過九點差分格式離散拉普拉斯算子。

顯式差分格式可滿足相場及溶質(zhì)場控制方程計算的穩(wěn)定性，而溫度場方程運用ADI交替顯隱式格式求解，其在二維模擬中無條件穩(wěn)定，因此僅考慮相場及溶質(zhì)場穩(wěn)定性條件，本文控制方程中增加的非線性項對時間步長要求更嚴(yán)格[23]：

(9)

式中 Δx為網(wǎng)格空間尺寸。

2.2 初始條件及邊界條件

結(jié)晶溫度總是低于融化溫度，即液體開始結(jié)晶時需要一個冷卻過程，同時液體中還必須有一定量的微小結(jié)晶核，結(jié)晶核的存在能夠縮短過冷卻過程[24]。晶核設(shè)置如圖1所示。本文設(shè)置液體區(qū)域初始無量綱溫度為-Δ，而無量綱融化溫度為0，從而形成過冷水區(qū)域，同時通過在特定計算區(qū)域設(shè)定固定半徑的圓來引入晶核并設(shè)置初始條件：

x2+y2≤r2:P=0,T=0,C=C0

(10)

x2+y2>r2:P=1,T=-Δ,C=C0

(11)

相場、溫度場及溶質(zhì)場均采用zero-Neumann邊界條件，即: ?P/?n=?T/?n=?C/?n=0。

圖1 晶核設(shè)置示意Fig.1 A schematic diagram of nucleation

2.3 計算優(yōu)化及參數(shù)設(shè)置

本文利用Python語言下的Scipy科學(xué)計算庫求解三對角矩陣，簡化了計算代碼。鑒于經(jīng)典的枝晶生長理論都是關(guān)于枝晶尖端的局部模型，進行定量分析時主要關(guān)注枝晶尖端特征數(shù)據(jù)，因此本文在計算區(qū)域底部中心設(shè)置形核并設(shè)定枝晶相場值沿生長主軸左右對稱，節(jié)省了計算時間。本文設(shè)定網(wǎng)格區(qū)域Nx×Ny為300×300，網(wǎng)格間距Δx及Δy均為0.03，時間步長dt為0.000 1 s，晶核半徑0.05，以上長度值已經(jīng)無量綱化，噪聲強度設(shè)為3。

3 數(shù)值模擬結(jié)果及對比

3.1 冰晶生長過程

在沒有外部擾動的冷水結(jié)冰時，冰晶粒與增長較快的水平軸平行，所以冰晶在水平方向上生長更快，由極小的冰球體發(fā)展到薄圓冰盤，冰盤會繼續(xù)發(fā)展成六角星形[6]，從圓盤到六角星形意味著冰快速生長的開始，六角冰晶彼此凍結(jié)在一起后會形成冰殼，厚度逐漸增加后便形成了一定體積的冰，完成整個結(jié)冰過程。冰晶生長過程見圖2。

圖2 冰晶生長過程示意圖(二維)Fig.2 The growth process of ice crystals (2-D)

過冷水滴凍結(jié)過程通常分為再輝和主凍結(jié)2個階段，再輝階段主要表現(xiàn)為枝晶生長，原則上對水滴凍結(jié)再輝過程模擬即為枝晶生長模擬，由圖3可以觀察初始時水滴與壁面接觸點形成晶核，晶核迅速生長成枝晶形狀而布滿整個水滴，形成糊狀[25]。對枝晶生長的準(zhǔn)確模擬有助于理解常被忽略的水滴凍結(jié)再輝過程。

圖3 水滴凍結(jié)再輝過程實驗觀察結(jié)果[26]Fig.3 Experimental Results of the recalescence process of water droplet freezing crystals[26]

陶樂仁利用低溫顯微鏡觀察了冰晶由晶核生長為類似雪花形狀的過程[26]，實驗結(jié)果見圖4。本文在給定上述基本參數(shù)基礎(chǔ)上，分別設(shè)置時間迭代步數(shù)nt為5、50、200和2 000。由圖5可見，隨著時間迭代步數(shù)的增加，初始晶核開始長大，首先由圓形生長為六邊形，接著六邊形的6個角開始突出生長，進而演化為6個分支，至此形成了六角冰晶的雛形。模擬結(jié)果與實驗觀察結(jié)果相一致。

圖4 低溫顯微鏡下觀察的冰晶生長過程[26]Fig.4 The growth process of crystal observed under microscope[26]

圖5 模擬冰晶生長過程示意Fig.5 A schematic diagram for simulating the growth process of ice crystals

3.2 模擬冰晶與自然界雪晶對比

加州理工學(xué)院Kenneth G. Libbrecht教授在《Field Guide to Snowflakes》一書中展現(xiàn)了其團隊拍攝到的形式多樣的雪花照片[27]，呈現(xiàn)了大自然界中雪花(見圖6(d))。本文模型可以模擬雪花的形貌，保持上述其他參數(shù)不變，劃分網(wǎng)格區(qū)域為800×800，無量綱時間為0.65時，由圖6可見，本文模擬的六角冰晶形貌與大自然中真實存在的雪花是相似的，兩圖中各主支均生長出側(cè)分支，甚至二次分支，有效說明了本文相場模型的真實性。本文結(jié)果與雪花的對比只是象征意義上的，相場模型是研究過冷熔體結(jié)晶過程，而雪花是水蒸氣凝華而成，因此相場模型并不描述雪花形成規(guī)律，可作研究參考。

圖6 數(shù)值模擬結(jié)果與真實雪花對比示意Fig.6 Comparison between simulation results and snowflakes

3.3 尖端特征數(shù)據(jù)的獲取

為了更好地描述枝晶生長行為，凝固學(xué)提出了一些枝晶生長參數(shù)，如尖端生長速度、半徑、溫度等。本文基于每一時間步的模擬結(jié)果，采用數(shù)學(xué)方法獲得上述數(shù)據(jù)，以便定量分析枝晶生長規(guī)律。

定義φ(x,y,t)=0.5為固液界面，通過線性插值獲得t時刻的尖端位置(xtip,ytip)，跟蹤每個時間步的尖端位置，即可獲得尖端速度和尖端溫度，而尖端半徑的獲取則相對復(fù)雜，本文采用拉格朗日方法進行多項式插值，求得尖端半徑[20]：

(12)

圖7 枝晶尖端數(shù)據(jù)隨無量綱時間的變化Fig.7 The variation of dendrite tip data with time

4 參數(shù)敏感性分析

Wheeler模型中無量綱過冷度、各向異性強度、噪聲強度等參數(shù)會影響冰晶生長(包括尖端速度、尖端半徑等)和形貌(側(cè)枝大小、取向等)，本文將對模型主要參數(shù)進行敏感性分析，與真實冰晶形貌進行對比，進而提高冰晶生長模擬的準(zhǔn)確性。

4.1 無量綱過冷度Δ

過冷度對冰晶的形核、生長過程及最終形貌均有重要影響，本例保持其他參數(shù)不變，分別取Δ為0.50、0.60、0.70，對應(yīng)模擬及分析結(jié)果見圖8～9，由圖8可見，隨著Δ值增大，枝晶尖端速度變大，尖端溫度減小，這是因為在低Δ下，界面熱擴散層厚度比較大，會阻礙潛熱釋放，從而抑制了界面處噪聲擾動，導(dǎo)致生長緩慢；而在高過冷度下，界面處溫度梯度比較大，熱擴散層的厚度小，如此便有利于潛熱的釋放，從而加快生長。

圖9的模擬結(jié)果進一步驗證了上述結(jié)論：過冷度較小時，生成的枝晶形貌小，枝晶主干較粗，隨著過冷度增大，枝晶形貌逐步變大，主干變得更加細長，開始出現(xiàn)二次分枝，當(dāng)過冷度足夠大時甚至出現(xiàn)3次、4次分枝，這符合凝固理論中晶體連續(xù)生長機理。

圖8 無量綱尖端速度及溫度隨無量綱過冷度變化情況Fig.8 Changes of dendrite tip velocity and temperature under different dimensionless undercooling conditions

4.2 各向異性強度γ

各向異性強度γ表示界面表面張力、界面厚度及界面動力學(xué)各向異性的程度，γ對枝晶生長過程及結(jié)果有重要影響。本例保持其他參數(shù)不變，分別取γ為0.01、0.02、0.03、0.055、0.065、0.07，對應(yīng)模擬及分析結(jié)果見圖10、11，由圖10可見，隨著γ的增大，尖端速度變大，而尖端半徑及溫度不斷減小，這與微觀可解性(MSC)理論關(guān)于枝晶尖端穩(wěn)態(tài)行為結(jié)論一致[28]。由圖11可見，γ較小時枝晶生長各向異性不明顯，各分枝都比較粗大，主支生長優(yōu)勢不明顯，趨于各向同性。隨著γ增大，晶體沿主支方向生長優(yōu)勢開始凸顯，主支逐漸變得細長，γ=0.05開始出現(xiàn)側(cè)向分枝，這是由于γ增大時會使噪聲擾動更容易被放大，生長前沿變得更不穩(wěn)定，有益于枝晶分裂。但同時，隨著γ更大，大于0.06時有兩大主支被掩蓋，僅四主支生長凸顯，說明γ值不是越大越好，而是需取折中合理值。

圖10 不同各向異性強度對應(yīng)尖端特征數(shù)據(jù)變化Fig.10 Variation of corresponding characteristic data of different anisotropic strength

圖11 不同各向異性強度γ模擬結(jié)果示意Fig.11 Simulation results of different anisotropic strength γ

4.3 噪聲強度Ar

結(jié)冰過程擾動是枝晶側(cè)向分支產(chǎn)生的主要原因，本例保持其他參數(shù)不變，分別取Ar為0.0、0.5、1.0、3.0對應(yīng)結(jié)果見圖12～13，圖12給出了不同Ar下尖端速度及溫度變化情況，可見隨著Ar的增大，尖端數(shù)據(jù)出現(xiàn)波動現(xiàn)象，但仍在某一穩(wěn)定值附近波動。由圖13可見，隨著Ar增大，枝晶界面由光滑變?yōu)槌霈F(xiàn)側(cè)枝，因此Ar能促進生長競爭，進而形成具有發(fā)達分枝的枝晶形態(tài)[20]。

圖12 不同噪聲強度Ar對應(yīng)尖端特征數(shù)據(jù)變化Fig.12 Variation of corresponding characteristic data of different noise strength

圖13 不同噪聲強度Ar模擬結(jié)果示意Fig.13 Simulation results of noise strength Ar

5 結(jié)論

1)敘述了結(jié)冰原理，解釋了晶核—圓盤—六角冰晶—冰殼—厚冰過程，并通過引入晶核方式模擬了海水冰晶的生長，與顯微鏡實驗觀察結(jié)果及自然界真實存在的雪花進行了對比，驗證了本文相場模型的真實性；

2)借助Python語言中高效可靠的第三方科學(xué)計算庫對控制方程求解過程進行了優(yōu)化，有效提高了計算效率和準(zhǔn)確性，并從形核生長位置、對稱區(qū)域等方面節(jié)省了大量計算時間；

3)對無量綱過冷度Δ、各向異性強度γ、噪聲強度Ar等參數(shù)對枝晶尖端生長的影響進行了定量分析。結(jié)果表明，過冷度是海水結(jié)冰的驅(qū)動力，過冷度直接影響界面熱擴散層厚度及潛熱釋放，Δ越大，冰晶界面由固相向液相生長越快，且能促進側(cè)枝生長，形成分枝發(fā)達的枝晶形貌；各向異性強度對枝晶生長行為有重要影響，γ越大，枝晶生長界面處擾動更容易被放大，導(dǎo)致界面前沿不穩(wěn)定，促進了側(cè)枝生長，枝晶形貌變得更復(fù)雜，但γ增大到一定程度后，冰晶會不再對稱，形貌開始失真，因此各向異性強度γ應(yīng)選取適中值；噪聲強度Ar的大小表征生長界面處擾動強弱，Ar越大，枝晶界面生長越激烈，越容易出現(xiàn)側(cè)枝，但Ar過大時也會導(dǎo)致界面失穩(wěn)失真，不宜取過大值；

4)本文模擬結(jié)果與傳統(tǒng)枝晶生長結(jié)論一致，模擬過程中對不用過冷度下枝晶生長速度進行了預(yù)測，并且考慮了界面動力學(xué)效應(yīng)等因素，能為過冷水滴凍結(jié)再輝階段模擬提供依據(jù)，進而完善水滴尺度結(jié)冰過程研究。

本文模擬了海水六角冰晶的生長，揭示了冰晶形成機理。后續(xù)研究將進一步建立相場模擬枝晶生長介觀尺度與水滴凍結(jié)細觀尺度的聯(lián)系，建立結(jié)冰規(guī)律，開展多尺度的數(shù)理建模與預(yù)測手段研究，進而為極地船舶防/除冰設(shè)計及安全提供理論基礎(chǔ)。