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

1（成都理工大學(xué) 核技術(shù)與自動(dòng)化工程學(xué)院 成都 610059）

2（中國核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室 成都 610213）

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

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

1 模型和算法

1.1 控制方程

本文涉及的控制方程主要包括質(zhì)量守恒方程、動(dòng)量守恒方程、能量守恒方程，且全部控制方程均基于二維軸對(duì)稱圓柱體坐標(biāo)系導(dǎo)出。

1） 質(zhì)量守恒方程

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

2） 動(dòng)量守恒方程

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

3） 能量守恒方程

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

1.2 數(shù)值算法

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

1） 第一步：網(wǎng)格內(nèi)源相計(jì)算

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

圖1 第一步：網(wǎng)格內(nèi)源項(xiàng)計(jì)算數(shù)值算法框架體系Fig.1 Step 1: numerical algorithm framework for cell-inside source term calculation

2） 第二步：網(wǎng)格間對(duì)流計(jì)算

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

圖2 第二步：網(wǎng)格間對(duì)流計(jì)算數(shù)值算法框架體系Fig.2 Step 2: numerical algorithm framework for intercell convection calculation

3） 第三步：壓力迭代

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

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

4） 第四步：時(shí)間步末更新

重復(fù)第二步具體過程，但其網(wǎng)格壓強(qiáng)、對(duì)流項(xiàng)速度將采用第三步最新結(jié)果，最終更新宏觀密度、比內(nèi)能、速度，即。其具體算法框架與第二步極為相似，限于篇幅不再贅述，其中上標(biāo)n+ 1 表示當(dāng)前時(shí)間步值。

1.3 輔助模型

1） 兩相流流型模型

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

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

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

2） 相界面輸運(yùn)模型

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

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

3） 相間曳力模型

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

4） 虛擬質(zhì)量力模型

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

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

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

蒸汽混合物虛擬質(zhì)量力為：

5） 傳熱傳質(zhì)模型

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

a） 傳質(zhì)模型

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

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

b）傳熱模型

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

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

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

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

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

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

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

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

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

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

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

本文考慮液態(tài)水為可壓縮流體，式（30）給出了液態(tài)水溫度與其比內(nèi)能和壓強(qiáng)間的關(guān)系，其中：TS，L為液態(tài)水飽和溫度，采用液態(tài)水比內(nèi)能計(jì)算；?TL?p為液態(tài)水溫度對(duì)壓強(qiáng)偏導(dǎo)數(shù)，也使用液態(tài)水比內(nèi)能計(jì)算。

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

2 程序驗(yàn)證

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

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

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

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

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

圖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 壓力波二維軸對(duì)稱傳播特性研究

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

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

首先，針對(duì)該問題開展了網(wǎng)格敏感性分析，選用了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，共計(jì)13套不同網(wǎng)格方案，均考慮等間距劃分，并以8.0×10-4s時(shí)刻計(jì)算區(qū)域內(nèi)最大壓強(qiáng)值為觀測(cè)點(diǎn)，對(duì)網(wǎng)格敏感性進(jìn)行評(píng)價(jià)。圖8即為評(píng)價(jià)結(jié)果，很明顯，當(dāng)網(wǎng)格劃分方案大于100×100時(shí)，數(shù)值模擬結(jié)果已具較好穩(wěn)定性，據(jù)此下文模擬選用了120×120 網(wǎng)格劃分方案，開展壓力波二維傳播模擬與特性分析。

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

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

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

4 結(jié)語

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

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

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

作者貢獻(xiàn)聲明龔政宇實(shí)施研究，分析/解釋數(shù)據(jù)，起草撰寫文章；辜峙钘提供研究思路，指導(dǎo)論文寫作與修改，提供研究經(jīng)費(fèi)支持；潘麒文負(fù)責(zé)材料支持，提供相關(guān)文獻(xiàn)；張牧昊負(fù)責(zé)文章內(nèi)容的審閱指導(dǎo)；戴嘉寧參與文獻(xiàn)收集整理工作，數(shù)據(jù)整理。