壓力驅動的毛細管內氣體流動特性數值研究

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

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

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

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

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

1 控制方程及物理模型

1.1 控制方程

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

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

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

(1)

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

馬赫數Ma如下：

(2)

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

結合式(2)和理想氣體狀態方程p=ρRT，式(1)可以變為：

(3)

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

(4)

(5)

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

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

T=T0，p=pstg

(6)

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

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

p=pL，MaL≤1.0

(7)

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

1.2 阻力系數關系式

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

當雷諾數Re小于2 100時，有：

4f=64/Re

(8)

當雷諾數Re大于2 100時，有：

(9)

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

2 數值求解方法

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

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

(10)

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

(11)

將式(10)、(11)應用于圖1的位置1～4，可得到6個方程，其中p1和p5已知，根據文獻[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下的若干組可能的亞音速和臨界音速流動狀態參數(p2～p4，Ma1～Ma4和K23)。

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

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

(12)

(13)

Δx采用下式計算：

(15)

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

毛細管的總長度L為：

(16)

毛細管內質量流量m為：

(17)

這樣，在已知毛細管直徑和絕對粗糙度后，針對每組可能的流動狀態參數，均可以計算得到相對應的毛細管長度和毛細管質量流量。在此基礎上，采用毛細管實際長度進行插值，即可得到實際的質量流量及流動狀態參數。采用FORTRAN語言編制計算程序實現該數值求解方法，計算流程如圖2所示。

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

3 數值計算方法的驗證

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

3.1 實驗簡介

鹿兒島大學研究團隊采用的實驗裝置如圖3[20]所示。實驗采用的工質為氮氣，溫度為300 K。實驗時將氮氣注入毛細管上游容器從而在毛細管內部形成壓力驅動的流動，毛細管下游環境壓力為100 kPa。毛細管為光滑玻璃管，長度120 mm，直徑397 μm，相對粗糙度為0.002%[21]。實驗采用流量計測量通過毛細管的質量流量，為得到毛細管內部的壓力、溫度及馬赫數分布，分別在距離入口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所示，模型將毛細管上游和下游均處理為圓柱形容器，容器直徑和長度遠大于毛細管直徑，從而消除了容器對毛細管內部流動的影響。

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

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

3.3 對比分析

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

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

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

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

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

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

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

4 毛細管內徑及兩側壓差對流動的影響

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

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

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

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

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

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

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

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

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

5 結論

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

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

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