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

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

0 引言

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

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

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

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

1 數(shù)學模型

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

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

1.1 LBM模型

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

式中：

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

Δt——時間步長，s；

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

τf——無量綱松弛時間；

Fi——離散速度空間外力源項，N；

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

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

F——外力項，N，在本研究中具體為因體積分數(shù)差及溫差產生浮力作用。

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

式中：

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

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

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

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

C0——初始體積分數(shù)；

T0——初始溫度，K。

類似地，體積分數(shù)場及溫度場計算方程也可以轉換成對應的LBM演化方程，具體表達式為：

式中：

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

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

τD、τT——體積分數(shù)場及溫度場的松弛時間；

Gi、Hi——凝固過程中一個時間步長釋放的溶質和潛熱后增加的體積分數(shù)及溫度源項。

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

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

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

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

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

1.2 CA模型

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

式中：

C1*——界面元胞平衡體積分數(shù)；

C0——初始氧體積分數(shù)；

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

T1eq——初始氧體積分數(shù)所對應的液相線溫度，K；

m1——液相線斜率；

?！狦ibbs-Thomson系數(shù)，m…K；

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

δs——表面能各向異性強度；

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

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

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

由此可以得到式(4)和式(5)中體積分數(shù)源項Gi及溫度源項Hi：

式中：

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

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

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

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

1.3 CA與LBM的耦合

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

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

1.4 邊界條件

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

圖1 凝固場邊界示意圖

2 結果與分析

2.1 純擴散下氮氧晶粒的生長

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

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

圖2 氮氧晶體形貌及凝固場氧體積分數(shù)分布圖

圖3 不同時刻下y = 0.019 μm處氧體積分數(shù)曲線

2.2 初始氧體積分數(shù)的影響

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

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

圖4 不同初始氧體積分數(shù)下固相率隨時間變化圖

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

圖5 不同初始氧體積分數(shù)下y = 0.019 μm處 氧體積分數(shù)分布變化圖

2.3 初始過冷度的影響

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

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

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

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

2.4 自然對流的影響

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

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

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

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

圖9 純擴散和自然對流下x=0.019 μm處氧體積分數(shù)對比圖

3 結論

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

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

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