基于LBM-CA模型的自然對流下液氫中固空枝晶生長模擬

戴聞驍 雷 剛 鄭曉紅 梁文清 疏志勇 錢 華

(1 東南大學能源與環境學院 南京 210096)

(3 東南大學航天低溫推進劑國家重點實驗室合作基地 南京 210096)

(3 航天低溫推進劑國家重點實驗室 北京 100028)

1 引 言

液氫的儲存形式具有綜合密度高、方便使用、技術成熟的優點,尤其在航天用的火箭發動機上具有良好的應用前景。由于在液氫的加注,轉運等過程中不可避免滲入空氣,研究表明,液氫中積累的固空是富氧的[1],當固空中的氧濃度達到一定比例時,液氫系統易發生爆炸[2-4]。因此研究液氫系統中固空凝固過程對提高液氫系統安全性能有重要應用價值。

當液氫流動停止時,凝固釋放的潛熱和溶質導致局部密度不均將引發自然對流。相關研究表明[5],自然對流將對凝固過程中枝晶的生長產生影響。由于液氫的低溫特性,且固空晶體在系統中隨機位置形核生長,現有實驗設備難以直接清晰觀察到固空顆粒的微觀形成過程,更難對其成分進行定量測量分析。劉海生[2-3]通過理論計算和摻雜-復溫的實驗方法測得在液氫系統復溫過程中氮氧的濃度變化,得到固空晶體中心貧氧外表富氧的宏觀定性結論,并觀察到雪花狀固空晶體,但未能定量給出氮氧濃度分布情況。

鑒于難以進行微觀的實驗觀察和定量測量,本研究采用數值模擬的方法研究自然對流作用下固空的微觀結構。傳統的CFD 方法難以模擬凝固過程中晶體具體的微觀組織,為將微觀組織形成過程與流場、溫度場、濃度場相耦合,選用格子玻爾茲曼(LBM)模型耦合元胞自動機(CA)的模擬方法。孫東科[6]耦合LBM 與CA 方法模擬了在復雜流場下枝晶生長的特性,該方法具備了計算速度快、天然并行、能夠處理復雜邊界等優點。但由于四邊形網格的局限性,LBMCA模型多用于模擬鐵、銅、鋁[7-8]等具有面心立方和體心立方晶體結構的四重對稱微觀枝晶生長。周偉煜[9]采用LBM-CA 方法模擬了固空四重對稱枝晶生長的過程,雖然得到氧濃度在固空顆粒中的分布情況,但與劉海生觀察到的雪花狀結構不符,即仍未解決網格局限性問題。在自然對流條件下流場受枝晶結構影響較大,因此LBM-CA 方法用于模擬自然對流作用下固空晶體六重對稱生長時,網格各向異性將影響計算的準確性。本文模擬中引入減弱因子[10]并優化差分方法[11]減小正方形網格各向異性及離散各向異性的影響。

本研究將空氣簡化為N2-21%O2的氮氧混合物,以過冷度作為生長驅動力,采用改進的LBM-CA 方法研究在自然對流條件下,溶質析出與潛熱擴散對枝晶形貌的影響,探究固空顆粒中氮氧濃度的分布情況。

2 計算模型

本研究中使用LBM-CA模型將固空枝晶生長與流場的演化結合到一起,CA 方法模擬枝晶生長過程,LBM模型模擬流場的演化,兩種模型相互耦合模擬在復雜流場下固空的六重對稱枝晶生長過程。

2.1 CA模型

網格點的凝固程度用固相率分數fs表示,根據網格點的固相率分數fs,將網格分為液相網格(L),固相網格(S)及界面網格(S/L)。

對仍未完全凝固的網格,即界面網格,以過冷度作為枝晶生長的動力,界面網格的固相率不斷增加最終成為新的固相網格,并按照Moore 元胞類型,捕獲周圍液相網格使其成為新的界面網格(圖1)。

圖1 Moore 元胞捕獲原則Fig.1 Moore cell capture principle

界面處固相率的增長Δfs與枝晶生長速度vn有關:

式中:Δt與Δx分別為CA模型下的時間步長和空間步長。bred為形狀減弱因子,早期的CA 方法中,枝晶生長的各向異性依賴于網格,枝晶的生長方向還是沿著網格線的方向。早期這種人為的各向異性重現了枝晶生長存在的各向異性,但是當枝晶的擇優生長方向與網格不重合時,CA 方法存在的網格各向異性問題就凸顯出來。本研究引入了削弱因子(reduction factor)以減弱四邊形網格的影響[10]。

枝晶的生長速度vn為局部過冷度的函數:

其中:μk為各向異性系數:

式中:μk0為平均向異性系數,εμ為各向異性因子,θx-θ0表示界面法向量與擇優生長方向的夾角,(°);m表示不同的多重對稱性,在6 重對稱枝晶生長中m取6。

總過冷度ΔT由熱過冷Tt(K),成分過冷Tc(K)和曲率過冷TR(K)組成,即:

式中:Clx為當地液相濃度,%;C0為初始濃度,%;當液相線斜率ml為負時局部濃度越高越不利于凝固的進行。Γ為Gibbs-Thomson 系數,K為界面曲率,其可根據物理意義表示為:

2.2 LBM模型

格子玻爾茲曼方法(LBM)根據離散的Boltzmann方程進行演化,在演化過程中分為遷移和碰撞步,因其天然并然性大大加快了計算速度。

本文采用常用基于單松弛時間的D2Q9模型即LBGK[12]方法,其演化方程為:

式中:fi(x,t)為粒子分布函數為粒子平衡態分布函數,δt為格子時間步長,ei為離散速度,τf為松弛時間Fi為外力源項,可表示為:

式中:cs為格子聲速,等于。c,c為格子速度,c=δx/δt。δx、δt分別為格子步長和格子時間步長,F為自然對流的浮力,根據Boussinesp 近似,該項可表示:

式中:ρ0為流體的初始密度,kg/m3;βT,βC分別為溫度膨脹系數和濃度膨脹系數;g為重力加速度,m/s2;C0為初始濃度,%;T0為初始溫度,K。

離散速度ei,wi為權重系數,當i=0 時,wi=4/9,當i=1—4 時,wi=1/9,當i=5—8 時,wi=1/36。

類似地,溶質場及溫度場的演化也可用對應的格子Boltzmann 方程計算:

式中:ΔCl為凝固過程中由于潛熱釋放和溶質析出對周圍濃度,%;ΔTH為凝固過程中由于潛熱釋放和溶質析出對周圍濃度溫度,K。

式中:k為溶質分配系數;Lh為單位質量溶液所含凝固潛熱,kJ/kg;cp為等壓熱容,kJ/(kg·K)。

速度場,溫度場,濃度場的平衡態分布函數分別為:

流場、濃度場、溫度場的松弛時間為運動粘度ν、溶質擴散系數DL和熱擴散系數α的函數:

根據LBGK模型將粒子分布函數恢復到密度ρ,(kg/m3);流速u,(m/s);濃度C,%;溫度T,K;等宏觀物理量。

2.3 CA 與LBM 耦合

初始狀態下,在計算區域內放置一個晶核作為枝晶生長起點,根據上述捕獲規則,固相網格捕獲周圍液相網格,使其成為固液界面網格。通過LBM模型得到流場的速度、溫度、濃度分布,決定該時間步長內固液界面的濃度和溫度。根據界面的溫度及濃度,利用CA 法計算該時間步長內的枝晶生長。凝固釋放出的潛熱和溶質作為源項添加到LBM 演化方程中,在自然對流過程中,因濃度和溫度的不均衡,在重力作用下引發的自然對流也需考慮。模擬程序流程圖如圖2 所示。

圖2 LBM-CA模擬流程圖Fig.2 Flow chart for simulating dendritic growth

2.4 邊界條件

在計算區域充滿了一定過冷度和一定濃度的氮氧混合液,熱邊界均采用非平衡外推,凝固潛熱最終釋放到外界壞境中保證凝固過程一直進行,濃度場采用無擴散邊界條件。計算區域四周及因枝晶生長在計算區域內部形成的復雜固液界面均設為無滑移邊界,并在LBM 計算中采用標準反彈格式處理。由于相對于溶質在液相中的擴散系數,固相中擴散系數小3 個量級,所以忽略溶質在固相中的擴散,因而溶質無法以擴散的形式穿過固液界面。

3 模擬結果與討論

將計算區域劃分為300 ×300 的網格,網格尺寸Δx=0.5 μm。LBM模型中的松弛因子取τf=1,因此格子粘度為1/6。t=0 時刻,計算區域中心放置一半徑為5 格子的圓形晶核,初始過冷度為0.6 K,冷卻速率為5 K/s。模擬中所需的參數如表1 所示[9,13]。

表1 模擬時所用物性參數Table 1 Parameters used in present work

3.1 單個固空枝晶在自然對流下的生長

本部分研究了單個固空枝晶在自然對流下的生長情況。取溫度瑞利數Ra=5 000,濃度瑞利數Ra=5 000,擇優生長角度θ0與重力方向夾角為0°及90°時,自然對流與純擴散條件下枝晶生長情況對比如圖3、圖4 所示(枝晶形貌由黑色輪廓線表示)。在純擴散條件下,凝固析出的溶質及釋放的潛熱僅依靠擴散向四周傳遞,各枝晶界面濃度、溫度幾乎一致,呈對稱分布。由于界面枝晶生長各向異性,晶核沿著特定的方向生長,固-液界面由圓形發展為正六邊形。與純擴散條件相比,在t=0.1 s 時,自然對流條件下枝晶生長表現出的不對稱性已相當明顯,當擇優生長方向不同時,固空枝晶表現出的生長不對稱性也不盡相同,但均呈現出下側枝晶的生長受到促進,上側枝晶的生長普遍受到抑制。

圖3 擇優生長角θ0=90°,t=0.1 s 時枝晶形貌及氧濃度分布Fig.3 Dendritic morphology and oxygen concentration at t=0.1 s with prefered growth orientation θ0=90°

圖4 擇優生長角θ0=0°,t=0.1 s 時枝晶形貌及氧濃度分布Fig.4 Dendritic morphologies and oxygen concentration at distribution t=0.1 s with prefered growth orientation θ0=0°

由于枝晶生長過程中在界面處析出溶質和釋放潛熱,導致其附近流體密度與較遠處流體密度有較大區別,在重力的作用下引起自然對流,左右兩側形成對稱分布的兩個順時針方向的旋渦。局部形成了小的旋渦,在小漩渦處出現了貧氧區。本研究主要考慮濃度場的變化引起的成分過冷度不同,從圖3、圖4所示的自然對流條件下濃度分布可以看出,自然對流促使溶質再分配出現局部差異。由式(9)可知成分過冷與界面濃度有關,各生長方向界面濃度不一致,因而導致枝晶生長速度上的差異。在自然對流的作用下,下側枝晶處于迎流位置,析出溶質及凝固潛熱易被流場帶到枝晶上側,從而使流場上游的枝晶尖端界面處濃度更低,過冷度更大,促進了該部分枝晶的生長。相反地,流場下游即枝晶上側尖端,溶質富集,擴散層較厚,溶質不易向四周擴散,生長受到抑制。該結果與對流條件下的枝晶生長結論[14]一致。

進一步分析自然對流對濃度場的影響,圖3、圖4可以看出自然對流更利于溶質的擴散,使溶質分布更均勻。圖5 為擇優生長方向為0°,t=0.1 s 時,模擬區域中線(y=75 μm)處氧濃度的變化曲線。從圖5可以看出,固空顆粒呈現出外表富氧,中心貧氧的特性,與劉海生等[3]宏觀試驗結果一致。與劉海生等人的試驗相比,本文的模擬給出了微觀的固空枝晶結構及更為精確的濃度分布結果。在t=0.1 s 時,純擴散條件下的氧濃度峰值為23.71%,自然對流條件下的濃度峰值為23.14%且高氧濃度區域更少,LTTCHFIELD 等[15]通過槍擊試驗所得的氧濃度安全閾值為24%,可見自然對流將降低液氫系統中固空帶來的安全風險。

圖5 t=0.1 s 時y=75 μm 處不同條件下氧濃度分布Fig.5 Oxygen concentration distribution under different conditions with y=75 μm and at t=0.1 s

3.2 自然對流強度對固空晶體的影響

衡量對流強度的兩個無量綱參數是濃度瑞利數和溫度瑞利數,由于在枝晶生長過程中濃度瑞利數的影響較小,因此本文僅分析溫度瑞利數對固空枝晶的影響。圖6 給出了不同對流強度下枝晶生長情況。當溫度瑞利數Ra從3 000 增加至7 000,自然對流對下側枝晶生長的促進越來越明顯。t=0.1 s 時,下側枝晶的長度分別為:36.5 μm,41 μm,48 μm。同時,上側枝晶臂隨著流場強度的增加逐漸消失。

圖6 t=0.1 s 時不同對流強度下枝晶形貌及氧濃度分布Fig.6 Dendritic morphology and oxygen concentration distribution under different convection intensities at t=0.1 s

雖然上側枝晶臂在流場的作用下生長受到抑制,從圖7 不同瑞利數下固相率隨時間變化曲線可以看出,自然對流對枝晶生長的固相率凈增長仍起到促進作用。隨著自然對流強度的增加,固相率增長更快。在初始階段,由于流場較弱,枝晶生長與純擴散的條件差異較小,不同條件下固相率增速無較大區別。隨著溶質的不斷析出與潛熱釋放,區域的濃度差及溫度差差異增大,自然對流強度增強,流場速度更快,自然對流對枝晶生長的影響越來越明顯。自然對流使模擬區域的溶質分布更均勻,局部氧高濃度區減小,枝晶界面處的濃度擴散層更薄,促進了枝晶的生長。

圖7 不同對流強度下固相率隨時間變化趨勢Fig.7 Solid-phase fraction changes with time under different convection intensities

從圖8 不同條件下氧濃度峰值隨時間變化也可以看出,純擴散條件下的氧濃度峰值在0.02 s 后均高于存在自然對流,進一步說明自然對流能促進溶質的擴散,促進固相率的增長。注意到Ra=3 000 時,氧濃度峰值在t=0.01 s 時出現極大值達到22.58%,由于在凝固初期自然對流強度較小,溶質未充分向遠處低濃度區擴散,導致溶質富集在枝晶上側,隨著自然對流強度增強,該峰值濃度下降。Ra=5 000 時凝固后期溶質峰值高于Ra=3 000,這是因為自然對流提高了凝固速率,更多的溶質析出,然而其自然對流強度不及Ra=7 000 時,不能及時向周圍傳遞,導致濃度逐漸超過Ra=3 000 時的峰值濃度。

圖8 不同對流強度下濃度峰值隨時間變化Fig.8 Peak concentration changes with time under different convection intensities

3.3 多晶粒的影響

在中心晶粒四周放置4 個晶核,研究多晶粒生長時各晶粒之間的相互影響,模擬結果如圖9 所示。凝固初始階段,晶核均對稱生長,如圖9a。隨著凝固的進行,自然對流強度增強,在枝晶附近形成了大小不等的旋渦,溶質擴散更充分,氧濃度在整個計算區域內分布更均勻,各枝晶的生長逐漸發展。圖9b、圖9c可以看出,隨著枝晶生長,不斷有溶質析出,中心處枝晶析出的溶質由于被周圍晶粒阻擋,難以向周圍擴散,同時又由于周圍枝晶排出的溶質也向中心聚集,在中心區域形成了溶質富集區,抑制了中心晶粒的生長,向柱狀晶轉變。外部的枝晶生長受自然對流的影響較為顯著,下側枝晶的生長速度明顯快于上側枝晶,即自然對流將下游析出溶質帶到上游,促使了下游枝晶的生長。當t=0.15 s 時,枝晶的增長使得流動區域更為狹窄,形成了更多的小漩渦,加速了溶質的擴散,濃度分布更為均勻。

圖9 多個固空顆粒在自然對流作用下枝晶形貌隨時間變化Fig.9 Dendrite morphology of several solid particles changes with time under natural convection

4 結 論

(1)建立改進的LBM-CA模型,減小網格各向異性影響,模擬自然對流條件下液氫中固空六重對稱枝晶生長過程。結果進一步驗證了固空顆粒外表富氧中心貧氧的特點,自然對流下固空枝晶生長不再呈現出對稱性,上游枝晶生長得到促進,下游受溶質聚集的影響枝晶生長受到抑制。

(2)自然對流有利于析出溶質在區域內分布更均勻,氧濃度峰值較純擴散條件有所下降,這有利于降低液氫系統的安全風險。自然對流條件下枝晶界面濃度擴散層更薄,促進了凝固過程,且自然對流強度越強,對固空的生長更有利。

(3)對多晶粒在自然對流條件下的生長模擬結果表明,受周圍枝晶的影響,中心顆粒的枝晶生長受到明顯抑制。