液氫中固空枝晶生長凝固定量相場法研究

李超龍,文鍵,王磊,厲彥忠,,雷剛,屠基元

(1.西安交通大學能源與動力工程學院,710049,西安; 2.航天低溫推進劑技術(shù)國家重點實驗室,100028,北京; 3.皇家墨爾本理工大學工學院,VIC 3083,澳大利亞墨爾本;4.清華大學核能與新能源技術(shù)研究院,100084,北京)

面對日益嚴峻的氣候問題,清潔低碳的氫能是一種備受青睞的替代能源[1]。其中,液氫液氧組合的推進劑是大推力火箭發(fā)動機的理想選擇[2]。然而,由于液氫的低溫特性,暴露于液氫的空氣會迅速凍結(jié),包裹于液氫中的固空顆粒在微弱能量的激發(fā)下存在爆炸風險。此外,固空的積聚沉積容易堵塞細小的管道和閥門,使得火箭發(fā)動機燃燒不穩(wěn)定。Cassutt等[3]利用熔絲引爆300 g固空和4.73 L液氫的混合物,發(fā)現(xiàn)固空中的氧含量是引發(fā)爆炸的關(guān)鍵因素;以槍擊的形式引爆,能夠引起液氫爆炸的固空氧含量最低為24%;當與惰性氣體共存時,該值提升至40%。由于溶解度的差異,凍結(jié)的固體空氣中氧溶質(zhì)的分布并不均勻。劉海生等[4]等研究了固空在液氫中的沉積現(xiàn)象,發(fā)現(xiàn)固空呈現(xiàn)雪花狀。此外,通過測量升溫過程氧氮逸出質(zhì)量分數(shù)發(fā)現(xiàn)固空顆粒表面富氧,內(nèi)部貧氧的氧溶質(zhì)分布特征。然而,這些實驗未能得到固空枝晶的全過程生長圖像,也未能直接得到枝晶生長凝固過程中枝晶內(nèi)部和液相中氧溶質(zhì)的分布。

考慮到液氫的低溫及不穩(wěn)定性,通過實驗手段獲得固空枝晶微結(jié)構(gòu)演化和溶質(zhì)分布信息是相當具有挑戰(zhàn)的。幸運的是,數(shù)值模擬方法如元胞自動機(CA)、相場法等已成為復現(xiàn)枝晶凝固微觀結(jié)構(gòu)演化的有效手段[5-10]。課題組在先前的研究[11-12]中運用格子玻爾茲曼和元胞自動機耦合方法初步探索了固空枝晶生長凝固行為,但所得到的固空枝晶形態(tài)較為簡單,僅包含一次枝晶臂,并不能反映真實固空枝晶微觀結(jié)構(gòu)。相場模型不需要明確追蹤移動界面,彎曲枝晶界面對枝晶生長動力的貢獻即Gibbs-Thomson效應(yīng)可以被精確表達。得益于眾多學者的巨大努力,相場模型取得了顯著進展[13-16]。Karma等[15]基于薄界面限制定量模擬了純質(zhì)的等軸枝晶生長,在此基礎(chǔ)上,該模型陸續(xù)被拓展到二元和多元組分系統(tǒng)的凝固行為的定量模擬。Rojas等[17]采用相場結(jié)合格子玻爾茲曼法探究了對流存在下枝晶的生長凝固行為。Nabavizadeh等[18]使用相場法對微重力條件下琥珀腈溶液的枝晶生長的實驗數(shù)據(jù)進行了驗證,得到的枝晶尖端生長速度與實驗符合良好。由于具有六重對稱性的固空枝晶微觀結(jié)構(gòu)演化[19]及氧溶質(zhì)分布對于液氫系統(tǒng)的安全至關(guān)重要,在完善先前研究不足的動機下,我們致力于開發(fā)能夠復現(xiàn)固空枝晶真實形態(tài)的定量相場模型。

為了探究液氫中固空枝晶生長凝固特性,本文建立了二維定量相場模型,且為提高數(shù)值求解過程的穩(wěn)定性,對相場序參量進行了非線性預(yù)處理,研究了連續(xù)冷卻條件下固空單枝晶和多枝晶微觀結(jié)構(gòu)演化和氧溶質(zhì)分布規(guī)律,可為液氫系統(tǒng)的安全使用提供理論指導。

1 數(shù)值模型與求解

1.1 定量相場模型

本文將空氣簡化為僅含有氧氣和氮氣的二元組分,空氣降溫后形成成分均勻的二元溶液,氧溶質(zhì)的質(zhì)量分數(shù)為0.233。采用完善的定量相場模型探究固空枝晶微結(jié)構(gòu)演化過程和溶質(zhì)偏析現(xiàn)象,可以在固液界面厚度W0大于毛細長度d0的條件下保持定量特性。由于空氣溶液的熱擴散系數(shù)高出溶質(zhì)擴散系數(shù)兩個數(shù)量級,計算域中的溫度梯度忽略不計,并且不考慮凝固潛熱,因此具有恒定冷卻速率的空氣枝晶等溫凝固,其計算域內(nèi)溫度場表述為

T=Tint-RDt

(1)

式中:Tint是初始溫度;RD是冷卻速率,當RD=0時計算域保持初始溫度,即空氣溶液具有恒定過冷度;t是時間。

連續(xù)場φ被用于描述與空間和時間對應(yīng)的相,對于固相,φ=1,對于液相,φ=-1,在固液界面區(qū)域,φ在-1和1之間光滑地變化。為了使相場變量φ在較大網(wǎng)格下的時間演化中保持數(shù)值穩(wěn)定,相場φ的本構(gòu)方程通過非線性預(yù)處理的定量相場ψ重新定義

(2)

非線性預(yù)處理定量相場ψ的演化方程為

(3)

無量綱氧溶質(zhì)場U的擴散方程為

(4)

其中

(5)

(6)

(7)

(8)

當Ω=1時,方程(8)退化為方程(7)。as(n)是表征枝晶各向異性的函數(shù),對于具有六重對稱性的固空枝晶的二維計算,as(n)表達式[21]如下

as(n)=1+ε1·

(9)

as(n)=as(θ)=1+ε1cos[6(θ-θ0)]

(10)

式中:θ0是相對于x軸的晶體取向角度,當θ0=0°時,方程(10)和方程(9)等效。

1.2 模型數(shù)值求解

定量相場模型在x和y方向的空間離散均采用具有固定步長Δx的有限差分方法,其在時間域的演化采用歐拉顯式時間方案,時間步長為Δt。空間和時間分別以固液界面厚度W0和松弛時間τ0為單位進行縮放。對于具有六重對稱性的固空枝晶,當規(guī)則笛卡爾正交網(wǎng)格被用于定量相場模型(QPFM)的有限差分離散時,虛假網(wǎng)格各向異性的引入不可避免地影響了枝晶尖端的生長動力[22],對于具有不同生長取向的多枝晶生長凝固過程,這種影響尤其明顯。實際上,在無其他因素干擾時,不同取向的固空枝晶生長形態(tài)應(yīng)具有旋轉(zhuǎn)不變性。要保持枝晶生長的旋轉(zhuǎn)不變性,方程(3)和(4)的離散要求最大化地保持各向同性。為達到這一目的,本研究采用了Ji等[23]提出的一種基于兩個基礎(chǔ)離散基的線性組合來保持笛卡爾網(wǎng)格中各向同性化離散的有限差分法。相場φ通過預(yù)處理相場ψ的演化獲得,當|ψ|≥10.27時,|1-φ|≤10-6,此時認為處于單一均勻相。對于單枝晶計算,相場模型四周邊界采用零擴散Neumann邊界;對于多枝晶計算,則四周采用周期性邊界。計算初始在計算域中放置晶核并賦予枝晶取向角,對于多枝晶,不同枝晶的取向角可以是隨機的。對于具有不同枝晶生長取向角的多枝晶計算,在保持一個唯一相場φ的條件下,引入一個枝晶取向場P用于記錄當?shù)氐闹∠蚪?而當?shù)氐闹∠蚪抢^承與其距離最近的固空枝晶。采用的模型參數(shù)及物性參數(shù)見表1。

表1 模型參數(shù)及物性參數(shù)[24]

2 計算結(jié)果及分析

2.1 自由枝晶生長:與解析理論的對比

由于空氣枝晶的尺度微小和液氫的低溫特性,因此難以獲取固空枝晶生長的直接實驗數(shù)據(jù),考慮與解析理論對比來驗證相場模擬的可靠性。流行的枝晶生長解析理論一般通過形態(tài)穩(wěn)定性分析和溶解性理論得到枝晶穩(wěn)態(tài)生長模式與過冷度之間的關(guān)系。Ivantsov第一個提出了描述控制枝晶生長的過冷度ΔT與枝晶生長Pélect數(shù)Pc之間的關(guān)系的枝晶穩(wěn)態(tài)生長解析理論并定義了Ivantsov函數(shù)。Pc=vtipρtip/2D,vtip是尖端生長速度,ρtip是尖端曲率半徑。尖端速度vtip通過記錄尖端位置獲得,而尖端曲率半徑則通過尖端輪廓文件擬合得到。Alexandrov等[25]將六重對稱枝晶的各向異性考慮在內(nèi)推導了決定vtip和ρtip的選擇標準并給出了修正的Ivantsov函數(shù)

(11)

對于2D計算,方程(11)中j取2,過飽和度Ω與Pc的關(guān)系可定義如下

(12)

圖1展示了不同無量綱過冷度下過飽和度Ω與Pc之間的依賴關(guān)系,紅色虛線代表解析理論預(yù)測結(jié)果,而散點數(shù)據(jù)為當前QPFM模擬結(jié)果。可以看到,當前模擬與解析結(jié)果相對誤差范圍為-3.75%～4.33%,表明兩者取得了良好的匹配性。

圖1 過飽和度Ω與Pc的依賴關(guān)系:當前模型與解析理論的對比

2.2 單枝晶生長:初始過冷度的影響

本節(jié)研究初始過冷度對固空枝晶形態(tài)演化和氧溶質(zhì)分布的影響,并選擇兩條特殊位置曲線定量表述氧溶質(zhì)分布規(guī)律。對于所有案例,計算區(qū)域網(wǎng)格為設(shè)定為500×500,經(jīng)過無關(guān)性檢驗,網(wǎng)格步長選定為Δx=0.4W0。計算域初始氧溶質(zhì)質(zhì)量分數(shù)c0=0.233,半徑r=115d0的圓形晶核被置于計算域中心,晶核內(nèi)氧質(zhì)量分數(shù)為kc0,枝晶生長取向角θ0=0°。冷卻速率保持RD=10 K/s,時間步長Δt=0.001τ0,模擬持續(xù)20 000步。

圖2(a)～(d)左半邊展示了20 000Δt時刻不同初始過冷度(1.0、1.5、2.0、2.5 K)下固空枝晶的形態(tài)。可以看出,不同初始過冷度下固空枝晶均呈現(xiàn)六重對稱雪花形態(tài),但在初始過冷度為1.0 K時僅分化出少量的二次枝晶臂,枝晶二次枝晶臂隨初始過冷度增大而更加發(fā)達,枝晶形態(tài)復雜度增加。

(a)1.0 K (b)1.5 K

圖2中右半邊則展示了20 000Δt時刻不同初始過冷度下氧溶質(zhì)質(zhì)量分數(shù)分布情況,可以明顯看到,氧溶質(zhì)質(zhì)量分數(shù)分布呈現(xiàn)枝晶內(nèi)部貧氧外部富氧的分布特征。對于不同初始過冷度,氧溶質(zhì)質(zhì)量分數(shù)在固液界面附近達到最大并高于初始質(zhì)量分數(shù)(0.233),而固空枝晶內(nèi)部氧溶質(zhì)則低于初始質(zhì)量分數(shù),這與早前的實驗結(jié)果[4]相符合。隨著初始過冷度的增加,固空枝晶二次枝晶臂更為發(fā)達,二次枝晶臂的存在不利于氧溶質(zhì)及時的向外擴散,在枝晶臂的間隙位置氧溶質(zhì)積聚從而導致了更高的氧質(zhì)量分數(shù)。

圖3展示了枝晶尖端生長速度和固相率Fs(固體元胞占計算域元胞比例)隨時間的變化。固相率隨時間增長一直增加并且由于二次枝晶臂的發(fā)展固相率后期的增加速率變快,同一時刻,固相率隨初始過冷度增加而變大。在18 000Δt時刻固相率由初始過冷度1.0 K的0.222增加到2.5 K的0.301。初始時刻,不同初始過冷度下枝晶尖端均具有一個較高的生長速度,隨著凝固的進行枝晶尖端生長速度逐漸降低并達到穩(wěn)定階段。穩(wěn)定階段的枝晶生長速度隨過初始冷度增加而增大,初始過冷度由1.0 K增加到2.5 K時,固空枝晶尖端速度由0.364 mm/s增長到0.673 mm/s。

圖3 枝晶尖端速度和固相率隨時間的變化

由于溶解度的差異,固空枝晶生長過程中會發(fā)生氧溶質(zhì)再分配現(xiàn)象。固體空氣中氧溶質(zhì)的溶解度較低,因此固空枝晶凝固生長時部分氧溶質(zhì)會被釋放到液相中形成溶質(zhì)偏析現(xiàn)象。盡管課題組之前的研究建立了數(shù)值模型來研究固空枝晶的這種溶質(zhì)偏析現(xiàn)象,但固空枝晶的形態(tài)被簡化,模型的定量能力不足[11-12]。接下來,將選取兩條特殊曲線位置著重以定量的形式揭示QPFM得到的固空枝晶生長凝固過程中的氧溶質(zhì)分布規(guī)律。

圖4展示了20 000Δt時刻不同初始過冷度下沿計算域水平中線(L1線)的相場和氧溶質(zhì)質(zhì)量分數(shù)分布。可以看到,相場參數(shù)數(shù)值在中心區(qū)域保持1(固相),而在兩邊區(qū)域保持-1(液相)不變,數(shù)值突變區(qū)域即為固液糊狀區(qū),固液界面即固空枝晶尖端具體位置以φ=0等值線對應(yīng)位置確定,由于生長速度的差異,枝晶尖端到邊界的距離有所不同。在固空枝晶核心區(qū)域,氧溶質(zhì)質(zhì)量分數(shù)由內(nèi)向外逐漸升高,當氧溶質(zhì)質(zhì)量分數(shù)達到0.111時,固空枝晶開始從核心區(qū)域分化出一次枝晶臂。不同初始過冷度下沿與P1線重合的一次枝晶臂其晶軸氧溶質(zhì)質(zhì)量分數(shù)幾乎保持恒定且只表現(xiàn)出微弱的差別。由于氧溶質(zhì)的釋放,氧質(zhì)量分數(shù)在糊狀區(qū)陡然升高,糊狀區(qū)氧質(zhì)量分數(shù)峰值和初始過冷度呈現(xiàn)正相關(guān),值得注意的是,氧質(zhì)量分數(shù)峰值出現(xiàn)的位置并非固液界面處。當初始過冷度為1.0、1.5、2.0、2.5 K時,糊狀區(qū)氧溶質(zhì)質(zhì)量分數(shù)峰值Cm分別為0.402、0.407、0.425、0.469,而固液界面處氧溶質(zhì)質(zhì)量分數(shù)Ci分別為0.265、0.268、0.280、0.313,均高于初始質(zhì)量分數(shù)0.233。

圖4 不同初始過冷度下沿L1線相場和氧溶質(zhì)分布

由于沿一次枝晶臂和二次枝晶臂的軸線形成了氧溶質(zhì)質(zhì)量分數(shù)較低的脊柱(圖2)。為了全面描述氧溶質(zhì)分布規(guī)律,圖5給出了沿紅色虛線圓弧C1(以枝晶中心為圓心,半徑R=190Δx)的氧溶質(zhì)質(zhì)量分數(shù)分布,角度θ范圍為150°～210°。可以看到,不同初始過冷度下自一次晶軸方向(180°)向兩側(cè),氧溶質(zhì)質(zhì)量分數(shù)逐漸升高,并在固液界面附近達到峰值。當初始過冷度為1.0、1.5、2.0、2.5 K時,氧溶質(zhì)質(zhì)量分數(shù)峰值Cm分別為0.433、0.449、0.466、0.488。對于初始過冷度2.5 K,氧溶質(zhì)質(zhì)量分數(shù)分布曲線出現(xiàn)多個波谷,這是因為此時枝晶有著發(fā)達的二次枝晶臂,這些波谷即對應(yīng)二次枝晶臂軸線位置。

2.3 多枝晶生長:冷卻速率的影響

應(yīng)用QPFM對多枝晶在連續(xù)冷卻條件的生長演化進行了研究。初始過冷度為ΔT0=1.0 K,冷卻速率分布設(shè)置為2、5、10、15 K/s,其他條件與上節(jié)保持相同。

圖6展示了不同冷卻速率下固空枝晶形態(tài)隨時間的演化,冷卻速率分別為2、5、10、15 K/s,第一列到第4列對應(yīng)的時刻分別是4 000Δt、15 000Δt、30 000Δt、45 000Δt。初始階段,各枝晶之間距離較遠,相互之間幾乎沒有影響,各枝晶沿各自取向發(fā)展并形成一次枝晶臂。隨著時間的增長,一次枝晶臂逐漸長大并開始分化出二次枝晶臂,不同枝晶的枝晶臂開始接近并影響了枝晶臂的原有生長軌跡,位于中心的枝晶的生長受到明顯抑制。在連續(xù)冷卻條件下,二次枝晶臂隨時間增加持續(xù)續(xù)粗化、融合,對于冷卻速率RD=15 K/s,枝晶最后幾乎充滿了計算域。

圖7展示了不同冷卻速率下固相率隨時間的變化,可以看到固相率隨時間增長而增大,并且更大的冷卻速率對應(yīng)的固相率更高。連續(xù)冷卻條件下,在30 000Δt時刻之后固相率仍持續(xù)增長但增長率降低。在60 000Δt時刻,當冷卻速率分別為2、5、10、15 K/s時,固相率Fs分別為0.446 0、0.658 5、0.830 5、0.911 1。圖8展示了與圖6相對應(yīng)的連續(xù)冷卻條件下多枝晶氧溶質(zhì)質(zhì)量分數(shù)分布。同樣地,枝晶生長初期因各枝晶間距離較遠,氧質(zhì)量分數(shù)分布不受其他枝晶干擾。隨著時間的增長,枝晶臂開始接近并形成間隙,枝晶凝固生長過程排出到液相的氧溶質(zhì)積聚在枝晶臂間隙中,使得局部的飽和度增加從而抑制了當?shù)刂П鄣陌l(fā)展。圖9給出了不同冷卻速率下計算域內(nèi)氧溶質(zhì)質(zhì)量分數(shù)最大值隨時間的變化情況。可以看到,計算域內(nèi)氧溶質(zhì)質(zhì)量分數(shù)最大值隨時間幾乎呈現(xiàn)線性增長趨勢增長率隨冷卻速率增大而增大。在45 000Δt時刻,當冷卻速率分別為2、5、10、15 K/s時,氧溶質(zhì)質(zhì)量分數(shù)最大值Cm分別為0.387 6、0.491 9、0.668 4、0.840 6。

圖7 不同冷卻速率下固相率隨時間的變化

(a)2 K/s

圖9 不同冷卻速率下最大氧溶質(zhì)質(zhì)量分數(shù)隨時間的變化

3 結(jié) 論

本文建立了經(jīng)過非線性預(yù)處理的定量相場模型,著重對單枝晶和多枝晶在連續(xù)冷卻條件的生長凝固行為進行了研究;獲得了固空枝晶形態(tài)的瞬態(tài)演化和定量氧溶質(zhì)分布,增進了液氫中固空枝晶生長行為的理解,對液氫系統(tǒng)安全使用具有一定參考意義,主要結(jié)論如下。

(1)氧溶質(zhì)質(zhì)量分數(shù)總體呈現(xiàn)枝晶內(nèi)部含氧量低外部含氧量高的分布特征,在一次枝晶臂和二次枝晶臂軸線上形成低氧溶質(zhì)質(zhì)量分數(shù)脊柱。

(2)固空枝晶尖端生長速度隨初始過冷度的增加而增大,并且二次枝晶臂更為發(fā)達,形態(tài)變得復雜。當初始過冷度由1 K增加到2.5 K時,固空枝晶尖端速度由0.364 mm/s增長到0.673 mm/s。

(3)多枝晶生長時固空枝晶的對稱性被破環(huán);連續(xù)冷卻條件下,固空枝晶二次枝晶臂隨時間增長逐漸粗化、融合,最終充滿計算域,計算域內(nèi)氧溶質(zhì)質(zhì)量分數(shù)最大值隨凝固進行持續(xù)升高,氧溶質(zhì)質(zhì)量分數(shù)在晶界位置達到最大值。冷卻速率由2 K/s升高到15 K/s,氧溶質(zhì)質(zhì)量分數(shù)最大值由0.387 6升高到0.840 6。

考慮到固空枝晶生長過程可能受外界熱力環(huán)境影響,將包含凝固潛熱在內(nèi)的能量方程耦合到當前模型即進行非等溫相場模擬是未來進一步的研究方向。