液氫中固空枝晶生長凝固定量相場(chǎng)法研究

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

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

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

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

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

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

1.1 定量相場(chǎng)模型

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

T=Tint-RDt

(1)

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

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

(2)

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

(3)

無量綱氧溶質(zhì)場(chǎng)U的擴(kuò)散方程為

(4)

其中

(5)

(6)

(7)

(8)

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

as(n)=1+ε1·

(9)

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

(10)

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

1.2 模型數(shù)值求解

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

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

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

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

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

(11)

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

(12)

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

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

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

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

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

(a)1.0 K (b)1.5 K

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

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

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

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

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

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

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

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

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

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

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

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

(a)2 K/s

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

3 結(jié) 論

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

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

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

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

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