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

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

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

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

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

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

1.1 定量相場模型

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

T=Tint-RDt

(1)

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

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

(2)

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

(3)

無量綱氧溶質場U的擴散方程為

(4)

其中

(5)

(6)

(7)

(8)

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

as(n)=1+ε1·

(9)

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

(10)

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

1.2 模型數值求解

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

表1 模型參數及物性參數[24]

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

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

(11)

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

(12)

圖1展示了不同無量綱過冷度下過飽和度Ω與Pc之間的依賴關系,紅色虛線代表解析理論預測結果,而散點數據為當前QPFM模擬結果?？梢钥吹?當前模擬與解析結果相對誤差范圍為-3.75%～4.33%,表明兩者取得了良好的匹配性。

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

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

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

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

(a)1.0 K (b)1.5 K

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

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

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

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

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

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

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

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

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

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

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

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

(a)2 K/s

圖9 不同冷卻速率下最大氧溶質質量分數隨時間的變化

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

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

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

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

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