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

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

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

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

1 引 言

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

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

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

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

2 計算模型

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

2.1 CA模型

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

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

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

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

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

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

其中:μk為各向異性系數(shù):

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

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

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

2.2 LBM模型

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

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

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

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

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

離散速度ei,wi為權(quán)重系數(shù),當i=0 時,wi=4/9,當i=1—4 時,wi=1/9,當i=5—8 時,wi=1/36。

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

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

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

速度場,溫度場,濃度場的平衡態(tài)分布函數(shù)分別為:

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

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

2.3 CA 與LBM 耦合

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

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

2.4 邊界條件

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

3 模擬結(jié)果與討論

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

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

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

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

圖3 擇優(yōu)生長角θ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 擇優(yōu)生長角θ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°

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

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

圖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 自然對流強度對固空晶體的影響

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

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

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

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

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

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

3.3 多晶粒的影響

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

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

4 結(jié) 論

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

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

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