Nuclear Technology ISSN: 0029-5450 (Print) 1943-7471 (Online) Journal homepage: http://www.tandfonline.com/loi/unct20 On the Subcooled Critical Flow Model in RELAP5/ MOD3 W. S. Yeung & J. Shirkov To cite this article: W. S. Yeung & J. Shirkov (1996) On the Subcooled Critical Flow Model in RELAP5/MOD3, Nuclear Technology, 114:1, 141-145, DOI: 10.13182/NT96-A35230 To link to this article: http://dx.doi.org/10.13182/NT96-A35230 Published online: 13 May 2017. Submit your article to this journal View related articles Full Terms & Conditions of access and use can be found at http://www.tandfonline.com/action/journalInformation?journalCode=unct20 Download by: [University of Florida] Date: 25 October 2017, At: 21:41 ON THE SUBCOOLED CRITICAL FLOW MODEL IN RELAP5/MOD3 W. S. YEUNG and J. SHIRKOV Yankee Atomic Electric Company Bolton, Massachusetts 01740 HEAT TRANSHR AND FLUID f Lii,V TECHNICAL NOTE KEYWORDS: RELAP5/M0D3, choking, fluid flow Downloaded by [University of Florida] at 21:41 25 October 2017 Received December 12, 1994 Accepted for Publication July 12, 1995 regions have been compared to data from separate effect and integral test data. During the code development stage, An analysis of an anomaly in the subcooled critical flow several deficiencies in the choking model have been identimodel in the RELAP5/MOD3 computer code is presented. fied. These deficiencies were subsequently corrected.3 One Specifically, the code produces a discontinuity in going fromof the deficiencies concerns a discontinuity in single-phase unchoked subcooled liquid flow (i.e., subsonic flow) to sub- steam critical flow rate due to incorrect throat density calcooled choked flow (i.e., sonic flow). The same anomaly has culation. However, no mention is made of the discontinubeen reported elsewhere. The root cause for this anomaly ity regarding unchoked and choked flow in the subcooled has been analyzed, and it is found that the user-supplied region. One possible explanation is that almost all assessjunction loss coefficient and discharge coefficient play an ments are made against transient experiments in which subimportant role in the occurrence of this anomaly. The analy-cooled critical flow occurs at the start of the simulation, and sis is verified by assessment against a test problem simulat- unchoked subcooled liquid flow is seldom simulated. It is ing single-phase liquid flow through a convergent nozzle also for this reason that the observed discontinuity has no with a fixed upstream pressure and a varying downstream impact on loss-of-coolant accident-type calculations using pressure. A corrective measure to eliminate the discontinu- the RELAP5 series of codes. ity is suggested. The objective of this paper is to investigate the root cause of this discontinuity. First, a test problem that reveals the discontinuity is described. An analysis is then presented for the root cause of the discontinuity. Results from actual RELAP5/MOD3 calculations of the simple problem are I. INTRODUCTION used to substantiate the analysis. Finally, a modification in the subcooled choking criterion, which eliminates the disDuring the course of applying the RELAP5/MOD3 continuity, is suggested. Version 3.1 computer code to calculate the steady-state flow rate of subcooled liquid through a control valve, a disconII. DESCRIPTION OF TEST PROBLEM tinuity in the mass flow rate was discovered. The discontinuity occurred as the flow went from unchoked to choked flow as the downstream pressure was lowered. Other versions Consider the ideal (i.e., frictionless) flow of subcooled of the code (including a version of the RELAP5/MOD1 liquid through a converging nozzle as shown in Fig. 1. The code) have been used, but the same discontinuity was calupstream conditions are fixed at Pup and T. The downstream culated to exist for each version. It appears that this disconpressure is Pdaw„. As the downstream pressure decreases tinuity exists in the RELAP5 code series. from Pup, the liquid flow rate through the nozzle increases. The increase continues until the local pressure at the nozThe same anomaly has been reported by Petelin et al.1 zle exit reaches a value Pvap, at which the liquid begins to during their investigation of the International Standard flash. The flow becomes choked, and further decrease of Problem 27 (BETHSY small-break Test 9.1b). Those auPdown wiH not increase the flow rate through the nozzle. thors suggested that the anomaly is closely connected with The pressure Pvap, corresponding to the inception of the the increase in the average junction liquid density and that flashing, is usually less than the saturation pressure this parameter may influence the choking criterion in a specorresponding to the liquid temperature T. Figure 2 depicts cial way when the abrupt area change option is applied. the steady-state flow rate as a function of Pdow„, as physiThe RELAP5 choking model has been well documented in Vol. IV of the RELAP5/MOD3 code manual.2 The model cally expected. has been extensively assessed by the RELAP5 user commuFigure 3 shows the nodalization of the test problem. The nity and by the code developers. Calculated critical flow nozzle is represented by a single volume (203), which is conrates in the subcooled, two-phase, and single-phase steam nected on either side to a time-dependent volume (201 and 201 are fixed at Pup and T. The initial conditions for SV 203 are the same as those for TV 201. For TV 205, the conditions are specified at Pdown and X — 1.0. For each value of Pdown. the problem is run until steady state is obtained. The corresponding flow rate is recorded. Figure 4 shows the calculated results for the following input data: Aup = 0.1 m 2 , AT = 0.001 m 2 , L = \m,Pup = 6.0 MPa, T = 500 K. The abrupt area option is used for the valve junction, with zero user input loss coefficient. As can be seen, the flow increases gradually as Pdown decreases. At - 2 . 3 MPa, the flow becomes choked, and the flow rate suddenly jumps to a larger value and remains there for Pdown < 2.3 MPa. The cause for this discontinuity is discussed in Sec. III. Fig. 1. Flow of subcooled liquid through nozzle. Downloaded by [University of Florida] at 21:41 25 October 2017 III. ROOT CAUSE ANALYSIS OF DISCONTINUITY The theoretical aspects of the RELAP5 choking model have been given elsewhere.2 Here we shall focus on singlephase subcooled liquid. In essence, a critical choking velocity Vc is calculated at the throat based on the Bernoulli equation and an empirical correlation for Pvap. The result is given by 2 (PUn Pvap) 2, 2. V — ' c= V ' up+' (1) where Vup is the upstream velocity and p is the liquid density, which is assumed constant in this study. From continuity, KpAup = VCAT . (2) Hence, Eq. (1) can be rewritten as K2 P,down/P,up Fig. 2. Variation of steady-state flow rate with downstream pressure. TV 201 P=Pup,T A=Aup SJ 202 SV 203 Pup, T A=Aup L VJ 204 AT TV 205 P=Pdown X=1.0 A=Adown Fig. 3. Nodalization of test problem. 2 ( ^up p[ 1 - Pvap) (3) (AT/Aup)2] where Pvap is evaluated from the Jones-Alamgir-Lienhard correlation 4 and is a function of Vc. Equation (3) therefore represents an implicit equation in Vc and can be solved by iteration. If a discharge coefficient CD is specified, Eq. (3) will be modified to V} = 2 ( Pup Pvap) Co p[\ - (Ar/Aup)2] (4) If the flow is not choked, the junction velocity will be given by the Bernoulli equation applied across the junction. Refer to the nodalization diagram, Fig. 3, and apply the Bernoulli equation from volume 203 to volume 205: Pup + iPVlp - Pdown + \pVdown + K\pVj From continuity, Vup= VjAT/Aup and Vdown = Hence Eq. (5) becomes . (5) VjAT/Adow„. 2 (Pup ^down Vf = 205). The nozzle exit is modeled by a trip valve junction (204), and the nozzle inlet is modeled by a single junction (202). The flow area for TV 201, SJ 202, and SV 203 is the same and is denoted by Aup, the exit area is denoted by Ar, and the area for TV 205 is denoted by Adow„. The choking option is specified at the valve junction only. Wall friction is deactivated for frictionless flow. The conditions for TV K + (AT/Adown)2 - (AT/Aup)2 (6) In Eq. (6), K is the junction loss coefficient, either internally calculated (if the abrupt area option is used) or user specified, or both. The choking criterion is applied as follows. At each time step, the junction velocity from Eq. (6) is compared to the critical velocity from Eq. (4). If Vj < Vc, the flow is not choked, and the junction velocity will be given by Eq. (6) Downloaded by [University of Florida] at 21:41 25 October 2017 90 0-i 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Downstream Pressure (MPa) Fig. 4. Variation of steady-state flow rate against downstream pressure. and is a function of the downstream pressure. If Vj > Vc, the flow is choked, and the junction velocity will be set equal to Vc regardless of the downstream pressure. At the point of choking, therefore, Vj = Vc. Thus, equating Eqs. (4) and (6), one could solve for the downstream pressure necessary to choke the flow as Keff Cp ( PUp ~ Pvap) Pdown Pup 1- (AT/AUP)2 (7) where Ke/f is defined as y+(—y \A ) • \AUnJ (8) down III.A. Nature of Discontinuity The cause of the observed discontinuity in the subcooled liquid region can be best explained with Eq. (7), and depends on its solution. Define a parameter such that Z = KeffCo 1- (AT/Aup)2 (9) that the solution of Pdow„ significantly impacts the choking logic and hence the steady-state flow rate response. To illustrate the effect of Z on the steady-state flow rate, the test problem was run with various combinations of K and CD. The geometrical and boundary conditions were the same as before. The smooth area option was used with user input K values. Figure 5 summarizes the results for various Z. Each data point represents a steady-state run. Each mass flow rate result was normalized with respect to CDmcrit, where m cril is the critical flow rate corresponding to unit discharge coefficient. The effect of Z is evident. For Z > 1, the steady-state flow rate is discontinuous across the choked point. For Z < 1, the flow rate is generally smooth across the choked point. The following explanation is put forth. 1. Z < 1: Because Pup is always >Pvap, Eq. (10) indicates that the downstream pressure for choking is >P v a p . The steady-state flow rate is continuous across the choked point, as it should be. Note that for small values of Z, the choking downstream pressure can be relatively high. Indeed, if Z = 0, Pdown = Pup for choking, which is not physical. It may be concluded that Z cannot be zero for this test problem. 2. Z > 1: For this case, Pdown will be <Pvap for choking. For large values of Z, Pdovm may be significantly less than Pvap Pdown ~ Pvap = (1 — 2)(Pup — Pvap) • (10) and may even become negative. According to our analysis, the steady-state flow rate should still be continuous, as long Depending on the values of Z, the downstream pressure for as Pd0wn remains positive. In addition, if Pdown is negative choking may be greater or less than Pvap. It will be shown Equation (7) can then be rewritten as Yeung and shirkov 1.0 SUBCOOLED CRITICAL FLOW MODEL I ill I 0.9 0.8 0.7 0 1 0.6 Iu. •o 0.5 | (0 Downloaded by [University of Florida] at 21:41 25 October 2017 I z 0.4 0.3 0.2 0.1 0.0 1.5 Pdown/Psat Fig. 5. Variation of steady-state flow rate against downstream pressure for various Z. according to Eq. (10), then the flow should not be choked at all. However, the code does not behave in this manner. The discontinuity is a manifestation of choking being predicted earlier than our simple analysis would indicate. For instance, consider the Z = 2.62 curve. The RELAP5/MOD3 code predicted choking even though the junction flow rate from the normal momentum solution is - 3 0 % lower than the critical flow rate. It is possible that during the transient calculation, the junction is found choked due to numerical oscillations (particularly at the first time step). Once the junction is choked, the unchoking test is not sufficient to unchoke the junction. As a result, the junction remains choked as steady state is approached, even though the junction velocity is less than the subcooled critical velocity. This explanation forms the basis for a proposed modification, to be discussed in Sec. IV, to eliminate the discontinuity. The discontinuity is therefore seen as a result of the combined effect of code logic, the manner in which steady-state solution is reached (which depends on specification of boundary conditions), and user input data (i.e., the Z parameter). In particular, unphysical results may arise if inconsistent user input K factors, flow areas, and discharge coefficients are specified. For instance, if Z is specified such that Pdow„ is negative, then our analysis shows that the flow will never be choked. This is unphysical and is because of the inappropriate user input data. IV. PROPOSED MODIFICATION A modification has been suggested in Ref. 1. Here, we suggest another corrective measure to eliminate the discontinuity for any value of Z, based on the present analytical study. Since the discontinuity may not exist if realistic input data are used for the break junction, this modification is suggested for the sake of numerical smoothness. The root cause for the discontinuity is related to the inability of the unchoking test to unchoke the junction, which renders the flag CHOKE in subroutine JCHOKE being always true once it is set. In the subcooled choking test logic, the junction is choked if either the flag CHOKE is true or the junction velocity Vj is greater than or equal to the critical velocity Vc. The modification therefore removes the test on the flag CHOKE in the logic and bases the choking decision solely on Vj > Vc. Thus, the condition Vj > Vc serves as both choking and unchoking test for the junction. With this modification, the code will behave exactly the same way as the simple analysis of Sec. Ill indicates. Results for several cases of Z > 1 have been obtained using the modified code version. Figure 6 shows the steadystate flow rate against the downstream pressure. No discontinuity is calculated. In addition, the calculated flow is never choked for the entire range of Pdown. The impact of this modification on the code prediction of separate effect and integral tests must be investigated before adaptation. Downloaded by [University of Florida] at 21:41 25 October 2017 1.5 PdowrVPsat Fig. 6. Variation of steady-state flow rate against downstream pressure with code modification. V. CONCLUSIONS REFERENCES A root cause analysis of an observed discontinuity in the subcooled choked/unchoked flow has been presented. It has been found that the user input junction loss coefficient and discharge coefficient affect the occurrence of the discontinuity. The root cause is traced to the inability of the code to unchoke a junction even though the junction velocity is well below the critical velocity for the specified flow conditions. By modifying the choking test to be based only on the relative magnitudes of the junction and critical velocity, the discontinuity is eliminated. Recently, code developers at Idaho National Engineering Laboratory have completed modifying the unchoking tests to eliminate this anomaly. 5 The corresponding code updates will be available in the next code release, Version 3.2. 1. S. PETELIN, O. GORTNAR, and B. MAVKO, "RELAP5 Subcooled Critical Flow Model Verification," Trans. Am. Nucl. Soc., 69, 546 (1993). 2. K. E. CARLSON et al., "RELAP5/MOD3 Code Manual," NUREG/CR-5535, EGG-2596, Vol. 4, EG&G Idaho (June 1990). 3. W. L. WEAVER, "Improvements to the RELAP5/MOD3 Choking Model," EGG-EAST-8822, Idaho National Engineering Laboratory (Dec. 1989). 4. O. C. JONES, Jr., "Flashing Inception in Flowing Liquids," J. Heat Transfer, 102 (1980). 5. P. MURRAY, Idaho National Engineering Laboratory, Personal Communication (July 28, 1994). W. S. Yeung (BS, mechanical engineering, University of Lowell, 1976; PhD, mechanical engineering, University of California, 1979) is a principal engineer at Yankee Atomic Electric Company (YAEC). His current technical interest is in thermalhydraulic analysis using RELAP5/MOD3 and water hammer fluid transient analysis. J. Shirkov (MS, Moscow Institute of Energetics, 1987) is a nuclear engineer at YAEC. His current interests include computational methods in thermal hydraulics and heat exchanger analysis.