過冷度影響海水結(jié)冰形狀與速度的相場模擬

白 旭， 楊蘇杰

(江蘇科技大學(xué)，船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江，212003)

近年來，隨著全球氣候變暖，北極地區(qū)冰川消融，北極航道得到開發(fā).但是，由于北極地區(qū)特殊的地理環(huán)境，船舶在北極航行時極易受到結(jié)冰的影響[1].已有的針對船舶結(jié)冰的研究大多集中在宏觀層面，即不同氣候條件下，海洋結(jié)構(gòu)物各部位產(chǎn)生的結(jié)冰以及結(jié)冰造成的影響，結(jié)冰的原理也主要來自于經(jīng)驗公式和觀測結(jié)果總結(jié)，很少有對其產(chǎn)生原因和微觀機理的分析[2-4].通過研究海水結(jié)冰的微觀機理，用實際環(huán)境參數(shù)來模擬海水結(jié)冰過程，可以得到更加真實的冰晶生長情況，為過冷條件下海水飛沫結(jié)冰提供微觀理論基礎(chǔ).針對海水的微觀結(jié)冰過程，現(xiàn)在國內(nèi)相關(guān)研究主要集中在單個晶核的凝固現(xiàn)象，而在實際情況下，外界環(huán)境所導(dǎo)致的結(jié)冰幾乎都是多晶核條件下的結(jié)冰.因此，本研究不僅會分析過冷度對海水結(jié)冰的影響，同時會涉及多晶核條件下的海水凝固情況.研究水滴凝固的微觀機理的重點在于追蹤冰晶凝固時固液界面的移動，而通過引入相場法，建立溫度、相場等多個微分方程，可以將真實環(huán)境條件引入模擬計算，模擬固液擴散界面，從而得到冰晶的生成演化過程[5-6].

相場法最早是在20世紀(jì)80年代由Collins等提出[7].1993年，Kobayashi首先建立了關(guān)于相場和溫度的相場模型，在考慮各向異性的情況下，模擬了金屬凝固時的枝晶生長[8].Wheeler等[9]模擬了Ni的凝固過程，并提出了二元合金的相場方程，即WBM模型，并將相場模擬與實際物理量聯(lián)系起來，定量分析參數(shù)對凝固現(xiàn)象的影響.Fan 等[10]基于平均濃度場和熱、質(zhì)及動量輸運現(xiàn)象，模擬了蛋白質(zhì)溶液在冷凍過程中蛋白質(zhì)的分離與重分布.Yang等[11]使用冰模板，考慮石墨烯再分配和冰晶的各向異性，模擬了在水凍結(jié)時石墨烯氣凝膠的形成過程.

國內(nèi)學(xué)者對于相場模型的研究也有一定進展.楊燕等[12]利用Wheeler模型模擬了過冷條件下不同過冷度對冰晶生長的影響.朱昌盛等[13]基于KKS模型模擬了在等溫和不等溫兩種條件下潛熱對二元合金凝固的影響.于艷梅等[14]對單相和二元合金分別進行相場模擬，分析空間尺度、外界噪聲及各向異性等參數(shù)對過冷熔體凝固過程的影響.韓端鋒等[15-16]分別使用Kobayashi模型和Wheeler模型，模擬了單晶核條件下純水和海水的凝固過程，定量分析了外界參數(shù)對冰晶凝固的影響.

針對船舶結(jié)構(gòu)的宏觀結(jié)冰與微觀機理的研究方面，Ryerson等[17]利用船舶切割下的冰晶進行研究，將結(jié)冰的宏觀現(xiàn)象和微觀機理相聯(lián)系，得到了結(jié)構(gòu)結(jié)冰的物理性質(zhì)和結(jié)冰部位的微觀分析.Singh等[18]利用分子動力學(xué)定量研究了光滑粗糙表面上過冷水滴的行為，分析了在粗糙表面上凝固溫度和液滴黏附功與粗糙度之間的聯(lián)系.Sultana等[19]利用數(shù)值模擬分析了單個鹽水滴在撞擊極冷表面后的相變情況，探討了液相率、鹽度等物理性質(zhì)對液滴凍結(jié)的影響.

本文基于Wheeler相場模型，運用有限差分法，對不同過冷度下海水結(jié)冰情況進行模擬，研究過冷度對海水凝固時冰晶生長形狀和速度的影響，為船舶結(jié)冰預(yù)報提供理論和技術(shù)支持.

1 Wheeler相場模型理論

1.1 相場方程

相場法以金茲堡-道朗方程為研究基礎(chǔ)[20]，通過引入相場參數(shù)φ來表示固液界面的移動.當(dāng)φ=1時表示固相，φ=0時表示液相，0<φ<1時表示固液界面，通過數(shù)值的變化來實現(xiàn)對固液界面的跟蹤.

金茲堡-道朗自由能方程為

(1)

式中：f(φ,T)和β(T)分別為自由能密度函數(shù)和熱力學(xué)驅(qū)動力函數(shù)，β(T)<1/2；W為單位體積能量；p(φ)為固相分數(shù)，1-p(φ)為液相分數(shù)，p(φ)=φ3(10-15φ+6φ2).

相場方程為

(2)

其中：δ為特征長度；ω為長度參考尺寸；μ為界面遷移率；σ為表面能；TM為海水凝固點；K為熱擴散系數(shù)；L為單位體積潛熱；cp為比熱容；ΔT為溶液過冷度，ΔT=TM-T0，T0為溶液初始溫度.

溫度場方程為

(3)

式中：p(φ)的導(dǎo)數(shù)p′(φ)=30φ2(1-φ)2；無量綱溫度u=(T-TM)/(TM-T0).

Wheeler相場將實際物理量與相場參數(shù)相結(jié)合，因此在計算中，應(yīng)該使用實際物理參數(shù)進行計算.但是為了將界面厚度控制在可計算范圍內(nèi)，取ω為2.7×10-4cm，將其無量綱化，實際界面厚度等于無量綱界面厚度與長度參考尺寸的乘積，網(wǎng)格間距也按照同樣方法進行換算.

1.2 溶質(zhì)場方程

在模擬海水凝固時，將海水視為水和鹽的二元溶液，相比于一元溶液，二元溶液需要考慮溶質(zhì)的影響.溶質(zhì)對凝固的影響主要有兩種方式，一種是引入溶質(zhì)對過冷度的影響，另一種是在方程中添加濃度項[21].本文通過引入過冷度變化，來實現(xiàn)溶質(zhì)擴散對溶液凝固的影響.溶質(zhì)場方程為

(4)

(5)

式中：C為海水中鹽的摩爾分數(shù)，海水中鹽的質(zhì)量分數(shù)約為3.5%，則鹽的摩爾分數(shù)為0.011；k0為平衡分配系數(shù)；D為無量綱擴散系數(shù)；Ds為固相擴散率；Dl為液相擴散率.

模擬過程中，計算域內(nèi)的過冷度受到溶質(zhì)濃度影響，實時過冷度為

(6)

式中：mL為液相線斜率；C0為溶液初始質(zhì)量分數(shù).

1.3 參數(shù)取值

本文所取物性參數(shù)均按照質(zhì)量分數(shù)為3.5%的海水溶液選取[15]，具體取值如表1所示.

表1 海水溶液物性參數(shù)Tab.1 Physical properties of seawater solution

2 數(shù)值求解與可靠性驗證

在得到相場模型后，需要對其進行數(shù)值求解.有限差分法原理簡單，易于編程，因此本文采用有限差分法來求解控制方程.有限差分法是將要求的函數(shù)在時間和空間上離散，用差分代替微分，并用離散的值來表示函數(shù)在計算域內(nèi)的變化，結(jié)合初始條件和邊界條件，求出函數(shù)在計算域內(nèi)各網(wǎng)格節(jié)點的值.

2.1 有限差分法及離散方式

求解相場控制方程時，空間離散采用中心差分方法，時間離散采用顯式差分方法.另外，對于方程中出現(xiàn)的拉普拉斯算子，采用九點差分方法進行離散：

φi+1, j+1+φi-1, j-1+φi-1, j+1+φi+1, j-1]/

[3(Δx)2]

式中：Δx、Δy為網(wǎng)格間距；Δt為時間步長.對于溫度場和溶質(zhì)場方程，采用與相場方程同樣的差分格式.

在計算域內(nèi)，為了滿足計算精度，保證計算結(jié)果收斂，需要同時考慮相場、溫度場及溶質(zhì)場.計算域內(nèi)的時間步長和空間步長需滿足以下條件

(7)

(8)

2.2 初始條件和邊界條件

本文計算域設(shè)置為正方形，網(wǎng)格數(shù)設(shè)置為1 500×1 500，網(wǎng)格間距Δx和Δy在無量綱化后均為0.03，時間步長Δt=0.000 1 s.在計算域內(nèi)，按照正方形設(shè)置4個初始晶核，晶核半徑為R.初始晶核距離計算域邊界均為50R，晶核之間距離為50R.晶核分布如圖1所示.

圖1 初始晶核大小與位置Fig.1 Initial nuclei size and position

晶核外相場參數(shù)φ設(shè)為1，無量綱溫度u為-Δ，質(zhì)量分數(shù)C為C0，具體晶核位置和條件設(shè)置如下式：

(x-500)2+(y-500)2≤R2時，

φ=0,u=-Δ,C=C0

(9)

(x-1 000)2+(y-500)2≤R2時，

φ=0,u=-Δ,C=C0

(10)

(x-500)2+(y-1000)2≤R2時，

φ=0,u=-Δ,C=C0

(11)

(x-1 000)2+(y-1 000)2≤R2時，

φ=0,u=-Δ,C=C0

(12)

相場、溫度場和濃度場均采用Neumann邊界條件，即

(13)

2.3 模擬冰晶與自然冰晶的對比

為了得到海水結(jié)冰的過程及規(guī)律，首先對冰晶凝固過程進行模擬，并將其和實際自然條件下生成的冰晶進行對比.使用上文設(shè)置參數(shù)，網(wǎng)格劃分為500×500, 模擬時間為7 500Δt，得到結(jié)果如圖2所示.冰晶在凝固時，先由晶核發(fā)展為冰盤，隨后主要沿著6個方向發(fā)展，形成主枝，主枝之間的夾角約為60°，主枝隨后向兩側(cè)長出側(cè)枝和多次分枝，最終形成六角冰晶.通過比較模擬的冰晶形狀和現(xiàn)實的冰晶形狀[22]，發(fā)現(xiàn)實驗結(jié)果吻合良好，可以進行進一步模擬計算.

圖2 模擬冰晶與實際冰晶形狀對比Fig.2 Comparison of simulated ice crystal and actual ice crystal

3 過冷度對晶核生長形狀及速度的影響

3.1 過冷度對晶核生長形狀的影響

根據(jù)凝固理論，過冷度是液體凝固的驅(qū)動力，在冰晶的形核和生長中起重要作用.本文通過分析冰晶形貌、生長速度以及枝晶間距來描述過冷度對冰晶生長的影響.

選擇0.80、0.85、0.90、0.95及1.00這5個不同無量綱過冷度作為模擬溫度，各無量綱過冷度對應(yīng)實際過冷度分別為 -63.34、-67.30、-71.26、-75.22 及 -79.18 K.總模擬時間為 7 500Δt，模擬結(jié)果如圖3所示.

從圖中可以看出，當(dāng)Δ=0.80時，枝晶生長緩慢，僅有少量二次分枝，枝晶直到模擬結(jié)束未發(fā)生交匯，仍然是晶核各自發(fā)展.隨著無量綱過冷度的增大，枝晶生長速度加快，主枝變粗，分枝數(shù)量增加，但枝晶之間仍未發(fā)生交匯.當(dāng)Δ=0.90時，生長速度進一步增加，出現(xiàn)二級以上多次分枝，枝晶在模擬過程中出現(xiàn)交匯.交匯后主枝停止生長，分枝未受影響，繼續(xù)向外生長.Δ=1.00時，粗大的主枝占據(jù)了極大的空間，使得分枝數(shù)量減少.當(dāng)主枝交匯后，主枝逐漸變粗，分枝繼續(xù)發(fā)展，逐漸填滿晶核間空隙，幾乎形成冰盤.

圖4為各過冷度下分別模擬 10 000Δt得到結(jié)果.通過比較各過冷度凝固結(jié)果可以看出，在同時存在多個晶核時，過冷度越大，枝晶就越早接觸，枝晶內(nèi)部空隙就越小，以至于最終完全接觸，內(nèi)部不再存在空隙.過冷度較小時，橫向與縱向主枝有粗細差別，橫向主枝較粗，而過冷度較大時，橫向與縱向主枝粗細基本相同.

圖3 冰晶凝固過程Fig.3 Solidification process of ice crystals

圖4 10 000Δt各無量綱過冷度凝固結(jié)果Fig.4 Solidification results of each dimensionless undercooling at 10 000Δt

3.2 過冷度對晶核生長速度的影響

為了定量分析過冷度對海水凝固的影響，使用枝晶尖端生長速度vtip來描述冰晶生長情況，枝晶尖端生長速度與無量綱時間和無量綱過冷度之間關(guān)系如圖5所示.

圖5(a)為無量綱過冷度為0.8時，尖端速度隨無量綱時間變化情況.根據(jù)計算結(jié)果，枝晶在生長過程中較為穩(wěn)定，尖端速度在開始時較大，隨后慢慢減小，近似呈正弦函數(shù)波動.枝晶尖端速度隨無量綱過冷度的變化如圖5(b)所示，在無量綱過冷度從0.8增至1.0的過程中，枝晶尖端平均速度分別為0.288、0.360、0.432、0.504、0.576 cm/s，枝晶生長速度隨著無量綱過冷度的增加而明顯增加.無量綱過冷度對冰晶凝結(jié)的主要影響在于，無量綱過冷度越大，界面厚度就越薄，凝固時散發(fā)的潛熱就越多，從而使冰晶凝固速度更快.因此，代入實際物理量后可知，過冷度越低，冰晶凝固越快，且凝固速度與過冷度基本為線性關(guān)系.

圖5 枝晶尖端速度變化情況Fig.5 Changes of dendrite tip velocity

圖6 枝晶尖端距離及枝晶距離斜率變化Fig.6 Changes of dendrite tip distance and dendrite distance slope

3.3 過冷度對多晶核生長形狀的影響

本文在計算域內(nèi)設(shè)置4個晶核，在晶核初始間距相同的情況下，通過計算不同過冷度下各晶核主枝之間的距離變化來分析過冷度對海水結(jié)冰的影響情況.

不同過冷度下枝晶尖端距離Dtip隨無量綱時間的變化如圖6(a)所示.隨著無量綱時間的增加，枝晶尖端間距逐漸減小，且二者呈線性關(guān)系.當(dāng)無量綱過冷度達到0.90時，枝晶間距最終為0，代表不同晶核生長出的枝晶最終接觸.根據(jù)圖中趨勢，無量綱過冷度越大，接觸所需時間越短.圖6(b)為將枝晶尖端距離隨時間變化的曲線擬合為一次函數(shù)后的曲線斜率Ktip，其代表枝晶間距的變化快慢，即枝晶生長速度的大小.根據(jù)曲線所示，斜率與無量綱過冷度之間近似呈線性關(guān)系，說明了隨著無量綱過冷度的增大，枝晶生長速度的變化也會相應(yīng)地近似呈線性增大.

4 結(jié)論

基于Wheeler相場模型，使用有限差分法，將海水看作鹽-水二元溶液，把實際物理量轉(zhuǎn)換為方程中的無量綱數(shù)，提取出枝晶生長速度作為定量分析海水結(jié)冰的指標(biāo)，并通過不同過冷度下各晶核枝晶之間的間距變化來描述過冷度對海水結(jié)冰的影響.具體結(jié)論如下：

(1) 無量綱過冷度為0.90以下時，枝晶較細，沒有或僅有二級分枝，枝晶生長較慢.無量綱過冷度為0.90和0.95時，枝晶較粗，枝晶生長速度較快，出現(xiàn)二級以上分枝.無量綱過冷度為1.00時，主枝與分枝明顯變粗，生長速度明顯加快，分枝之間幾乎無空隙.

(2) 無量綱過冷度為0.80時，枝晶生長慢，形狀細，晶核內(nèi)部空隙大，最終形成各主枝組合成的枝狀冰晶.無量綱過冷度為0.90和0.95時，有多次分枝，分枝之間有一定空隙，且空隙隨著分枝發(fā)展而減少.無量綱過冷度1.00時，枝晶生長速度快，分枝之間空隙小，晶核內(nèi)部空隙被填滿，形成冰盤.

(3) 在各過冷度下，枝晶生長速度較為穩(wěn)定，尖端最大速度隨無量綱時間呈正弦變化.無量綱過冷度從0.80增至1.00時，枝晶尖端生長速度從0.12 cm/s增至0.24 cm/s，生長速度與過冷度之間為線性關(guān)系.

(4) 枝晶尖端間距隨著時間的增加而減小，且二者呈線性關(guān)系，間距曲線斜率與無量綱過冷度之間近似為線性關(guān)系，其反映了枝晶生長速度的與無量綱過冷度之間也近似為線性關(guān)系，與之前求得的枝晶尖端速度相吻合.