基于LBM-CA 的六角冰晶非平衡生長(zhǎng)模擬研究

藍(lán)碧蘭 郭 晴 全曉軍

(上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院 上海 200240)

1 引言

結(jié)冰現(xiàn)象在自然界和工業(yè)生產(chǎn)中廣泛存在,包括飛機(jī)機(jī)翼結(jié)冰、冷凍食品工藝及冰模板法制備多孔材料等。冰晶的生長(zhǎng)分為非平衡階段及平衡階段。六角枝晶也稱枝晶,是非平衡凝固過(guò)程中最常見(jiàn)的一種晶體生長(zhǎng)形式。枝晶的生長(zhǎng)提供平衡階段的初始生長(zhǎng)條件,深入探究枝晶生長(zhǎng)的機(jī)理有利于加強(qiáng)對(duì)冰晶形貌的調(diào)控。Zhang 等[1]基于冰模板法制備了具有可調(diào)孔徑和結(jié)構(gòu)的納米纖維素氣凝膠,通過(guò)改變冷凍溫度、溫度梯度證明了冷凍溫度和溫度梯度會(huì)對(duì)冰晶的形貌和生長(zhǎng)速率產(chǎn)生影響。Fukasawa 等[2]發(fā)現(xiàn)宏觀溫度梯度場(chǎng)會(huì)改變冰晶的形貌。研究表明冷凍工藝條件對(duì)冰晶的尺寸、形貌具有較大影響,但實(shí)驗(yàn)操作過(guò)程中往往由于結(jié)晶的隨機(jī)性以及生長(zhǎng)時(shí)間較短等因素,導(dǎo)致直接對(duì)冰晶進(jìn)行定量觀測(cè)和微觀機(jī)理研究存在局限性。

為了進(jìn)一步探究冰晶生長(zhǎng)的微觀機(jī)理,打破實(shí)驗(yàn)探究的局限性,一系列數(shù)值模擬方法被廣泛用于結(jié)晶過(guò)程的研究,包括蒙特卡洛法、頂點(diǎn)法、相場(chǎng)法,元胞自動(dòng)機(jī)法等。蒙特卡洛法簡(jiǎn)潔方便、易于編程,但存在模擬時(shí)間與實(shí)際時(shí)間無(wú)法對(duì)應(yīng)的問(wèn)題。頂點(diǎn)法高效性只體現(xiàn)于頂點(diǎn)的計(jì)算,使用受限。相場(chǎng)法(Phase Field,簡(jiǎn)稱PF)和元胞自動(dòng)機(jī)法是目前在晶體生長(zhǎng)過(guò)程中應(yīng)用較為廣泛的兩種方法。韓端鋒等[3]采用經(jīng)典Kobayashi 相場(chǎng)模型討論了界面寬度均值、融化溫度等多個(gè)參數(shù)對(duì)冰晶生長(zhǎng)及形貌的影響,結(jié)果表明冰晶形貌對(duì)參數(shù)的依賴性較高。白旭等[4]采用Wheeler 相場(chǎng)模型,分析不同過(guò)冷度對(duì)海水冰晶生長(zhǎng)的影響,結(jié)果表明冰晶生長(zhǎng)速率隨過(guò)冷度增加而增加。然而該模擬受計(jì)算精度及收斂性要求的限制,網(wǎng)格數(shù)量需求較多,對(duì)算力要求相對(duì)較高。相場(chǎng)法普遍存在界面厚度較難確定、結(jié)果對(duì)參數(shù)較為敏感的局限性。并且相場(chǎng)法計(jì)算量較大,而當(dāng)前計(jì)算機(jī)數(shù)據(jù)處理能力有限,因此相場(chǎng)法僅適用于微小區(qū)域的模擬。

與相場(chǎng)法相比,元胞自動(dòng)機(jī)(cellular automaton,CA)方法計(jì)算速度較快,能夠適應(yīng)較大范圍內(nèi)的凝固模擬。LBM(Lattice Boltzmann Method,LBM)是一種區(qū)別于傳統(tǒng)傳輸模型的模擬手段。一方面LBM 避免了求解非線性方程組,簡(jiǎn)化了計(jì)算過(guò)程;另一方面,LBM 方法具有先天并行性、物理意義清晰的優(yōu)點(diǎn)。耦合CA 方法與LBM 方法的模型具備并行性高、計(jì)算速度快、計(jì)算范圍廣的優(yōu)點(diǎn)。周靖超等[5]基于LBM-CA 方法模擬強(qiáng)制對(duì)流作用下單相合金的非對(duì)稱生長(zhǎng)行為,模擬結(jié)果表明對(duì)流速度越大、初始過(guò)冷度越小、冷卻速率越小會(huì)使枝晶非對(duì)稱生長(zhǎng)行為趨向更強(qiáng)。張勇等[6]基于LBM 方法研究了熱障陶瓷涂層凝固生長(zhǎng)過(guò)程,等離子噴涂涂層在正溫度梯度方向呈柱狀晶生長(zhǎng),溫度梯度越大,柱狀晶越發(fā)達(dá)。上述單相合金固溶體及等離子噴涂涂層晶體均呈四角晶形貌,而過(guò)冷水結(jié)冰多為六角冰晶。與四角晶相比,六角冰晶存在網(wǎng)格離散導(dǎo)致的各向異性的問(wèn)題,是模擬中需要重點(diǎn)解決的問(wèn)題。Guo 等[7]基于LBM-CA 方法研究自然對(duì)流及強(qiáng)制對(duì)流對(duì)冰晶的生長(zhǎng)速率及形貌的影響。模型引入減弱因子削弱了各向異性,但該文章并未研究溫度梯度及冷卻速率等冷凍工藝對(duì)冰晶生長(zhǎng)的影響。

目前采用LBM-CA 方法模擬冰晶生長(zhǎng)的研究仍處于起步階段。鮮少有學(xué)者研究溫度梯度與冷卻速率對(duì)枝晶生長(zhǎng)的影響。因此,本文基于GPU 加速,建立LBM-CA 模型,探究過(guò)冷度、溫度梯度及冷卻速率等因素對(duì)枝晶形貌及生長(zhǎng)速率的影響,并引入減弱因子和優(yōu)化網(wǎng)格差分方法以削弱網(wǎng)格各向異性及離散各向異性對(duì)結(jié)果的影響[7]。

2 數(shù)值模型

2.1 CA 模型

過(guò)冷度ΔT是結(jié)晶過(guò)程的動(dòng)力,主要包括熱過(guò)冷(ΔTt)、曲率過(guò)冷(ΔTr)及成分過(guò)冷(ΔTc),本文是純水結(jié)晶過(guò)程的模擬,因此ΔTc=0。

式中:Γ0為平均吉布斯-湯姆遜系數(shù),εΓ為由各向異性界面張力引起的各向異性系數(shù),K(x)為固液界面處的界面曲率,即:

不同的網(wǎng)格格點(diǎn)可以根據(jù)固相率分為固相元胞(S)、界面元胞(S/L)、液相元胞(L)。

對(duì)于尚未完全凝固的界面元胞,固相率將繼續(xù)增加直至為1,成為固相網(wǎng)格,并且按照元胞自動(dòng)機(jī)Moore 類型進(jìn)行界面元胞的捕獲。

界面處的法向角可以用式(6)計(jì)算得到:

根據(jù)尖銳界面模型,枝晶尖端生長(zhǎng)速率是過(guò)冷度(ΔT)的函數(shù):

式中:μk(θ)為各向異性動(dòng)力學(xué)系數(shù),可以通過(guò)式(8)計(jì)算:

式中:εμ為關(guān)于生長(zhǎng)速率波動(dòng)的函數(shù),θp為優(yōu)先生長(zhǎng)方向。

2.2 格子玻爾茲曼模型(LBM)

溫度場(chǎng)碰撞遷移格式為:

式中:τΓ為溫度場(chǎng)的松弛函數(shù),τΓ=3α+0.5,α為熱擴(kuò)散系數(shù);φ(x)為溫度場(chǎng)的源項(xiàng),主要受凝固過(guò)程放熱的影響:

溫度場(chǎng)的平衡分布函數(shù)表達(dá)式為:

2.3 各向異性縮減處理

2.3.1 網(wǎng)格各向異性

為了降低網(wǎng)格各向異性,引入一個(gè)假想存在的各向同性標(biāo)量場(chǎng)[8],在該標(biāo)量場(chǎng)中標(biāo)量的定義如下:

擴(kuò)散方程為:

式中:D=0.24。

引入縮減因子bred(x):

經(jīng)過(guò)比較,使用虛擬各向同性場(chǎng)?(x)能夠在界面處得到一個(gè)更為準(zhǔn)確且穩(wěn)定的法向角,因此法向角的計(jì)算公式調(diào)整為式(15):

2.3.2 離散各向異性

離散各向異性對(duì)界面曲率的大小有較大的影響,通過(guò)調(diào)整固相元胞周圍8 個(gè)相鄰網(wǎng)格的權(quán)重以及添加縮減因子可以使得界面曲率的取值進(jìn)一步貼近實(shí)際情況。其中:

同時(shí)在界面曲率K(x)的分母項(xiàng)中添加各向異性削減因子,表示如下:

2.4 CA 模型與LBM 模型的耦合

CA 模型主要進(jìn)行凝固過(guò)程冰晶生長(zhǎng)過(guò)程的模擬,而LBM 模型主要進(jìn)行溫度場(chǎng)的模擬。通過(guò)CA模型計(jì)算的尖端生長(zhǎng)速度可以計(jì)算出每一次迭代的固相率增量:

為了保證LBM 模型與CA 模型的相關(guān)性,LBM模型中無(wú)量綱擴(kuò)散系數(shù)及動(dòng)力粘度的選取如式(19):

式中:LBM 中δx及δt均為1,CA 模型中Δt及Δx均為7.02 ×10-7。

上述各式中部分常數(shù)取值見(jiàn)表1。

表1 水的特性參數(shù)Table 1 Physical parameters of water

3 過(guò)冷度對(duì)枝晶生長(zhǎng)的影響

3.1 過(guò)冷度對(duì)枝晶生長(zhǎng)形貌的影響

過(guò)冷度是凝固的主要驅(qū)動(dòng)力,在枝晶生長(zhǎng)中起重要作用。不同初始過(guò)冷度下,枝晶形貌顯著不同。如圖2 所示,枝晶的尖銳程度和枝狀分化程度隨過(guò)冷度增加而增加,主要原因是初始過(guò)冷度增大會(huì)導(dǎo)致生長(zhǎng)驅(qū)動(dòng)力增大、熱擴(kuò)散層減薄,從而提高枝晶生長(zhǎng)速率。如圖3 所示,低過(guò)冷度下枝晶尖端被更厚的熱擴(kuò)散層包圍,熱量擴(kuò)散受阻礙。高過(guò)冷度下表現(xiàn)為界面層減薄,熱量擴(kuò)散速率增加,枝晶分化更明顯。

圖1 Moore 型格點(diǎn)捕獲規(guī)則示意圖Fig.1 Schematic of the Moore type capture rules

圖2 t=1.20 s 不同過(guò)冷度條件下枝晶生長(zhǎng)形貌圖Fig.2 Dendrite growth morphology under different subcooling degree in same time

圖3 不同過(guò)冷度枝晶溫度云圖對(duì)比圖Fig.3 Contrast diagram of dendrite temperature cloud of different undercooling degrees

如圖4a 所示,過(guò)冷度越大,枝晶臂越發(fā)達(dá),尖端越深入過(guò)冷水。圖4b 顯示枝晶臂間出現(xiàn)不同程度的“縮頸”效應(yīng)。過(guò)冷度越大,“縮頸”效應(yīng)越明顯。主要原因是過(guò)冷度越大,尖端生長(zhǎng)越快,從而潛熱釋放加快。但枝晶臂間空間狹小,熱量擴(kuò)散困難。因此過(guò)冷度越大,臂間熱積聚效應(yīng)越明顯,尖端生長(zhǎng)受更大程度地抑制。

圖4 t=1.26 s 不同過(guò)冷度枝晶形貌對(duì)比圖Fig.4 Contrast diagram of dendrite morphology at different undercooling degrees in same time

3.2 過(guò)冷度對(duì)枝晶生長(zhǎng)速度的影響

Ivantsov[9]基于熱量擴(kuò)散的假設(shè)給出了過(guò)冷度與枝晶生長(zhǎng)Peclet 數(shù)(Pe)的關(guān)系:

式中:E1(Pe)表達(dá)式為:

Pe表示對(duì)流速率與擴(kuò)散速率之比,表達(dá)式為:

Pe既是尖端速率也是尖端半徑的函數(shù),僅由式(21)無(wú)法求解。

引入Langer 和Muller-Krumbhaar[10-11]提出的臨界穩(wěn)定性理論。LMK 臨界穩(wěn)定性原理指出枝晶尖端穩(wěn)定性取決于一個(gè)無(wú)量綱數(shù)σ:

定義σ?為由實(shí)驗(yàn)及數(shù)值計(jì)算確定的穩(wěn)定性因子,σ?通常取0.02 或0.025。當(dāng)σ>σ?時(shí),將會(huì)產(chǎn)生側(cè)向分枝不穩(wěn)定性使σ減小;當(dāng)σ>σ?時(shí),又將導(dǎo)致尖端不穩(wěn)定性使σ增大。因此,σ趨向于穩(wěn)定在σ?附近。因此,根據(jù)式(21)—式(24)即可求得生長(zhǎng)速率隨過(guò)冷度的變化曲線。

如圖5 所示,初始階段枝晶生長(zhǎng)速率較大,而后生長(zhǎng)速率逐漸降低直至趨于穩(wěn)定。主要原因是枝晶生長(zhǎng)初始階段,枝晶尖端仍保持較高的過(guò)冷度,生長(zhǎng)速率較大。而后隨著界面元胞不斷釋放潛熱,界面前沿溫度不斷上升直至潛熱釋放與熱擴(kuò)散達(dá)到動(dòng)態(tài)平衡時(shí),溫度基本不變,生長(zhǎng)速率趨于穩(wěn)定。此外,圖5對(duì)比了過(guò)冷度為0.65 K 時(shí)LMK 理論及本文模擬結(jié)果的生長(zhǎng)速率變化。LMK 理論中當(dāng)σ?給定,枝晶尖端始終保持相同速度生長(zhǎng)。而本文模擬結(jié)果除在初始階段受初始過(guò)冷度影響而高于LMK 曲線外,整體與LMK 曲線較為一致。

圖5 不同過(guò)冷度條件下尖端生長(zhǎng)速率隨時(shí)間變化圖Fig.5 Diagram of tip growth rate changing with time under different degree of undercooling

如圖6 所示,本模擬結(jié)果尖端生長(zhǎng)速率與初始過(guò)冷度呈近似線性關(guān)系,主要原因是當(dāng)Pe較小時(shí),對(duì)流換熱量與熱傳導(dǎo)相比處于較低水平,相變釋放的潛熱主要依靠熱傳導(dǎo)向周圍傳遞。過(guò)冷水與枝晶取相同的導(dǎo)熱系數(shù),熱傳導(dǎo)各向同性。因此,界面元胞過(guò)冷度主要受初始過(guò)冷度影響。由式(7)可知,當(dāng)尖端μk(θ)變化較小時(shí),生長(zhǎng)速率與初始過(guò)冷度呈近似直線關(guān)系。

圖6 t=0.84 s 時(shí)尖端生長(zhǎng)速率隨過(guò)冷度變化關(guān)系圖Fig.6 Growth rate of tip varies with degree of undercooling at t=0.84 s

令LMK 理論中σ?分別為0.02 及0.025 計(jì)算出不同過(guò)冷度下枝晶尖端的生長(zhǎng)速率并與本模擬結(jié)果進(jìn)行對(duì)比。如圖6 所示,二者在小過(guò)冷度范圍內(nèi)生長(zhǎng)速率均隨過(guò)冷度增加而增加,趨勢(shì)一致。當(dāng)過(guò)冷度過(guò)大或過(guò)小時(shí),受網(wǎng)格尺寸及算力的限制理論分析結(jié)果與模擬結(jié)果出現(xiàn)一定的偏差?？傮w上,本模擬結(jié)果與LMK 理論分析結(jié)果較為吻合。

4 溫度梯度

4.1 溫度梯度對(duì)枝晶形貌的影響

在計(jì)算域范圍內(nèi)施加沿X軸負(fù)方向的溫度梯度。不同的溫度梯度下枝晶呈現(xiàn)顯著不同的形貌。如圖7 所示,負(fù)溫度梯度方向的枝晶臂均受到不同程度的促進(jìn),其中與溫度梯度反向平行方向的枝晶臂受到最大程度的促進(jìn)。這是因?yàn)榕c溫度梯度反向平行的方向是溫度下降最快的方向,尖端過(guò)冷度最大,促進(jìn)作用最顯著。其余兩個(gè)枝晶臂因?yàn)榕cX軸正方向存在一定夾角,所以促進(jìn)效應(yīng)減弱。然而,正溫度梯度方向枝晶臂的生長(zhǎng)受到了不同程度的抑制,其中與溫度梯度平行方向枝晶的生長(zhǎng)受到最大程度的抑制。主要原因是與溫度梯度平行的方向是溫度升高最快的方向,過(guò)冷度最小,抑制作用最顯著。但由于周圍流體處于過(guò)冷狀態(tài),所以枝晶仍具備生長(zhǎng)動(dòng)力。其余兩個(gè)枝晶臂因?yàn)榕cX軸負(fù)方向存在一定夾角,抑制作用減弱。因此,在宏觀溫度梯度的作用下,枝晶形貌得到控制,枝晶實(shí)現(xiàn)了定向凝固。

圖7 t=1.24 s 不同溫度梯作用下枝晶形貌圖Fig.7 Dendrite morphology under different temperature gradient at t=0.84 s

從圖8 可看到,溫度梯度越大,正溫度梯度方向受抑制越嚴(yán)重,負(fù)溫度梯度方向受促進(jìn)越顯著。主要原因是溫度梯度越大,正溫度梯度方向溫度上升越快,生長(zhǎng)速率越小,抑制作用越嚴(yán)重;而負(fù)溫度梯度方向溫度下降越快,生長(zhǎng)速率越大,促進(jìn)作用越顯著。因此,合理選擇溫度梯度是對(duì)冰晶形貌進(jìn)行調(diào)控的有效方法。

圖8 t=0.84 s 不同溫度梯度作用下枝晶形貌對(duì)比圖Fig.8 Contrast diagram of dendrite morphology under different temperature gradient at t=0.84 s

4.2 溫度梯度對(duì)枝晶尖端生長(zhǎng)速率的影響

如圖9 所示,負(fù)溫度梯度方向枝晶臂生長(zhǎng)速率不斷上升,而正溫度梯度方向枝晶臂生長(zhǎng)速率不斷下降。該現(xiàn)象與圖5 等溫枝晶后期生長(zhǎng)速率趨于穩(wěn)定的現(xiàn)象不同。主要原因是負(fù)溫度梯度方向枝晶尖端可以不斷深入溫度更低的過(guò)冷水中,相變產(chǎn)生的熱量能夠迅速地向周圍過(guò)冷液體擴(kuò)散。相反地,正溫度梯度方向枝晶尖端周圍過(guò)冷水溫度提升,凝固推動(dòng)力減弱,尖端生長(zhǎng)速率不斷下降。

圖9 不同溫度梯度下枝晶尖端生長(zhǎng)速率隨時(shí)間變化圖Fig.9 Diagram of tip growth rate changing with time under different temperature gradient

圖10 是t=0.82 s 時(shí)尖端生長(zhǎng)速率隨溫度梯度的變化圖。在負(fù)溫度梯度方向枝晶尖端生長(zhǎng)速率隨溫度梯度的增大而增大;而正溫度梯度方向,枝晶尖端生長(zhǎng)速率隨溫度梯度的增加反而下降。原因是在溫度梯度作用下枝晶生長(zhǎng)偏離對(duì)稱等溫生長(zhǎng),促進(jìn)作用與抑制作用均隨溫度梯度的增大而增大。

圖10 t=0.82 s 時(shí)的尖端生長(zhǎng)速率隨著溫度梯度的變化圖Fig.10 Growth rate of tip varies with t temperature gradient at t=0.84 s

值得注意的是,此時(shí)生長(zhǎng)速率與溫度梯度并不為線性關(guān)系。兩條曲線的斜率均隨溫度梯度的增大而增大。主要原因是溫度梯度的作用具有正反饋特性。溫度場(chǎng)的差異導(dǎo)致生長(zhǎng)速率的差異,而生長(zhǎng)速率的差異又通過(guò)影響界面所處位置而作用于生長(zhǎng)速率本身。例如負(fù)溫度梯度方向,溫度梯度越大,枝晶生長(zhǎng)速率越大,尖端越深入低溫過(guò)冷水,生長(zhǎng)速率將進(jìn)一步提升。

5 冷卻速率

如圖11 所示,恒定冷卻速率作用下枝晶表現(xiàn)出高度對(duì)稱性。因?yàn)槭┘雍愣ɡ鋮s速率后冰晶尖端過(guò)冷度仍保持相同,6 個(gè)枝晶臂的生長(zhǎng)速率一致,枝晶對(duì)稱生長(zhǎng)。值得注意的是,枝晶的尺寸隨冷卻速率的增加而增大。冷卻速率越大,枝晶生長(zhǎng)速率越快,相同時(shí)刻枝晶越大;而冷卻速率越小,枝晶生長(zhǎng)速率越小,相同時(shí)間內(nèi)枝晶將越小。因此可以通過(guò)合理地調(diào)控冷卻速率達(dá)到控制枝晶大小的目的。

圖11 不同冷卻速率下枝晶生長(zhǎng)形貌對(duì)比圖Fig.11 Contrast diagram of dendrite morphology under different cooling rate

6 結(jié)論

基于格子玻爾茲曼方法及元胞自動(dòng)機(jī)方法,利用尖銳界面模型,探究過(guò)冷度、溫度梯度及冷卻速率對(duì)枝晶形貌的影響,結(jié)論如下:

(1)基于GPU 加速,采用LBM-CA 算法,引入網(wǎng)格各向異性削弱因子,模擬了冰晶在過(guò)冷水中的生長(zhǎng),結(jié)果與實(shí)際情況吻合,各向異性削弱因子表現(xiàn)良好。

(2)過(guò)冷度影響枝晶生長(zhǎng)速率從而影響枝晶形貌及大小,枝晶尖端生長(zhǎng)速率隨過(guò)冷度的增大而增大,呈現(xiàn)近似直線的關(guān)系。模擬結(jié)果與LMK 臨界穩(wěn)定性理論吻合,驗(yàn)證了模擬結(jié)果的正確性。

(3)通過(guò)施加不同溫度梯度實(shí)現(xiàn)冰晶的定向凝固。與溫度梯度相反的方向,枝晶生長(zhǎng)受到促進(jìn),其中與溫度反向平行方向受到最大程度的促進(jìn)。與溫度梯度相同的方向,枝晶生長(zhǎng)受到抑制,其中平行方向受到最大程度的抑制。

(4)冷卻速率影響枝晶的尺寸。施加恒定的冷卻速率可以得到對(duì)稱的六角冰晶,但枝晶的大小與冷卻速率呈現(xiàn)正相關(guān)關(guān)系。因此,合適的冷卻速率是制備優(yōu)異性能材料過(guò)程中需要重要考慮的因素。