李 偉 蘇光輝 秋穗正 田文喜 陳 平 李垣明

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

2（西安交通大學 核科學與技術(shù)學院 西安 710049）

一種穩(wěn)定性增強及高精度數(shù)值方法在RELAP5中的實現(xiàn)與評價

李 偉1蘇光輝2秋穗正2田文喜2陳 平1李垣明1

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

2（西安交通大學 核科學與技術(shù)學院 西安 710049）

在計算穩(wěn)定性的改進方面，修正RELAP5程序中的虛擬質(zhì)量力(Virtual mass force)形式，同時添加了新的界面壓力項(Interface pressure)；在計算精度的改進方面，采用具有總變差減小(Total variation diminish, TVD)特點的高精度通量限制器(Flux limiter)方法取代RELAP5程序原本的一階迎風方法來離散質(zhì)量和能量守恒方程中的對流項。在模擬水平管道內(nèi)空泡份額微擾隨時間發(fā)展的數(shù)值實驗中，相比改進前的RELAP5，改進后的RELAP5計算得到的微擾幅度并未增長；在模擬液相沉降的數(shù)值實驗中，改進前的RELAP5程序計算得到了不真實的空泡份額分布，而改進后的RELAP5在不同的網(wǎng)格數(shù)量下能夠得到收斂的穩(wěn)定解。對Ransom 水龍頭數(shù)值實驗和Marviken CFT 15大破口噴放實驗的計算表明，改進后的具有二階TVD格式的RELAP5程序能夠得到更接近實驗數(shù)據(jù)的計算結(jié)果。

RELAP5，穩(wěn)定性，精確性，雙曲性，虛擬質(zhì)量力，界面壓力項

RELAP5程序長期以來在反應堆系統(tǒng)的熱工水力及安全分析領(lǐng)域都有著廣泛的應用，包括乏燃料池(Spent fuel pool)的安全設(shè)計[1]。但RELAP5的基本偏微分方程組包含虛的特征根，因此不具備雙曲性，這種兩流體模型不具備適定性[2-3]。在實際應用中可能導致由不適定性引起的數(shù)值解發(fā)散，出現(xiàn)不真實的計算結(jié)果。在數(shù)值離散方面，RELAP5作為系統(tǒng)程序采用了稀疏的交錯網(wǎng)格，動量守恒方程中的對流項采用一階迎風方法進行離散。這有利于增強程序的數(shù)值穩(wěn)定性和提高計算效率，然而也使得數(shù)值擴散過大，影響計算精度。

本文基于RELAP5/MOD3程序，從改善計算穩(wěn)定性和提高數(shù)值精度方面對其進行性能改進并且通過實驗驗證。在計算穩(wěn)定性改進方面，分析了動量守恒方程中不同的虛擬質(zhì)量力形式以及界面壓力項對RELAP5六方程特征根的影響，修正了RELAP5原本的虛擬質(zhì)量力形式和增加界面壓力項使其特征根全部為實數(shù)來增強模型的雙曲性。設(shè)計一個空泡份額微擾隨時間變化的數(shù)值實驗來驗證改進后的RELAP5程序在計算穩(wěn)定性方面的性能；通過一個相沉降數(shù)值實驗來驗證改進后的RELAP5程序預測結(jié)果在不同網(wǎng)格劃分方案下的行為。在數(shù)值精度改進方面，Wang等[4]研究了二階總變差減小(Total variation diminish, TVD)格式在系統(tǒng)程序TRACE中的實現(xiàn)。TRACE是系統(tǒng)程序TRAC和RELAP5組合后的版本。發(fā)現(xiàn)二階TVD格式比一階迎風格式的計算更準確，且與一階迎風格式一樣高效和穩(wěn)定。張小英等[5]開發(fā)了二階精度的一維兩流體模型程序TFIT。該程序未考慮能量方程，用于絕熱不可壓縮問題。通過對Ransom水龍頭問題進行對比計算，二階TVD格式比一階精度格式的計算準確度明顯提高，而且同樣保持良好的穩(wěn)定性。本文采用結(jié)合了Minmod flux limiter的二階TVD格式代替RELAP5程序中對流項的一階迎風離散格式，用于實際工程計算。基于Ransom水龍頭數(shù)值實驗和實際工程試驗Marviken CFT 15對改進后的程序進行驗證，評估對數(shù)值擴散問題的改善。

1 計算穩(wěn)定性改進

1.1RELAP5兩流體模型的非雙曲性

為分析RELAP5兩流體模型的非雙曲性，將其6個守恒方程改寫為如下的矩陣方程形式：

表示RELAP5的直接求解變量，分別為壓力、空泡份額、氣相速度、液相速度、氣相內(nèi)能和液相內(nèi)能；A、B和C為包含未知量的矩陣[6×6]。設(shè)為式(1)的特征根，它們滿足以下特征方程：

只有以下條件成立時[5]，特征根才全部為實數(shù)：

或：

1.2虛擬質(zhì)量力的影響

虛擬質(zhì)量力[6]在多相流問題中刻畫了離散相流體在連續(xù)相流體中運動時，帶動其周圍的連續(xù)相流體運動所需要的額外力。RELAP5采用如下形式的虛擬質(zhì)量力：

式中：ma具有氣液混合物等效音速的意義。如果ε在數(shù)量級上足夠小，即氣、液相的速度大小相對值遠小于混合物音速時，經(jīng)分析須滿足如下的條件才能保證所有的特征根均為實數(shù)：

在大部分情況下，ε都可以被認為是數(shù)量級很小的量。當氣、液相速度差與音速相當時，Tiselj等[7]認為無論虛擬質(zhì)量力系數(shù)如何取值也無法使得方程組具備雙曲性。本文在采用式(10)的虛擬質(zhì)量力系數(shù)的基礎(chǔ)上，進一步研究界面壓力項對兩流體模型適定性的影響。

1.3界面壓力項的影響

在RELAP5的模型假設(shè)中，氣、液和相間界面的壓力相等，因此也稱為單壓力模型，這造成了守恒方程中某些重要信息的缺失，是導致模型不具備適定性的重要因素。在另一個著名的反應堆熱工水力系統(tǒng)分析程序CATHARE中，通過在動量守恒方程中引入相間界面壓力項[8]來考慮各相壓力的差異，即：

式中：iPΔ為界面壓力項系數(shù)。界面壓力項在動量方程中引入了空泡份額的空間偏微分項。同時存在虛擬質(zhì)量力和界面壓力項時，方程組的特征方程化解為：

式(14)判別式非負，求得界面壓力項系數(shù)須滿足如下關(guān)系式：

從物理量綱來看iPΔ具有動壓的意義。當虛擬質(zhì)量力的系數(shù)為零時，式(15)與CATHARE程序中的形式一致。當虛擬質(zhì)量力的系數(shù)滿足式(10)時，即時，無需界面壓力項方程組也能滿足雙曲性要求，否則界面壓力項能夠?qū)μ摂M質(zhì)量力的效果起到補充作用（例如在水平分層流動中虛擬質(zhì)量力系數(shù)為零的情況下）。雖然虛擬質(zhì)量力和界面壓力項的目的都與改善RELAP5兩流體模型的適定性有關(guān)，但是二者之間存在一定的區(qū)別：虛擬質(zhì)量力在兩相流中是具有對應物理意義的，而界面壓力修正項則缺乏堅實的物理根據(jù)。它雖然考慮到單壓力模型中界面壓力與各相平均壓力之間并不相等，但是真正的界面壓力非常復雜，因此采用經(jīng)驗關(guān)系式近似。不可把式(11)和(12)看作是真實物理過程中的界面壓力，可以理解為使兩流體模型的方程組具有雙曲性而人工添加的微分項，所以為防止過分強調(diào)其影響，一般在合適的取值范圍內(nèi)取最小的值。在程序的實際實施過程中，還通常對iPΔ乘以一個稍大于1的系數(shù)ξ。此外，RELAP5在水平分層流動時考慮到了界面壓力與各相平均壓力的不同而使用了這種壓力修正項，CATHARE則在所有流型中予以使用，以改善基本守恒方程的雙曲性。

2 計算精度改進

RELAP5采用了一階精度的數(shù)值方案進行對流項的離散，大量實驗和工程應用證明這對大部分系統(tǒng)瞬態(tài)行為的預測能夠達到相應的精度要求，同時也簡化了程序的編程難度，提高了計算效率（特別是在進行不確定性分析時）。RELAP5采取的是交錯網(wǎng)格方案，質(zhì)量和能量守恒方程中的對流項采用一階迎風格式[9]進行離散：

采用Waterson[11]關(guān)于Flux limiter的通用概念，對流項的離散格式可表達為低階（一階迎風）和高階格式的組合：

式中：KxΔ(LxΔ)為標量控制體單元的長度；/xφ??為標量控制體單元界面上的參數(shù)梯度；()rψ為Flux limiter函數(shù)衡量控制體單元界面j上參數(shù)的平滑程度。Flux limiter函數(shù)的具體形式將會對程序的計算精度和穩(wěn)定性產(chǎn)生重要影響。合適的Flux limiter函數(shù)形式應該符合Flux limiter和TVD格式的初衷，即提高程序計算精度的同時，保證程序具有合理的穩(wěn)定性。Flux limiter的形式也不宜過于復雜，以致降低整個程序的運行效率。出于實際計算中的穩(wěn)定性和精度考慮，本文采用Minmod flux limiter函數(shù)

3 改進措施的驗證

3.1計算穩(wěn)定性改進措施驗證

數(shù)值實驗在RELAP5程序的開發(fā)與驗證過程中起到重要作用。采用類似文獻[8]中的空泡份額微擾實驗對穩(wěn)定性改進措施進行驗證。在長度5 m的水平管道中，流動著0.1 MPa的飽和蒸汽和飽和水，速度分別保持0.1 m·s-1和1.0 m·s-1。除了距離入口2.0 m處的初始空泡份額為0.5008，其他位置的初始空泡份額均為0.5。進行數(shù)值模擬時，整個管道劃分為99個控制體，不計壁面摩擦。在這種情況下，由于數(shù)值擴散的作用，空泡份額微擾將隨時間而衰減。但是對于不具備雙曲性的模型，擾動可能迅速增長，最終計算結(jié)果發(fā)散。在本例中虛擬質(zhì)量力的效應相當微弱，因此對于改進前的RELAP5，式(7)可能無法得到滿足，如圖1(a)所示，改進前的RELAP5計算得到的空泡份額值迅速發(fā)散。如圖1(b)所示，對于改進后的RELAP5在虛擬質(zhì)量力和界面壓力項的共同作用下模型保持了雙曲性，微擾的增長受到抑制，保持了計算的穩(wěn)定性。圖2為改進后的RELAP5在不考慮虛擬質(zhì)量力時計算得到的空泡份額分布隨時間的變化，可見微擾迅速衰減。本例說明虛擬質(zhì)量力以及界面壓力項的作用在程序中得到了實現(xiàn)，改進后的RELAP5相比改進前的RELAP5計算穩(wěn)定性明顯提高。

圖1 空泡份額分布隨時間變化(a) 改進前的RELAP5，(b) 改進后的RELAP5Fig.1 Time-dependent distribution of the void fraction.(a) Default RELAP5, (b) Improved RELAP5

圖2 改進后的RELAP5不考慮虛擬質(zhì)量力時空泡份額分布隨時間的變化Fig.2 Time-dependent distribution of the void fraction without virtual mass term of improved RELAP5.

考慮更復雜的情況，采用如圖3所示的液相沉降數(shù)值實驗進行驗證。利用改進后的RELAP5分別計算了控制體劃分數(shù)目(Control Volume, CV)為25、49和99的情況，并且與改進前的RELAP5進行比較。圖4給出了改進前/后的RELAP5計算得到控制體1中的空泡份額計算值隨時間的變化。隨著控制體數(shù)目增多，控制體單元尺寸減小，改進前的RELAP5模型非雙曲性、不適定的缺點表現(xiàn)得更明顯。在精細的網(wǎng)格下更容易出現(xiàn)不穩(wěn)定性，計算過程中也出現(xiàn)了不真實的結(jié)果。而改進后的RELAP5模型適定性得以改善，即使計算網(wǎng)格加密計算結(jié)果也不容易出現(xiàn)不穩(wěn)定的行為，達到了改進的目的。

圖3 沉降數(shù)值實驗示意圖Fig.3 Sketch of the phase sedimentation numerical test.

圖4 控制體1內(nèi)的空泡份額隨時間變化 (a) 改進前的RELAP5，(b) 改進后的RELAP5Fig.4 Variation of the void fraction in CV #1 with time. (a) Default RELAP5, (b) Improved RELAP5

3.2計算精度改進措施驗證

采用Ransom[12]提出的水龍頭(Water faucet)數(shù)值實驗來對本文所作的精度改進進行驗證。Water faucet問題常被用于驗證兩相流求解器的數(shù)值性能。如圖5所示，在長度12 m、直徑1 m的豎直圓管中初始時刻充滿常溫常壓靜止的空氣，初始時刻圓管中還存在一股射流，其速度處處為10 m·s-1，且軸向的空泡份額均為0.2。管道入口空泡份額和液體速度分別保持為0.2 m·s-1和10 m·s-1，出口保持常壓。計算開始后，水柱各部分在重力的作用下加速運動。不計管道壁面摩擦，整個過程絕熱。Ransom在假設(shè)液相不可壓縮的基礎(chǔ)上推導了該問題的精確解：

式中：0gα為初始空泡份額；gα為空泡份額；z為距管道入口的距離，m；0fν為液相初始速度，m·s-1；H為管道總長度。

圖5 Ransom水龍頭數(shù)值實驗示意圖Fig.5 Sketch of the Ransom water faucet numerical test.

分別劃分96和120個控制體時，改進前/后計算得到的空泡份額分布與分析解的對比如圖6所示。與改進前的RELAP5一階迎風格式(Upwind Difference, UD)相比，隨著控制體數(shù)目的增多，二階TVD格式(Minmod)更接近精確解(Analytic)。一階迎風格式在控制體單元尺寸過大的情況下，數(shù)值擴散更加明顯。從本次數(shù)值實驗可以看出，改進后的RELAP5二階精度格式體現(xiàn)了程序捕捉流場間斷的優(yōu)勢。

圖6 空泡份額沿流動方向上的分布(a) 96個控制體，(b) 120個控制體Fig.6 Variation of the void fraction along the flow. (a) 96 CVs are employed, (b) 120 CVs are employed

采用Marviken CFT 15全尺寸大型臨界流實驗來進一步驗證改進措施[13]。實驗中存放氣液混合物的壓力容器高24.55 m，平均內(nèi)徑5.2 m，排出管長6.3 m，內(nèi)徑0.72 m，噴嘴直徑0.5 m，長徑比為3.6。采用的控制體劃分方案如圖7所示，采用Henry-Fauske臨界流噴放模型。

圖7 Marviken CFT 15控制體劃分方案Fig.7 Nodalization scheme for the Marviken CFT 15 experiment.

圖8給出了改進前/后的RELAP5計算得到的CFT 15破口流量與實驗數(shù)據(jù)的比較。從圖8中可見，改進后的RELAP5計算得到的破口流量計算值與實驗數(shù)據(jù)符合得很好，只在兩相噴放階段出現(xiàn)了一定的偏差。如圖8(b)所見，改進后的RELAP5由于精度提高，計算穩(wěn)定性增強，能夠更好地預測從破口噴放的流量隨時間的變化。

圖8 改進前的RELAP5 (a)和改進后的RELAP5 (b)計算得到的破口流量Fig.8 Break flow rate calculated by the default RELAP5 (a) and improved RELAP5 (b).

4 結(jié)語

為改善RELAP5兩流體模型非雙曲性的缺點，通過分析特征根得到虛擬質(zhì)量力和界面壓力項，用于保證模型應滿足的條件；為改善RELAP5兩流體模型僅具有一階精度的不足，將二階TVD格式應用到程序中，同時為兼顧計算穩(wěn)定性選擇了能夠獲得二階精度、形式簡單的Minmod 函數(shù)作為Flux limiter替代RELAP5的一階迎風格式。實驗驗證結(jié)果表明，對RELAP5的改進措施是可行、合理的，達到了改進程序性能的目的，為更好地發(fā)揮系統(tǒng)程序在安全分析領(lǐng)域的作用具有積極意義，對改善其他最佳估算程序的性能也具有借鑒意義。

1 Ge L, Wang H T, Zhang G L, et al. Thermal-hydraulic design and transient analysis of passive cooling system for CPR1000 spent fuel storage pool[J]. Nuclear Science and Techniques, 2016,27(1): 8. DOI: 10.1007/s41365-016-0017-6

2 Lyczkowski R W, Gidaspow D, Solbrig C W, et al. Characteristics and stability analyses of transient one-dimensional two-phase flow equations and their finite difference approximations[J]. Nuclear Science and Engineering, 1978,66(3): 378-396. DOI: 10.13182/ NSE78-4

3 Lax P D. On the nation of hyperbolicity[J]. Communications on Pure and Applied Mathematics, 1980,33(3): 395-397. DOI: 10.1002/cpa.3160330309

4 Wang D, Mahaffy J H, Staudenmeier J, et al. Implementation and assessment of high-resolution numerical methods in TRACE[J]. Nuclear Engineering and Design, 2013,263(0): 327-341. DOI: 10.1016/ j.nucengdes.2013.05.015

5 張小英, 丁斐, 陳佳躍. 反應堆一維兩流體模型二階精度數(shù)值解法研究[J]. 核動力工程, 2013,34(4): 27-32. DOI: 10.3969/j.issn.0258-0926.2013.04.007

ZHANG Xiaoying, DING Fei, CHEN Jiayue. Study on numerical solution of the two dimensional two fluid model of the reactor[J]. Nuclear Power Engineering, 2013,34(4): 27-32. DOI: 10.3969/j.issn.0258-0926.2013.04.007

6 Drew D, Cheng L, Lahey R. The analysis of virtual mass effects in two-phase flow[J]. International Journal of Multiphase Flow, 1979,5(4): 233-242. DOI: 10.1016/ 0301-9322(79)90023-5

7 Tiselj I, Petelin S. Modelling of two-phase flow with second-order accurate scheme[J]. Journal of Computational Physics, 1997,136(2): 503-521. DOI: 10.1006/jcph.1997.5778

8 Hogon L, Lee U, Kim K, et al. Application of hyperbolic two-fluids equations to reactor safety code[J]. Nuclear Engineering and Technology, 2003,35(1): 45-54

9 RELAP5 Code Development Team. RELAP5/MOD3 code manual Vol I: code structure, system models, and solution methods[R]. Idaho Falls: Idaho National Laboratory, 2001

10 Versteeg H K, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method[M]. Essex, UK: Pearson Education, 2007

11 Waterson N P, Deconinck H. Design principles for bounded higher-order convection schemes - a unified approach[J]. Journal of Computational Physics, 2007,224(1): 182-207. DOI: 10.1016/j.jcp.2007.01.021

12 Ransom V. Numerical benchmark tests[A]. Hewitt G, Delhay J, Zuber N. Multiphase science and technology[M]. Washington DC: Hemisphere Publishing Corporation, 1987

13 Kim K, Kim H J. Assessment of RELAP5/MOD2 critical flow model using Marviken test data 15 and 24[R]. Washington, DC: Office of Nuclear Regulatory Research, U.S. Nuclear Regulatory Commission, 1992

Implementation and assessment of a stability-enhancing and high-resolution numerical scheme in RELAP5

LI Wei1SU Guanghui2QIU Suizheng2TIAN Wenxi2CHEN Ping1LI Yuanming1
1(Science Technology on Reactor System Design Technology Laboratory,Nuclear Power Institute of China,Chengdu 610213,China)
2(School of Nuclear Science and Technology,Xi’an Jiaotong University,Xi’an 710049,China)

Background：RELAP5 has been widely used in the field of nuclear thermal-hydraulics and safety analysis. However, improvement issues in terms of numerical stability and accuracy exist due to its basic single-pressure two-fluid model.Purpose：This study aims to enhance the numerical stability and numerical accuracy of RELAP5 by modifying the mathematical structure of RELAP5 basic model.Methods：Firstly, the default mathematical expression of the virtual mass force was modified in the RELAP5 code. Secondly, an interface pressure term was incorporated in the code, and a total variation diminish (TVD)-type flux limiter method was implemented in the RELAP5 code in place of the default 1st-order upwind scheme for the convective terms in mass and energy conversation equations.Results：In a numerical test where the evolvements of void fraction perturbation were concerned, the modified RELAP5 predicted no growth of the perturbation amplitude in contrast to a rapidly developing divergence for the default RELAP5. In addition, for the phase segmentation problem, convergence was accomplished on finer mesh by the modified RELAP5, while non-physical distribution of the void fraction was rendered by the default RELAP5 especially for finer mesh. When applied to the Ransom water faucet numerical test and Marviken CFT 15 experiments, the improved RELAP5 with 2nd-order of accuracy was more reliable in that thepredictions were closer to the experimental results.Conclusion：Compared with the original RELAP5, the modified or improved RELAP5 showed evident numerical stability and accuracy during the assessment in this work.

RELAP5, Stability, Accuracy, Hyperbolicity, Virtual mass, Interface pressure

TL33

10.11889/j.0253-3219.2016.hjs.39.110601

李偉，男，1986年出生，2015年于西安交通大學核科學與技術(shù)專業(yè)獲博士學位，主要從事核燃料設(shè)計、反應堆熱工水力及安全分析First author: LI Wei, male, born in 1986, graduated from Xi’an Jiaotong University with a doctoral degree in 2015, major in nuclear science and technology, mainly focusing on nuclear fuel research and design, reactor thermal-hydraulic and safety analysis

2016-05-30，

2016-09-01