壓力容器內超臨界二氧化碳噴放特性分析及模型驗證

明 楊,方華偉,劉 凱,易經緯,劉玖松,趙富龍,*,田瑞峰

(1.哈爾濱工程大學 核科學與技術學院,黑龍江 哈爾濱 150001;2.哈爾濱工程大學 黑龍江省核動力裝置性能與設備重點實驗室,黑龍江 哈爾濱 150001;3.中國核動力研究設計院,四川 成都 610213)

超臨界二氧化碳能量轉換系統具有循環效率高、系統緊湊、靈活性高、應用場景廣泛等優點。近年來,與第4代反應堆結合的超臨界二氧化碳布雷頓循環已成為能源行業的研究熱點方向之一[1-2]。然而,超臨界二氧化碳的獨特物性不僅對循環設備的設計與制造帶來了特殊要求,同時也為系統出現泄漏或破口事故下的安全分析帶來了一系列的未知與挑戰。

超臨界二氧化碳流體噴放過程隨著上游滯止狀態的不同,可能出現單相、兩相,甚至多相臨界流,其噴放進程較為復雜[3]。此外,噴放過程會導致循環壓力的快速降低,二氧化碳的密度、比熱容、導熱系數等參數將出現大幅變化,這將對系統的穩定安全運行帶來較大的威脅。因此,研究超臨界二氧化碳的破口噴放特性,明確循環內冷卻劑的流失速率,了解壓力容器內的工質狀態極為重要。

在破口事故的預防與處理方面,諸多學者進行了大量的研究。早在1963年,Faltetti與Moulton[4]開展了蒸汽與水兩相混合物的臨界流實驗研究,并建立了兩相均勻平衡模型。Moody[5]分別從動量守恒與能量守恒方程出發,推導出臨界流動的理論模型。之后,Henry與Fauske[6]同樣基于能量守恒方程,獲得了更能更適用于計算短管噴放的Henry-Fauske模型。在此之后,眾多學者針對以亞臨界水工質的臨界流模型進行了相關的實驗與模擬研究,并對已有模型進行修正與優化。

近年來,隨著超臨界二氧化能量轉換系統的興起,超臨界流體的破口噴放現象與臨界流特性研究逐步成為關注重點之一。在實驗方面,德國漢堡科技大學的Gebbeken等[7]進行了超臨界二氧化碳噴放的實驗研究,觀察到了容器內超臨界二氧化碳的閃蒸現象以及噴放過程中壓力曲線的劇烈彎曲,他們認為發生閃蒸的壓力主要取決于初始條件。Hébrard等、Ahmad等與Guo等[8-10]分別進行了中等規模尺寸容器與大型管道內超臨界二氧化碳噴放過程的相關研究,在實驗進程中均觀測到了飽和狀態下二氧化碳的長時間穩定噴放。Fan等[11]研究了超臨界二氧化碳的臨界流動行為,結果表明臨界流量隨上游溫度的增大而減小,隨上游壓力的增大而增大。Zhou等[12]在容積為50 L的壓力容器上進行了超臨界二氧化碳噴放實驗研究,結果表明噴放持續的時間主要與破口直徑和初始壓力有關。在模擬方面,Wang等[13]建立了超臨界二氧化碳壓力容器的噴放模型,并采用均相模型和相分離模型分析了容器內的相分布。Min等[14]基于Henry-Fauske模型計算了超臨界二氧化碳通過噴嘴的臨界流動,并與實驗裝置得到的實驗數據進行了對比。章靜等[15]建立了簡單容器噴放的數學物理模型,并將計算結果與實驗數據進行了比較,驗證了模型的合理性。李偉卿等[16]基于超臨界二氧化碳臨界流實驗的測量結果,驗證了臨界流熱平衡通用模型的準確性。

總體而言,超臨界二氧化碳噴放過程的機理尚未明確,利用現有的實驗數據得到的經驗關系式適用的幾何參數與工況范圍較為有限,噴放過程的特性分析與流量預測等方面的工作亟待開展。因此本文基于Moody與Henry-Fauske提出的單組分混合物兩相臨界流分析方法,采用Modelica語言開發一套瞬態分析程序(SCRSAP-LOCA),針對超臨界二氧化碳噴放過程中的臨界流現象進行模擬,與現有實驗數據進行對比驗證,并分析容器內不同初始條件對噴放過程造成的影響。

1 計算模型

1.1 簡單容器模型

本研究借鑒壓水堆與超臨界水堆LOCA事故的研究方法,采取簡單容器模型研究超臨界二氧化碳熱力循環中壓力容器發生破口時的噴放進程。簡單容器模型基于等熵假設,認為噴放過程中流體處于熱力平衡狀態,忽略容器與外界的熱量交換、壓力梯度和流動阻力,并采用可壓縮流體模型[13 , 15]。簡化后的瞬態過程由質量守恒方程、能量守恒方程及狀態方程共同描述[13]:

(1)

(2)

h=h(ρ,s)

(3)

式中:V為容器的體積,m3;p為容器內的壓力,Pa;Gout為破口處的流量,kg/s;ρ、s、T、h分別為工質的密度、熵、溫度與焓;Q為總換熱量,J。根據質量守恒與能量守恒方程可以計算得到容器內工質的密度與熵,并由狀態方程得到焓、溫度等參數。

1.2 容器內工質的相模型

根據已有的實驗研究,噴放初始階段壓力容器內的壓力較高,二氧化碳處于超臨界狀態或液相。而伴隨噴放的進行,容器內壓力將快速降至二氧化碳的飽和點,氣化的二氧化碳將因密度差遷移至容器頂部,容器內可以觀測到明顯的氣液界面[7,17]。

如圖1所示,在初始的單相噴放階段,容器內的二氧化碳并未產生相變,即破口處與容器內部的二氧化碳均處于超臨界態或單相液態,可采取均勻流模型:

圖1 簡單容器模型及噴放進程Fig.1 Simple vessel model and ejection process

hout=h0

(4)

式中,h0為容器內工質的平均焓,J/kg。

而隨著噴放過程的進行,氣液界面出現后的相分離效應不可忽略,破口處的工質性質應使用分相流模型表征,當氣液界面高于破口高度時,認為噴放出的工質為當前壓力下的飽和液體,當氣液界面低于破口高度時,則認為是飽和氣體[13],具體的方程如下所示:

hout=hf(p)Z>Ze

hout=hg(p)Z

(5)

式中:hf為飽和液焓;hg為飽和氣焓;α為空泡份額;Z為自由液面高度,m;Ze為破口高度;Zall為容器的總高度。

1.3 破口噴放流量模型

1) 非臨界流動模型

在噴放初期,對于二氧化碳超臨界區的非臨界噴放,由于流體的壓縮性較差,該過程可視為不可壓縮流體的等熵流動,一般選用Bernoulli方程計算破口處的質量流量[11]:

p0>pc,T0>Tc,x0≤0

(6)

式中:p0為上游壓力;pb為背壓;pc與Tc分別為二氧化碳的臨界壓力與溫度;x0為含氣率;Cd為流量系數,其定義為真實流量與理想流量的比值,取值與破口的幾何形狀及長徑比L/d相關,在建模過程中需根據計算結果進行修正。

現有的超臨界二氧化碳噴放實驗中測定的Cd的取值范圍為0.61～0.9[11,16,18]。本文在進行多組工況的計算與對比后,最終發現Cd選取保守值0.61時,獲得的非臨界流動流量計算結果與Wang等[13]的實驗結果最為接近。

2) 單相過冷臨界流動模型

當容器內工質的壓力高于臨界壓力,但溫度低于臨界溫度時,容器內工質由超臨界區進入密相區,破口處為過冷臨界噴放狀態。本文采取Henry-Fauske模型計算該階段的破口流量[19-20],臨界流量表達式為:

p0>pt,x0>0

(7)

式中:x為含氣率;v為比容,m3/kg ;cp為比定壓熱容,J/(kg·K);下標0代表容器內的滯止參數,下標g與f分別代表氣相與液相參數;N為表征破口處發生部分相變的參數,當含氣率小于0.14時,N=x/0.14,而含氣率大于0.14時,N=1;γ為等熵膨脹系數,γ=cp/cv;n的定義為熱力平衡多變指數,其表達式為:

(8)

臨界壓比η的表達式為:

(9)

式中,β、α0與αt的表達式可在文獻[19-20]中查詢。而當容器內的二氧化碳處于過冷狀態時,臨界流量表達式可進一步近似為:

p0>pt,x0≤0

(10)

在容器內滯止參數p0與x0已知的情況下,可將臨界流量計算公式與臨界壓比進行迭代求解,具體過程如下:首先假設破口處壓力等于上游滯止壓力,即臨界壓比η=1。同時調用pt下對應的物性參數,代入式(9),可求得新的臨界壓比。若兩次壓力比數值不同,則將下游壓力減掉一個小量(0.001p0),并求得新的臨界壓比,直到兩次壓力比數值相同。此時得到的壓力比為臨界壓力比,求得的下游臨界壓力pt=ηp0,進而計算得到臨界流量。

3) 兩相臨界流動模型

隨著容器內工質含氣率的增長,破口處逐步達到穩定的兩相臨界流動狀態,現有較為常用的兩相臨界流動模型為Henry-Fauske模型與Moody模型,二者均基于等熵假設。

本文在建模時同時考慮了上述兩種模型,以驗證其對二氧化碳工質的適用性。其中,Henry-Fauske模型的主要方程如式(7)～(10)所述,Moody模型的原理如方程(11)所述,該模型由能量守恒出發,認為氣液兩相處于熱力平衡狀態并擁有各自的流動速度,其主要方程[21]為:

p0>pt,x0>0

(11)

式中:hfg=hg-hf;sfg=sg-sf;h0為容器內的滯止焓;S為滑速比,S=(vg/vf)1/3。上式意味著,當容器內的滯止參數已知時,質量流量Gc可由滑速比S和壓力p計算。此外滑速比S表達式中的液體比容vf和氣體比容vg由壓力p決定。因此Moody模型最終的計算過程如下。

首先給定上游的工況條件,確定容器內兩相混合物的滯止焓與滯止壓力。之后給下游出口端一個初始的壓力,根據動量方程原理,破口處壓力應小于滯止壓力,令pt略低于p0并調用此時下游壓力下對應的熱力參數,并代入式(11)求得對應的質量流量。之后逐步降低下游壓力,求出一系列的質量流量。選取所有pt對應的流量并查詢極大值,即為兩相臨界流流量,其對應的破口壓力為臨界壓力。

4) 單相氣臨界流動模型

隨著容器內工質含氣率的進一步增長,當自由液面高度在破口高度以下時,噴放出的工質將由兩相混合物轉變為單相二氧化碳氣體。認為此時的二氧化碳氣體為理想氣體,臨界壓力與臨界流量由下式[15]計算:

(12)

(13)

2 瞬態代碼的實現

本文的超臨界二氧化碳噴放瞬態分析程序(SCRSAP-LOCA)采用Modelica語言編譯,同時基于模塊化思想,根據計算功能分為簡單容器模塊、相判斷模塊、破口噴放計算模塊、物性查詢模塊等子模塊。程序主體算法與計算流程如圖2所示,采取有限差分方法求常微分方程,并將單相過冷臨界流模型與兩相臨界流模型中涉及到的迭代求解編譯成函數調用的方式,時間步長設置在0.01～0.1 s之間。物性數據來源于NIST RefProp軟件[22],并對二氧化碳在臨界點附近的物性數據采取局部加密方法,最終實現的物性查詢最大相對誤差在5%以內。

圖2 程序主體算法與計算流程Fig.2 Program main algorithm and calculation flow

3 結果與討論

3.1 計算結果與實驗數據對比

首先將瞬態計算結果與Wang等[13]的實驗數據進行了比對,實驗裝置主要由壓力容器、加熱帶、保溫層、管道、閥門、噴嘴、采集系統等構成。壓力容器的總體積為0.05 m3,噴嘴長度為5 mm,直徑為1.01 mm,噴嘴高度位于壓力容器總高度的1/2處。實驗的初始工況列于表1[13],程序設置的初始參數與實驗保持一致。

表1 實驗壓力容器內的初始參數Table 1 Initial parameter in experimental pressure vessel

不同初始條件下壓力容器破口流量與壓力瞬變如圖3所示,噴放過程可大致分為3個階段,破口流量與容器壓力呈階梯式下降。

圖3 不同初始條件下流量與壓力計算結果與實驗數據[13]Fig.3 Calculation results and experimental data of mass flow rate and pressure with different initial conditions[13]

第1階段為高質量流量快速泄壓階段,流量與壓力曲線的斜率最大。破口流量的范圍為0.045～0.06 kg/s,容器內壓力迅速降至6.5～7 MPa,同時質量流量與工質密度迅速減小。初始階段容器中二氧化碳為超臨界態或液態,其密度較高且可壓縮性差,破口處于非臨界流動或單相過冷臨界流動狀態。第1階段持續時長為100 s左右,階段結束的標志為壓力與流量曲線斜率突然降低后出現了拐點。

第2階段為中質量流量穩定泄壓階段,流量與壓力的曲線斜率最小。破口流量的范圍在0.025～0.028 kg/s之間,呈階梯平臺狀緩慢下降,而容器內壓力接近線性下降。此時容器內二氧化碳為兩相態,破口處于兩相臨界流動狀態,容器內工質含氣率逐漸增大,但此時自由液面高度仍高于破口高度。第2階段持續時長為350 s左右,階段結束的標志為壓力與流量曲線出現第2個拐點,并且斜率出現一定程度的增加。

第3階段為低質量流量勻速泄壓階段,流量與壓力曲線均呈凹狀,斜率先增加后逐步減小。破口流量的范圍為0.005～0.025 kg/s,容器內壓力于1 600 s左右降至3 MPa左右。此時容器內二氧化碳仍為兩相態,但隨著容器內工質含氣率的增大,自由液面高度低于破口高度,因此破口處工質變為單相氣體,質量流量明顯降低。

圖4為初始條件為10.0 MPa、308 K時容器內二氧化碳的P-D(壓力-密度)與P-T(壓力-溫度)曲線。由圖4可知,噴放第1階段開始時破口處為單相噴放,容器內二氧化碳處于超臨界狀態,其密度為719.51 kg/m3。當噴放第2階段開始時,容器內二氧化碳的密度降低至609.73 kg/m3,P-D曲線與二氧化碳飽和液線相交,容器內發生相變并出現氣液界面。而在噴放第3階段內,P-D曲線逐步接近飽和氣體密度曲線,P-T曲線也愈發向氣相區靠近。當自由液面降至破口以下時,噴放出二氧化碳飽和氣體,因此容器內工質密度的降低速率逐步降低。

圖4 初始條件為10 MPa、308 K的P-D曲線與P-T曲線Fig.4 P-D and P-T curves with initial condition of 10 MPa and 308 K

總體而言,結果表明程序計算的結果與實驗數據吻合較好,流量與壓力的計算誤差均不超過23%。在兩相臨界流動模型方面,Moody模型計算的破口質量流量與實驗最為接近,而Henry-Fauske模型計算得到的流量相對較高,這導致流量曲線第2拐點的較早出現。程序的壓力計算結果與實驗數據在噴放的第1、2階段較為接近,但在第3階段略低于實驗數據,這是由于程序忽略了容器內的壓力梯度,計算所得的平均壓力可能與實驗中采集點得到的局部壓力存在一定誤差。

3.2 初始壓力對噴放過程的影響

圖5示出了溫度為297 K,壓力分別為10、15、20 MPa的不同初始條件下的主要參數計算結果。如圖5a與圖5c所示,在初始溫度相同時,噴放第1階段的初始流量與初始壓力呈正相關,這是由于更高的初始壓力使流體具有高的初始密度。同時,隨著初始壓力的增加,容器內壓力與流體密度降低的速率也逐步增加,圖5c中P-D曲線與飽和曲線的交點即為噴放第2階段的起始點。如圖5a與圖5b所示,隨著初始壓力的增大,噴放第2階段開始時容器內的壓力變低,但兩相臨界噴放的質量流量幾乎不隨初始壓力的變化而變化。此外,如圖5d所示,更高的初始壓力使自由液面高度下降的速率更慢,這導致噴放更晚進入第3階段,因此兩相臨界噴放持續時間將更長。

圖5 不同壓力初始條件下主要參數計算結果Fig.5 Main parameter calculation results under different initial pressure conditions

3.3 初始溫度對噴放過程的影響

圖6示出了壓力為15 MPa,溫度分別為280、297、330、350 K的主要參數計算結果,可看到在不同的初始溫度條件下,容器內二氧化碳的噴放行為將出現較大的差別。

圖6 不同溫度初始條件下主要參數計算結果Fig.6 Main parameter calculation results under different initial temperature conditions

當初始溫度低于臨界溫度(304.13 K)時,容器內的二氧化碳為液態,更低的初始溫度使流體具有更高的初始密度與初始噴放流量,對應噴放第2階段的壓力更低,兩相臨界噴放的持續時間更長。而當初始溫度高于臨界溫度時,容器內二氧化碳處于超臨界狀態。當初始溫度為330 K時,流體初始密度為635.51 kg/m3,噴放過程中自由液面下降的速度顯著加快,這使得兩相臨界噴放的時間大幅縮短。值得注意的是,由圖6c與6d所示,當初始溫度進一步升高至350 K時,流體密度為449.20 kg/m3,超臨界二氧化碳的物性更接近于氣體。由圖6c可知,P-D曲線由二氧化碳的超臨界區進入到氣相區,隨后與飽和線相交,這表明容器內的二氧化碳發生了相變。液化的二氧化碳在密度差的作用下沉積至容器底部,導致了圖6d中的相對液面高度曲線在200 s左右升高,這與Tian等[17]的實驗測量結果相符。由于破口處將持續維持較長時間的單相氣體臨界噴放狀態,導致流量與壓力曲線變得平滑,拐點幾乎消失。

4 結論

因此本文基于Modelica語言,開發了一套針對超臨界二氧化碳壓力容器破口噴放的瞬態分析程序(SCRSAP-LOCA),并對噴放過程中不同的流動現象進行了建模與分析,選取了多組工況的計算結果與現有實驗數據進行了對比驗證,分析了容器內不同初始條件對噴放過程造成的影響,得到的主要結論如下。

1) 程序的計算結果與實驗數據吻合較好,流量與壓力的計算誤差均不超過23%。二氧化碳壓力容器的噴放可大致分為3個階段,破口流量與容器壓力曲線呈階梯式下降。

2) 程序采用的兩相臨界流動模型中,Moody模型計算的破口質量流量與實驗值最為接近,而Henry-Fauske模型計算得到的流量與實驗值相比高15%～25%左右。

3) 初始溫度為297 K但初始壓力增加時,第1階段的初始噴放流量更大,第2階段開始時容器內的壓力更低,兩相臨界噴放持續時間更長,但破口處的兩相臨界噴放流量幾乎不隨初始壓力的變化而變化。

4) 不同的初始溫度使容器內二氧化碳的噴放行為具有很大的差別。當壓力容器內的二氧化碳為超臨界態時,初始溫度的升高會使流體密度大幅降低,自由液面下降的速度顯著加快,兩相臨界噴放的時間大幅縮短。當初始溫度進一步升高,超臨界二氧化碳愈發趨近于類氣區域時,流量與壓力曲線將變得平滑,此外噴放過程中容器內會發生相變,液化的二氧化碳會在密度差的作用下沉積至容器底部。

本文研究結果可為壓力容器內超臨界二氧化碳噴放特性分析提供參考,并對超臨界二氧化碳動力轉換系統的泄漏事故進程研究與安全分析評價提供支撐。