固空沉積的數(shù)值模擬

周偉煜，梁文清*，錢華，雷剛，荀其寧

（1-東南大學(xué)能源與環(huán)境學(xué)院，江蘇南京 210096；2-航天低溫推進(jìn)劑技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100028； 3-國防科技工業(yè)應(yīng)用化學(xué)一級計(jì)量站，山東濟(jì)南 250031）

0 引言

隨著航天事業(yè)的發(fā)展，液氫作為一種高比沖特性的環(huán)境友好性推進(jìn)劑組元，其應(yīng)用水平和應(yīng)用規(guī)模有了極大的提高。而在推進(jìn)液氫工業(yè)化過程中，液氫系統(tǒng)安全性問題顯得尤為重要。液氫系統(tǒng)在液氫加注、氣氫排放、系統(tǒng)吹除置換中，都極易造成少量空氣的滲入。氮、氧在液氫中溶解度低，液氫中固空的積累不可避免。而固空中氮氧比例分布是影響固空安全性的重要因素[1]。

國外開展的固空爆炸試驗(yàn)研究表明，當(dāng)固空中氧體積分?jǐn)?shù)大于空氣中氧體積分?jǐn)?shù)時(shí)，會有爆炸甚至爆轟危險(xiǎn)[2]。國內(nèi)劉海生等[3-4]開展固空理論和實(shí)驗(yàn)研究，通過摻雜-復(fù)溫-氣相色譜檢測的方式得到固空顆粒表面富氧、內(nèi)部貧氧的結(jié)論。在理論研究方面，通過對比微觀下氮氧分子間吸引能得到與試驗(yàn)相一致的結(jié)論。這些研究僅能對固空中氧分布進(jìn)行定性描述，而未給定量結(jié)果。由于實(shí)驗(yàn)研究難以對氮氧成分進(jìn)行定量探究，故采用凝固模擬來解決固空沉積氮氧分布問題。

目前已有較多關(guān)于相變過程的模擬及實(shí)驗(yàn)研究。齊雋楠等[5]研究了疏水和親水豎直圓柱面上的滴狀冷凝特性以及膜狀冷凝特性。郭亞等[6]建立三分螺旋折流板冷凝器模型，并模擬對比了單、雙頭變角度三分螺旋折流板與變間距弓形折流板立式冷凝器的凝結(jié)換熱性能。胡海濤等[7]建立金屬平片冷表面濕空氣凝水過程熱質(zhì)傳遞特性的數(shù)值模型。黃飛等[8]提出了外科手術(shù)冷凍過程的動(dòng)態(tài)模型。而凝固領(lǐng)域較常采用主要有相場（PF）模型[9-10]、元胞自動(dòng)機(jī)（CA）模型[11-12]。鄒陽等[13]采用相場法模擬制冷劑中冰晶的生長過程，結(jié)果與實(shí)驗(yàn)現(xiàn)象較為吻合。與相場模型相比，元胞模型計(jì)算效率更高，并在ZHU等[14]的改進(jìn)下已具備定量模擬能力。對于流場部分，格子玻爾茲曼（LBM）因其效率高、并行能力強(qiáng)、易處理復(fù)雜邊界的特點(diǎn)，更適用于介觀尺度下的流場計(jì)算。康利云等[15]通過建立LBM傳熱模型，計(jì)算泡沫多孔材料有效導(dǎo)熱系數(shù)，具有較高精度。

本文采用格子玻耳茲曼和元胞自動(dòng)機(jī)相結(jié)合的模型，并結(jié)合氮氧凝固相圖研究固空凝結(jié)過程中具體的氮氧體積分?jǐn)?shù)分布，討論過冷度、初始氧體積分?jǐn)?shù)和自然對流對凝固的影響。

1 數(shù)學(xué)模型

實(shí)際情況下，空氣進(jìn)入液氫后轉(zhuǎn)變?yōu)楣炭盏倪^程可能有兩種，一種是先冷凝后凝固，另一種是凝華。為簡化分析，僅考慮氮氧混合液在一定過冷度或冷卻速率下的凝固過程。凝固過程同時(shí)涉及溶質(zhì)場、溫度場的變化，氮氧混合液物理參數(shù)也會相應(yīng)地變化。若全面考慮凝固對物理參數(shù)的影響，將極大地增加計(jì)算量，因此在模擬中對相關(guān)參數(shù)進(jìn)行簡化處理。

同時(shí)為保證流場計(jì)算的穩(wěn)定和效率，作如下假設(shè)：1）液相為不可壓牛頓流體；2）固液兩相熱擴(kuò)散系數(shù)相等且為常數(shù)；3）忽略溫度和體積分?jǐn)?shù)變化對液相的影響；4）忽略晶體在液相中的移動(dòng)。

1.1 LBM模型

本文采用基于LBGK近似的LBM模型用于計(jì)算流場，其演化方程為：

式中：

fi(x,t)——粒子分布函數(shù)；

Δt——時(shí)間步長，s；

ei——離散速度，m/s；

τf——無量綱松弛時(shí)間；

Fi——離散速度空間外力源項(xiàng)，N；

cs——格子聲速，m/s，即

c——格子速度，其值為1 m/s；

F——外力項(xiàng)，N，在本研究中具體為因體積分?jǐn)?shù)差及溫差產(chǎn)生浮力作用。

根據(jù)Boussinesq近似，浮力項(xiàng)可表示為[16]：

式中：

ρ0——初始密度，kg/m3；

βT——溫度膨脹系數(shù)，K-1；

βC——體積分?jǐn)?shù)膨脹系數(shù)，vol-1；

g——重力加速度，m/s2；

C0——初始體積分?jǐn)?shù)；

T0——初始溫度，K。

類似地，體積分?jǐn)?shù)場及溫度場計(jì)算方程也可以轉(zhuǎn)換成對應(yīng)的LBM演化方程，具體表達(dá)式為：

式中：

gi(x,t)——體積分?jǐn)?shù)分布函數(shù)；

hi(x,t)——溫度函數(shù)；

τD、τT——體積分?jǐn)?shù)場及溫度場的松弛時(shí)間；

Gi、Hi——凝固過程中一個(gè)時(shí)間步長釋放的溶質(zhì)和潛熱后增加的體積分?jǐn)?shù)及溫度源項(xiàng)。

流場、體積分?jǐn)?shù)場、溫度場的松弛時(shí)間可分別由動(dòng)力黏度ν、氧的擴(kuò)散系數(shù)D1、熱擴(kuò)散系數(shù)α求得，即下式：

該模型的宏觀密度ρ、速度u、體積分?jǐn)?shù)Cl、溫度T定義如下：

本文采用D2Q9離散速度模型，其速度配置為：

由此可確定粒子、體積分?jǐn)?shù)及溫度平衡態(tài)分布函數(shù)為：

式中，wi為權(quán)系數(shù)，與離散速度方向相對應(yīng)的值為w0=4/9，w1-4=1/9，w5-8=1/36。由上述(1)～(16)可耦合流場、體積分?jǐn)?shù)場、溫度場。而研究中采用被動(dòng)標(biāo)量法，通過速度場u的改變繼而更新式(15)和式(16)，從而體現(xiàn)流場對體積分?jǐn)?shù)及溫度場的影響。

1.2 CA模型

根據(jù)ZS模型[17]，凝固生長的驅(qū)動(dòng)力來源于界面處實(shí)際體積分?jǐn)?shù)與平衡體積分?jǐn)?shù)之差。在界面處過冷度由熱過冷、成分過冷、曲率過冷組成，并滿足熱力學(xué)平衡，由此得到界面處平衡體積分?jǐn)?shù)：

式中：

C1*——界面元胞平衡體積分?jǐn)?shù)；

C0——初始氧體積分?jǐn)?shù)；

T1*——界面元胞溫度，K；

T1eq——初始氧體積分?jǐn)?shù)所對應(yīng)的液相線溫度，K；

m1——液相線斜率；

Γ——Gibbs-Thomson系數(shù)，m…K；

d(θ)——界面能各向異性函數(shù)，此處考慮到氮氧晶體在低溫下呈現(xiàn)立方晶體形狀；

δs——表面能各向異性強(qiáng)度；

θ0——擇優(yōu)生長方向與x軸的夾角，°；

K——界面曲率，m-1，根據(jù)其物理意義表示為：

由式(19)可得到ZS模型中界面前沿新增固相分?jǐn)?shù)，同時(shí)控制固相分?jǐn)?shù)不大于1，得到單位時(shí)間新增固相分?jǐn)?shù)sdf計(jì)算式(20)和式(21)：

由此可以得到式(4)和式(5)中體積分?jǐn)?shù)源項(xiàng)Gi及溫度源項(xiàng)Hi：

式中：

k——溶質(zhì)再分配系數(shù)；

C1——氮氧混合液中氧份數(shù)；

L——氮氧混合液凝固潛熱，kJ/kg；

Cp——氮氧混合液定壓比熱容，kJ/(kg…K)。

1.3 CA與LBM的耦合

在凝固初始階段，在凝固場中放置一個(gè)氮氧晶粒，固相中氧體積分?jǐn)?shù)為kC0，其余區(qū)域充滿過冷度為ΔT0的過冷氮氧混合液。由于氧在液氮中擴(kuò)散系數(shù)是其在固氮中擴(kuò)散系數(shù)的103數(shù)量級，因此忽略氧在固態(tài)晶體中的擴(kuò)散過程。隨著凝固前沿不斷推進(jìn)，界面元胞析出溶質(zhì)氧至前沿，同時(shí)釋放潛熱產(chǎn)生溫度增量，造成體積分?jǐn)?shù)和溫度分布不均，最終導(dǎo)致凝固場密度不均。低密度液相上浮，高密度下降，產(chǎn)生自然對流。而對流對晶體生長的影響主要通過速度u改變體積分?jǐn)?shù)、溫度分布函數(shù)，進(jìn)而影響晶體生長。具體物理參數(shù)的設(shè)置見表1，為簡化計(jì)算，忽略溫度等對相關(guān)物理量的影響而設(shè)為定值。

表1 模擬所用物性參數(shù)[18]

1.4 邊界條件

氮氧凝固場邊界條件見圖1，四周為恒壁溫邊界，溫度恒定為T0，采用非平衡態(tài)外推格式實(shí)現(xiàn)；體積分?jǐn)?shù)設(shè)為無擴(kuò)散邊界，即流場設(shè)為無滑移邊界條件，即ux=uy=0，也通過非平衡態(tài)外推格式計(jì)算。當(dāng)晶體逐漸生長呈復(fù)雜形貌時(shí)，流場與晶體間的固液(S/L)界面設(shè)為無滑移邊界，并用反彈格式處理。忽略氧在固相中的擴(kuò)散作用，氧分子無法以擴(kuò)散形式通過S/L界面，故同樣采用反彈格式處理界面體積分?jǐn)?shù)邊界。

圖1 凝固場邊界示意圖

2 結(jié)果與分析

2.1 純擴(kuò)散下氮氧晶粒的生長

為研究空氣比例的氮氧混合物在凝固過程中體積分?jǐn)?shù)分布，模擬時(shí)取氮氧混合液中初始氧體積分?jǐn)?shù)為C0= 0.21，溫度過冷度ΔT0= 0.8 K，界面能各向異性系數(shù)δs= 0.2，Gibbs-Thomson 系數(shù)Γ = 1.0×10-7m…K，Δx = 3×10-7m，整個(gè)凝固場劃分300×300個(gè)均分網(wǎng)格。圖2為固相率為0.04時(shí)氮氧晶體形貌及氧體積分?jǐn)?shù)分布，可見氮氧晶體在形貌和體積分?jǐn)?shù)分布上均呈現(xiàn)四方對稱性。隨著凝固的進(jìn)行，界面處不斷析出溶質(zhì)，而氧在液相中的擴(kuò)散系數(shù)較小，因而無法及時(shí)擴(kuò)散，在外沿處形成一層氧的高體積分?jǐn)?shù)區(qū)。由于研究中將界面體積分?jǐn)?shù)與平衡體積分?jǐn)?shù)差C1-C1eq作為凝固驅(qū)動(dòng)勢，凝固生長出的晶體臂阻礙中間區(qū)域氧的擴(kuò)散，因而氧在根部富集形成“頸縮”現(xiàn)象[19]。相反地，晶體臂前端析出的氧更易與大空間低體積分?jǐn)?shù)區(qū)接觸，故生長較為迅速。

圖3為不同凝固時(shí)刻下y=0.019 μm處氧體積分?jǐn)?shù)分布圖。圖中所示的曲線圖呈現(xiàn)中間低兩邊高的“凹”字形特征，即中間貧氧、外部富氧。分析可知兩側(cè)外側(cè)氧體積分?jǐn)?shù)恒定為0.21，此處為未發(fā)生凝固的液相。隨著凝固時(shí)間延長，S/L界面推移，高體積分?jǐn)?shù)所對應(yīng)的位置向兩側(cè)移動(dòng)，外側(cè)恒定體積分?jǐn)?shù)區(qū)域縮短。從外側(cè)環(huán)境體積分?jǐn)?shù)至體積分?jǐn)?shù)峰值之間形成一定的體積分?jǐn)?shù)梯度，由氧的有限擴(kuò)散能力所致。當(dāng)內(nèi)側(cè)經(jīng)過高濃度界面處后，體積分?jǐn)?shù)驟然下降，這是由于晶體界面厚度一般為納米級，實(shí)際工程中可以忽略。采用的CA模型為尖銳界面模型，忽略界面厚度，故內(nèi)側(cè)低體積分?jǐn)?shù)即為固態(tài)晶體中氧體積分?jǐn)?shù)值。此外，y = 0.019 μm處氧體積分?jǐn)?shù)最高為0.244，而在整個(gè)凝固場體積分?jǐn)?shù)最高值為0.273，均已超過LITCHFIELD等[20]通過槍擊試驗(yàn)所得到氧體積分?jǐn)?shù)安全閾值0.24，存在爆炸風(fēng)險(xiǎn)。凝固前沿與低體積分?jǐn)?shù)液相接觸面積大，且體積分?jǐn)?shù)勢差較大、擴(kuò)散快，因而凝固時(shí)間對于體積分?jǐn)?shù)峰值及固相中氧體積分?jǐn)?shù)影響均不大。

圖2 氮氧晶體形貌及凝固場氧體積分?jǐn)?shù)分布圖

圖3 不同時(shí)刻下y = 0.019 μm處氧體積分?jǐn)?shù)曲線

2.2 初始氧體積分?jǐn)?shù)的影響

為研究不同初始氧體積分?jǐn)?shù)對于氮氧晶體形成的影響，模擬中選取3個(gè)不同初始體積分?jǐn)?shù)C0分別為0.20、0.25和0.30，控制初始過冷度為ΔT0=0.8 K，并以凝固速率以及氮氧體積分?jǐn)?shù)分布為研究對象。

由圖4可見，固相率隨時(shí)間以近似正比關(guān)系增長。固相率隨時(shí)間增長斜率可視為凝固速率。C0=0.20、0.25、0.30下凝固速率對應(yīng)為0.135 K/s、0.042 K/s、0.023 K/s。隨著初始氧體積分?jǐn)?shù)的增加，凝固速率卻下降。當(dāng)凝固時(shí)間為0.44 s時(shí)，C0=0.20對應(yīng)的固相率為0.044，而此時(shí)C0=0.30所對應(yīng)的固相率為0.011，兩者相差3倍。這是由于凝固時(shí)析出的氧在晶體表面富集，初始體積分?jǐn)?shù)越大意味著晶體界面附近體積分?jǐn)?shù)值越高。在擴(kuò)散系數(shù)一定的條件下，初始體積分?jǐn)?shù)越高，生長越緩慢。由式(20)可知，固相率增量與初始體積分?jǐn)?shù)值成反比關(guān)系，因而凝固至相同固相率，初始體積分?jǐn)?shù)越大所需的時(shí)間越久。

圖4 不同初始氧體積分?jǐn)?shù)下固相率隨時(shí)間變化圖

圖5中，C0=0.20、0.25、0.30分別對應(yīng)的體積分?jǐn)?shù)峰值為0.23、0.28、0.34，固相體積分?jǐn)?shù)為0.17、0.20、0.23。由此可見，初始體積分?jǐn)?shù)對體積分?jǐn)?shù)峰值及固相體積分?jǐn)?shù)的影響程度幾乎一致，同時(shí)可推測氧體積分?jǐn)?shù)峰值Cmax正比于C0，若固空中氧體積分?jǐn)?shù)安全閾值為0.40，則混入的氮氧混合物雜質(zhì)中氧體積分?jǐn)?shù)應(yīng)小于0.34。

圖5 不同初始氧體積分?jǐn)?shù)下y = 0.019 μm處 氧體積分?jǐn)?shù)分布變化圖

2.3 初始過冷度的影響

為研究不同初始過冷度對于氮氧晶體形成的影響，模擬中選取3個(gè)不同過冷度ΔT0=0.5 K、0.7 K、0.8 K，控制初始氧體積分?jǐn)?shù)為C0=0.20，并以凝固速率以及氮氧體積分?jǐn)?shù)分布為研究對象。

由圖6可見，隨著初始過冷度的增加，凝固速率迅速上升。相同凝固時(shí)間為0.44 s，ΔT0=0.8 K時(shí)固相率為0.044，而ΔT0=0.5 K時(shí)固相率僅為0.006，兩者相差約7倍。ΔT0=0.5 K、0.7 K、0.8 K凝固速率分別為0.017 K/s、0.054 K/s、0.135 K/s。由此可見，過冷度相比于初始體積分?jǐn)?shù)對于凝固速率的影響大得多。從式(17)可以看出，溫度對于凝固的影響體現(xiàn)在ΔT/ml，而體積分?jǐn)?shù)影響主要體現(xiàn)在C1-C1eq上，而液相線斜率ml小于1，故數(shù)值上溫度變化影響較初始體積分?jǐn)?shù)更大些。

圖6 不同過冷度下固相率隨時(shí)間變化圖

從圖7可見，初始過冷度增加，凝固生成的固相中氧體積分?jǐn)?shù)也增加，而對于凝固場體積分?jǐn)?shù)峰值，ΔT0=0.5 K時(shí)峰值最小，ΔT0=0.7 K和ΔT0=0.8 K時(shí)峰值幾乎一致。由此可以推測相同初始成分氮氧混合物，過冷度增加時(shí)，晶體外沿的氧體積分?jǐn)?shù)增加，此時(shí)安全風(fēng)險(xiǎn)越大。

2.4 自然對流的影響

當(dāng)將晶體凝固過程中引發(fā)的自然對流考慮在內(nèi)時(shí)，從圖8可見，晶體呈非對稱性生長。凝固初期，計(jì)算域內(nèi)體積分?jǐn)?shù)和溫度差異不大，因而自然對流強(qiáng)度較弱，區(qū)域內(nèi)流速較小。從圖8(a)中流線上看，左右側(cè)分布有兩個(gè)渦，流線相對稱。此時(shí)氮氧晶體近似于純擴(kuò)散時(shí)生長，在4個(gè)方向上生長程度相近。經(jīng)過一段時(shí)間后，明顯下側(cè)方向生長較上側(cè)快，上側(cè)幾乎不生長。這是由于自然對流由密度差造成，而密度差的產(chǎn)生源于凝固過程中氧體積分?jǐn)?shù)差及溫度差。凝固中不斷釋放的熱量和析出的溶質(zhì)使界面溫度和體積分?jǐn)?shù)持續(xù)上升，促進(jìn)了自然對流的增強(qiáng)。對流也影響體積分?jǐn)?shù)和溫度的分布，使高氧體積分?jǐn)?shù)、高溫度的流體上浮至流場上游，抑制晶體上端部分生長。而下游的氧體積分?jǐn)?shù)較小溫度較低，更利于晶體的生長。從圖8(d)流線可看出，隨著晶體左右兩側(cè)的生長，流場將在晶體臂的阻隔作用下由原來的2個(gè)渦轉(zhuǎn)變?yōu)?個(gè)渦。因而自然對流對氮氧晶體凝固形態(tài)產(chǎn)生重要影響，也強(qiáng)化了溶質(zhì)和熱量的傳輸，導(dǎo)致豎直方向體積分?jǐn)?shù)梯度較水平方向的梯度大得多。下面以x=0.019 μm處的體積分?jǐn)?shù)分布研究自然對流對氧體積分?jǐn)?shù)峰值及氮氧分布的影響。

圖7 不同過冷度下y=0.019 μm處氧體積分?jǐn)?shù)分布變化圖

圖8 不同時(shí)刻下自然對流凝固場氮氧晶體形貌氧體積分?jǐn)?shù)分布圖

從圖9可以看出，自然對流下固相范圍比純擴(kuò)散時(shí)小。這是因?yàn)榫w凝固生長驅(qū)動(dòng)力主要來源于成分過冷和熱過冷，而在氮氧凝固時(shí)熱過冷相比于成分過冷，數(shù)值上更具主導(dǎo)性。在純擴(kuò)散中將計(jì)算區(qū)域溫度視為恒定過冷度，而實(shí)際凝固時(shí)釋放潛熱產(chǎn)生溫度增量，在冷卻速率有限時(shí)，將引起全場平均溫度升高，抑制晶體平均生長速率。因此盡管自然對流促進(jìn)氧分子在液相中的擴(kuò)散，界面體積分?jǐn)?shù)下降，但成分過冷對生長促進(jìn)程度不足以彌補(bǔ)溫度升高而產(chǎn)生的抑制作用。另一方面，下游較上游固相區(qū)域更大，這與此前分析結(jié)果一致。在氧體積分?jǐn)?shù)分布方面，自然對流降低界面氧體積分?jǐn)?shù)梯度，促使整個(gè)計(jì)算液相中氧趨于一致。最高氧體積分?jǐn)?shù)位于上游界面處，約為0.23，界面上下游峰值體積分?jǐn)?shù)差約0.01。

圖9 純擴(kuò)散和自然對流下x=0.019 μm處氧體積分?jǐn)?shù)對比圖

3 結(jié)論

1）本文建立CA-LBM模型，模擬單個(gè)氮氧晶體在一定過冷度下生長過程，并研究晶體生長形貌以及凝固場氧體積分?jǐn)?shù)分布情況。在氮氧混合液氧體積分?jǐn)?shù)為空氣中比例時(shí)，新生成的晶體呈內(nèi)部貧氧、外部富氧，這與劉海生等[3]所得實(shí)驗(yàn)結(jié)果一致，而晶體界面處局部體積分?jǐn)?shù)峰值已達(dá)0.244，存在一定安全風(fēng)險(xiǎn)。

2）模擬單個(gè)氮氧晶體在不同過冷度、不同初始氧體積分?jǐn)?shù)條件下的生長，分析得到凝固速率及凝固場氧體積分?jǐn)?shù)峰值隨過冷度增大而增大；初始氧體積分?jǐn)?shù)越大，凝固速率越小，氧體積分?jǐn)?shù)峰值越大。

3）將自然對流考慮在凝固過程中，發(fā)現(xiàn)晶體呈非對稱性生長，上游較下游生長緩慢。自然對流使晶體凝固速率下降，同時(shí)內(nèi)部流場的存在使得整個(gè)凝固場的氧體積分?jǐn)?shù)趨向均勻。一定過冷度下，與純擴(kuò)散相比，自然對流下晶體內(nèi)部氧體積分?jǐn)?shù)減小，界面處氧體積分?jǐn)?shù)峰值也下降，但外部流場平均體積分?jǐn)?shù)上升。