Simulation of the thermal field of submerged supercritical water jets at near-critical pressures

J. of Supercritical Fluids 75 (2013) 128–137

Contents lists available at SciVerse ScienceDirect

The Journal of Supercritical Fluids journal homepage: www.elsevier.com/locate/supflu

Simulation of the thermal field of submerged supercritical water jets at near-critical pressures Martin J. Schuler, Tobias Rothenfluh, Philipp Rudolf von Rohr ∗ ETH Zurich, Institute of Process Engineering, Sonneggstrasse 3, CH-8092 Zurich, Switzerland

a r t i c l e

i n f o

Article history: Received 12 October 2012 Received in revised form 12 December 2012 Accepted 13 December 2012 Keywords: Supercritical water jet Heat transfer Pseudo-critical point Near-critical pressure Turbulent Prandtl number Numerical model

a b s t r a c t For the investigation of fluid flow in the supercritical state, computational fluid dynamics (CFD) have become increasingly important to gain information that cannot be obtained experimentally. Especially at near-critical pressures, the application of CFD methods is still a challenging task due to the strongly varying thermo-physical properties that complicate the solution process of the Favre-averaged conservation equations. The focus of the present work is on a numerical model (ANSYS FLUENT® ) simulating flow and temperature fields of supercritical water jets discharged in a subcritical water bath at near-critical pressures. However, the used turbulence model tended to overestimate the values of the turbulent thermal conductivity near the pseudo-critical point of water. Various approaches based on a variable turbulent Prandtl number are presented to deal with these difficulties. Additional experiments were performed in a high-pressure vessel to record axial temperature profiles of the supercritical water jet and validate the numerical model. © 2012 Elsevier B.V. All rights reserved.

1. Introduction – background and motivation Several experimental and theoretical investigations dealing with submerged fluid jets under supercritical conditions are available in literature, especially in the field of combustion science. The efficiency of combustion processes can often be increased by atomization and break-up of the fluid jets (e.g. fuel) at elevated pressure conditions in the combustion chamber. A common way of investigating jet break-up over a broad range of thermodynamic conditions is the injection of a single, non-reacting fluid (e.g. nitrogen) at sub- (T < Tcp ) or supercritical (T > Tcp ) temperatures into an environment at various pressures ranging from subcritical (p < pcp ) to supercritical (p > pcp ). The temperature in the expansion chamber can be adjusted to sub- or supercritical values [1–5]. In jet break-up under subcritical pressures, ligaments and drops are formed during jet disintegration. Approaching the critical point (p = pcp , T = Tcp ), surface tension and latent heat vanish and thus the jet evolution is dominated by turbulent mixing [6]. For the experimental characterization of jet break-up in the supercritical state, Raman methods and shadowgraphy are frequently used [2].

∗ Corresponding author. Tel.: +41 44 632 92 63; fax: +41 44 632 13 25. E-mail addresses: [email protected] (M.J. Schuler), rothenfl[email protected] (T. Rothenfluh), [email protected] (P. Rudolf von Rohr). URL: http://www.ltr.ethz.ch (P. Rudolf von Rohr). 0896-8446/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.supflu.2012.12.023

Generally, the thermo-physical properties of supercritical fluids are in between those of gases and liquids. At the critical point (p = pcp , T = Tcp , thermodynamic singularity) the isobaric heat capacity (cp ), the isentropic compressibility and the thermal conductivity () are infinite [6]. By increasing the pressure even further beyond pcp the peak of the isobaric heat capacity is lowered, broadened and shifted toward higher temperatures. So, for each supercritical pressure there is a corresponding temperature yielding the maximum value of cp at the so-called pseudo-critical point (PCP). This temperature is called the pseudo-critical temperature (PCT). The line connecting all these pseudo-critical points is named pseudo-critical line (PCL) [7]. In the vicinity of this PCL all thermo-physical properties are very sensitive to temperature variations by undergoing sharp changes [6–9]. Performing experimental studies on supercritical jets are often troublesome, since the experimental setup has to be exposed to elevated temperatures and pressures and access for measurement equipment is often limited. Therefore, computational fluid dynamics (CFD) are frequently applied [6,8]. Due to the strongly varying thermo-physical properties around the pseudo-critical point solving the coupled conservation equations (mass, momentum, and energy) is a challenging task that has to be undertaken with care [6]. The applicability and validity of ordinary turbulence models originally developed for constant property flows under supercritical conditions is discussed in literature [6,7]. At supercritical pressures close to the pseudo-critical point (PCP), the molecular Prandtl number (Pr = cp /) varies significantly and the use of a constant turbulent Prandtl number (Prt ) in the turbulence models is often

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

inappropriate [7]. Obtaining numerical stability and convergence for such flows is ambitious due to the thermodynamic nonidealities and transport anomalies [6–8]. In spite of these difficulties, ordinary turbulence models have been used extensively to simulate supercritical flows. The Reynolds-averaged Navier–Stokes (RANS) equations, closed with the widely established k–ε turbulence models, have been applied to gain deeper insight into the highly dynamic flow [3–5,10]. Additionally, the more sophisticated, but computationally expensive, Large-Eddy Simulation approach was successfully applied to supercritical flow [1,11]. Generally, turbulent momentum, energy and mass transfer occur at different rates in the flow. The ratio between turbulent momentum and energy transfer is captured by the turbulent Prandtl number (Prt ) in common RANS turbulence models [12,13]. By applying the gradient transport hypothesis for turbulent momentum and energy transfer, Prt can be defined as the ratio of momentum eddy diffusivity to thermal eddy diffusivity Prt = t /˛t = t /˛t  [14]. For a wide variety of flows, Prt is treated as a model constant with a default value of 0.85 in ANSYS FLUENT® [13]. For homogenous turbulent flow fields a constant value between 0.9 and 1 is usually assigned to Prt [15,16]. For round turbulent jets, the use of different values for the jet’s near- and far-field are suggested in literature. In the near-field of the jet, a turbulent Prandtl number in the range of 0.4–0.5 is proposed, whereas values between 0.7 and 0.8 are recommended for the self-similar far-field of the jet [17,18]. The transport equation for the total specific energy (e) in ANSYS FLUENT® is: ∂ ∂ ∂ (e) + [ui (e + p)] = ∂t ∂xi ∂xj



eff

∂T + ui (ij )eff xj



(1)

where eff is the effective thermal conductivity, p represents the pressure, and (ij )eff is the deviatoric stress tensor (defined later). For the realizable k–ε turbulence model [19], eff is defined as the sum of the molecular thermal conductivity () and the turbulent thermal conductivity (t ) [12]: eff =  + t =  +

cp t Prt

(2)

t accounts for heat transfer due to turbulent fluctuations. The isobaric heat capacity (cp ), the turbulent viscosity (t ) and the turbulent Prandtl number (Prt ) are used to evaluate t as shown in Eq. (2) above. In fully turbulent flow fields  is negligible in comparison to t . Viscous heating effects can be included in the model by the deviatoric stress tensor (ij )eff [12]:



(ij )eff = eff

∂uj ∂xi

+

∂ui ∂xj



−

∂uk 2 ı  3 ij eff ∂xk

(3)

where uj is the mean velocity in xj direction, eff the effective viscosity and ıij represents the Kronecker delta. Since Prt was found to vary across the different regions of shear flow [17,18], several authors chose an approach with a variable Prt number in their models. Aouissi et al. [16] successfully applied such an approach to heated jets and diffusion flames. In RANS modeling of aero-propulsive flows, variable Prt numbers have also proven a successful concept to improve numerical predictions [14,15]. In contrast to other components such as nitrogen, water has only been rarely used in supercritical jet studies. This is mainly due to the considerably higher pressure and temperature values water = = 220.64 bar, Tcp required to reach the critical point (pwater cp ◦ 373.946 C) [20]. To the authors’ knowledge only two experimental investigations [21,22] regarding submerged supercritical water jets have been published up to date. In the comprehensive experimental study of Rothenfluh et al. [22], submerged supercritical water jets in a slowly co-flowing annular cooling water stream

129

were investigated over a wide range of conditions. Four different nozzle diameters were applied in the experimental series. Entrainment effects and turbulent mixing of the hot supercritical water with the cold surrounding cooling water are quantified by determining the length of the supercritical plume. Optical Schlieren techniques were used to detect the supercritical plume length. These lengths are found to be equal to the injector’s nozzle diameters and almost independent of the jet’s nozzle exit temperature and the jet’s supercritical water mass flow rates under almost all experimental conditions. A theoretical work by Sierra-Pallares et al. [23] presents a numerical study focusing on mixing at the macro and micro scale of supercritical water jets. Similar to the present study, also numerically determined axial temperature profiles are compared to experimental values [21]. The presented work is motivated by the idea of using submerged supercritical water jets in an alternative drilling method named “Hydrothermal Spallation Drilling” (HSD) [21]. HSD is a promising alternative drilling technology that could prove to be economically advantageous over conventional rotary techniques for drilling deep wells in hard rock formations [24,25]. This drilling technique uses the properties of certain rock types to disintegrate into small fragments when heated up rapidly by a highly energetic impinging jet [26,27]. In water (resp. water-based drilling fluid) filled boreholes below 2 km depth, water exceeds its critical pressure and thus hydrothermal flames or supercritical water jets are favored to provide the required heat for thermal rock-fragmentation. A possible realization of a spallation drilling head consists of a combustion chamber fed by water, fuel and an oxidant. The water inflow is used to control momentum flux and nozzle exit temperature of the discharged hot water jet. Fuel and oxidant are preheated and subsequently ignited to form a hydrothermal flame in the aqueous environment of the combustion chamber. The water in the chamber is thus heated up to high, supercritical temperatures and ejected through a nozzle together with the combustion products. This highly energetic supercritical water jet is finally directed toward the rock surface to induce thermal fragmentation. The thermal field of these jets is of major interest for HSD and is therefore investigated in the present study. In the experiments of the present study the thermal field of submerged supercritical water jets was investigated for a wide range of operating conditions. Electrically preheated supercritical water (SCW) was injected into a vessel containing subcritical water at 224 bar. A nozzle diameter (d0 ) of 3 mm was used to discharge the supercritical water jet into the subcritical water bath in gravitational direction. A fine thermocouple was moved along the jet’s axis to record temperature profiles for validation of the numerical flow simulations [22]. In the theoretical part of this work numerical models were developed and implemented in the commercial CFD code ANSYS FLUENT® (12.1.4) [13] to simulate flow and temperature fields of submerged supercritical water jets. The model is based on conservation of momentum, mass and energy, and includes the variation of the water’s physical properties according to the National Institute of Standards and Technology (NIST) Standard Reference Database [28]. For the simulation of round jets at high Reynolds numbers, the use of the realizable k–ε turbulence model [19] is suggested in literature [12,19]. This model was also applied in the present study, where temperatures along the jet axis were extracted from the simulations and compared to the experimental data. Difficulties in the simulations occurred around the pseudocritical point of water and became manifest in the axial temperature profiles. Different approaches based on a variable turbulent Prandtl number (Prt ) were applied to deal with these problems. First, the scalar-variance model according to Brinckman et al. [14] was implemented to calculate a local Prt number. This model has been proved to have a wide applicability ranging from low-speed and

130

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

low-temperature flows to high-speed and hot, reacting flows. For this study, the model of Brinckman et al. was combined with the realizable k–ε turbulence model [12,19]. In addition to the application of this scalar-variance model, also different Prt functions were applied around the pseudo-critical point of water to account for the strong variations of the molecular Prandtl number. The present work also highlights additional challenges related to the strongly varying thermo-physical properties around the PCP when simulating supercritical water systems with commercially available CFD tools. 2. Methods 2.1. Experiments The deionized water for the supercritical jet was electrically preheated to supercritical temperatures. The jet was discharged into a bath of subcritical water at 224 bar in gravitational direction. A discharge nozzle with a diameter of 3 mm was used. The high-velocity jet was surrounded by a co-annular, low velocity flow of deionized cooling water (CW) at 20 ◦ C. Two symmetrical outlets were located at the top of the vessel for the exhaust water consisting of the mixed streams SCW and CW. The mixing temperature at the outlets was subcritical (60–140 ◦ C). Optical access to the jet is given by four sapphire glass windows mounted on two opposite sides of the vessel. A more detailed description of the experimental setup is provided by Rothenfluh et al. [22]. To gain information about the jet’s thermal field, temperature profiles were recorded along the jet axis. For this purpose a fine sheath thermocouple (∅ 0.25 mm, K-type, tolerance class 1, Thermocoax (France)) was inserted from a port at the bottom of the vessel to be moved along the jet’s axis. The thermocouple was fixed in a ceramic rod (∅ 0.4 mm × ∅ 0.8 mm, Friatec (Germany)) for stabilization with its tip protruding about 3 mm from the rod. Vertical movement of the temperature probe was realized by a high precision linear actuator (Schneider Electric, Germany) situated outside the vessel. A CCD camera was used to control the probe’s position through the windows of the pressure vessel. Axial temperature profiles in the experiments were recorded starting from 20 mm below the nozzle and moving gradually to a position flush with the nozzle orifice. Close to the nozzle orifice a spatial resolution (step size of thermocouple) of 0.3 mm was chosen for the axial temperature profiles. Beyond 6 mm from the nozzle outlet, however, the axial step size was increased to 0.5 mm, since spatial temperature gradients become more moderate there. 2.2. Numerical model 2.2.1. Governing equations and turbulence model The conservation equations for momentum and mass were solved as Favre-averaged Navier–Stokes equations [12,29,30]. Also the energy conservation equation is needed for the calculation of the thermal field in a non-isothermal jet. To close these equations, further model equations and assumptions are needed to determine the unknown turbulent Reynolds stresses and the turbulent heat flux. In the realizable k–ε turbulence model, an additional equation for closure is obtained by introducing the Boussinesq approximation. Also, two additional transport equations for the turbulent kinetic energy (k) and its dissipation rate (ε) are solved [12,19,29,30]. The realizable k–ε turbulence [12,19] model avoids the “round-jet plane-jet anomaly”. This term describes the effect that two-equation turbulence models predict spreading rates of round jets rather poorly as opposed to plane jets, where predictions are satisfactory [12,19,31]. One of the main differences between the realizable k–ε turbulence model and other approaches of the k–ε

family is the definition of the turbulent viscosity (t ): k2 1 with Sij = ε 2

t = C (Sij , ˝ij , k, ε)

and ˝ij =

1 2



∂uj ∂ui − ∂xi ∂xj



∂uj ∂xi



+

∂ui ∂xj



(4)

whereas C is treated as a constant in most models (e.g. C = 0.09 in the standard k–ε model [12]), its value in the realizable k–ε turbulence model depends on the mean rates of strain (Sij ), mean rate of rotation (˝ij ) and on k and ε. In detail C is evaluated as follows [12]:

 1.0 Sij Sij + ˝ij ˝ij (5) with U ∗ = ∗ A0 + AS (kU /ε) √ where A0 = 4.04 and AS = 6 cos( ). The quantity is defined as [12]: C =

=

√ 1 cos−1 (W 6) with W = 3

Sij Sjk Ski



(

Sij Sij )

(6)

3

A variable C helps in making the realizable k–ε model more consistent with the underlying turbulent physics when calculating the Reynolds stresses [12,19]. The values of C in the investigated turbulent jet regions of interest vary between ∼0.015 to values up to ∼0.25. For most constants in the chosen turbulence model, the default values suggested by ANSYS FLUENT® were applied in the present work. For the turbulent Prandtl number (Prt ), however, a value of 0.7 was used instead of the default value equal to 0.85. This modification was introduced, because the major part of the investigated round turbulent jet is expected to be in the self-similar region, for which Prt values between 0.7 and 0.8 are recommended in literature [17,18]. Additionally, the scalar-variance model proposed by Brinckman et al. [14] was implemented to calculate local values of Prt based on the turbulent velocity (k/ε) and turbulent thermal (ke /εe ) time scale of the flow. Thus two additional transport equations for the energy variance ke (Eq. (7)) and its dissipation rate εe (Eq. (8)) were incorporated. These equations generally have a similar form as the k–ε turbulence model equations. Brinckman et al. [14] extended and generalized the model proposed by Sommer et al. [32], which is based on a temperature variance formulation for a variable turbulent Prandtl number in compressible flows. The final result is an internal energy variance method for non-reacting and reacting aeropropulsive flows. For highly nonuniform thermo-physical fluid properties (especially cp ) in single-component and multi-species flows, the step mentioned above is recommended [14]. The supplementary equations for ke and εe were implemented by user-defined scalars (UDS) in ANSYS FLUENT® [13,33]. ∂ ∂ ∂ (uj ke ) = (ke ) + ∂t ∂xj ∂xj

  

˛+

˛t ıke



∂ke ∂xj



 + 2˛t

∂e ∂xj

2 − 2εe

(7)

uj is the mean velocity in xj direction, ˛ = /cp  defines the thermal diffusivity and ıke = ıεe = 1.0. ∂ ∂ ∂ (εe ) + (uj εe ) = ∂t ∂xj ∂xj + ˛t





Cd1

− εe Cd4

εe ε + Cd2 ke k

εe ε + Cd5 ke k

 

 ˛+

 ∂e 2



∂xj + WDF

˛t ıεe

∂ε 

+ Cd3 Pk

e

∂xj εe k (8)

The term Pk incorporates production of turbulent kinetic energy due to mean velocity gradients (Gk ), generation of turbulence due

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

to buoyancy (Gb ) and compressibility effects on turbulence by the dilatation dissipation term (Ym ) [12,13]: Pk = Gk + Gb − Ym = t S 2 + ˇgi with S =



t ∂T k − 2ε 2 Prt ∂xi a

2Sij Sij

(9)

ˇ is the coefficient of thermal expansion [28] and gi is the component of the gravitational vector in xi direction. The near-wall damping function WDF according to Sommer et al. [32] is implemented 2



WDF = e−(Ret /80)  (Cd4 − 4)

εˆ e (εe ∗ )2 εˆ εe + Cd5 εe − ke k ke



(10)

 √ 2 2 where the terms εˆ e = εe − ˛(∂ ke /∂y) , εˆ = ε − (2/)(∂ k/∂y) ∗ 2 and εe = εe − (˛ke /y ) are multiplied with the damping term 2

e−(Ret /80) . All constants (Cd1 = 2, Cd2 = 0, Cd3 = 0.72, Cd4 = 2.2, and Cd5 = 0.8) are unchanged with respect to the recommendation in [14]. Based on the turbulent velocity (k/ε) and turbulent thermal (ke /εe ) time scale, the local thermal eddy diffusivity is:



˛t = C f k

k ke ε εe

(11)

with C = 0.14. The damping function, f [14,32], is: f = e−(Ret /80)

2

C1 Re0.25 t

+ (1 − e−y

∗ /A+

)

2

(12)

where the turbulent Reynolds number, Ret = k2 /ε; C1 = 0.1; A+ = 45; and the dimensionless wall distance, y* = uε y/, is defined by the Kolmogorov velocity scale, uε = (ε/)1/4 , and the wall normal distance, y. Near the wall, ke equals zero and the εe value is



2

given by ˛wall [(∂ ke /∂y) ]wall . Finally, the local turbulent Prandtl number is determined with: t Prt = C2 (13) ˛t The additional constant C2 = 0.5 does not appear in the expression proposed by Brinckman et al. [14]. t is evaluated by the realizable k–ε turbulence model [12,19]. 2.2.2. Geometry and computational domain In the experiments, two symmetrical outlets were located on top of the vessel [22]. This is not exactly consistent with the simulation domain of the present study, where the outlet was located at the bottom (see Fig. 1) due to simplicity and symmetry reasons. In the work provided by Rothenfluh et al., penetration lengths of

131

supercritical water jets were determined by a nonintrusive optical Schlieren method. These supercritical plume lengths, however, were not influenced by the choice of the outlet position (top or bottom) over the whole range of applied experimental conditions [22]. Moreover, simulations with bottom and top outlet position were executed and their axial temperature profiles were compared to each other. The temperature field up to 20 mm below the nozzle exit was not affected by the location of the outlet. The axisymmetric simulation domain depicted in Fig. 1 only included the core part of the high pressure vessel (radius 18.2 mm, length 123.5 mm) and excluded the outer parts, where velocities were approximately zero. The simulation of the injection system for supercritical water and cooling water (left part of Fig. 1) exactly reflected the experimental conditions and also accounted for heat transfer in the solid material. All solid parts of the injection system were made of stainless steel 1.4435 (316 L). An insulating, annular air gap inside the inner part of the injector between the central, hot water stream (SCW) and the co-annular cooling water (CW) flow minimized heat losses of the hot water stream. Constant material properties were assigned to all solid steel parts according to the data sheet of the supplier. Also the air gap insulation of the injector was treated as a solid with its physical properties changing as a function of temperature [28]. In addition to the computational domain, Fig. 1 also shows a typical temperature field of the investigated jet flow. A more detailed description of the injection system and the high-pressure vessel can be found in [22]. GAMBIT® (2.2.30) was used to generate the computational domain. In the simulations for a nozzle diameter of 3 mm, a structured mesh containing approximately 875,000 cells was used to guarantee a mesh-independent solution. Enhanced wall treatment including the required mesh refinement criteria in wall normal direction was applied in the whole domain. Mass-flow-inlet boundary conditions for the SCW and CW streams were applied. A pressure-outlet was chosen as boundary condition for the effluent flow [12,13]. 2.2.3. Thermo-physical fluid properties Strong variations of thermo-physical properties, especially close to the pseudo-critical point for a given pressure, require a careful implementation of these properties. The NIST Standard Reference Database 23 (REFPROP® , Version 8.0) was used to evaluate the thermo-physical properties of water. This database represents the currently most accurate database available from NIST [28]. The water’s physical properties (density (), thermal conductivity (), dynamic viscosity (), isobaric heat capacity (cp ) and speed of sound (a)) at the respective temperatures and pressures were calculated by this REFPROP® 8.0 database [28]. The uncertainties of

Fig. 1. Contour plot of the temperature field in the 2D-axisymmetric computational domain including the injection system. The flow was simulated for a nozzle diameter of 3 mm at a pressure of 224 bar. The thermo-physical properties of water were evaluated as a function of temperature and pressure according to the REFPROP® 8.0 database. Mass flow rates of 4 g/s for SCW and 65 g/s for CW were applied. The realizable k–ε model in ANSYS FLUENT® was used.

132

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

these property values for the applied temperature and pressure ranges can be found in [20,28]. All fluid properties were embedded in ANSYS FLUENT® via user defined functions (UDF) [33]. If a constant pressure in the whole computational domain is assumed, isobaric material properties can be applied. For isobaric simulations the entire property data was tabulated within the temperature range of interest using a resolution of 0.1 K. For temperatures falling in between the tabulated values linear interpolation was applied. However, the assumption of a constant pressure in the entire flow field does not exactly reflect the real situation in the experiments. The dynamics of the flow cause small variations in pressure across the flow field. To investigate also the pressure sensitivity of the jet flow, all routines of the REFPROP® 8.0 database [28] for the calculation of the fluid properties were directly implemented in ANSYS FLUENT® . These routines were used to evaluate the corresponding fluid properties in each cell according to the local pressure and temperature assigned to it. For the operating pressure of 224 bar used in the experiments, the pseudo-critical temperature of water is found at 375.21 ◦ C. All thermo-physical properties undergo strong changes in the vicinity of this PCT, which particularly applies to the isobaric heat capacity (cp ) reaching extraordinarily high values at this point. Around the PCT, the so-called “pseudo phase change” (PPC) occurs and thus water passes from a liquid-like state into a gas-like state with increasing temperature. Hence, the thermo-physical properties also change from a liquid-like to a gas-like behavior [20,28]. Fig. 2 shows the dimensionless variations of the thermo-physical properties of water at 224 bar as a function of the fluid temperature. All properties were normalized by their respective values at 25 ◦ C and 224 bar. For the solution process it is important that the water’s specific enthalpy (h) and isobaric heat capacity (cp ) are implemented in a consistent manner, i.e. cp = dh/dT and dh = cp dT, respectively. The reason for this is the way ANSYS FLUENT® calculates the water’s specific enthalpy by integrating the cp values provided by the user with respect to temperature [12,33]. 2.2.4. Solver settings and solution procedure The ANSYS FLUENT® pressure-based coupled solver was applied. Steady state simulations including gravity (9.81 m/s2 ) and viscous heating were performed. Second-order scheme for pressure and second-order upwind scheme for all other variables were chosen for spatial discretization. Gradients were evaluated by the least squares cell-based method [12,13]. The default under-relaxation factors recommended by ANSYS FLUENT® had to be reduced in order to reach a convergent solution. A reduction in between 25% and 75% compared to the default values in the solver settings was required to guarantee numerical stability. The simulations were

Dimensionless property [-]

2.5 density dynamic viscosity thermal conductivity cp

2.0 1.5 1.0 0.5 0.0 25

100

200

300

400

500

600

Temperature [°C] Fig. 2. Thermo-physical properties of water at 224 bar as a function of temperature and normalized by the value at 25 ◦ C and 224 bar [28].

parallelized on the Brutus Cluster at ETH Zurich. The Dynamic Gradient Adaption Approach in ANSYS FLUENT® [12,13] was used for mesh refinement to reach a mesh-independent solution. An expanded mesh convergence study was performed using the experimental operating conditions that yielded the highest gradients for all variables in the flow field. The mesh study primarily showed that the computational grid had to be vigorously refined in the free shear layer of the jet, where temperatures were close to the pseudo-critical point and fluid property variations were more pronounced. Selected variables and parameters characterizing the jet flow (e.g. potential core length) were recorded after each iteration step of the simulations to check, if a iteration-independent and mesh-independent (independent of number of cells) solution was reached.

3. Results and discussion 3.1. Axial temperature profile In Fig. 3(a) an experimental temperature profile recorded along the jet axis at 224 bar is shown and compared to the corresponding temperature profiles extracted from three different flow simulations. However, the profiles of the three different simulations cannot be distinguished because the deviations were negligibly small. The differences in these flow simulations refer to the calculation of the individual thermo-physical fluid properties of water: The first two simulations both used the REFPROP® 8.0 database to calculate the water’s thermo-physical properties in each computational cell. Whereas isobaric conditions (constant pressure of 224 bar) were applied to the case denoted “Simulation REFPROP isobaric”, the calculation of property data also accounted for small variations in the absolute pressure in the simulation denoted “Simulation REFPROP”. The third simulation denoted “NIST Real Gas Model” used the ANSYS FLUENT® NIST Real Gas Model [12,13] to evaluate the fluid properties as a function of temperature and pressure. In the NIST Real Gas Model, an earlier version of the REFPROP® database (Version 7.0) is applied with respect to the two other approaches explained above. In all mentioned cases Prt was set to 0.7. The zoom of Fig. 3(b) showing a small section of Fig. 3(a) visualizes the small differences between all calculations. The simulations using temperature and pressure dependent fluid properties (“Simulation REFPROP” and “NIST Real Gas Model”), however, were extremely computationally expensive, since fluid properties were evaluated iteratively in every cell of the domain. For this reason isobaric simulations (“Simulation REFPROP isobaric”) were favored in this study. For the experiment and all simulations, a nozzle diameter of 3 mm with a SCW mass flow rate of 4 g/s was used. The mass flow rate of the annular co-flowing cold water stream was 65 g/s with an inlet temperature of 20 ◦ C. The thermal inlet boundary condition of the SCW stream (see Fig. 1 left) had to be adjusted to guarantee a nozzle exit temperature (T0 ) of 460 ◦ C. The flow at the matched nozzle exit conditions was clearly in the turbulent regime, because a nozzle exit Reynolds number (Re0 ) of roughly 60,000 was reached. In the near-field of the jet ranging from the nozzle outlet position (temperature T0 ) to the position where the pseudo-critical temperature of water (375.21 ◦ C) was reached, an acceptable agreement between experiment and simulations is observed in Fig. 3(a). The cooling of the jet, however, is slightly over predicted. In a narrow temperature band around the PCP, the simulation-based axial temperature profiles undergo considerable flattening. This behavior does not reflect the underlying physics of supercritical submerged jet flow, but was rather found to be a peculiarity of the underlying turbulence model. The sharp peak in the water’s specific heat (cp ) at the PCP tends to lead to a local overestimation of the turbulent

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

400

Axial temperature [°C]

Axial temperature [°C]

400

Experiment Simulation REFPROP isobaric Simulation REFPROP NIST Real Gas Model

500

300 200 100

(a) 0 0

5

10

15

20

133

25

398 396 394 392

(b) 390 1.90

Axial distance from nozzle [mm]

1.95

2.00

2.05

2.10

2.15

2.20

Axial distance from nozzle [mm]

Fig. 3. (a) Comparison of calculated axial temperature profiles at 224 bar with the experiment. The profiles of the three different simulations cannot be distinguished because the deviations were negligibly small. A nozzle diameter of 3 mm with a SCW mass flow rate of 4 g/s was used. The mass flow rate of the annular co-flowing cold water stream was 65 g/s (inlet temperature: 20 ◦ C). All simulated profiles were obtained with the realizable k–ε model in ANSYS FLUENT® using a turbulent Prandtl number of 0.7. Three different ways of calculating the thermo-physical properties of water were used: For “Simulation REFPROP isobaric” fluid properties were evaluated as function of temperature at a constant pressure of 224 bar. “Simulation REFPROP” calculates pressure and temperature dependent fluid properties. Both simulations use the REFPROP® 8.0 database. In the simulation denoted “NIST Real Gas Model” the ANSYS FLUENT® NIST Real Gas Model was implemented. (b) Zoom of a small section of (a) illustrating the differences between the calculated temperature profiles.

Simulation 224bar Simulation 250bar Simulation 300bar Simulation 350bar Simulation 400bar

Axial temperature [°C]

500 400 300 200 100 0

0

5

10

15

20

25

Axial distance from nozzle [mm] Fig. 4. Comparison of calculated axial temperature profiles for different pressure levels. A nozzle diameter of 3 mm, a SCW mass flow rate of 4 g/s and a CW mass flow rate of 65 g/s (20 ◦ C) were used. Prt was set to 0.7 in the applied realizable k–ε turbulence model.

higher pressures. This applies also to all other thermo-physical properties of water. In the implemented realizable k–ε turbulence model, it is exclusively the quantity t that controls turbulent heat transfer (see Eqs. (1) and (2)). Since turbulent heat transfer is a crucial mechanism in the investigated jet flow, the turbulent thermal conductivity (t ) assumes great importance in the calculation of the presented temperature profiles. The values of t as a function of the jet’s axial coordinate are given in Fig. 5 for all simulations of Fig. 4. For supercritical pressures close to pcp (e.g. 224 bar), the values of t are locally overestimated at the pseudo-critical point. As a result of these high t -values, the temperature profile for 224 bar exhibits a significant flattening at the PCP, which does not reflect the experimental findings (see Fig. 3a). This temperature plateau develops, because heat is conducted too rapidly across the water zones close to the PCP in the simulations, where the high cp values are encountered. Another indication for the overestimation of t is the fact that the simulated supercritical jet of Fig. 3a is cooled down too fast after exiting the nozzle when compared to the experiment. This is another hint that turbulent heat transfer (represented by t ) is overestimated close to the PCT. The reason for the overestimation of the turbulent thermal conductivity (t ) can be found in the way it is calculated in the realizable k–ε turbulence model used here: t =

cp t Prt

λ t turb. therm. conductivity [W/mK]

thermal conductivity (t ) as discussed later. Even a very vigorous mesh refinement around the PCP did not have the desired effect of removing the rather unphysical temperature plateau appearing in all profiles of the simulations. This temperature plateau around the PCP is responsible for the subsequent deviations between simulated and experimental temperature values in the jet’s far-field. However, the development of the experimental and calculated curves is still similar. Since the spatial temperature gradients are predicted correctly, the temperature profiles of the simulation and the experiment run parallel to each other in the jet’s far-field. The development of the calculated temperatures along the jet’s axis is illustrated in Fig. 4 for five different supercritical pressure levels. Isobaric conditions were assumed for the calculation of the thermo-physical properties. When the operating pressure is increased, the temperature plateau at the PCP in the profiles is shifted toward higher temperatures and becomes less pronounced. For a pressure above 350 bar, the flattening of the temperature profile finally disappears. The simultaneous shifting and weakening of the temperature plateau in these profiles is primarily a result of the PCP-related heat capacity (and consequently Pr) peak behavior becoming less pronounced and being shifted to higher temperatures as pressure increases. But not only the development of the cp values across the pseudo-critical point are smoothed out with

(14)

10000 Simulation 224bar Simulation 250bar Simulation 300bar Simulation 350bar Simulation 400bar

8000 6000 4000 2000 0

0

5

10

15

20

25

Axial distance from nozzle [mm] Fig. 5. Calculated axial profiles of the turbulent thermal conductivity for different supercritical pressure levels. In the simulations for a nozzle diameter of 3 mm, the realizable k–ε model was used in ANSYS FLUENT® . A SCW mass flow rate of 4 g/s and a CW mass flow rate of 65 g/s at 20 ◦ C were applied. Prt was set to 0.7.

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

This term appears in the total energy transport Eq. (1). As a result of this definition, the development of t as seen in Fig. 5 directly reflects the corresponding development of the water’s isobaric heat capacity (cp ). It is mainly the assumption of a constant turbulent Prandtl number (0.7 in the simulations above) that leads to a local over prediction of t at the pseudo-critical point. However, the assumption of a constant turbulent Prandtl number is questionable, especially for flows at near-critical conditions, and can lead to large errors of numerical results [7]. There is evidence that Prt should not be treated as a constant, but somehow linked to the molecular Prandtl number (Pr) [34–36].

(15)

The above expression was implemented in the following way: Prt exhibits a predefined, maximum value Pˆ rt at the PCT to account for the locally overestimated value of t as discussed in the previous section. On both sides of the PCT the turbulent Prandtl number decreases rapidly and takes the default value of 0.7 for temperatures lower than T1 (below the PCT) and higher than T2 (above the PCT). Hence, function (15) above only comes into play in a pre-defined, narrow temperature interval [T1 , T2 ] around the PCT, whereas the turbulent Prandtl number for all other temperatures is at 0.7. To fulfill these restriction a Gaussian distribution function was chosen for (15): 2

(16)

with a1 = 1.75, b1 = 648.4 and c1 = 15.67. The resulting function (16) is defined by the three points (633.357 K/0.7), (648.357 K/1.75) and (663.357 K/0.7) inside a temperature band of 30 K. Therefore turbulent Prandtl numbers were only adapted within a small section of the jet, where t reaches the elevated values. A physically more meaningful type of function shows a dependency of Prt on the molecular Prandtl number (Pr = cp /) instead of temperature. For this purpose the introduction of the turbulent Peclet number is convenient: Pet = Ret Prt =

t k2 Pr with Ret =  ε

(17)

The final expression applied in the numerical calculations is given by Prt = k0

1

 Pe

Cratio

Ret

t

Sim. Prt = 0.7 Sim. Prt = f(T) eq. (16) Sim. Prt = Pe t/Re t eq. (18) Sim. Scalar variance model

2.5 2.0 1.5 1.0 0.5 0.0 2

4

6

8

10

12

Axial distance from nozzle [mm]

In this section, various ideas to implement a variable turbulent Prandtl number are presented and discussed. The isobaric heat capacity (cp ), and consequently t , peaks around the PCT. cp was implemented as a temperature depended quantity in the model and thus one possibility could be to formulate Prt also as a function of the fluid temperature (T):

Prt = a1 e−((T −b1 )/c1 )

3.0

0

3.2. Axial temperature profile with a variable turbulent Prandtl number

Prt = f (T )

Prt turbulent Prandtl number [-]

134

(18)

where the constants k0 = 6.75 and Cratio = 3 are explained later in the text. The call of Eq. (18) during simulations has to be controlled by certain criteria. For the definition of such criteria, ratios between turbulent and molecular transport properties were introduced ((t /) and (t /)). Setting the turbulent quantities (t and t ) in relation to the molecular quantities ( and ) makes only sense when a certain level of turbulence is present in the flow. Therefore Eq. (18) comes only into play when Ret has a local value exceeding 2500. In the computational domain Ret values up to roughly 43,000 were reached. The ratio (t /)/(t /) has proven very sensitive around the PCP of water. Values up to approximately 60 were detected there for Prt = 0.7 (without the use of Eq. (18)). Therefore, this ratio was finally applied in a second criterion to

Fig. 6. Axial profiles of Prt for the different model approaches discussed in Section 3.2. In all simulations the realizable k–ε model was used. A SCW mass flow rate of 4 g/s, a CW mass flow rate of 65 g/s at 20 ◦ C and a nozzle diameter of 3 mm were applied (224 bar).

control the use of Eq. (18). If (t /)/(t /) exceeded a set threshold value of Cratio = 3, Eq. (18) was called to reduce (t /) by increasing Prt around the PCP (see also Eq. (14)). After this correction step, the maximum value of (t /)/(t /) around the PCP in the entire domain was in the range of 45. The appropriate value for Cratio has to be determined individually for each turbulent flow field. Cratio should be chosen in a way, that Eq. (18) only acts around the PCP of water. k0 was introduced to ensure a smooth intervention of Eq. (18). Here, too, expression (18) was only used in a small region around the PCP, if the following conditions were fulfilled: (t /) > Cratio (t /)

for Ret > 2500 with Cratio = 3

(19)

The value of Prt in all other regions was set to 0.7. The last approach used the scalar-variance model according to Brinckman et al. [14], where two additional transport equations for ke and εe were solved to determine a local turbulent thermal time scale (ke /εe ). This time scale allowed for the calculation of the turbulent thermal eddy diffusivity (˛t ) and the local value of Prt (see Eqs. (11) and (13)). To favor numerical stability when using this model, a simplified computational domain was used compared to the domain depicted in Fig. 1. In Fig. 6 the axial profiles of Prt are illustrated for all approaches discussed above. In the simulation named “Prt = 0.7”, a constant value for the turbulent Prandtl number equal to 0.7 was applied. The approaches “Prt = f(T)” and “Prt = Pet /Ret ” were using the Eqs. (16) and (18), respectively, resulted both in a sharp increase of Prt around the PCT. The scalar-variance model [14] showed a decrease of Prt to a minimum value of 0.2 followed by a significant increase of Prt around the PCP along the axis of the jet. Elevated Prt values were detected for temperatures far beyond the PCT. This means that the scalar-variance model is less sensitive around the PCP compared to the other two approaches (see additionally Figs. 9–11). However, both in the potential core (Prt = 0.4 − 0.5) and the selfsimilar far-field of the jet (Prt = 0.7 − 0.8), the predicted Prt values are consistent with values found in literature [17,18]. Only for the scalar-variance model, the peak value of Prt did not appear at the PCT and was slightly shifted toward lower temperatures (see Fig. 11). The corresponding axial profiles of t for all approaches are presented in Fig. 7(a) and (b). The increase of Prt around the PCP for all models with a variable turbulent Prandtl number yielded maximal values of t significantly reduced compared to the case of a constant turbulent Prandtl number of 0.7. Also the width of the t peak was influenced by the varying Prt values. A comparatively narrow peak of t was detected for the formulation of Prt as a function of Pet and Ret (see Fig. 7(b)). For this approach only a small amount of water was at temperatures near the PCT undergoing high values of

3

80 60 40 20

(a)

0 0

2

4

8

6

10

135

50

3

Sim. Prt = 0.7 Sim. Prt = f(T) eq. (16) Sim. Prt = Pe t/Re t eq. (18) Sim. Scalar variance model

100

λ t turb. therm. cond. [ *10 W/mK]

λ t turb. therm. cond. [ *10 W/mK]

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

12

40 30 20 10

(b)

0 1

2

3

4

5

Axial distance from nozzle [mm]

Axial distance from nozzle [mm]

Fig. 7. (a) Axial profiles of the turbulent thermal conductivity for the different approaches discussed in Section 3.2. In all simulations the realizable k–ε model was applied. A SCW mass flow rate of 4 g/s, a CW mass flow rate of 65 g/s at 20 ◦ C and a nozzle diameter of 3 mm were applied (224 bar). (b) Zoom of a small section of (a) illustrating the differences between the calculated t profiles.

400

Axial temperature [°C]

Axial temperature [°C]

Sim. Prt = 0.7 Sim. Prt = f(T) eq. (16) Sim. Prt = Pe t/Re t eq. (18) Sim. Scalar variance model

(a)

500

300 200 100

(b)

500

400

300

200

0 0

5

10

15

20

25

0

1

Axial distance from nozzle [mm]

2

3

4

5

6

Axial distance from nozzle [mm]

Fig. 8. (a) Comparison of the calculated axial temperature profiles for the different approaches with a measured one (filled squares). In all simulations the realizable k–ε model was used. A SCW mass flow rate of 4 g/s and a CW mass flow rate of 65 g/s at 20 ◦ C and a nozzle diameter of 3 mm were applied (224 bar). (b) Zoom of a small section of (a) emphasizing the differences between the selected temperature profiles.

2.0 1.5 1.0 0.5 0.0

400

600 Pr t cp λt T

400

300 200

200

100 0

0 0

1

2

3

4

5

6

50 40 30 20 10 0

λ t turb. therm. cond. [kW/mK]

2.5

500

800

T axial temperature [°C]

3.0

c p isobaric heat capacity [kJ/kgK]

The results of all applied approaches using a variable turbulent Prandtl number are compared to each other (see Figs. 6–8). However, it is still unclear, how the crucial quantities (t , Prt , cp , T) interact with each other around the PCP for each individual approach. In order to illustrate these interactions, the axial profiles of cp , t , Prt and temperature (T) are given for the three considered ideas. Fig. 9 shows the behavior of the crucial variables for the approach denoted “Prt = f(T)” (see Eq. (16)), Fig. 10, on the other hand, illustrates the case “Prt = Pet /Ret ” (see Eq. (18)), whereas Fig. 11 visualizes the variables for the scalar-variance model [14]. There is evidence that a variable turbulent Prandtl number is a possible way to improve results of the thermal field when simulating the supercritical water jets investigated here. This is mainly

Prt turb. Prandtl number [-]

cp and t . This behavior has proven beneficial to the mitigation of the plateau in the temperature profiles discussed in Section 3.1 and illustrated in Fig. 8. For all applied approaches of a varying turbulent Prandtl number, t was decreased by the increase of Prt around the PCP in the near-field of the jet (see Fig. 7(a)). Hence the cool down of the jets after the nozzle exit was delayed and the sharp temperature drop in the supercritical range was shifted further downstream compared to the default case of Prt = 0.7 (see Fig. 8(b)). The spatial temperature gradients predicted by the scalar-variance model were too high in the near-field of the jet. The variable turbulent Prandtl number also influenced the flattening of the axial temperature profile. In contrast to the pronounced temperature plateau for the simulation with Prt = 0.7, all considered approaches achieve a certain reduction of the corresponding temperature plateaus. Therefore the agreement between experiment and simulations in the thermal far-field of the jet improved (see Fig. 8(a)), particularly for the scalar-variance model. This model by Brinckman et al. matched up the thermal far-field of the jet with the experimental data much better than any other approach applied in this study. This is mostly a result of the over-predicted temperature decay in the near-field of the jet. The idea of implementing a Prt number as a function of temperature (“Prt = f(T)”) was able to reduce the flattening in the temperature profile with respect to Prt = 0.7 (Fig. 8(b)), but elongated the supercritical region of the jet compared to the experiment. In the thermal far-field of the jet, the performance of both approaches (“Prt = 0.7” and “Prt = f(T)”) was similar (Fig. 8(a)). For the case “Prt = Pet /Ret ”, the flattening in the axial temperature profile almost vanishes, but the length of the supercritical jet region was overestimated compared to the experiment.

Axial distance from nozzle [mm]

Fig. 9. Axial profiles of cp , t , Prt and temperature (T) for the approach of calculating a variable turbulent Prandtl number denoted “Prt = f(T)” according to Eq. (16). These profiles illustrate how these quantities interact for calculating the turbulent heat transfer around the PCP of water.

2.5 2.0 1.5 1.0 0.5 0.0

500

800

400

600 Pr t cp λt T

400

300 200

200

100

50

T axial temperature [°C]

3.0

0

0 0

1

2

3

4

5

40 30 20 10 0

6

λ t turb. therm. cond. [kW/mK]

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137

c p isobaric heat capacity [kJ/kgK]

Prt turb. Prandtl number [-]

136

Axial distance from nozzle [mm]

2.0 1.5 1.0 0.5 0.0

500

Pr t cp λt T

600

400 300

400 200 200

100 0

0 0

1

2

3

4

5

6

50 40 30 20 10 0

λ t turb. therm. cond. [kW/mK]

2.5

800

T axial temperature [°C]

3.0

c p isobaric heat capacity [kJ/kgK]

Prt turb. Prandtl number [-]

Fig. 10. Axial profiles of cp , t , Prt and temperature (T) for the approach of calculating a variable turbulent Prandtl number denoted “Prt = Pet /Ret ” according to Eq. (18). These profiles illustrate how these quantities interact for calculating the turbulent heat transfer around the PCP of water.

Axial distance from nozzle [mm]

Fig. 11. Axial profiles of cp , t , Prt and temperature (T) for the approach of calculating a variable turbulent Prandtl number denoted “Scalar-variance model” [14]. These profiles illustrate how these quantities interact for calculating the turbulent heat transfer around the PCP of water.

due to the calculated turbulent heat transfer, which plays a crucial role in the cooling of such jets. Modeling heat transfer around the PCP at near-critical pressures is a challenging task, but a variable turbulent Prandtl number has proven to be an appropriate tool to improve numerical predictions. However, all approaches presented here were not able to completely remove the unphysical behavior of the axial temperature profiles around the PCP (temperature plateau). Deviations between the numerical results and the experimental values were still present, especially in the thermal far-field of the jet. 4. Conclusions Flow calculations of fluids in the supercritical state considerably above the critical point are feasible due to moderate fluid property variations. However, CFD simulations at near-critical conditions involve the fluid’s transition from sub- to supercritical and are therefore far more troublesome due to the strongly varying properties around the PCP. This work clearly shows that experimental validation of such flow simulations is an indispensable step when trying to develop an appropriate numerical model. The present study investigated round supercritical water jets discharged in a bath of subcritical water above the critical pressure. When temperature and pressure dependent thermo-physical properties (NIST database) and a constant turbulent Prandtl number were used in the simulations, the turbulent thermal conductivity t (turbulent heat transfer) was locally overestimated in all computational cells having temperatures close to the PCT. This leads to an unphysical flattening of the axial temperature profiles around the PCP, which was not evident in the experimental measurements. Several model approaches are presented and discussed to tackle these difficulties. The basic concept is the implementation of a

variable turbulent Prandtl number that increases around the PCP. Therefore, all considered approaches are able to reduce the peak value of t and hence elongate the hot supercritical region in the near-field of the jet compared to the simulation with a constant Prt value (0.7). Furthermore, the plateau in the calculated temperature profiles was reduced improving the agreement with the experiment in jet’s far-field. All presented model approaches (including Prt = 0.7) are able to predict the length of the hot supercritical region of the jet (in the near-field) satisfactorily. This means that the overall turbulent heat transfer between water in the gas-like state (supercritical water plume) and the surrounding cooling water in the liquid-like state is predicted reasonable well. Finally, it can be concluded that a locally varying turbulent Prandtl number addresses the difficulties of the underlying turbulence model in predicting the turbulent heat transfer around the PCP. However, all approaches presented here are not able to resolve the mentioned difficulties completely. Efforts are ongoing in our research group to further improve the used approaches for the application of flows at near-critical conditions. Due to its generic nature, the scalar-variance model and the formulation of Prt as a function of Pet /Ret are favored for further improvement and future investigations. Acknowledgements This work was funded by Swisselectric Research and the Swiss National Science Foundation (SNF). We gratefully acknowledge this financial support. Moreover, the authors would like to thank P. Feusi, B. Kramer and M. Meuli for their technical support. References [1] Y.M. Yin, X.Y. Lu, Effects of injection temperature on the jet evolution under supercritical conditions, Chinese Science Bulletin 54 (22) (2009) 4197–4204. [2] M. Oschwald, et al., Injection of fluids into supercritical environments, Combustion Science and Technology 178 (1–3) (2006) 49–100. [3] R. Branam, W. Mayer, Characterization of cryogenic injection at supercritical pressure, Journal of Propulsion and Power 19 (3) (2003) 342–355. [4] J. Sierra-Pallares, et al., Numerical analysis of high-pressure fluid jets: application to RTD prediction in supercritical reactors, Journal of Supercritical Fluids 49 (2) (2009) 249–255. [5] J.M.M. Barata, I. Gokalp, A.R.R. Silva, Numerical study of cryogenic jets under supercritical conditions, Journal of Propulsion and Power 19 (1) (2003) 142–147. [6] J. Bellan, Supercritical (and subcritical) fluid behavior and modeling: drops, streams, shear and mixing layers, jets and sprays, Progress in Energy and Combustion Science 26 (4–6) (2000) 329–366. [7] X. Cheng, T. Schulenberg, Heat transfer at supercritical pressures - Literature review and application to an HPLWR (FZKA 6609), Institut für Kern- und Energietechnik, Forschungszentrum Karlsruhe GmbH, Karlsruhe, 2001. [8] V. Yang, Modeling of supercritical vaporization, mixing, and combustion processes in liquid-fueled propulsion systems, Proceedings of the Combustion Institute 28 (2000) 925–942. [9] S. Mokry, et al., Development of supercritical water heat-transfer correlation for vertical bare tubes, Nuclear Engineering and Design 241 (4) (2011) 1126–1136. [10] T. Kim, Y. Kim, S.K. Kim, Numerical study of cryogenic liquid nitrogen jets at supercritical pressures, Journal of Supercritical Fluids 56 (2) (2011) 152–163. [11] N. Zong, et al., A numerical study of cryogenic fluid injection and mixing under supercritical conditions, Physics of Fluids 16 (12) (2004) 4248–4261. [12] ANSYS FLUENT Theory Guide, 2011. [13] ANSYS FLUENT User’s Guide, 2011. [14] K.W. Brinckman, W.H. Calhoon, S.M. Dash, Scalar fluctuation modeling for highspeed aeropropulsive flows, AIAA Journal 45 (5) (2007) 1036–1046. [15] N. Chidambaram, S.M. Dash, D.C. Kenzakowski, Scalar variance transport in the turbulence modeling of propulsive jets, Journal of Propulsion and Power 17 (1) (2001) 79–84. [16] M. Aouissi, A. Bounif, K. Bensayah, Scalar turbulence model investigation with variable turbulent Prandtl number in heated jets and diffusion flames, Heat and Mass Transfer 44 (9) (2008) 1065–1077. [17] A.N. Secundov, S.F. Birch, P.G. Tucker, Propulsive jets and their acoustics, Philosophical Transactions of the Royal Society A: Mathematical Physical and Engineering Sciences 365 (1859) (2007) 2443–2467. [18] N.J. Georgiadis, J.R. DeBonis, Navier–Stokes analysis methods for turbulent jet flows with application to aircraft exhaust nozzles, Progress in Aerospace Sciences 42 (5/6) (2006) 377–418.

M.J. Schuler et al. / J. of Supercritical Fluids 75 (2013) 128–137 [19] T.H. Shih, et al., A new kappa–epsilon eddy viscosity model for high reynoldsnumber turbulent flows, Computers and Fluids 24 (3) (1995) 227–238. [20] W. Wagner, A. Pruss, The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use, Journal of Physical and Chemical Reference Data 31 (2) (2002) 387–535. [21] C.R. Augustine, Hydrothermal spallation drilling and advanced energy conversion technologies for engineered geothermal systems, PhD Thesis, Massachusetts Institute of Technology, 2009. [22] T. Rothenfluh, M.J. Schuler, P.R. von Rohr, Penetration length studies of supercritical water jets submerged in a subcritical water environment using a novel optical Schlieren method, The Journal of Supercritical Fluids 57 (2) (2011) 175–182. [23] J. Sierra-Pallares, P. Santiago-Casado, F. Castro, Numerical modelling of supercritical submerged water jets in a subcritical co-flow, The Journal of Supercritical Fluids 65 (0) (2012) 45–53. [24] J. North, Improvements in or relating to drilling with gas liquid swirl generator hydrocyclone separation combustion thermal jet spallation., International PCT Patent, WO/1996/003566, 1996. [25] R.M. Potter, J.W. Tester, Drilling of Vertical Boreholes by Thermal Processes: Including Rock Spallation and Fusion, U.S. Patent No. 5,771,984, 1998. [26] F.W. Preston, H.E. White, Observations on spalling, Journal of the American Ceramic Society 17 (1934) 137–144.

137

[27] R.E. Williams, The thermal spallation drilling process, Geothermics 15 (1) (1986) 17–22. [28] E.W. Lemmon, M.L. Huber, M.O. McLinden, NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 8.0, Standard Reference Data Program, National Institute of Standards and Technology, Gaithersburg, 2007. [29] S.B. Pope, Turbulent Flows, Cambridge University Press, New York, 2000. [30] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics, 2nd edition, Pearson Education Limited, Edinburgh, 2007. [31] S.B. Pope, Explanation of turbulent round-jet-plane-jet anomaly, AIAA Journal 16 (3) (1978) 279–281. [32] T.P. Sommer, R.M.C. So, H.S. Zhang, Near-wall variable-Prandlt number turbulence model for compressible flows, AIAA Journal 31 (1) (1993) 27–35. [33] ANSYS FLUENT UDF Manual, 2011. [34] D.M. McEligot, M.F. Taylor, The turbulent Prandtl number in the near-wall region for low-Prandtl-number gas mixtures, International Journal of Heat and Mass Transfer 39 (6) (1996) 1287–1295. [35] A. Malhotra, S.S. Kang, Turbulent Prandtl number in circular pipes, International Journal of Heat and Mass Transfer 27 (11) (1984) 2158–2161. [36] W.M. Kays, Turbulent Prandtl number – where are we, Journal of Heat Transfer – Transactions of the ASME 116 (2) (1994) 284–295.