基于STORM試驗的GASFLOW氣溶膠再懸浮模型評估與改進

黃 挺，江 斌，陳 煉

(國核華清(北京)核電技術(shù)研發(fā)中心有限公司，北京 102209)

核反應(yīng)堆發(fā)生嚴重事故期間，在氫氣燃燒核安全殼泄壓等特定場景下，安全殼內(nèi)氣流擾動劇烈，會導(dǎo)致大量已沉寂的氣溶膠再懸浮，安全殼空間內(nèi)放射性源項增加，最終可能導(dǎo)致釋放到環(huán)境的源項增加[1]。為更加準確地對事故晚期源項進行評估，需開發(fā)模擬氣溶膠再懸浮行為的計算模型，并利用相關(guān)試驗進行評估及驗證。目前已有部分程序或模型經(jīng)過了試驗驗證，并在反應(yīng)堆嚴重事故分析程序中獲得了一定應(yīng)用[2-4]。

GASFLOW是由美國洛斯阿拉莫斯國家實驗室(LANL)和德國卡爾斯魯厄研究中心(FzK)共同開發(fā)的三維計算流體力學程序[5]，除可模擬氫氣的擴散、混合分布、分層及燃燒等三維流動現(xiàn)象外，還針對氣溶膠的傳輸、沉降及再懸浮等建立了相應(yīng)的分析模型[6]。GASFLOW中氣溶膠再懸浮模型基于力學平衡原理，并根據(jù)Caberjos和Klingzing的試驗研究結(jié)果進行了修正[7]。該試驗氣溶膠工質(zhì)的粒徑范圍為10～1 000 μm，相對較大，對于核反應(yīng)堆而言，嚴重事故下氣溶膠的粒徑大多集中在0.1～1 μm范圍內(nèi)。因此，GASFLOW中再懸浮模型是否適用于該粒徑范圍需進行進一步評估及驗證。

STORM試驗中的氣溶膠再懸浮試驗被廣泛應(yīng)用于氣溶膠再懸浮模型的驗證中，其試驗工況涵蓋了氣溶膠粒徑小于10 μm的條件[2]，且與大型非能動先進壓水堆在大破口失水事故后安全殼內(nèi)氣溶膠粒徑分布較為一致[8]。國際標準例題ISP40報告中，有近10種程序或模型選取STORM試驗中的SR11工況作為基準進行了對比計算[9]。國內(nèi)近些年的模型評估及驗證工作，也基本使用了STORM試驗的研究結(jié)果[10-12]。

本文以GASFLOW為基礎(chǔ)，通過與STORM試驗中SR11工況的試驗數(shù)據(jù)進行對比研究，對程序中再懸浮模型的適用性進行評估，并基于試驗數(shù)據(jù)對模型進行改進及驗證。

1 模型描述

GASFLOW中通過獲取沉積壁面的固體顆粒在平行壁面氣流中的最小拾取速度來判定其是否發(fā)生再懸浮[7]。圖1示出穩(wěn)定且完全發(fā)展的湍流場中，靜止于壁面的單一球形顆粒的受力分析，所受各力如下。

圖1 穩(wěn)定且完全發(fā)展的湍流場中靜止于壁面的單一球形顆粒的受力Fig.1 Force acting on single sphere at rest on wall with steady, fully developed turbulent flow

重力Fg為：

(1)

浮力Fb為：

(2)

黏附力Fa[13]為：

(3)

拖曳力Fd為：

(4)

升力Fl[14]為：

(5)

摩擦力Ff為：

Ff=fsFn=fs(Fg+Fa-Fb-Fl)

(6)

當拖曳力等于摩擦力時，氣溶膠顆粒開始運動并離開沉積表面，此時：

Fd=fs(Fg+Fa-Fb-Fl)

(7)

將式(1)～(5)代入式(7)可得：

(8)

式(8)是以Ugcp為未知數(shù)的方程，通過求解該方程可得出單一顆粒產(chǎn)生再懸浮的最小氣體流速，將其用Ugpu0表示，即當Ugcp大于該速度時，氣溶膠粒子就可懸浮。

對Ugpu0的修正參考了Caberjos和Klingzing的試驗研究結(jié)果，其修正[15]如下：

Ugpu=fArUgpu0

(9)

(10)

其中，Ar為阿基米德數(shù)。

(11)

2 幾何建模及計算條件

2.1 幾何建模及假設(shè)

STORM試驗分為兩個獨立階段：第1階段研究氣溶膠沉積現(xiàn)象，以此獲取再懸浮試驗所需的沉積邊界條件；第2階段研究不同氣體流速下的氣溶膠再懸浮現(xiàn)象。STORM試驗的試驗段是一內(nèi)徑63 mm、長度5.005 5 m的細長管道。在再懸浮試驗階段，高溫氣體從一側(cè)管口進入，攜帶再懸浮的氣溶膠顆粒從另一側(cè)管口流出，根據(jù)管內(nèi)沉積氣溶膠質(zhì)量的減少來獲取不同氣體流速下的再懸浮率[8]。

管內(nèi)的氣體流速是影響氣溶膠再懸浮的最主要因素之一。由于管內(nèi)氣體流動方向與管壁平行，因此y和z方向的節(jié)點劃分對網(wǎng)格內(nèi)氣體流速的影響幾乎可忽略不計。x方向網(wǎng)格節(jié)點的劃分雖然對氣流速度場的分布有一定影響，但對每時間段內(nèi)的最高氣體流速的影響較為有限，對再懸浮率計算結(jié)果的影響較小。由于SR11工況中氣體最高流速約為130 m/s，過度加密x方向節(jié)點會在增加網(wǎng)格數(shù)量的同時減小時間步長，從而導(dǎo)致計算時間過長。因此在使用GASFLOW進行建模時，綜合考慮網(wǎng)格規(guī)模及計算時間的可接受性，對STORM試驗段模型進行簡化，將其假設(shè)為一正方形管口的管道，管口邊長為55.83 mm，使其流動面積與設(shè)計參數(shù)保持一致。在x、y和z方向分別設(shè)置25、9和9個節(jié)點，即1 536(24×8×8)個有限元網(wǎng)格，如圖2所示，假設(shè)氣溶膠均勻沉積在管道底部。

圖2 STORM試驗段網(wǎng)格劃分示意圖Fig.2 Schematic diagram of STORM test section mesh

2.2 計算條件

表1中列出SR11工況的初始條件和邊界條件[9]。

表1 SR11工況的初始條件和邊界條件Table 1 Initial condition and boundary condition of SR11 case

續(xù)表1

針對表1中的溫度范圍，在計算過程中取最低溫度和最高溫度的均值。除上述試驗參數(shù)外，表2列出GASFLOW計算所需的初始條件。

表2 GASFLOW計算初始條件Table 2 Initial condition for GASFLOW calculation

2.3 模型選取

GASFLOW中有兩個湍流模型可供選擇，即代數(shù)模型和k-ε模型。本文研究中選用計算速度相對較快的代數(shù)模型。在計算氣溶膠顆粒輸運時，GASFLOW具有單向動量耦合和雙向動量耦合兩個模型。其中，單向動量耦合模型不考慮離散粒子相對連續(xù)流場相的影響。考慮到本文模型中的氣溶膠粒子無論在體積份額還是質(zhì)量份額方面與流場中氣體相比都很小，對流場的影響可忽略不計，因此選擇單向動量耦合模型進行計算。

3 計算結(jié)果及分析

3.1 模型評估

GASFLOW再懸浮率計算結(jié)果與STORM試驗SR11工況結(jié)果對比如圖3所示。由圖3a可見：GASFLOW計算得出的再懸浮率隨時間的增長而不斷增大，計算最終得出的再懸浮率為51.2%，與試驗最終的74.1%有較為明顯的差距。由圖3b可見：當氣體流速為106 m/s時，GASFLOW計算結(jié)果與SR11工況結(jié)果最為接近；當氣體流速小于106 m/s時，計算結(jié)果與SR11工況結(jié)果相比偏高；當氣體流速大于106 m/s時，計算結(jié)果與SR11工況結(jié)果相比偏低。

管內(nèi)壁表面粗糙度是影響再懸浮率的主要參數(shù)之一，在GASFLOW再懸浮模型中表現(xiàn)為拖曳力系數(shù)粗糙度因子(cdrf)和顆粒與沉積表面的平衡間距粗糙度因子(sdzrf)。敏感性分析結(jié)果表明[16]，再懸浮率隨cdrf和sdzrf均呈單調(diào)變化，cdrf和sdzrf增大，則再懸浮率增大，反之則減小。通過調(diào)節(jié)敏感性參數(shù)cdrf和sdzrf并不能改變圖3中GASFLOW計算結(jié)果的變化趨勢，使其與試驗數(shù)據(jù)更加接近。

圖3 GASFLOW計算結(jié)果與SR11工況結(jié)果對比Fig.3 Comparison between calculation result of GASFLOW and SR11 case result

綜上可知，GASFLOW中的氣溶膠再懸浮模型在應(yīng)用于粒徑小于10 μm的條件時，其分析結(jié)果與試驗數(shù)據(jù)存在較為明顯的差距。為進一步提升GASFLOW中的氣溶膠再懸浮模擬的準確性，利用SR11工況的試驗結(jié)果對模型進行改進。

3.2 模型改進

根據(jù)ISP40報告中所記錄的SR11工況試驗數(shù)據(jù)可知，隨氣體流速的增加，仍沉積在試驗段管壁的氣溶膠顆粒的幾何平均粒徑不斷減小[8]。這說明在特定氣體流速下，氣溶膠顆粒的粒徑是決定其是否發(fā)生再懸浮的重要因素之一，且隨氣體流速的增大，可發(fā)生再懸浮的氣溶膠粒徑減小。依照上述分析，可在試驗數(shù)據(jù)基礎(chǔ)上利用氣體流速與發(fā)生再懸浮氣溶膠粒徑的相關(guān)性，對GASFLOW中的氣溶膠再懸浮模型進行改進，使改進后的模型分析結(jié)果與試驗數(shù)據(jù)更為吻合。模型改進方法如下。

1) 假設(shè)特定氣體流速Ug下存在再懸浮氣溶膠的最小粒徑dmin，即粒徑大于dmin的顆粒均發(fā)生再懸浮。再懸浮率可等價于粒徑大于等于dmin的累積概率。

圖4 不同氣體流速下的再懸浮最小粒徑Fig.4 Minimum resuspension particle size at different air velocities

2) 根據(jù)氣溶膠粒徑的累積概率分布曲線與圖3a中不同階段的再懸浮率，得出不同再懸浮率所對應(yīng)的dmin，再利用圖3b中再懸浮率與Ug的對應(yīng)關(guān)系，可得出Ug與dmin的關(guān)系曲線，如圖4所示。

3) 利用GASFLOW計算氣體流速Ug條件下粒徑為dmin的氣溶膠顆粒所對應(yīng)的Ugcp和Ugpu0。在該條件下，粒徑為dmin的氣溶膠顆粒處于發(fā)生再懸浮的臨界狀態(tài)，由此可得Ugpu=Ugcp，式(9)中的修正系數(shù)可以表示為fAr=Ugcp/Ugpu0。

4) 根據(jù)粒徑dmin和式(11)計算出對應(yīng)的Ar，之后根據(jù)fAr和Ar關(guān)系曲線(圖5)擬合出新的關(guān)系式并應(yīng)用于模型中。

圖5 fAr與Ar的關(guān)系曲線Fig.5 Relation curve between fAr and Ar

由圖5可見，fAr與Ar在雙對數(shù)坐標下基本呈線性關(guān)系，可將fAr表示為Ar的冪函數(shù)，通過對圖5中數(shù)據(jù)進行擬合可得：

fAr=34.02Ar-0.168 8

(12)

3.3 模型驗證

將GASFLOW中氣溶膠再懸浮模型中的式(10)替換為式(12)，重新對SR11的再懸浮工況進行模擬計算。改進后模型的計算結(jié)果與SR11工況的對比如圖6所示。

由圖6可見，與GASFLOW中的原模型相比，改進后模型的再懸浮率計算結(jié)果無論是從變化趨勢還是數(shù)值上均與試驗結(jié)果更加接近。如計算最終得出的再懸浮率為75.12%，與原來的51.2%相比，與試驗值74.1%的誤差更小，其變化趨勢也與試驗數(shù)據(jù)更加吻合。

需要明確的是，氣溶膠再懸浮現(xiàn)象發(fā)生的機理十分復(fù)雜，粒徑與氣體流速只是其中較為重要的兩個影響因素，除此之外還受到其他多種因素的影響(如顆粒黏附力的分布形式，氣溶膠的沉積形態(tài)，氣溶膠顆粒之間的相互作用，氣流方向與沉積平面的角度等)，無論試驗還是數(shù)值分析，均存在很大的不確定性。上述對GASFLOW中氣溶膠再懸浮模型的改進，雖然計算結(jié)果在數(shù)值上更加接近試驗數(shù)據(jù)，但仍有進一步改進的空間，這需要對氣溶膠再懸浮的機理有更加深入的認識。

圖6 改進后GASFLOW的計算分析結(jié)果與SR11工況結(jié)果對比Fig.6 Comparison between calculation result of improved GASFLOW and SR11 case result

此外，STORM試驗的試驗段為細長管，其氣體流速也始終平行于沉積面，這與嚴重事故下安全殼內(nèi)的環(huán)境條件存在顯著差距。模型分析的準確性離不開試驗的驗證，因此開展更為接近實際狀態(tài)的大容器內(nèi)氣溶膠再懸浮試驗是十分必要的。

4 結(jié)論及建議

本文利用ISP40中STORM試驗的SR11工況對GASFLOW中的再懸浮模型進行了適用性評估，并根據(jù)試驗數(shù)據(jù)所反映出的氣體流速與再懸浮最小粒徑的關(guān)系，對GASFLOW中的再懸浮模型進行了改進，所得結(jié)論及建議如下：

1) GASFLOW中基于Caberjos和Klingzing的試驗研究結(jié)果所修正的再懸浮模型，適用粒徑范圍為10～1 000 μm，對于小于10 μm的情況，該模型分析結(jié)果與ISP40中試驗結(jié)果的差距較為明顯；

2) 利用ISP40中試驗數(shù)據(jù)擬合得出修正系數(shù)fAr的新關(guān)系式，GASFLOW中的再懸浮模型經(jīng)改進后，其計算結(jié)果無論在變化趨勢還是數(shù)值上，均與試驗結(jié)果更加吻合；

3) 建議開展更為接近嚴重事故下安全殼內(nèi)實際狀態(tài)的氣溶膠再懸浮試驗，并基于試驗結(jié)果對分析模型進行驗證及改進，擴展現(xiàn)有再懸浮模型的適用范圍，提升模型分析的準確性。