龔政宇 辜峙钘,2 潘麒文 張牧昊 戴嘉寧

1（成都理工大學 核技術與自動化工程學院 成都 610059）

2（中國核動力研究設計院 核反應堆系統(tǒng)設計技術重點實驗室 成都 610213）

壓力波傳播特性是反應堆事故安全分析的關注點之一，反應堆系統(tǒng)內壓力波傳播將引發(fā)水力學載荷問題，獲得壓力波的傳播數據是開展結構載荷分析、評價的必要條件［1］。核能系統(tǒng)反應堆內的壓力波傳播現(xiàn)象一般由瞬態(tài)事故引發(fā)，例如，蒸汽輪機/主蒸汽管閥突開、突閉，或諸如冷卻劑喪失（Loss Of Coolant Accident，LOCA），蒸汽發(fā)生器傳熱管破裂（Steam Generator Tube Rupture accident，SGTR）［2］等事故。大部分壓力波傳播數值模擬采用可壓縮歐拉流體動力學模型［3］。而對于兩相流動、傳熱、傳質條件下的壓力波傳播問題，通常選用兩流體六方程模型，被系統(tǒng)分析程序廣泛使用，如RELAP5，TRACE，WAHA［4-5］。Sokolowski等［6］針對典型系統(tǒng)分析程序RELAP5 模擬壓力波傳播的能力開展了分析與評價，包括單相氣體介質、單相蒸汽介質及汽水兩相介質條件下壓力波傳播問題。Barna和Serre使用蒸汽冷凝引發(fā)的水錘效應及壓力波傳播案例，對WAHA3 和CATHARE 程序也進行了分析與評價［7-8］。另外，兩流體六方程模型一致采用平衡態(tài)壓力模型，即氣液兩相的動量守恒方程共用同一壓力，針對該問題，相關研究機構已開始研發(fā)新一代反應堆系統(tǒng)分析程序，如RELAP7，其核心便是雙壓力七方程模型。Delchini等［9］基于單相、兩相激波管基準例題，分析、評價了RELAP7程序捕捉壓力波傳播行為的能力，結果表明：在壓力波傳播模擬方面，該程序的確優(yōu)于RELAP5程序。

壓力波傳播行為通常存在多維特性，幾何結構越復雜，該特征越顯著，例如SGTR事故下壓力波在傳熱管束范圍內的復雜傳播、反射及疊加行為，而前述系統(tǒng)分析程序通常屬于一維程序，無法處理壓力波的多維傳播問題，因此開展壓力波多維傳播模擬技術、方法研究具有重要意義。針對該問題，本文就汽、水兩相流中存在的壓力波二維軸對稱傳播行為，開展了數值模型與算法研究，控制方程采用二維軸對稱圓柱坐標系，考慮了汽、水兩相非平衡態(tài)相變傳熱、傳質行為，控制方程求解使用了時間步分離四步算法，開發(fā)了程序代碼TPFPWPC-2D （Two Dimensional Two Phase Flow Pressure Wave Propagation Code）。采用典型汽、水兩相激波管基準例題開展了程序驗證，同時為考查該程序模擬壓力波二維傳播特性的能力，對圓柱體空間內壓力波的軸對稱二維傳播特性進行了數值模擬與分析研究。

1 模型和算法

1.1 控制方程

本文涉及的控制方程主要包括質量守恒方程、動量守恒方程、能量守恒方程，且全部控制方程均基于二維軸對稱圓柱體坐標系導出。

1） 質量守恒方程

就質量組分而言，考慮三種不同質量組分，即液態(tài)水、水蒸氣、不可凝氣體，其質量守恒方程如式（1）～（3）所示。其中：ΓC、ΓV為相變質量源項；αL為液態(tài)水體積份額；αG為蒸汽混合物（水蒸氣和不可凝氣體）體積份額；ρˉL、ρˉGh2o、ρˉGgas分別是液態(tài)水、水蒸氣、不可凝氣體宏觀密度；ρL、ρGh2o、ρGgas為相應理論密度；uL、uG分別是液態(tài)水、蒸汽混合物徑向速度；νL、νG為相應軸向速度?？梢钥闯鏊魵夂筒豢赡Y性氣體共享同一速度場。

2） 動量守恒方程

本文考慮了兩類速度場組分，即液態(tài)水速度場和蒸汽混合物速度場，其徑向動量守恒方程由式（4）、（5）描述，軸向動量守恒方程如式（6）、（7）所示。其中：p為網格壓強；FDR，G，R、FDR，L，R表示氣、液兩相徑向曳力；FDR，G，Z、FDR，L，Z表示氣、液兩相軸向拽曳力；FVM，G，R、FVM，L，R為 氣、液 兩 相 徑 向 虛 擬 質 量 力；FVM，G，Z、FVM，L，Z表示氣、液兩相軸向虛擬質量力；ρˉG為蒸汽混合物宏觀密度，表示為ρˉG=ρˉGh2o+ρˉGgas。

3） 能量守恒方程

考慮兩種能量組分，即液態(tài)水、蒸汽混合物，其能量守恒方程如式（8）、（9）所示。其中：QΓ，G、QΓ，L為非平衡態(tài)相變傳熱；QG、QL為平衡態(tài)溫差傳熱；a為接觸面積；h為換熱系數；ΔT為溫差。由于本文模擬對象涉及強傳熱、傳質，而且對流效應也非常顯著，因此忽略了網格間熱傳導的貢獻，也忽略了相間相互作用力做功的貢獻。

1.2 數值算法

本文模擬對象具有強傳熱、傳質，源項變化劇烈，壓力波傳播急劇等典型特征，同時其“速度-壓強”耦合尤為緊密，基于完全統(tǒng)一、全隱式求解算法對本模擬對象而言不太現(xiàn)實，主要是收斂性難以保障。針對此問題，本文采用了時間步分割四步算法技術，以提高其數值計算收斂性與穩(wěn)定性。

1） 第一步：網格內源相計算

首先，聯(lián)立迭代求解水蒸氣質量守恒方程，液態(tài)水、蒸汽混合物能量守恒方程，僅考慮傳熱、傳質源項，其核心是基于汽、液兩相狀態(tài)方程（Equation Of State，EOS）的非平衡態(tài)傳熱、傳質迭代計算，迭代收斂后，進而更新各組分宏觀密度、比內能、溫度及傳質速率。然后，根據傳質速率單獨求解各組分動量守恒方程，僅考慮其相變源項，由于當前“速度-壓強”并不匹配，因此相變源項中速度采用顯式格式，最終更新速度場。上述求解過程整體邏輯框架體系如圖1所示，其中上標1 表示第一步值，上標n表示前時間步值。

圖1 第一步：網格內源項計算數值算法框架體系Fig.1 Step 1: numerical algorithm framework for cell-inside source term calculation

2） 第二步：網格間對流計算

首先，單獨求解全部質量守恒方程，只考慮對流貢獻，且采用顯式格式；然后，單獨求解全部能量守恒方程，只考慮對流貢獻，且也采用顯式格式；最后，聯(lián)立迭代求解動量守恒方程，僅忽略相變源項，其中壓力項、對流項均采用顯式格式，拽曳力、虛擬質量力均采用隱式格式。由于在第2步中，“速度-壓強”依然不匹配，因此其對流項、壓強均采用顯式格式。上述求解過程整體邏輯框架體系如圖2 所示，其中上標2 表示第二步值，上標0 表示前時間步值。

圖2 第二步：網格間對流計算數值算法框架體系Fig.2 Step 2: numerical algorithm framework for intercell convection calculation

3） 第三步：壓力迭代

首先，確定殘差方程，本文共考慮4個，第1個是壓力殘差（EOS 壓力和網格壓力之差），第2、3 分別是液態(tài)水、蒸汽混合物質量守恒方程殘差，其對流項采用半隱格式，第4 個是蒸汽混合物能量守恒方程殘差，其比內能對流項采用顯式格式（與膨脹功相比，忽略對流效應影響），空泡份額對流項采用半隱格式。然后，將上述守恒方程殘差對p、pe、pw、pn、ps、ρˉL、ρˉG、TG進行泰勒展開，即可獲得迭代方程組。最后，充分利用狀態(tài)方程EOS 模型，迭代求解即可更新網格壓強，進而更新兩相速度、各組分宏觀密度、各組分比內能等其他狀態(tài)參數。上述求解過程整體邏輯框架體系如圖3 所示，其中上標3 表示第三步值。

圖3 第三步：壓力迭代算法框架體系Fig.3 Step three: pressure iterative algorithm framework

4） 第四步：時間步末更新

重復第二步具體過程，但其網格壓強、對流項速度將采用第三步最新結果，最終更新宏觀密度、比內能、速度，即。其具體算法框架與第二步極為相似，限于篇幅不再贅述，其中上標n+ 1 表示當前時間步值。

1.3 輔助模型

1） 兩相流流型模型

兩相流流型決定了各相組分間的接觸方式，繼而影響交界面面積、曳力系數、傳熱系數及傳熱、傳質計算。文中考慮三種不同流型，包括泡狀流、彌散流和過渡流，具體流型由蒸汽相體積份額判斷。當αG≤0.3時，視為純泡狀流，當αG≥0.7時，視為純彌散流，當0.3 ＜αG＜0.7 時，視為過渡流。對于控制體網格而言，可統(tǒng)一視為由泡狀流區(qū)域、彌散流區(qū)域疊加構成，并采用泡狀流區(qū)域、彌散流區(qū)域體積份額fB、fD進行具體表征，其中fB、fD按式（10）計算，其中αB、αD分別是泡狀流、彌散流特征值，分別默認取0.3、0.7。xB、xD為組分在泡狀流、彌散流區(qū)域內的分配系數，具體按式（11）計算。則泡狀流區(qū)域內液態(tài)水（連續(xù)相）體積份額可表示為αL，B=αLxB，彌散流區(qū)域內液態(tài)水（液滴）體積份額可表示為αL，D=αLxD。蒸汽混合物在泡狀流區(qū)域、彌散流區(qū)域內的體積份額按式（12）、（13）計算。

泡狀流區(qū)域內水蒸氣（氣泡）體積份額為：

彌散流區(qū)域內水蒸氣（連續(xù)相）體積份額為：

2） 相界面輸運模型

為捕捉在流動，傳熱、傳質等過程中相界面演化行為，考慮了相界面輸運方程，具體包括泡狀流區(qū)域內蒸汽泡和彌散流區(qū)域內液滴的相界面輸運行為，其輸運方程分別如下，其中：aG，B、aD，D分別是泡狀流區(qū)域內氣泡相界面積、彌散流區(qū)域內液滴相界面積。

就本文模擬而言，泡狀流區(qū)域內蒸汽泡相界面輸運考慮了對流、膨脹/壓縮、成核及破裂行為；彌散流區(qū)域內液滴僅考慮了對流效應。其中，氣泡成核即式（14）中源項SG，B，N，具體模型參考了先進流體動力學程序AFDM［10］，具體采用平衡態(tài)模型，其平衡態(tài)參數aeG，B，N、時間常數τN及源項SG，B，N按式（16）計算。蒸汽泡破裂模型采用了韋伯破裂模型，也使用平衡態(tài)方法計算，其平衡態(tài)參數aeG，B，W、時間常數τW及源項SG，B，W按式（17）計算。

3） 相間曳力模型

在泡狀流區(qū)域，“液態(tài)水-水蒸氣”接觸方式為連續(xù)液體包裹氣泡，在彌散流區(qū)域，其接觸方式為連續(xù)蒸汽包裹液滴。因此，總體上“液態(tài)水-水蒸氣”接觸方式可按“連續(xù)相-彌散相”處理。下文以連續(xù)相對彌散項施加的徑向曳力進行闡述，正如式（18）所示，其中：KDR，D，R為“連續(xù)相-彌散相”曳力系數；aCD為交界面積；uC、uD分別表示連續(xù)相和彌散相徑向速度；μC為連續(xù)相動力黏度；ρC為連續(xù)相密度；rD為彌散相半徑。系數CD參考Ishii發(fā)展的相關模型［11］。

4） 虛擬質量力模型

在兩相流中，當彌散相（如氣泡、液滴）相對連續(xù)相加速時，將產生“虛擬質量效應”，由于本文需要模擬“汽-水”兩相流壓力波傳播，因此考慮了虛擬質量力模型。具體參考了Fullmer和Bertodano報道的相關模型［12］，系統(tǒng)分析程序RELAP5/MOD3.3 也使用了該模型。

液態(tài)水所受虛擬質量力為：

式中：αL為液態(tài)水體積份額；ρL為液態(tài)水理論密度；αG為蒸汽混合物體積份額；注意到上述虛擬質量力采用的是矢量形式，CVM為常數，通常取0.5。

蒸汽混合物虛擬質量力為：

5） 傳熱傳質模型

就本模擬而言，考慮了“汽-水”兩相傳熱、傳質模型，尤其是交界面相變傳熱模型。

a） 傳質模型

傳質考慮的是交界面相變傳質，“液態(tài)水-水蒸氣”交界面凈傳熱速率可表示為：

式中：aL，G表示“液態(tài)水-蒸汽混合物”交界面積；RL，Gh2o是“液態(tài)水-水蒸氣”接觸比例；hL為交界面之液態(tài)水側換熱系數；hG為交界面之蒸汽混合物側換熱系數；TLi，Gh2o為“液態(tài)水-水蒸氣”交界面溫度，取飽和溫度與平衡態(tài)傳熱交界面溫度之較大者。qiL，Gh2o＞0表示交界面存在凈熱量傳出，則交界面發(fā)生了水蒸氣凝結現(xiàn)象，凝結速率按式（22）計算；qiL，Gh2o＜0表示交界面存在凈熱量傳入，則交界面發(fā)生了液態(tài)水汽化現(xiàn)象，汽化速率按式（23）計算；qiL，Gh2o= 0 表示交界面處于傳熱平衡態(tài)，無相變發(fā)生。

b）傳熱模型

就文中模擬對象而言，包括非平衡態(tài)相變傳熱和平衡態(tài)溫差傳熱。液態(tài)水、水蒸氣的非平衡態(tài)相變傳熱分別按式（24）、（25）計算。

液態(tài)水、水蒸氣的平衡態(tài)溫差傳熱分別按式（26）、（27）計算。

式中：RL，Ggas表示“液態(tài)水-不可凝氣體”接觸比例；TiL，Ggas表示“液態(tài)水-不可凝氣體”交界面溫度，取其平衡態(tài)溫度。

6）狀態(tài)方程EOS模型

狀態(tài)方程EOS 是封閉可壓縮兩相流守恒方程的必備條件，同時也是本文開展傳熱、傳質數值迭代，壓力迭代的重要保障。本文模擬涉及液態(tài)水、水蒸氣、不可凝氣體三種物質組分的狀態(tài)方程。

a）水蒸氣狀態(tài)方程

水蒸氣狀態(tài)方程采用了MRK（Modified Redlich-Kwong）模型，下式給出了水蒸氣分壓與其溫度和比體積之間的關系。

其中：x(TG)按式（29）計算；Th2o，crt為水臨界溫度。

b）不可凝氣體狀態(tài)方程

不可凝氣體采用理想氣體狀態(tài)方程模型。

c）液態(tài)水狀態(tài)方程

本文考慮液態(tài)水為可壓縮流體，式（30）給出了液態(tài)水溫度與其比內能和壓強間的關系，其中：TS，L為液態(tài)水飽和溫度，采用液態(tài)水比內能計算；?TL?p為液態(tài)水溫度對壓強偏導數，也使用液態(tài)水比內能計算。

類似地，式（31）給出了液態(tài)水比體積與其比內能和壓強之間的關系，形式與式（30）類似。

2 程序驗證

為驗證本文開發(fā)的汽-水兩相二維軸對稱壓力波傳播模擬程序TPFPWPC-2D，選用了典型的汽-水兩相激波管基準例題，如圖4所示，雖然該激波管屬于一維問題，但可以驗證該程序的關鍵模型與算法。該激波管長100 m，中心插入一隔板，將激波管劃分為兩個不同區(qū)域，左側為15 MPa的汽水混合物高壓區(qū)，右側為10 MPa 的汽水混合物低壓區(qū)，具體初始參數如圖4 所示。由于該程序具有二維軸對稱特性，恰可建立上述激波管，具體在徑向劃分3 個網格，軸向劃分200 個網格，激波管幾何模型豎直放置，則需忽略重力與基準例題保持一致。另外，建模時需要考慮到徑向均勻性，即壓力波僅沿軸向傳播。激波管上、下兩端采用絕對反射邊界條件。很明顯，汽-水兩相流激波管基準例題并無解析解，因此采用了程序對比方式進行驗證，具體選用了典型系統(tǒng)分析程序RELAP5 公開發(fā)布的模擬結果作為對比基準［6］。

圖4 兩相激波管示意圖Fig.4 Schematic diagram of two-phase shock tube

TPFPWPC-2D 程序在方程求解上采用了時間步分割四步算法，即人為將控制方程中傳熱、傳質、相變所致動量交換分離處置，這極可能造成時間步長敏感性問題，針對此問題，首先開展了時間步長敏感性分析，共選用15 個時間步長，即5 000 μs、2 500 μs、1 000 μs、750 μs、500 μs、250 μs、100 μs、75 μs、50 μs、25 μs、10 μs、7.5 μs、5 μs、2.5 μs、1 μs，對此基波管進行模擬。圖5給出了不同時間步長下此激波管中心位置0.081 s時刻的壓強值，很明顯時間步長在1～100 μs范圍內具有良好的數值穩(wěn)定性。同時，測試表明若時間步長大于5 ms，計算則發(fā)散（庫朗限制）。

圖5 汽水兩相激波管TPFPWPC-2D時間步長敏感性分析Fig.5 Time-step sensitivity analysis of TPFPWPC-2D based on water-vapor two-phase shock tube

根據上述時間步長敏感性分析結果，最終選用了10-5s 作為該驗證用時間步長，且上述文獻中RELAP5 也是采用了此步長。圖6 分別給出了0.081 s 時刻壓力、汽水兩相溫度、速度、空泡體積份額的分布情況。從圖6（a）可知，隨著隔膜突然撤離，汽-水兩相高壓區(qū)內形成沖擊波并傳入汽-水兩相低壓區(qū)，致使低壓區(qū)內壓力上升；相反汽-水兩相低壓區(qū)內形成膨脹波并傳入汽-水兩相高壓區(qū)，致使高壓區(qū)內壓力急劇下降。從圖6（b）可知，汽-水兩相溫度分布與壓力分布具有一定相似性。從圖6（c）可知，激波管中心兩側區(qū)域出現(xiàn)了汽-水兩相高速運動，最大速度達到約48 m·s-1。從圖6（d）可知，由于汽-水兩相膨脹、汽化效應，高壓區(qū)偏中心區(qū)域空泡份額出現(xiàn)了急劇上升，可見明顯極大值點；相反，由于汽-水兩相壓縮、凝結效應，低壓區(qū)偏中心位置空泡份額出現(xiàn)了急劇下降，可見明顯極小值點。從壓力、溫度、速度及空泡份額分布特性比較可以看出，TPFPWPC-2D程序與系統(tǒng)分析程序RELAP5模擬結果整體符合較好，但TPFPWPC-2D 與RELAP5 依然存在局部差異，很可能是由兩程序所考慮的模型精細度或個別模型差異造成的，也可能與數值算法有關系。

圖6 兩相基波管壓力分布(a)、溫度分布(b)、速度分布(c)和空泡份額分布(d)Fig.6 Pressure profile (a), temperature profile (b), velocity profile (c) and void fraction profile (d) of the two-phase shock tube

3 壓力波二維軸對稱傳播特性研究

為研究圓柱體區(qū)域內汽-水兩相流動、傳熱、傳質條件下壓力波二維軸對稱傳播特性，使用文中開發(fā)的TPFPWPC-2D 程序開展了數值模擬研究。數值模擬幾何模型如圖7 所示，建模考慮了徑向和軸向（即對稱軸）兩方向，實際模擬的幾何區(qū)域為一個高1.0 m，半徑為1.0 m 的圓柱體。徑向0～0.2 m、軸向0.1～ 0.3 m 為高壓區(qū)域，其余為低壓區(qū)域，且初始時刻兩區(qū)域處于相互隔離狀態(tài)，數值模擬開始后，兩區(qū)域流體瞬間接觸。

圖7 二維壓力波傳播實例示意圖Fig.7 Schematic of 2D pressure wave case

首先，針對該問題開展了網格敏感性分析，選用了20×20、40×40、50×50、60×60、70×70、80×80、90×90、100×100、110×110、120×120、130×130、140×140、150×150，共計13套不同網格方案，均考慮等間距劃分，并以8.0×10-4s時刻計算區(qū)域內最大壓強值為觀測點，對網格敏感性進行評價。圖8即為評價結果，很明顯，當網格劃分方案大于100×100時，數值模擬結果已具較好穩(wěn)定性，據此下文模擬選用了120×120 網格劃分方案，開展壓力波二維傳播模擬與特性分析。

圖8 二維壓力波傳播網格敏感性分析Fig.8 Mesh sensitivity analysis of 2D pressure wave propagation

圖9 給出了壓力波的二維傳播過程，注意該圖已進行對稱性拓撲。從數值模擬結果可以看出，隨著兩個不同區(qū)域內流體被解除隔離，來自高壓區(qū)域的流體向低壓區(qū)域射流，進而產生了向外快速傳播的沖擊波（shock wave），致使外部低壓區(qū)域內壓力劇增。同時，來自低壓區(qū)域內的膨脹波（rarefaction wave）也隨之產生，并傳播進入高壓區(qū)域，致使其內部壓力降低。隨著時間推移，壓力波不斷向邊界傳播，約0.4 ms 后，壓力波已抵達底部邊界，被邊界壁面完全反射，進而產生反向沖擊波，反向沖擊波與迎面沖擊波隨即相遇，繼而出現(xiàn)了疊加效應，致使下邊界附近出現(xiàn)了局部高壓區(qū)。在下一個階段，隨著壓力波不斷向外擴散，沖擊波逐漸減弱。

圖9 二維壓力波傳播實例數值模擬結果Fig.9 The simulation results of 2D pressure wave propagation case

4 結語

本文基于汽-水兩相歐拉流體動力學理論，開發(fā)了一種基于二維軸對稱圓柱坐標系的壓力波傳播模擬二維程序TPFPWPC-2D，考慮了汽-水兩相非平衡態(tài)相變傳熱、傳質，偏微分方程求解使用了時間步分離四步算法。并采用了典型的汽水兩相激波管基準例題開展了程序的驗證工作及時間步長敏感性分析，數值模擬結果與典型的系統(tǒng)分析程序RELAP5進行了對比。最后為考查該程序模擬壓力波二維傳播的能力，對圓柱體空間區(qū)域內的壓力波傳播行為進行了二維軸對稱數值模擬，同時還開展了網格敏感性分析。本文研究得到的結論如下：

1）針對汽-水兩相激波管基準例題，TPFPWPC-2D 程序計算獲得的壓力、溫度、速度及空泡份額分布數據與典型的系統(tǒng)分析程序RELAP5 符合較好，程序總體上可以很好地模擬強傳熱傳質兩相條件下壓力波傳播行為，進而從整體上驗證了該程序的結構模型及數值算法的正確性。同時，從時間步長敏感性分析結果得知，本文采用的時間分割四步算法具有較為顯著的時間步長敏感性。

2）針對圓柱體空間內壓力波二維軸對稱傳播算例，采用TPFPWPC-2D 程序開展了數值模擬，獲得壓力波在圓柱體內的二維軸對稱傳播特性，尤其是捕捉到了壓力波在二維傳播過程中反射、疊加行為，進而從一定層面上確認了該程序模擬壓力波二維傳播問題的能力。

作者貢獻聲明龔政宇實施研究，分析/解釋數據，起草撰寫文章；辜峙钘提供研究思路，指導論文寫作與修改，提供研究經費支持；潘麒文負責材料支持，提供相關文獻；張牧昊負責文章內容的審閱指導；戴嘉寧參與文獻收集整理工作，數據整理。