壓力驅(qū)動的毛細(xì)管內(nèi)氣體流動特性數(shù)值研究

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

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

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

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

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

1 控制方程及物理模型

1.1 控制方程

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

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

假定毛細(xì)管與外界絕熱，當(dāng)流動的克努森數(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)組成了毛細(xì)管內(nèi)一維絕熱流動的控制方程組。

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

T=T0，p=pstg

(6)

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

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

p=pL，MaL≤1.0

(7)

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

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

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

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

4f=64/Re

(8)

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

(9)

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

2 數(shù)值求解方法

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

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

(10)

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

(11)

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

當(dāng)流動為亞音速和臨界音速時，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)。

當(dāng)流動為壅塞流時，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é)點(diǎn)敏感性分析確定)，每一份即為微分方程(3)～(5)的有限差分單元。第i個差分單元有如下方程：

(12)

(13)

Δx采用下式計(jì)算：

(15)

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

毛細(xì)管的總長度L為：

(16)

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

(17)

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

圖2 計(jì)算流程圖Fig.2 Computational flow diagram

3 數(shù)值計(jì)算方法的驗(yàn)證

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

3.1 實(shí)驗(yàn)簡介

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

圖3 毛細(xì)管流動特性實(shí)驗(yàn)裝置Fig.3 Experimental setup for flow characteristic through capillary

3.2 CFD模型

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

圖4 毛細(xì)管流動特性實(shí)驗(yàn)裝置的CAD模型Fig.4 CAD model for experimental setup for flow characteristic through capillary

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

3.3 對比分析

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

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

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

圖8 計(jì)算結(jié)果與實(shí)驗(yàn)值的對比Fig.8 Contrast of calculated and experimental results

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

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

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

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

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

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

對毛細(xì)管內(nèi)徑較大的情況，隨壓差的增大，毛細(xì)管內(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)，滯止壓力提高導(dǎo)致氣體流速增大，流速增大又造成阻力系數(shù)增大，因此泄漏流量只有小幅增長或基本不變，當(dāng)跨過平臺區(qū)后，泄漏流量隨滯止壓力的提高而逐漸增大。進(jìn)一步分析如圖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，原因在于毛細(xì)管內(nèi)徑越大，相同滯止壓力下的雷諾數(shù)越大，因此層流到湍流的過渡區(qū)越早出現(xiàn)。當(dāng)毛細(xì)管內(nèi)徑增大到一定程度時，質(zhì)量流量平臺幾乎消失，如圖8a所示。

5 結(jié)論

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

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

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