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

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

(3 東南大學(xué)航天低溫推進(jìn)劑國家重點(diǎn)實(shí)驗(yàn)室合作基地 南京 210096)

(3 航天低溫推進(jìn)劑國家重點(diǎn)實(shí)驗(yàn)室 北京 100028)

1 引 言

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

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

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

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

2 計(jì)算模型

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

2.1 CA模型

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

對(duì)仍未完全凝固的網(wǎng)格,即界面網(wǎng)格,以過冷度作為枝晶生長的動(dò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模型下的時(shí)間步長和空間步長。bred為形狀減弱因子,早期的CA 方法中,枝晶生長的各向異性依賴于網(wǎng)格,枝晶的生長方向還是沿著網(wǎng)格線的方向。早期這種人為的各向異性重現(xiàn)了枝晶生長存在的各向異性,但是當(dāng)枝晶的擇優(yōu)生長方向與網(wǎng)格不重合時(shí),CA 方法存在的網(wǎng)格各向異性問題就凸顯出來。本研究引入了削弱因子(reduction factor)以減弱四邊形網(wǎng)格的影響[10]。

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

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

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

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

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

2.2 LBM模型

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

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

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

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

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

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

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

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

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

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

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

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

2.3 CA 與LBM 耦合

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

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

2.4 邊界條件

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

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

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

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

3.1 單個(gè)固空枝晶在自然對(duì)流下的生長

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

圖3 擇優(yōu)生長角θ0=90°,t=0.1 s 時(shí)枝晶形貌及氧濃度分布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 時(shí)枝晶形貌及氧濃度分布Fig.4 Dendritic morphologies and oxygen concentration at distribution t=0.1 s with prefered growth orientation θ0=0°

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

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

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

3.2 自然對(duì)流強(qiáng)度對(duì)固空晶體的影響

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

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

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

圖7 不同對(duì)流強(qiáng)度下固相率隨時(shí)間變化趨勢Fig.7 Solid-phase fraction changes with time under different convection intensities

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

圖8 不同對(duì)流強(qiáng)度下濃度峰值隨時(shí)間變化Fig.8 Peak concentration changes with time under different convection intensities

3.3 多晶粒的影響

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

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

4 結(jié) 論

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

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

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