ML12145A438

Licensee Handout on Calibration and Benchmarking of Single and Two-Phase Jet Cfd Models

ML12145A438
Person / Time
Site:	South Texas
Issue date:	05/18/2012
From:	Ballew D, Gurecky W, Schnieder E University of Texas at Austin
To:	Balwant Singal Plant Licensing Branch IV
Singal B
Shared Package
ML12145A232	List: ML12145A312 ML12145A438 ML12145A466
References
TAC ME7735, TAC ME7736
Download: ML12145A438 (40)

Text

THE UNIVERSITY OF TEXAS AT AUSTIN Calibration and Benchmarking of Single-and Two-Phase Jet CFD Models May 18, 2012 Erich Schnieder, Davis Ballew, William Gurecky 5/18/2012

Table of Contents Introduction .................................................................................................................................................. 4

Background:

Jet Phenomena and Terminology ........................................................................................ 4 Critical Flow ........................................................................................................................................... 6 Compressible Effects ............................................................................................................................. 7 Modeling ....................................................................................................................................................... 9 Mesh ..................................................................................................................................................... 9 Solver .................................................................................................................................................... 9 Turbulence .......................................................................................................................................... 10 Single Phase Modeling ............................................................................................................................ 10 Boundary Conditions........................................................................................................................... 11 Single Phase Solution Methods and Controls ..................................................................................... 11 Two Phase Modeling ............................................................................................................................... 11 Mixture Model .................................................................................................................................... 11 Critical Mass Flux ................................................................................................................................ 12 Benchmark Results...................................................................................................................................... 14 Single Phase Benchmarks........................................................................................................................ 14 Mesh ................................................................................................................................................... 14 Results ................................................................................................................................................. 15 Sensitivity to Turbulent Intensity Rate .............................................. Error! Bookmark not defined.23 Sources of Error ................................................................................. Error! Bookmark not defined.23 Future Single Phase Scope .................................................................................................................. 24 References .......................................................................................................................................... 24 Multiphase Benchmarks ......................................................................................................................... 24 ANSI / Kastner Benchmark .................................................................................................................. 25 NPARC Benchmark .............................................................................................................................. 31 NED13 (C.F. Forrest) Benchmark ........................................................................................................ 32 Discussion............................................................................................................................................ 34 Computation Time ............................................................................. Error! Bookmark not defined.34 Critical Flow Model Benchmark and discussion.................................................................................. 35

References for Multiphase Section ..................................................................................................... 37 Conclusion and Future Scope...................................................................................................................... 39

Introduction In the event of a Loss of Cooling Accident (LOCA) triggered by a pipe rupture, the resultant high pressure two phase jet emanating from the break is energetic enough to cause damage to surrounding components and insulation. The debris generated from the jet impingement on nearby structures can potentially cripple the emergency containment spray system by the clogging sump screens located on the containment floor. Understanding the damaging potential of the jet via numerical and experimental study will aid in quantifying the expected damage seen from such an event.

The Fluent jet model in development is intended to fill gaps in experimental data. The results will be used to improve the fidelity of damaging Zone of Influence (ZOI) predictions in the currently-used ANSI/ANS standard [2]. Pressure and velocity fields will be had from the simulation results.

Compressible flow and turbulent features may have a strong impact on the dissipative qualities of the jet. The ability to model such phenomena is a key advantage of the CFD approach over the simplified treatment found in the ANSI model.

To assess the validity of Fluent as the simulation platform over the flow domain of interest, initial results will be compared to measurements and previously published model results. Measured data (e.g., pitot and impingement pressures) is only available for a small fraction of the possible upstream conditions and impingement geometries seen during a LOCA in a PWR or BWR. The ultimate objective of the CFD simulations is thus to provide pressure, temperature and velocity field calculations that, once calibrated against the sparse experimental data, will allow ZOI extents to be correlated to upstream and geometric conditions.

The available experimental data primarily focuses on impingement pressure measurements with relatively few data sets from pitot tube stagnation pressure measurements. This document focuses on the latter data set for single- and two-phase freely-expanding jets. Subsequent work will expand the validation exercise to include impinging jets.

Background:

Jet Phenomena and Terminology Shown in Figure 1, the ANSI jet model assumes the existence of a water core defined as a 45 degree liquid cone emitted from the inlet. The asymptotic plane is explained to be the last axial location where the jet does not interact with the environment and contains the same momentum flux as at the break plane. In addition the jet static pressure is at ambient pressure at the asymptotic plane.

Thermodynamic properties are calculated and used to determine steam quality. Jet density and expansion rates can then be obtained. A momentum balance approach allows one to obtain a thrust coefficient for the jet and provide accurate results for forces reacted out by the supporting structure [2].

Fiigure 1. Key ANSI terminology Issues arise when using the ANSI mod del alone to predict the destructive ZOI for the jet. Its derivation does not account for interp phase slip. Instead, a smeared phase approach is taaken to simplify the description of the jet beh havior. This model is not capable of resolving any ccompressible flow structures that may act to bring the t jet out of the damaging regime more quickly th han the ANSI model predicts. The bounding jet geo ometry is largely based on photographic evidence. The true jet structure is obfuscated by condensed d water droplets, making such evidence difficult to u utilize -

especially for a physical jet structure estimate. Furthermore, dissipative viscous turbuleent flow effects are not considered. The autho ors were conscious of the conservatism of their assu umptions in and justified them as necessary given n the simplifying modeling approach they pursued.

Figure e 2. Multiphase free jet illustration.

Figure 2 depicts the geometry of a flashing jet undergoing choked flow at its outlet, while Table 1 defines the symbols appearing in the figure. The key features of the jet can be summarized as follows. As the subcooled liquid crosses the exit plane, a sudden pressure drop induces violent phase change. The jet initially expands in the radial direction. A water core is formed jutting out just beyond the exit plane where pressures remain near in-pipe conditions. Vapor nucleation rapidly occurs just outside the water core. The mixture accelerates to satisfy continuity as its density falls. Further from the break plane, the jet continues to expand in the radial direction but less rapidly. Frictional shear forces slow the jet at large L/Ds and ultimately bring the jet out of the destructive regime.

The subcooled liquid core length is a function of the local sound speed ( , m/s), the speed v (m/s) at which the fluid is traveling at the break plane, and the pipe's inner diameter D (m) [9].

(1)

The formation of a liquid core is a consequence of sonic flow at the break plane. As pressure information can only travel at the local sound speed, the fluid emanating from the break plane near the pipe wall receives the ambient pressure information in less time than the centerline flow. In that time, the centerline flow travels further downstream, in this case beyond the break plane, without receiving downstream pressure information.

Table 1. Frequently referenced terms.

Symbol Unit Description

[Kg/s] Critical mass flow rate

[Kg/m^2*s] Critical mass flux

, [MPa], [C] or [K] Upstream pressure and temperature

[MPa] Break plane (throat) static pressure

- Critical pressure ratio

[m] Liquid core axial extent r [m] Radial distance from jet centerline L [m] Axial distance from break plane Critical Flow Across the flow regime of study, upstream conditions exist which give rise to flashing choked flow of an initially subcooled or saturated liquid. Similar phenomena must considered in the design of high pressure relief valves and in the design of flash evaporative distillation equipment, as such, a substantial amount of literature exists on the subject [1]. Single phase methods for determining the choked flow mass flux cannot be applied in the multiphase case - one may not necessarily assume a single phase exists at the break plane.

Critical flow models can be broken into two categories: those that incorporate non-equilibrium effects and those that do not. The non-equilibrium flow models attempt to predict vapor nucleation rates and interphase slip effects near the break location. These models are dependent on the pipe length, pipe roughness and break geometry to be used to their full potential. They have the benefit

of being significantly less computationally intensive than multiphase CFD simulation of the section of piping leading to the break.

Equilibrium flow models take a simplified thermodynamic approach to the critical flow problem. In many cases, the flow is assumed to be isentropic. The critical mass flux can be found by searching for a throat pressure and temperature, given upstream conditions, which force the average flow velocity at the choking location to reach the local speed of sound. All thermodynamic quantities are calculated at equilibrium conditions. Phase change kinetics are not considered or greatly simplified in many critical flow models. The goal is to consolidate the myriad physical processes which govern choked flow rate of an initially subcooled or saturated liquid discharge.

Obtaining a value for the critical mass flux is necessary if a mass flow rate inlet boundary condition is to be used for numerical simulation of the jet. Choosing and validating a critical flow model in the flow regime is an important task if this boundary condition will be employed in the CFD simulation. A brief overview of the critical flow models utilized in the simulations presented in this report is given in the Critical Mass Flux section on pg. 12.

Compressible Effects Figure 3 Expected flow geometry for underexpanded supersonic flow. Taken from Ref [4].

The free jets in this study develop the shock diamond pattern seen in Figure 3. A series of expansion fans and oblique shockwaves act upon the jet as it passes from the break plane to ambient conditions. Prandtl-Meyer and oblique shock relations can be used to determine the pressure, Mach number, and the flow angle relative to centerline at each state in the expansion-compression cycle. The behavior of a multiphase jet may approach the behavior of a single phase steam jet after a substantial amount of the liquid phase has transferred to vapor.

In the single phase case, upon exit from the nozzle, the high-pressure flow begins to expand into its surroundings as it crosses the expansion fans created by the supersonic flow. However, the flow also creates compression fans, which direct the flow back inwards toward the axis of flow. This compresses the flow, creating once again the underexpanded scenario described upon exit from the nozzle, leading to a continuing cycle of expansion and compression. Upon crossing the initial shock, the static pressure increases, while the total pressure and velocity drop very quickly. The thinness of the shock region means that the large gradients of total pressure and velocity occur over a very small distance, leading to numerical challenges when converging solutions near these regions. In an idealized situation, this expansion and compression would continue forever. However, frictional losses and turbulence cause the flow to lose significant amounts of pressure after each shock, until the barrel structure disappears completely as energy dissipates into the surrounding air.

Modeling Mesh The mesh was developed using ANSYS meshing program built into the Workbench productivity suite. A structured mesh was imposed so that we could finely control mesh density. It is possible to refine the mesh mid-run along high-pressure gradients. This is similar to Zerkle's [10] mesh adaptation scheme which increased the resolution of the standing shocks present in a supersonic steam jet. The mesh moves from a high density at the break plane to a relatively lower density moving axially and radially away from the jet core.

Figure 4. Meshing adaptation example. The blue regions indicate areas where the density of the mesh will be reduced. The Red regions indicate areas of step pressure gradients where the mesh density will be increased.

Solver The pressure based coupled solver was chosen as the overarching solver method. The coupled solver solves the momentum and pressure based continuity equations simultaneously [Fluent help]. This solver may have some advantages over sequential methods depending on the case. It is known to reduce convergence times especially in high velocity flow scenarios. This is an implicit method where the equations are solved via an iterative technique. This method is not compatible with the Eulerian multiphase model. Future simulations employing the Eulerian model will need to consider an uncoupled solution technique.

In the cases under study, a steady state solution is to be obtained for a boundary value problem governed by a system of nonlinear, coupled PDEs. The pseudo transient solution method recently implemented in Fluent (v13) is employed throughout the duration of the calculations. This solver method is especially beneficial when attempting to maintain stability of the numerical solution when sharp gradients are present in one or more field variables.

The pseudo transient method is a form of implicit under-relaxation for steady-state cases. Under relaxation is analogous to a negative feedback loop. Under relaxation incorporates a portion of the

previous iterations results into the current iteration's calculation to dampen numerical oscillation.

The maximum possible pseudo time step is governed by the stiffness of the system. A system is said to be stiff when one or more components of the solution respond on vastly different time or length scales than the others. The resulting algebraic system, written in matrix form, is difficult to invert.

To ensure stability, small pseudo time steps and a fine grid are essential. Fluent is capable of detecting divergence and will automatically reduce the pseudo time step when necessary. In some cases, however, the first iterations involve rapid divergence of temperature, pressure, and velocity fields. When this occurs it is possible to manually set the pseudo time step to a suitably small value to stabilize the initial solution. Once stabilized it is typically possible to allow Fluent to dynamically determine the pseudo time step necessary to maintain a stable solution.

Turbulence Turbulent mixing occurs at the interface between the quickly moving primary jet flow and the stagnant ambient surroundings. This region is characterized by highly rotational flow pockets where the average curl of the velocity field is extraordinarily high. The frictional shear present in this region gradually dissipates the jet's kinetic energy. A transient solution procedure coupled with a very fine mesh would allow one to gain more insight into the turbulence phenomena at hand.

The realizable K-Epsilon turbulence model is currently used due its ease of implementation. The turbulent viscosity ratio limit is left as default in the solution controls. The K-Epsilon model may not be the most suitable choice for the flow regimes arising within the jet. For example, it is known not to perform well when applied with adverse static pressure gradients, which do arise at locations in the jet. Other turbulence models may be better suited for modeling high sheer rates and may resolve small-scale turbulent effects despite a relatively large cell volume compared to these structures.

The solution is initialized with Fluent's hybrid initialization feature which solves for the potential flow field. This accelerates solution convergence versus initializing the flow field from inlet conditions alone.

Single Phase Modeling Although Fluent documentation recommends a density-based solver in problems involving compressible flow, the density-based solver is only applicable to single-phase flow, meaning it is not suited for use on the two-phase problems of interest later. Therefore, a pressure-based solver was used. All calculations were done using the steady solver to simulate steady-state and steady-flow conditions that dont take the passage of time into account across a control volume. An axisymmetric geometry was used in order to describe the flow on a plane crossing the centerline axis of the pipe. The energy model was turned on to account for the temperature and kinetic energy of the flow, as was a viscous model, to account for turbulence. Because Fluent documentation recommends the k- turbulence model as the first choice for scoping analyses, results were generated using this model. The k- turbulence model is recommended due to its robustness with both wall-bounded internal flow and free-shear flows. However, this model has only proven to give accurate results in cases involving small pressure gradients. Finally, water

vapor and air with their respective default properties were used in their respective benchmarking cases. To simplify the treatment of superheated steam and air in choked flow conditions, the fluids were modeled as an ideal gas, as they were in [2].

Boundary Conditions Inlet conditions for each simulation were taken from their respective benchmarking conditions specified in Table 1. A turbulent intensity and length method was used to model the turbulence at the inlet, using a length of 0.021 meters and a turbulent intensity of 2% as determined by the equation for turbulent intensity provided by FLUENT. The outlet zones were set at atmospheric pressure and room temperature.

Single Phase Solution Methods and Controls A coupled solver was used so that all of the governing equations were solved simultaneously, as opposed to the segregated solvers that solve the equations sequentially. This decision was supported by Fluent documentation, which recommends the coupled solver for flows involving compressible effects. A first order upwind spatial interpolation scheme offers stability as opposed to the more accurate but less stable second and third order interpolation schemes. As such, a first order interpolation scheme was selected in the initial stages of the calculations, and a second order interpolation was adopted after a sufficient number of iterations had been performed to ensure a converging solution.

Two Phase Modeling Mixture Model Fluent's mixture multiphase model with evaporation-condensation mass transfer is used to determine liquid to vapor phase change rates. This model can be used to initialize a more complex Eulerian multiphase run. The immediate goal is to produce realistic multiphase flow using the mixture model. The mixture model simulates n phases that can move at different velocities in the flow using the concept of slip velocities. One gas phase may be defined as compressible. The mixture model solves the continuity equation for the mixture, the momentum equation for the mixture, the energy equation for the mixture, the volume fraction equation for the secondary phases, and equations describing the relative velocities between phases.

Fluent's evaporation-condensation model handles liquid to vapor mass transfer. The governing mass transfer equation for this model when is displayed in Error! Reference source not found.Equation 2.

(2) has units of [kg/s/m^3]

is volume liquid phase fraction

is liquid density is a tunable relaxation constant can be found with knowledge of latent heat, accommodation coefficient (describes portion of vapor absorbed by liquid interface), and vapor bubble diameter. The determination of vapor bubble diameter can be difficult, as this value is dependent on the temperature, pressure, and velocity field. Thus, must be fine tuned to fit experimental data.

Classical approaches to vapor nucleation may not be adequate when attempting to model flashing flow. More robust equations of state can be included in an overarching mass transfer model. Van-Der Wall's EOS (1894) and later Redlich-Kwong EOS (1949) were developed to incorporate different modes of molecular energy storage into the equation of state. Rather than assuming the

``billiard ball model" indicative of an ideal gas these equations of state can approximate specific heat nonlinearity due to vibration and angular momentum molecular energy storage modes in addition to the primary kinetic mode.

Piecewise curve fits are drawn through steam table data and incorporated into user defined C functions that can be interpreted or compiled by Fluent at runtime. The functions return saturation temperature given pressure. These UDFs are supplied to the evaporation-condensation mass transfer method to obtain accurate results when compared to the default constant saturation temperature provided by Fluent. It is unclear why more advanced saturation temperature prediction methods are not available in this software package by default.

Critical Mass Flux The multiphase simulations presented in this paper utilize a mass flux inlet boundary condition.

Setting the inlet boundary condition to a specified mass flux ensures the phenomena of choked flow can be captured. In previous calculations a pressure inlet boundary condition was specified, however, FLUENT cannot simulate choking phenomena correctly without first simulating a section of pipe upstream of the break location. In critical flow upstream information drives the mass flow rate. The task of finding the critical mass flow rate given upstream conditions is delegated to the critical flow models in the current cases. In any case, the choked flow mass flow rate represents an upper discharge limit denoted as the critical mass flow rate ( ). The value of the critical mass flow rate given an initial upstream temperature and pressure must be applied as the inlet mass flux boundary condition.

The determination of the critical mass flux ( ) in the event of a pipe rupture has traditionally been handed by either the HEM or the Henry - Fauske (HF) model, each being slightly more accurate depending on the flow regime. The HEM and HF models were suggested for the ANSI model [2].

Currently, Leung's - method [6] is employed to calculate the critical mass flux for all subcooled initial conditions. At saturated conditions the short pipe lossless Henry (1970) [5] model is used.

The two models effectively handle all upstream conditions that may be encountered in BWR and PWR coolant discharge flows.

2 (3)

When highly subcooled initial conditions exist, the critical mass flux used as the inlet boundary condition was determined by the Bernoulli equation shown in Error! Reference source not found.Equation 3, recommended by Leung for initially subcooled liquid flow. This equation gives essentially the same result for the critical mass flux as the homogeneous equilibrium model for highly subcooled conditions.

When approaching saturation, Leung's - method transitions from a simple Bernoulli equation to a transcendental equation solved via iterative techniques for the critical mass flux. The correlation attempts to simplify two phase flow phenomena that develop in the pipe when the fluid is very near saturation. This equation predicts higher mass flow rates than the Bernoulli equation predicts at the same conditions, however critical mass flow rates calculated by this relation are still significantly lower than shown in some experimental data sets. [5, 7]

In addition to providing the critical mass flow rate, the critical flow models also provide the critical pressure ratio ( ) defined as the ratio between the break plane static pressure and the upstream stagnation pressure. The throat static pressure can be used to determine the saturation temperature of the fluid at the break plane. The result is not insignificant because it provides a sanity check for Fluents simulated break static pressure and mixture density.

Benchmark Results Single Phase Benchmarks The goal of modeling single phase supersonic free jet expansions is to replicate results from benchmarks on a simpler level than that of two-phase flow, allowing us to define modeling choices that also apply to the more complex two-phase flow simulation.

A list of benchmark sources and results is provided in Table 2. First, the methodology used in developing the model will be described in detail, noting the many user inputs that were selected in setting up and initializing the simulation. Then, the results gathered via these simulations will be compared to those of several benchmarks.

Table 2. Single Phase Benchmarks Benchmark Case Inlet Gauge Total Inlet Total Temperature Benchmarks Available Pressure CDI - [1] 4.39 MPa 949 K

• Centerline static Steam pressure plot - Fig. 9

• Centerline total pressure plot - Fig. 10 Zerkle - [2] 4.39 MPa 949 K

• Static pressure Steam contours - Fig. 3 & 4

• Total pressure contours - Fig. 5 & 6

• Centerline static pressure plot - Fig. 9

• Centerline total pressure plot - Fig. 10 Zerkle - [2] 4.26MPa 894 K

• Total pressure Air contours - Fig. 7 & 8 Masuda - [3] 4.02 MPa 548 K

• Radial Pitot pressure Steam profiles - Fig. 11-14 Mesh A mesh was created to describe the geometry of the flow region using the ANSYS Workbench suites meshing tools. All meshes were designed to match the geometries of the benchmarks, extending to 20 L/D in the x-direction and approximately 5 L/D in the y-direction. A bias is introduced in the y-direction to create a higher resolution among cells where large gradients are expected, specifically regions in close proximity to the break plane. Results in [2] imply that a quadrilateral mesh, while conducive to converging solutions, will result in a less accurate description of the flow, so a triangular mesh was used. Additionally, adaptive meshing was used to create a higher resolution in those cells with larger gradients. The mesh used is shown in Figure 5, in which the pipe break

plane is visible in the bottom left of the mesh. Important values pertaining to the mesh are shown in Table 3. Note that orthogonal quality values can range from 0 to 1, where values close to 0 correspond to low quality.

Table 3. Mesh Parameters.

Type Triangular Cells 27750 Minimum Orthogonal Quality .575 Figure 5 Mesh used in modeling single phase jet.

Results In [2], Zerkle aimed to correct results obtained in [1]. Simulations in [2] were run using version 5 of Fluent, while results in [1] were obtained using the NPARC code. Zerkle noted that results in [1]

show large, unexplained differences in total pressure downstream of the break plane for air and steam, and aimed to either correct this or explain the discrepancy. However, Zerkles results show that using Fluent to simulate this flow may give rise to nonphysical features in the flow due to the use of the control volume solver.

All results shown were carried out at the conditions for their respective benchmarking cases stated above. All contours and plots were designed so that their format would match those provided in benchmark cases. Initial conditions were set forth in [1] such that the steam and air would have similar Reynolds numbers [1]. The inlet was chosen to be a round hole with a diameter of 299 mm.

Static pressure fields shown in Figure 6 show that current results show a strong qualitatiive resemblance to results from [2]. The classic shock diamond structure shown in the resu ults of [2] is distinctly evident, matching the expeccted geometry shown in Figure 3. Similar observattions can be made from Figures 6-8 with both steaam and air. However, current contours fail to matcch the shock locations shown in [2], an issue that willw be addressed in the following section.

Figure 6. Comparison of current static prressure field (top) for water vapor with inlet conditions P P=4.39MPa; T=949K with [2] (bottom).

Figure 7. Comparison of current total prressure field (top) for water vapor with inlet conditions P=

=4.39 MPa; T=949K with [2] (bottom).

Figure 8. Comparison of current total presssure field (top) for air with inlet conditions P=4.26 MPa; T T=894K with

[2] (bottom).

Current centerline static and total preessures are compared to those of both [1] and [2] in n Figure 9 and Figure 10. Shock front locations in i current results generally correspond with the lo ocations in

[1]. However, the shock occurs appro oximately 2 L/D upstream from that of [2]. Meanw while, the

recovery of both static and total pressure after the initial shock is closer in magnitude to that of [2]

than [1]. Similar results can be drawn from Figure 8. This figure also offers information about the total pressure along the centerline prior to crossing the first shock boundary. While the total pressure in [1] tends to dissipate gradually, current total pressures and total pressures in [2] are mostly sustained from the break plane to the first shock boundary. Given that the results from [1]

and [2] were generated using the same boundary conditions, it is difficult to determine the accuracy of current results with respect to prior attempts to model the flow computationally. Because the modeling choices of [1] and [2] are not known, a root cause analysis of the differences in these results is not possible. Additionally, given the sparseness of experimental data for single phase freely expanding jets with similar upstream conditions, it is difficult to determine the accuracy of current results on a general level. However, it is clear that the current results qualitatively match certain aspects of those of both [1] and [2], indicating some degree of correlation between results.

0.7 0.6 0.5 0.4 Zerkle p/p0 0.3 CDI 0.2 Current 0.1 0

0 2 4 6 8 10 12

-0.1 x/D Figure 9. Comparison of current centerline static pressure plot with those of CDI [1] and Zerkle [2].

700 600 500 400 Psia Zerkle 300 CDI Current 200 100 0

0 5 10 15 20 x/D Figure 10. Comparison of current centerline total pressure plot with those of CDI [1] and Zerkle [2].

Current results of total pressure as a function of distance from the axis of rotation at various distances from the inlet are compared to experimentally obtained results from [3]. The pressures in [3] were determined with pitot tubes placed at different radial distances from the centerline of the flow at several distances from the break plane. Current results match the shape of those from

[3], capturing the dip in pressure in the core of the flow relative to the outer regions caused by the shocks. The dissipating viscous effects are also clear here, as the shock diamond structure breaks down with increasing distance from the outlet.

The pitot pressure measurements of [3] were obtained from pitot tubes oriented parallel to the jet centerline. The stagnation pressure ready by the tubes would therefore be expected to be lower than the total pressure predicted by Fluent, as the flow in the jet is not oriented normally to the pitot tube inlets except at the centerline. Even so, current results appear to be too conservative with respect to experimentally obtained results, as they indicate that the current flow maintains a high total pressure much further downstream than in [3].

0.6 L/D From Inlet 0.8 0.7 0.6 0.5 P/P0 0.4 Current 0.3 From [3]

0.2 0.1 0

0 1 2 3 4 5 6 Distance from Axis (r/D)

Figure 11. Comparison of radial total pressure plot at 0.6 L/D away from the inlet with [3]

3.3 L/D From Inlet 0.25 0.2 0.15 P/P0 0.1 Current From [3]

0.05 0

0 1 2 3 4 5 6 Distance from Inlet (L/D)

Figure 12. Comparison of radial total pressure plot at 3.3 L/D away from the inlet with [3]

4.3 L/D From Inlet 0.25 0.2 0.15 P/P0 0.1 Current From [3]

0.05 0

0 1 2 3 4 5 6 Distance from Inlet (L/D)

Figure 13. Comparison of radial total pressure plot at 4.3 L/D away from the inlet with [3]

13 L/D From Inlet 0.1 0.09 0.08 0.07 0.06 P/P0 0.05 0.04 Current 0.03 From [3]

0.02 0.01 0

0 1 2 3 4 5 6 Distance from Inlet (L/D)

Figure 14. Comparison of radial total pressure plot at 13 L/D away from the inlet with [3]

Figure 15 shows current total pressure plots for water vapor at the boundary condition specified by

[1] and [2], in which the turbulent intensity rate was varied from 2% (the value used in benchmarks) to 50%.It is clear that changing the level of turbulent intensity affects the amount of total pressure dissipated prior to shock, the amount of recovery after the shock, and the location of downstream shocks.

5.0E+06 4.0E+06 Total Pressure (Mpa) 3.0E+06 2% Intensity 5% Intensity 2.0E+06 10% Intensity 25% Intensity 50% Intensity 1.0E+06 0.0E+00 0 1 2 3 4 5 6 7 Distance From Inlet (m)

Figure 15. Effect of Turbulent Intensity Rate on Total Pressure along Centerline Discrepancies between individual results can arise from meshing choices: both polyhedron shape and cell number must be studied further to determine their effects on solutions. As mentioned before, [1] and [2] used different methodologies in determining centerline pressures. The NPARC code was used in [1], which is a time-marching solver that solves for time-dependent variables in conservation equations. Meanwhile, Fluents control-volume based solver was used in both [2] and current simulations. This fundamental difference in approaching the problem could be a significant source of error in results. Finally, Masuka used experimental results with pitot pressure tubes placed within the flow to collect data in [3], a practice that could lead to significant changes in what would otherwise be a free-jet. This placement of a physical object within the flow could account for the low pitot pressures obtained in [3] relative to current results. The orientation of the pitot tube relative to the flow direction impacts the pressure read from the device. Hence Fluent pressure field results do not necessarily correlate one-to-one with pitot pressure data.

Corrective measures should be taken if purely stagnation or static pressure readings are desired from a pitot tube not precisely parallel to the velocity vector at the point of interest. The tangential and axial components of the velocity are necessary to estimate the effect misalignment of the pitot tube has on the recorded pressures. A relation between the static and stagnation pressure for

compressible flow is used in conjunction with the measured pressures, misalignment angle and velocity field data to yield purely static or only the stagnation pressure. This corrective process is only possible if adequate field data is available.

1 (4)

The stagnation pressure ( ) of a compressible fluid is related to the static pressure via Error!

Reference source not found.Equation 4. The stagnation pressure for a compressible fluid is the static pressure the fluid would retain if brought to rest isentropically.

Future Single Phase Scope The methodology developed for the single phase simulations is being applied to models involving two phase flow. However, benchmarking of single-phase flow will continue until a higher degree of agreement is achieved against accepted results, specifically, the location of shocks and the small pressure drops in current results when compared to available experimental data.

Future work with single phase flow will include a mesh sensitivity sweep to determine the sensitivity of field results to mesh shapes and face orientation. Additional model parameter studies, such as the one performed with turbulent intensity in section 3.1.4, will be performed to isolate the effects of individual input parameters on results. An upstream conditions sweep will also be performed to generate a wide range of results over various upstream temperatures and pressures in an effort to determine a semi-empirical relationship between flow properties and input conditions.

References

[1] D. Zerkle, NRC Update #1: Prediction of Supersonic Expansions Relevant to Debris Generations, LA-UR-00-2247.

[2] M. E. Teske, A. H. Boschitsch, T.B. Curbishly, "Zone of Influence as Defined by Computational Fluid Dynamics", C.D.I. Report No. 96-0 l (in "Utility Resolution Guide for ECCS Suction Strainer Blockage, Vol3", NED0-32686-A).

[3] F. Masuda, T. Nakatogawa, K. Kawanishi, M. Isono, Experimental Study on an Impingement High-Pressure Steam Jet, Nuclear Engineering and Design 67 (1981) 273-286.

[4] Scott, J. (2005, April 17). Shock diamonds and mach disks. Retrieved from http://www.aerospaceweb.org/question/propulsion/q0224.shtml Multiphase Benchmarks

Table 4 presents the benchmark data sources considered in this section. The ANSI jet overview presented in Appendix I [2] provides a suite of impingement jet data. The ANSI pressure results presented in the ANSI report are obtained by isentropically bringing the flow to rest on a fictitious impingement surface. The ANSI result considers an infinite series of impingement surfaces on which impingement pressure distributions are taken rather than singular point pressure measurements, comparing ANSI pressure results to Fluent field data can be misleading. We do not expect to see good agreement, though this does not discount our simulation results, differences simply highlight the conservative and unphysical nature of the ANSI model.

Table 4. Benchmark overview.

Study Type Upstream Stagnation Subcool Benchmarks Pressure [MPa] [C] Available ANSI Empirical and 7.19 20 Impinging pressure thermodynamic and quality contours.

equilibrium Centerline pressure.

Kastner Empirical 7.19 20 Centerline pressure.

NPARC Simulation 7.5 Unk. Centerline Stagnation pressure.

13NED Experimental 6.3 30 Centerline Static (C.F. pressure.

Forrest)

The NPARC simulation results were obtained from reference [9]. The lack of information regarding upstream temperature conditions make this data set difficult to benchmark against. C. F. Forrest conducted a series of multiphase blowdown experiments. Pitot tubes were used to record the static pressure at various L/D (ref [3]).

ANSI / Kastner Benchmark Boundary conditions and solver settings for this case are given in Tables 5 through 7. The upstream stagnation pressure and temperature are given to be 1044 [psi] and 516F. These conditions represent an 18K subcool. , the primary user input to the evaporation-condensation phase change model, was set by guess and check methods constrained by convergence. A 2.5 currently causes convergence issues. A supplemental or entirely new phase change model will likely produce better results for flashing phase change rates.

Table 5. Boundary Conditions.

Po [MPa] To [K] Critical Mass flux

[Kg/(m^2*s)]

Kastner 7.198 542 ANSI 7.198 542 Fluent - 542 52676

Table 6. Solver Settings.

Setting Detailed Settingg Space Axis-Symmetric Time Steady Multiphase Mixture (Evap. - = 1.5 Cond)

Turbulence K-Epsilon realizable Solver Pseudo Transient Matrix Conditioner ILU Multigrid Settings Selective AMG solver RPM stabilized Momenttum and Energy Table 7. Mesh Parameters.

Type Quad Only Initial Cells 48000 Average Orthogonal Quality 1.0 Figure 16. Centerline sta agnation pressure comparison. Reproduced from [2].

The discrepancy seen in the initial staagnation pressure as compared to the ANSI results in Figure 16 arises from specifying a mass flow inllet boundary rather than a pressure inlet. The consservative nature of the ANSI model is shown by y the large assumed constant pressure region exten nding to 2 L/D - indicating a large distended liqu uid core length. The Fluent simulation result closeely resembles the Kastner curve though retaining r a larger stagnation pressure at 2 - 5 L/D.

Figu ure 17. Centerline static pressure.

Figu ure 18. Centerline Steam Quality.

Comparing Figure 18 to Figure 19 rev veals a large discrepancy in the quality predicted att the centerline. It is possible that Fluents evaporation-condensation mass transfer model is not capable of describing flashing liquid to vapor mass transfer accurately. Finding a substitute or complementary phase change model in the near future is a primary goal. If phase changge occurs at a faster rate the simulated jet geomettry as well as the location of the first shock will cha nge.

Figures 20 through 25 provide additio onal field data for this case.

Figure 19. ANSI steam quality contour. [2]

Figure 20.. ANSI stagnation pressure contour. [2]

Figure 21. The local sound speed in the co ompressible vapor phase. As the density of the jet falls th he local speed of sound decreases.

Figure 22. The local Mach number of the mixture. Rapid expansion near the break plane contrributes to extraordinarily high velocities surroundin ng the jet core. The Mach number at the break plane is very near unity, indicatiing our choked flow condition is met.

Figure 23. Phase change mass transfer ratte. The largest mass transfer rates can be seen near the in nner radius of the pipe break.

Figure 24. Turbulent kinetic energy [m^

^2/s^2]. The jet dissipates energy as it mixes with the env vironment at large L/Ds.

Figure 25. Static prressure field exhibiting shock diamond pattern.

F Figure 26. Total pressure field.

NPARC Benchmark The boundary conditions used to prod duce the following results were identical to those u used in the ANSI/Kastner benchmark above (Tab ble 5).

Figure 27. Stagnattion pressure centerline plot from reference [9].

The NPARC code result presented in figure f Figure 27 predicts a region of sustained staggnation pressure far into the open environmeent. Temperature boundary conditions are not fullyy specified in the NPARC data source so it is unclear if a high degree of subcooling or a fault in the computational model is at issue.

Pressure recovery via the reflected co ompression fan present after the first expansion is well defined. The location of the first shocck predicted by Fluent is fairly consistent with the N NPARC code results. The current Fluent results dissplay non-physical smoothing of pressure gradientts leading to a delay in the formation of the first shhock. The NPARC code does not predict a strong sh hock effects past the first pressure depressed locaation where the Fluent results indicate that a shock k diamond pattern will develop. A shock diamon nd pattern may aid in dissipating the jet energy whiile also acting to refocus the jet downstream the initial expansion allowing for some recovering of the jet's damaging potential intermittently. One would expect this pattern if the flow primarily i s in the vapor state and is at or approaching thet local sound speed.

NED13 (C.F. Forrest) Benchmark Benchmark results are not available at a this time. The calculation is currently in progresss. The boundary conditions and modeling op ptions are presented in Tables 8-10; the data beingg benchmarked is shown in Figure 28.

Table 8 : Boundary Co onditions for C.F. Forrest experimental benchmark.

Po [MPa] Subbcool [K] T @ break [K] Critical Mass flu ux

[Kg/(m^2*s)]

Experiment 6.3 30 -

Fluent - - 512 63783

Table 9. Solver settings.

Setting Detailed Settting Space Axiis-symmetric Time Steady Multiphase Mixxture (Evap. - Cond.) = 1.5 Turbulence K-E Epsilon realizeable Solver Pseeudo Transient Matrix Conditioner ILU U

Multigrid Settings Seleective AMG RPM Stabilized Momentum and Energy Table 10. Meshing Information.

Type Tri Only Initial Cells 44000 Average Orthogonal Quality 0.59 Figure 28. Experimenttal centerline static pressure. Reproduced from [3].

The steady state blowdown results ob btained by C.F. Forrest indicated that upstream preessure had little effect on the location of the first sub-ambient pressure region. The static pressure at the break plane is seen to be very near the liquid saturation pressure at the local break temperatu ure and pressure. We expect to see static presssures similar to the saturation pressure at the breeak plane in our Fluent simulations.

Discussion The smoothed pressure gradients seen in all the multiphase Fluent results are evidence of numerical diffusion. Adaptive meshing can reduce numerical diffusion in steep pressure gradient areas. A triangular mesh also will have a marked effect on reducing numerical diffusion and increasing shock resolution as seen in Zerkles [10] results. The overall convergence criteria may need to be set at a more demanding value. Currently, a solution is deemed to be converged when the L1 norm residual for the continuity equation reaches a valued of 1E-7.

In all the simulations presented in this report use a critical mass flux inlet boundary condition. In the December 2011 progress report, a (total) pressure inlet condition was used. When specifying an inlet pressure we found that the simulated mass flux was unrealistically large. In the current simulations we choose to enforce the choking phenomena at the break plane by using a critical flow model to determine the choked flow mass flow rate. Supplying this value to Fluent ensures that the choked condition is met. We have encountered difficulties matching the simulated static pressure at the break to the expected static pressure value at the break plane however. This error is also manifested in the simulated total pressure results. Therefore, currently, we may only match the inlet pressure or the critical mass flux accurately but not both at the same time. It was found that decreasing the incompressible liquid phase density had a significant increased the brake plane static pressure. A non-physical liquid phase density (<700 [kg/m^3]) can be supplied to drive the static pressure to acceptable values.

The multiphase model is currently under-predicting the rate of liquid to vapor mass transfer. The evaporation-condensation model provided by Fluent may not be able to faithfully capture flashing phase change rates. In a flashing jet the rate of liquid to vapor mass transfer is aided by turbulent break-up of the liquid droplets. When the liquid is in a superheated state, the driving force behind phase transformation is the amount of superheat. Once the fluid quickly releases excess thermal energy, through flashing vaporization, evaporation occurs at constant temperature. This constant temperature boiling is the thermodynamic region the evaporation-condensation model was designed to handle.

Fluent does not include a flashing phase change model under the blanketing mixture model by default, however, one may be implemented by a UDF. Fluent allows for more complex phase change models when using the Eularian Multiphase model. Investigating the Eularian multiphase model is a goal in the coming months.

A multiphase run with adaptive meshing schemes present (limiting the maximum number of cells to 2E5) and utilizing second order interpolation methods after satisfactory stability is seen, can take up to two weeks on a four core machine with 16 GB of ram. The total number of iterations typically necessary to converge a solution is on the order of 5E5. In the near future additional computational resources will be available totaling in 24 cores, each with a dedicated 1 GB of ram.

Figure 29. Fluent multiphase run r near completion. Total run time: approximately 1 wee ek.

Fluctuations seen in the residual plot leads one to believe that stability of the solution iss delicate even at a pseudo time step on the ord der of 1e-5 seconds. The spike in the residuals seen n at near the 150000 iteration mark indicates an adaptive meshing event. These events are scheduleed to happen at set intervals during the solution prrocess.

Critical Flow Model Benchmark and d discussion As the fluid nears saturation, experim mental evidence suggests that the critical mass flow w rate does not continue to degenerate as significcantly as is predicted by critical flow models shown n in Figure

30. Initially subcooled critical flow. Reproduced R from [7]. A pressure ratio of 1 indicatees saturated conditions.. Leung suggests that this could be due to nonequilibrium phase change and interphase slip velocity effects which are not takeen into account when treating the phases as a hom mogenized fluid.

Figure 30. Initially subcooled critical fllow. Reproduced from [7]. A pressure ratio of 1 indicatess saturated conditions.

A critical flow model developed by F.JJ. Moody was employed to predict the critical mass flux in a circular nozzle in the NED105 (Kastneer) case [7].

The current MATLAB code used to callculate the critical mass flux employs a different cr itical flow model depending on the upstream conditions. Leungs -model [6] is used for subcoolled conditions and Henrys 1970 critical flow f model [5] is employed for initially saturated c onditions.

A point of inflection can be seen near a pressure ratio of 0.9 where Leung's model transiitions from a modified Bernoulli equation at higherr levels of subcooling to a two phase critical flow co orrelation.

For a pressure ratio equivalent to uniity, Henry's HEM short pipe discharge model is utilized by the code. The thermodynamic data suppliied to the MATLAB code was obtained from referen nce [8].

Slight over prediction in the highly saaturated region can be seen when comparing the cu urrent code to past model results. For a majority of the flow regime plotted, at significant levels of subcooling, the Bernoulli equation (Error! Referen nce source not found.Equation 3) was used to calcu ulate the critical mass flow rate per Leungs reccommendation.

Figure 31. Initiallly saturated critical flow. Reproduced from [5].

Mentioned previously, for initially satturated conditions Henrys 1970 saturated short piipe discharge model is employed. Validattion of the current code is presented in Figure 31 aabove. The model does not include a frictional pip pe loss correction. An entrance loss coefficient, ho wever, is present in the Henry model. This parrticular model is said to be valid for a pipe which haas a length less than 13 break diameters [5].

References

[1] Adel K. El-Fiqi, N.H Ali, H.T. El-Desssouky, H.S. Faith, M.A. El-Hefni. Flash evaporation in a superheated water liquid jet. Desalina ation, 206:311-321, 2007.

[2] ANSI/ANS Jet Model. Appendix I, pages p I1-I52.

[3] C.F. Forrest, K.S. Shin. Measuremeents of impact loads and expansion of flashing wateer jets.

Nuclear Engineering and Design, 99:53-61, 1987.

[4] D. Tomasko G.G. Weigand, S.L. Tho ompson. Two-phase jet loads. Sandia National Lab boratories, NUREG/CR-2913(R4), 1983.

[5] Henry, R.E. The two-phase critical discharge of initially saturated or subcooled liqui d. Nuclear Science and Engineering, 41:336-342, 1970.

[6] M.A. Grolmes J.C. Leung. A generaal correlation for flashing choked flow of initially su ubcooled liquid. AIChE, 34(4):688-691, 1988.

[7] R. Rippel, W. Kastner. Jet impingem ment forces on structures. Nuclear Engineering and d Design, 105:269-284, 1988.

[8] P. Schmidt, D. Baker, O. Ezokoye, J. Howell. Thermodynamics an Integrated Learning System.

Wiley. 2006.

[9] Thomas B. Curbishley, Milton E. Teske, Alexander H. Boschitsh. Zone of influence as defined by computational fluid dynamics. Continuum Dynamics Inc., 96-01, 1996.

[10] Zerkle, NRC Update #1: Prediction of Supersonic Expansions Relevant to Debris Generations, LA-UR-00-2247.

Conclusion and Future Scope Benchmarking, calibration and model sensitivity studies will continue with the aim of achieving agreement against published numerical and experimental free jet results. Simulations of an impinging two phase jet must be conducted once the simpler free jet case is fully understood. Much of the experimental LOCA blowdown results are given by impingement pressure distributions on a flat plate or structure of interest. Validation against these results will ultimately determine the quality of Fluent two phase jet simulation predictions. A grid convergence study is to be conducted after accurate and consistent results can be obtained.

An effort in the near term to identify a suitable liquid to vapor mass transfer model will be made.

Currently the rate of mass transfer is significantly less than what the ANSI model predicts. Though the ANSI case is not an appropriate model to validate a phase change model against, phase change rates should result in a quality distribution similar to that seen in the ANSI quality contour.

Experimental density and/or quality results will be sought after in the literature.

The root cause of disagreement in break plane static pressures between must be explained and corrected. The density of the liquid phase, the inlet temperature and the rate of phase change at the liquid core interface all impact the local static pressure. The calculation of the static pressure is largely a thermodynamic problem. From a critical flow model the critical pressure ratio serves as an estimate to the expected pressure at the break plane. Fluent does not know what the break plane pressure should be when specifying a mass flow inlet, and specifying mass flow, temperature, density and pressure over constrains the problem. Proper settings in the phase change model and liquid density parameter should lead to correct static pressure results.

Achieving a solution with a minimal amount of numerical diffusion incurred along the way is an ongoing effort which will be aided by additional computational resources allowing for a finer mesh size (especially after refinement) and perhaps more memory and processor hungry stabilization schemes.

The ANSI model is quite conservative in its estimation of the ZOI volume at any given upstream condition. By building a more complete library of flashing jet behavior one can fill gaps in the understanding of the jets behavior left after experimental results are evaluated. With this knowledge it is possible to relax some of the conservative assumptions made in the development of the ANSI model. Decisions regarding relaxing these conservative constraints must be grounded in convincing evidence provided by simulated or experimental means. But experimental results do not provide pressure and velocity field data, only measurements at a handful of isolated locations.

The field data is necessary to predict the effect the jet will have on complex impingement surfaces.

Developing an improved ZOI model hinges on being able to accurately estimate the jets damaging potential at all points in the flow field. Data from a combination of our simulations and others experimental results will be condensed into a reduced order jet model that is simple to implement in Casa Grande simulation of a LOCA.

The formulation of this reduced order model will begin this summer. The first task will be to identify the family of dimensionless parameters that describe the jets behavior. An initial literature review has shown that two distinct correlations will likely be needed, one that describes field behavior close to the break plane versus another that characterizes far field, asymptotic jet behavior. A further challenge will arise in devising a means for interpolating between the near- and far-field correlations.