王 輝，邢 繼，孫中寧，谷海峰，孫曉暉，王一博

(1.哈爾濱工程大學 核安全與仿真技術(shù)國防重點學科實驗室，黑龍江 哈爾濱 150001；2.中國核電工程有限公司，北京 100840)

壓水堆核電廠嚴重事故期間堆芯熔化，釋放大量放射性裂變產(chǎn)物進入安全殼，這些裂變產(chǎn)物主要以氣溶膠的形態(tài)在安全殼內(nèi)輸運。當安全殼失效時，放射性氣溶膠釋放到外界環(huán)境，威脅公眾安全。即使安全殼未被破壞，這些氣溶膠也可能隨氣流經(jīng)由混凝土安全殼縫隙、貫穿件縫隙等窄縫通道進入外界，給周邊人員及環(huán)境帶來影響[1]。

實驗和理論研究[2-4]表明，縫隙內(nèi)氣溶膠存在顯著的滯留，在一定情況下，氣溶膠甚至會堵塞通道，從而阻止放射性的進一步釋放?？p隙內(nèi)氣溶膠的滯留受多種機理影響，如重力沉降、布朗擴散、湍流擴散和湍流碰撞等，與縫隙內(nèi)的氣體流動狀態(tài)密切相關(guān)。研究人員在實驗和理論分析基礎(chǔ)上總結(jié)得到了一系列基于流動參數(shù)的縫隙內(nèi)氣溶膠穿透系數(shù)關(guān)系式[5-8]，在使用這些關(guān)系式時首先需要確定縫隙內(nèi)氣體泄漏率及流動特性，例如，在使用Fuchs[5]針對水平通道層流提出的氣溶膠穿透系數(shù)關(guān)系式時，需要已知通道內(nèi)流體的平均速度。對于毛細管類型的縫隙，工程上多采用基于層流理論推導得到的哈根-泊肅葉公式計算泄漏率[9]，但核電廠嚴重事故下安全殼內(nèi)、外兩側(cè)壓差較大，縫隙內(nèi)流速較高，馬赫數(shù)(Ma)較大，甚至會達到壅塞狀態(tài)，哈根-泊肅葉公式存在局限。Anderson等[10]分別針分子流、層流和壅塞流提出了毛細管的氣體泄漏率計算公式，與實驗數(shù)據(jù)的對比表明分子流和層流公式符合較好，壅塞流公式誤差較大。計算流體力學(CFD)方法以其對局部流動細節(jié)的模擬能力，被研究人員廣泛用于模擬壓力驅(qū)動的微小通道內(nèi)的氣體流動[11-13]，但工程應用較為不便。

本文以規(guī)則幾何結(jié)構(gòu)的毛細管作為縫隙代表，基于一維可壓縮絕熱流動方程，提出一種壓力驅(qū)動的毛細管內(nèi)氣體流動計算方法，該方法能夠模擬亞音速、臨界音速及壅塞流動，并利用該方法研究毛細管兩側(cè)壓差和內(nèi)徑對氣體流動的影響。

1 控制方程及物理模型

1.1 控制方程

壓力驅(qū)動的毛細管內(nèi)氣體流動如圖1所示，數(shù)字1～5表示不同位置。其中1為毛細管上游容器處；5為毛細管下游環(huán)境，上游容器壓力高于下游環(huán)境，驅(qū)動氣體在毛細管內(nèi)部形成流動；2為毛細管入口處，從1到2流動有入口形阻壓降；3和4為毛細管出口，從3到4流動有出口形阻壓降；毛細管內(nèi)部從2到3流動有沿程摩擦阻力壓降。

圖1 壓力驅(qū)動的毛細管氣體流動示意圖Fig.1 Sketch of pressure-driven gas flow in capillary

假定毛細管與外界絕熱，當流動的克努森數(shù)Kn小于0.001時，流體為連續(xù)介質(zhì)，可采用一維等截面通道內(nèi)有摩擦阻力的絕熱流動方程描述：

(1)

式中：p為壓力，Pa；f為范寧阻力系數(shù)；DH為水力直徑，m；ρ為流體密度，kg/m3；v為流體速度，m/s；x為一維方向上的距離，m。

馬赫數(shù)Ma如下：

(2)

式中：γ為氣體絕熱指數(shù)；R為氣體常數(shù)，N·m/(kg·K)。

結(jié)合式(2)和理想氣體狀態(tài)方程p=ρRT，式(1)可以變?yōu)椋?/p>

(3)

p和溫度T也可寫為Ma的函數(shù)，如下：

(4)

(5)

式(3)～(5)組成了毛細管內(nèi)一維絕熱流動的控制方程組。

在毛細管入口(x=0)處，有：

T=T0，p=pstg

(6)

式中：T0為滯止溫度，K；pstg為滯止壓力，Pa。

在毛細管出口(x=L)處，有：

p=pL，MaL≤1.0

(7)

式中：pL為出口處的壓力，Pa；MaL為出口處的馬赫數(shù)。

1.2 阻力系數(shù)關(guān)系式

對于毛細管等微小通道內(nèi)的流動阻力，氣體稀薄性減弱了通道的阻力系數(shù)，而氣體密度變化導致的流動加速則增大了阻力系數(shù)，這兩種效應的耦合作用使得阻力系數(shù)實驗結(jié)果偏離度較大。本文考慮毛細管內(nèi)氣體為連續(xù)流動的情形，未包括過渡流和分子流，因此計算層流阻力系數(shù)時選擇了經(jīng)典的層流阻力關(guān)系式進行計算；在計算湍流阻力系數(shù)時，文獻[14]開展的實驗數(shù)據(jù)對比顯示光滑管情況下Blasius關(guān)系式符合較好，而粗糙管時則Colebrook關(guān)系式符合較好。在層流向湍流轉(zhuǎn)捩方面，文獻[15]根據(jù)實驗結(jié)果認為微小通道內(nèi)的轉(zhuǎn)捩點(以雷諾數(shù)Re表征)存在提前的情況，某些情況下阻力系數(shù)轉(zhuǎn)捩點甚至可以提前到500～600；但文獻[16-18]則指出在各自的實驗中并未發(fā)現(xiàn)轉(zhuǎn)捩點提前，本文計算時考慮阻力系數(shù)轉(zhuǎn)捩點仍為2 100。這樣，本文選用的阻力系數(shù)關(guān)系式如下。

當雷諾數(shù)Re小于2 100時，有：

4f=64/Re

(8)

當雷諾數(shù)Re大于2 100時，有：

(9)

式中：D為毛細管直徑，m；ε為絕對粗糙度，m。

2 數(shù)值求解方法

上述控制方程中，阻力系數(shù)與雷諾數(shù)有關(guān)，進而與流體溫度有關(guān)，流體溫度又取決于流動狀態(tài)，因此無法直接求解。本文采用如下方法進行求解。

對式(3)、(4)從位置i到位置j進行積分，有：

(10)

式中，Kij為從位置i到位置j的阻力系數(shù)。

(11)

將式(10)、(11)應用于圖1的位置1～4，可得到6個方程，其中p1和p5已知，根據(jù)文獻[19]，K12可取0.5，K34可取1.0。

當流動為亞音速和臨界音速時，p4=p5，Ma4≤1.0。初始假定Ma4=0.001，則p2、p3、Ma1、Ma2、Ma3、K23未知，方程組封閉，可以求解。依次以較小的幅度增大Ma4直至1.0，從而得到已知p1和p5下的若干組可能的亞音速和臨界音速流動狀態(tài)參數(shù)(p2～p4，Ma1～Ma4和K23)。

當流動為壅塞流時，p4>p5，Ma4=1.0。初始假定p4為稍大于p5的數(shù)值，此時未知的物理量仍為p2、p3、Ma1、Ma2、Ma3、K23，可以求解。依次以較小的幅度增大p4直至K23=0.0，同樣可以得到已知p1和p5下的若干組可能的壅塞流動狀態(tài)參數(shù)(p2～p4，Ma1～Ma4和K23)。

針對每組可能的流動狀態(tài)參數(shù)，將Ma3-Ma2均分成N份(N由節(jié)點敏感性分析確定)，每一份即為微分方程(3)～(5)的有限差分單元。第i個差分單元有如下方程：

(12)

(13)

Δx采用下式計算：

(15)

式中，Δx為差分單元的長度，m。

毛細管的總長度L為：

(16)

毛細管內(nèi)質(zhì)量流量m為：

(17)

這樣，在已知毛細管直徑和絕對粗糙度后，針對每組可能的流動狀態(tài)參數(shù)，均可以計算得到相對應的毛細管長度和毛細管質(zhì)量流量。在此基礎(chǔ)上，采用毛細管實際長度進行插值，即可得到實際的質(zhì)量流量及流動狀態(tài)參數(shù)。采用FORTRAN語言編制計算程序?qū)崿F(xiàn)該數(shù)值求解方法，計算流程如圖2所示。

圖2 計算流程圖Fig.2 Computational flow diagram

3 數(shù)值計算方法的驗證

為測量微小通道的阻力系數(shù)，日本鹿兒島大學開展了毛細管流動特性實驗[20-21]，實驗涵蓋了低壓和高壓工況，具有較好的代表性，因此采用該實驗驗證數(shù)值計算方法。另外，考慮到實驗監(jiān)測點有限，無法得到毛細管內(nèi)所有位置的流動信息，因此，本文也采用三維CFD方法對實驗進行數(shù)值模擬，模擬結(jié)果亦用于數(shù)值方法的驗證。

3.1 實驗簡介

鹿兒島大學研究團隊采用的實驗裝置如圖3[20]所示。實驗采用的工質(zhì)為氮氣，溫度為300 K。實驗時將氮氣注入毛細管上游容器從而在毛細管內(nèi)部形成壓力驅(qū)動的流動，毛細管下游環(huán)境壓力為100 kPa。毛細管為光滑玻璃管，長度120 mm，直徑397 μm，相對粗糙度為0.002%[21]。實驗采用流量計測量通過毛細管的質(zhì)量流量，為得到毛細管內(nèi)部的壓力、溫度及馬赫數(shù)分布，分別在距離入口0.09、0.10、0.11 m處鉆小孔進行測量。

圖3 毛細管流動特性實驗裝置Fig.3 Experimental setup for flow characteristic through capillary

3.2 CFD模型

采用ANSYS Fluent 19.2對鹿兒島大學的毛細管流動特性實驗建模計算，CAD模型如圖4所示，模型將毛細管上游和下游均處理為圓柱形容器，容器直徑和長度遠大于毛細管直徑，從而消除了容器對毛細管內(nèi)部流動的影響。

圖4 毛細管流動特性實驗裝置的CAD模型Fig.4 CAD model for experimental setup for flow characteristic through capillary

采用ICEM CFD軟件對CAD模型進行網(wǎng)格劃分，生成結(jié)構(gòu)化網(wǎng)格。經(jīng)網(wǎng)格敏感性分析，最終選取了網(wǎng)格數(shù)量約為170萬的劃分方案。在前處理中分別為模型設置壓力入口與出口，如圖4所示，壓力入口處的壓力即為容器內(nèi)部的滯止壓力；壓力出口取外界環(huán)境壓力1個大氣壓。采用理想氣體狀態(tài)方程計算氮氣密度。湍流模型選取了SSTk-Ω模型，與標準的k-ε模型相比，該模型直接求解黏性層，不需要壁面函數(shù)，但需在邊界層內(nèi)布置一定密度的高質(zhì)量網(wǎng)格。經(jīng)計算，本文的邊界層網(wǎng)格Y+約為2，滿足模型適用要求。

3.3 對比分析

利用本文提出的一維數(shù)值計算方法和CFD模型對上游滯止壓力為200 kPa的實驗工況進行計算。圖5示出了一維方法、CFD模型計算得到的壓力、溫度及馬赫數(shù)沿管長分布與實驗的對比?？梢钥闯?，沿著管長方向，在摩擦阻力作用下，氣體壓力下降，密度減小，流動加速，溫度降低。圖6、7中毛細管不同位置的速度和溫度云圖更直觀地示出了流動加速和溫度降低的過程。由圖5可看出，在x=0.09，0.10，0.11 m位置處，兩種方法的計算結(jié)果與實驗較為接近，其中一維方法得到的壓力、溫度和馬赫數(shù)與實驗值的最大偏差分別為5.1%、1.4%和10.0%，CFD模型得到的壓力、溫度和馬赫數(shù)與實驗值的最大偏差分別為4.7%、0.6%和12.0%，一維方法與CFD方法精度相當。在管長方向上，一維方法與CFD模型變化趨勢一致，偏差在4.0%以內(nèi)。另外，該工況下實驗測得的質(zhì)量流量為1.966 96×10-5kg/s，一維方法和CFD模型對應的質(zhì)量流量分別為1.834 42×10-5kg/s和1.703 00×10-5kg/s，3種方法的質(zhì)量流量比較接近，但CFD模型與實驗數(shù)據(jù)的相對偏差較大。

圖5 滯止壓力200 kPa時的壓力、溫度與馬赫數(shù)沿管長的分布Fig.5 Distribution of pressure, temperature and Mach number along capillary length at stagnation pressure of 200 kPa

圖6 毛細管內(nèi)部不同位置處的速度云圖Fig.6 Velocity contour at different locations in capillary

圖7 毛細管內(nèi)部不同位置處的溫度云圖Fig.7 Temperature contour at different locations in capillary

采用一維方法對更多實驗工況進行計算，圖8示出不同滯止壓力下毛細管質(zhì)量流量、壓力及馬赫數(shù)與實驗的對比。可看出，一維方法很好預測了質(zhì)量流量及毛細管內(nèi)部壓力分布，其中質(zhì)量流量計算結(jié)果與實驗值在滯止壓力為150 kPa時偏差為11.16%，偏差較大，其余偏差均在3%以內(nèi)；壓力偏差均小于3%，馬赫數(shù)分布的計算結(jié)果與實驗值的偏差和質(zhì)量流量與壓力的偏差相比較大，但絕大部分偏差小于10%。產(chǎn)生這種偏差可能的原因在于，本文計算時考慮為絕熱條件，文獻未詳細說明毛細管試驗件的溫度控制措施及結(jié)果，溫度偏差導致了密度偏差，進而造成了馬赫數(shù)的偏差。

圖8 計算結(jié)果與實驗值的對比Fig.8 Contrast of calculated and experimental results

基于鹿兒島大學開展的毛細管流動實驗以及本文建立的CFD模型的驗證結(jié)果表明，本文開發(fā)的一維計算方法可在大的壓差范圍內(nèi)分析計算不同內(nèi)徑的毛細管內(nèi)壓力驅(qū)動的氣體流動。與CFD方法相比，一維方法計算效率很高，1個工況的計算時間約為5 s，而CFD模型在使用24核達到收斂解時所需的時間約6 h。

4 毛細管內(nèi)徑及兩側(cè)壓差對流動的影響

鹿兒島大學開展的毛細管流動特性實驗雖然涵蓋了低壓和高壓情形，但毛細管內(nèi)徑單一，本文采用經(jīng)驗證的一維計算方法開展了不同內(nèi)徑和壓差下的毛細管流動計算，以研究流動的變化規(guī)律。圖9示出了不同內(nèi)徑的毛細管內(nèi)壓力驅(qū)動的氣體質(zhì)量流量隨上游滯止壓力的變化，并結(jié)合圖8進行分析。

由圖8c可看出，在下游背壓不變的情況下，隨著上游滯止壓力的提高，兩側(cè)壓差增大，管長各處的馬赫數(shù)逐漸增大。但當滯止壓力超過某一數(shù)值時(圖8c約為400 kPa)，管長各處的馬赫數(shù)基本不再變化，表明流動處于壅塞狀態(tài)。雖然x=11 cm處距出口僅為1 cm，但受該1 cm段摩擦損失及出口形阻損失作用，馬赫數(shù)穩(wěn)定在0.54左右，不能達到1.0。另外，在流動達到壅塞狀態(tài)后，雖然馬赫數(shù)不再隨滯止壓力的提高而增大，但滯止壓力的提高增大了氣體密度，因而通過毛細管的質(zhì)量流量隨滯止壓力的提高持續(xù)增大，如圖8a和圖9所示。

進一步分析圖9，當毛細管內(nèi)徑小于100 μm時，隨壓差的增大氣體質(zhì)量流量近似呈二次方增長，而當內(nèi)徑大于100 μm時，壓差超過某一數(shù)值后質(zhì)量流量增長趨勢向線性化過渡。其原因在于，在本文計算的壓差范圍內(nèi)，小內(nèi)徑毛細管處于層流狀態(tài)且氣體壓縮性不顯著(以25 μm毛細管為例，雷諾數(shù)介于8～20，出口馬赫數(shù)介于0.004～0.06)，這樣沿程氣體溫度降幅較小，經(jīng)典的哈根-泊肅葉公式完全適用，因此泄漏流量隨壓差的增加近似呈二次方增長。

圖9 不同內(nèi)徑的毛細管內(nèi)質(zhì)量流量隨滯止壓力的變化Fig.9 Variation of mass flow rate with stagnation pressure in capillaries with different inner diameters

對毛細管內(nèi)徑較大的情況，隨壓差的增大，毛細管內(nèi)部流態(tài)會從小壓差的層流到大壓差的湍流轉(zhuǎn)變，馬赫數(shù)也逐漸增大，氣體壓縮性顯著，哈根-泊肅葉關(guān)系式不再適用。

注意到，對于150 μm和200 μm的情形，在泄漏流量由二次方增長變?yōu)榫€性增長時，均出現(xiàn)了流量平臺，平臺范圍內(nèi)隨壓差的增大流量只有小幅增大或基本不變。以150 μm為例，質(zhì)量流量平臺出現(xiàn)在滯止壓力350～450 kPa。圖10示出沿程摩擦阻力系數(shù)K23隨滯止壓力的變化，由圖10可見，隨著滯止壓力的提高，阻力系數(shù)首先下降，在350～450 kPa范圍內(nèi)，阻力系數(shù)增大，之后又呈下降趨勢。根據(jù)這種阻力變化趨勢，可推斷隨著滯止壓力的提高，泄漏流量首先增大，在平臺區(qū)內(nèi)，滯止壓力提高導致氣體流速增大，流速增大又造成阻力系數(shù)增大，因此泄漏流量只有小幅增長或基本不變，當跨過平臺區(qū)后，泄漏流量隨滯止壓力的提高而逐漸增大。進一步分析如圖11所示的平臺區(qū)雷諾數(shù)，可看到隨滯止壓力的增大，雷諾數(shù)增大，但均位于層流到湍流的過渡區(qū)。根據(jù)經(jīng)典的阻力系數(shù)莫迪圖[19]，在過渡區(qū)阻力系數(shù)隨雷諾數(shù)的增大而增大，證明了前述分析的合理性。

圖10 阻力系數(shù)隨滯止壓力的變化Fig.10 Variation of friction coefficient with stagnation pressure

圖11 不同滯止壓力下雷諾數(shù)沿管長的分布Fig.11 Distribution of Reynolds number along pipe length under different stagnation pressures

比較圖9c和d，可發(fā)現(xiàn)，200 μm情形下流量平臺要早于150 μm，原因在于毛細管內(nèi)徑越大，相同滯止壓力下的雷諾數(shù)越大，因此層流到湍流的過渡區(qū)越早出現(xiàn)。當毛細管內(nèi)徑增大到一定程度時，質(zhì)量流量平臺幾乎消失，如圖8a所示。

5 結(jié)論

1) 針對壓力驅(qū)動的毛細管內(nèi)氣體流動，本文基于一維可壓縮絕熱方程開發(fā)了一種數(shù)值計算方法，能夠計算壓力驅(qū)動的毛細管內(nèi)亞音速、臨界音速與壅塞流動。經(jīng)過與鹿兒島大學開展的毛細管流動特性實驗和CFD數(shù)值模擬結(jié)果對比，流動趨勢一致，偏差較小，驗證了一維數(shù)值計算方法的準確性與可靠性。與CFD數(shù)值模擬方法相比，該方法大幅提高了計算效率。

2) 采用本文開發(fā)的數(shù)值方法分析了毛細管兩側(cè)壓差與內(nèi)徑對泄漏特性的影響，結(jié)果表明小內(nèi)徑毛細管的泄漏流量隨壓差的增大呈二次方增長，而當內(nèi)徑大于100 μm時，壓差超過某一數(shù)值后質(zhì)量流量增長趨勢向線性化過渡。當壓差使得流動處于層流和湍流過渡區(qū)時，受阻力系數(shù)變化趨勢影響，泄漏流量會出現(xiàn)一個平臺，內(nèi)徑越大，則平臺出現(xiàn)越早。

3) 本文數(shù)值計算方法的開發(fā)為基于氣溶膠穿透系數(shù)關(guān)系式的氣溶膠泄漏與滯留分析奠定了技術(shù)基礎(chǔ)。針對嚴重事故后安全殼內(nèi)高溫、高濕的熱工環(huán)境，后續(xù)將拓展該方法至壓力驅(qū)動的毛細管內(nèi)多組分氣體流動的計算分析。