鈉滴燃燒仿真模型研究

王 超，張智剛，彭康瑋，孫樹斌

(哈爾濱工程大學 核安全與仿真技術(shù)國防重點學科實驗室，黑龍江 哈爾濱 150001)

鈉泄漏和鈉火事故與鈉冷快堆發(fā)展密切相關(guān)，鈉泄漏和鈉火事故一直受到國際上以及中國實驗快堆安全分析的重視。而單一鈉滴的燃燒作為鈉火，特別是噴霧鈉火研究的基礎，其燃燒機理更應受到重視。在鈉滴傳統(tǒng)計算模型中使用的是球狀模型，而鈉滴在燃燒過程中呈現(xiàn)非球形態(tài)，因此仿真模擬精度有待提高。鈉滴在燃燒過程中表面生成錐狀氧化物[1]，故本文采用更接近實際燃燒特性的帶有錐狀的新模型，模擬計算在不同鈉滴初始溫度和氧氣濃度下鈉滴燃燒溫度隨時間的變化規(guī)律，并與相關(guān)實驗數(shù)據(jù)進行對比分析研究。

1 鈉滴燃燒模型的確定

1.1 幾何模型

在傳統(tǒng)鈉滴計算研究中，均選用球狀鈉滴模型用于分析計算，如圖1a所示。由于在鈉滴燃燒過程中表面會形成氧化物突起，便于與氧氣接觸及熱量傳遞，從相關(guān)鈉滴實驗可觀察到氧化物的形狀接近于錐狀，因此本文使用表面帶有錐狀突起的新模型(簡稱錐狀模型)替代球狀模型，如圖1b所示。

在錐狀模型中，假設如下：錐狀模型與球狀模型的質(zhì)量相等，燃燒過程中生成n個大小相同的錐狀氧化物，鈉滴直徑由D減小為D1；圓錐的底面直徑為aD1(a=0.05)，高為bD1(b=0.1)；考慮到各圓錐底面之間存在一定間隙，假設邊長為aD1的正方形分布于鈉滴表面，因此每個正方形對應一個圓錐底面。由圓錐狀凸起的底面積和球形鈉滴的表面積相互之間的關(guān)系可得：

β4π(D1/2)2=n(aD1)2

(1)

式中，β為圓錐底面積的總和占鈉滴表面積的百分比。

由錐狀模型與球狀模型的質(zhì)量相等可得：

π

(2)

根據(jù)鈉滴燃燒實驗現(xiàn)象中錐狀氧化物的大小，估算確定錐底面直徑和高，假定β=0.9，根據(jù)式(1)、(2)計算錐狀氧化物的個數(shù)n及生成氧化物后鈉滴的直徑D1。

由于表面積與體積之比在燃燒過程中對燃燒速率影響很大，表面積與體積之比越大，燃燒速率越大。以表面積與體積之比作為評價鈉滴燃燒特性的主要因素，因球狀模型在數(shù)學建模和數(shù)值計算上更有優(yōu)勢，對錐狀模型進行球殼狀模型虛擬變換，如圖1c所示。根據(jù)球殼狀模型與錐狀模型表面積相等，可得球殼狀模型的外徑Do，在質(zhì)量及密度相同的條件下，可得球殼狀模型的內(nèi)徑Di。

a——球狀模型；b——錐狀模型；c——球殼狀模型

1.2 單一鈉滴模型

分析鈉滴燃燒的實驗現(xiàn)象[1]可知，鈉與氧的反應一般分為3個階段：表面氧化階段、預燃階段和燃燒階段。根據(jù)鈉滴在不同燃燒階段是否產(chǎn)生錐狀氧化物，分別使用球狀模型和錐狀模型。

實驗證明，鈉滴燃燒階段的直徑變化遵循“D2”定律[2]，由于鈉滴直徑在單位時間內(nèi)的變化極小，為簡化計算，假設鈉滴在空氣中的反應過程中直徑變化始終遵循“D2”定律，即鈉滴直徑的平方隨時間線性遞減[3]：

(3)

式中：D為t時刻鈉滴的直徑，m；D0為鈉滴初始直徑(D0=0.005 m)；Kv為汽化系數(shù)[2,4]，即：

(4)

式中：λ為氣體導熱率，W/(m·K)；cp為氣體比定壓熱容，J/(kg·K)；B為傳遞系數(shù)，即：

(5)

式中：Hfg為汽化潛熱，J/kg；Tg為氣體平均溫度，K；Ts為鈉滴溫度，K；Hc為鈉反應熱，J/kg；y為氧氣摩爾份額；i為化學計量比。

靜態(tài)鈉滴的質(zhì)量變化速率可由汽化系數(shù)Kv表示，為：

(6)

式中：在表面氧化階段，ms為鈉滴反應速率，kg/s；在預燃和燃燒階段，ms為鈉滴氣化速率，kg/s。

鈉滴在表面氧化階段，反應熱一部分傳遞至鈉滴使其溫度升高，一部分傳遞給周圍環(huán)境，包括對流傳熱及熱輻射。由于鈉滴質(zhì)量較小和較高的熱傳導率，根據(jù)集總參數(shù)法可認為鈉滴各處溫度相等。建立能量守恒方程[5]為：

Φdt=mt1cpNadTs

(7)

式中：dt為時間步長，s；Nu為鈉滴表面到周圍氣體的Nusselt數(shù)；k為導熱系數(shù)，W/(m·K)；S為鈉滴表面積，m2；Φ為單位時間內(nèi)發(fā)出的熱輻射熱量，W；cpNa為鈉的比定壓熱容，J/(kg·K)；dTs為時間步長內(nèi)鈉滴溫度的變化，K。左邊第1項為鈉反應釋放的熱量，第2項為通過對流傳熱傳遞給周圍氣體的熱量，第3項為通過熱輻射傳遞給周圍氣體的熱量，等號右邊為鈉滴溫度升高所需的熱量。鈉滴在單位時間內(nèi)發(fā)出的熱輻射熱量Φ為：

(8)

式中，εNa為鈉的發(fā)射率。

鈉滴在表面氧化階段，氧氣不足，根據(jù)實驗現(xiàn)象并無錐狀氧化物生成，因此在本階段，鈉滴燃燒溫度的計算采用球狀模型，則：

πD3ρNa

(9)

式中，mt1為使用球狀模型階段未燃燒的鈉滴質(zhì)量，kg。

鈉滴在預燃階段，隨鈉滴表面溫度的升高，其表面發(fā)生氣化，反應熱一部分傳遞至鈉滴使鈉氣化，因此式(7)可變?yōu)椋?/p>

Φdt=mt2cpNadTs

(10)

考慮到部分發(fā)生氣化的鈉蒸氣并未與氧氣反應，而是直接釋放到環(huán)境中，因此發(fā)生氧化反應的鈉反應速率mp可表示為：

mp=γms

(11)

式中，γ為單位時間參與氧化反應的鈉與發(fā)生氣化的鈉的質(zhì)量比。

在預燃階段，鈉滴表面發(fā)生氣化，且表面產(chǎn)生許多錐狀氧化物，在該階段，適用錐狀模型，并采用球殼狀模型進行求解，則：

πρNa

(12)

式中，mt2為使用錐狀模型階段未燃燒的鈉滴質(zhì)量，kg。

鈉滴在燃燒階段，由于在鈉滴的周圍有一層擴散火焰鋒面，火焰鋒面發(fā)出的熱量同時向內(nèi)和向外導熱[6]，由于火焰鋒面的溫度極高，因此，反應熱中的一部分還需用來支持火焰鋒面的高溫，因此式(10)可變?yōu)椋?/p>

Hcmpdt-mscpNa(Tf-Ts)dt-Hfgmsdt-

(13)

式中，Tf為火焰鋒面溫度，K。

在燃燒階段，鈉滴表面溫度的升高，增加了鈉滴中氧的溶解度，在預燃階段生成的氧化物溶解于鈉滴中，鈉滴表面又變得光滑，在本階段，采用球狀模型進行計算。

此外，假設在鈉滴反應至5～10 s內(nèi)，停止通氧氣，氧氣濃度隨時間開始下降，假設氧氣濃度隨時間線性下降，若初始氧氣濃度為20%，則10 s時下降至16%；若初始為16%，則10 s時下降至12%；若初始為12%，則10 s時下降至10%；若初始為8%和4%，鈉燃燒氧氣濃度基本不變。此階段錐狀氧化物重新產(chǎn)生，在此階段使用球殼狀模型進行計算。為了提高模擬的準確度，需考慮不同氧氣摩爾份額、不同溫度下空氣和鈉的物性參數(shù)(表1)。

2 計算結(jié)果與分析

由于本文在鈉滴的不同反應階段使用的模型不同，因此需在不同鈉滴初始溫度和氧氣濃度的情況下確立各階段反應的開始和結(jié)束時間，根據(jù)鈉滴燃燒的3個階段所表現(xiàn)的特征現(xiàn)象[1]，即表面氧化階段——白色的氧化層覆蓋液滴表面；預燃階段——開始有氧化物產(chǎn)生，并伴隨反應產(chǎn)物的氣溶膠出現(xiàn)和燃燒階段——氧化物表面被橘黃色火焰包圍，隨后氧化物溶入鈉滴中，統(tǒng)計界定鈉滴燃燒的實驗數(shù)據(jù)中各反應階段的開始和結(jié)束時間[1]，如表2和3所列，并將其代入程序計算，得到仿真計算結(jié)果。表2中，初始溫度為300 ℃；表3中，氧氣濃度為20%。

表1 空氣和鈉的物性參數(shù)[7]

表2 不同氧氣份額下鈉滴各反應階段的開始和結(jié)束時間[1]

表3 不同鈉滴初始溫度下鈉滴各反應階段的開始和結(jié)束時間[1]

圖2示出鈉滴初始溫度為300 ℃時不同氧氣濃度下鈉滴溫度隨時間的變化。不同氧氣濃度下鈉滴溫度的變化趨勢均與實驗數(shù)據(jù)相近，但相對于實驗數(shù)據(jù)，鈉滴溫度的模擬結(jié)果在各反應階段的起始時刻均存在明顯的轉(zhuǎn)折，尤其在氧氣濃度為20%、16%和12%時，這是由于在不同反應階段采用了不同模型，使鈉滴在不同階段具有不同的溫度變化速率。在鈉滴預燃和燃燒階段，部分發(fā)生氣化的鈉蒸氣并未與氧氣反應，而是直接釋放到環(huán)境中，如何確定式(11)中的γ，是計算結(jié)果與實驗數(shù)據(jù)能否符合的關(guān)鍵。在預燃階段，由于鈉滴溫度的上升趨勢為非線性，無法準確計算γ，因此本文采用迭代法，將不同的γ代入式(11)，最終得到與實驗數(shù)據(jù)相近的鈉滴溫度變化趨勢。由計算結(jié)果可知，在鈉滴初始溫度為300 ℃及氧氣濃度為20%、16%、12%和8%時，鈉滴預燃階段的γ平均值分別為0.66、0.69、0.72和0.72。

圖2 不同氧氣濃度下鈉滴溫度隨時間的變化

由式(11)可知，γ受氣化作用和氧化作用劇烈程度影響，初始溫度相同，氧氣濃度越高氧化作用越強，溫度上升越快，隨溫度的升高，氣化作用增強，未參與氧化作用的鈉蒸氣份額也會增大，γ就會減小，雖然溫度升高也會增強氧化作用，但氧氣份額具有上限值，因此氧化作用的增強有限，氣化作用的增強為主要影響因素，因此初始溫度一定時，γ隨氧氣濃度的減小而增大。

在鈉滴燃燒階段，由實驗數(shù)據(jù)可知鈉滴的溫度不變，說明燃燒產(chǎn)生的熱與鈉蒸氣溫度升高所需熱量、汽化潛熱及釋放到周圍環(huán)境中的熱量平衡，即式(13)右邊為零。因此在燃燒階段，通過式(13)可計算得到γ。由計算結(jié)果可知，在鈉滴初始溫度為300 ℃及氧氣份額為20%和16%時，鈉滴燃燒階段的γ分別為0.56和0.44。在氧氣濃度隨時間下降的階段，不同氧氣濃度下鈉滴溫度的變化趨勢不同，且均為非線性，因此同樣采用迭代法確定γ，由計算結(jié)果可知，在氧氣濃度為20%、16%、12%和8%時，氧氣濃度隨時間下降的階段，γ分別為0.66、0.60、0.76和0.83。

氧氣濃度為20%時，不同初始溫度下鈉滴溫度隨時間的變化如圖3a所示，鈉滴溫度的變化趨勢與實驗數(shù)據(jù)基本符合。在鈉滴預燃階段，初始溫度為500、400、300和200 ℃時，通過計算可知，γ分別為0.85、0.73、0.66和0.62；在鈉滴燃燒階段，γ分別為0.63、0.59、0.56和0.55；在氧氣濃度隨時間下降階段，γ分別為0.76、0.70、0.66和0.63。

在各反應階段，γ隨初始溫度的升高而增大，這是由于溫度越高，氧化作用越強，因此會有更多的鈉蒸氣參與反應。預燃階段的γ均高于燃燒階段的，這是由于在燃燒階段的反應均發(fā)生在火焰鋒面上，由于反應劇烈，火焰鋒面周圍的氧氣濃度低于20%，使部分鈉蒸氣未燃燒，而是直接釋放到環(huán)境中；預燃階段鈉滴溫度低于燃燒階段，氣化作用較弱，但由于氧氣充足，因此參與反應的鈉較多，即γ較高。氧氣濃度隨時間下降階段的γ受兩種作用共同影響，初始溫度較高時，γ介于預燃階段和燃燒階段之間，初始溫度較低時，γ高于兩個階段。

在氧氣濃度為4%時，鈉滴溫度的上升速率遠高于實驗數(shù)據(jù)，如圖3b所示。由于氧氣不充足，鈉滴一直處于表面氧化階段，鈉滴表面被氧化層覆蓋，使鈉滴與氧氣氧化反應的實際接觸面積小于仿真計算時所使用的鈉滴表面積，鈉滴初始溫度分別為500、400、300和200 ℃時，若達到與實驗數(shù)據(jù)相吻合的鈉滴溫度，則ms分別減小至原來的60%、50%、40%和30%，由式(6)可知，ms正比于D，即正比于表面積的平方根，則鈉滴與氧氣的實際接觸面積占鈉滴表面積的百分比為36%、25%、16%和9%，初始溫度越高，實際接觸面積越大，這是由于溫度高促進了表面氧化物溶解于鈉中，使更多的氧氣和鈉參與了反應。

本文對使用了傳統(tǒng)球狀模型與錐狀模型的計算結(jié)果和實驗數(shù)據(jù)進行了比較分析，如圖4所示。通過對比發(fā)現(xiàn)，在預燃階段，使用球狀模型計算得到的鈉滴溫度略低于實驗數(shù)據(jù)，這是由于在此階段鈉滴表面產(chǎn)生的錐狀氧化物增大了鈉滴表面積，而球狀模型表面積小，導熱能力小于帶有錐狀物的鈉滴，因此溫度上升速率略低于實驗結(jié)果；在氧氣濃度下降階段，實驗結(jié)果中鈉滴的溫度接近恒定，并有下降趨勢，而使用球狀模型的計算結(jié)果中鈉滴溫度繼續(xù)上升，與實際不符。對比可得出使用了錐狀模型的計算結(jié)果與實驗數(shù)據(jù)基本符合，而使用球狀模型的計算結(jié)果與實驗數(shù)據(jù)有一定的差別，說明在鈉滴燃燒計算中錐狀模型具有一定的優(yōu)越性。

圖3 不同初始溫度下鈉滴溫度隨時間的變化

圖4 球狀模型、錐狀模型的計算結(jié)果與實驗數(shù)據(jù)的對比

3 結(jié)論

本文對鈉滴燃燒進行仿真研究，針對表面氧化、預燃及燃燒3個階段分別建立物理數(shù)學模型，通過引入?yún)⑴c氧化反應與發(fā)生氣化的鈉的質(zhì)量比γ，與實驗數(shù)據(jù)進行比較，分析影響鈉滴燃燒的因素和機理，可得如下結(jié)論。

1) 初始溫度一定時，預燃階段的γ隨氧氣濃度的增大而減??；燃燒階段在氧氣濃度為20%和16%時，γ分別為0.56和0.44；氧氣濃度隨時間下降的階段γ受氧化作用和氣化作用共同影響。

2) 氧氣濃度一定，不同初始溫度下，在各反應階段，γ隨初始溫度的升高而增大；預燃階段的γ均高于燃燒階段的；氧氣濃度隨時間下降的階段，初始溫度較高時，γ介于預燃階段和燃燒階段之間，初始溫度較低時，γ高于兩個階段。

3) 在氧氣不充足的條件下，鈉滴與氧氣的接觸面積小于鈉滴的表面積，鈉滴初始溫度為500、400、300和200 ℃時，鈉滴與氧氣實際的接觸面積占鈉滴表面積的百分比約為36%、25%、16%和9%，表明初始溫度越高，接觸面積越大。

4) 采用表面帶有錐狀突起的鈉滴新模型，通過與傳統(tǒng)球狀模型的仿真計算結(jié)果對比可知，錐狀模型的計算結(jié)果與實驗數(shù)據(jù)符合更好，因此錐狀模型在有關(guān)鈉滴燃燒模擬計算中具有優(yōu)越性，還需進一步驗證。

參考文獻：

[1] OTAKE S, NISHIMURA M, SUGIYAMA K, et al. Oxidation and combustion in liquid sodium droplet[C]∥17th International Conference on Nuclear Engineering(ICONE17). America: ASME, 2009: 411-415.

[2] 喻宏，徐銤，金德圭. 鈉冷快堆中噴霧鈉火的分析計算[J]. 核動力工程，2002，23(2)：21-26.

YU Hong, XU Mi, JIN Degui. Sodium spray fire analysis of sodium-cooled fast reactor[J]. Nuclear Power Engineering, 2002, 23(2): 21-26(in Chinese).

[3] 張斌，朱繼洲，韓浪. 鈉冷快堆鈉霧火事故三維數(shù)值模擬[J]. 核動力工程，2005，26(2)：105-109.

ZHANG Bin, ZHU Jizhou, HAN Lang. 3-dimensional numerical simulation of sodium spray fire accidents in LMFBRs[J]. Nuclear Power Engineering, 2005, 26(2): 106-109(in Chinese).

[4] TAKATA T, YAMAGUCHI A, MAEKAWA I. Numerical investigation of multi-dimensional characteristics in sodium combustion[J]. Nuclear Engineering and Design, 2003, 220(1): 37-50.

[5] 王學容，駱純珊，單建強，等. 鈉冷快增殖堆鈉霧火分析計算[J]. 核科學與工程，2001，21(3)：218-222.

WANG Xuerong, LUO Chunshan, SHAN Jianqiang, et al. Analysis and assessment of sodium spray fire in LMFBRs[J]. Chinese Journal of Nuclear Science and Engineering, 2001, 21(3): 218-222(in Chinese).

[6] 徐通模. 燃燒學[M]. 北京：機械工業(yè)出版社，2010.

[7] 楊世銘，陶文銓. 傳熱學[M]. 北京：高等教育出版社，2006.