Numerical Investigation of Two-Phase Shock Waves in CO₂ Flows Using a Modified Hertz–Knudsen Model

In recent years, the utilization of supercritical CO₂ (sCO₂) as a working fluid in several industrial applications has gained significant attention due to its environmental advantages and unique thermodynamic properties [1–5]. Supersonic two-phase flows are encountered in various components, including separators [6,7], ejectors [8–10], thermal vapor compression systems [11–14], and sCO₂ compressors [15,16]. Accurately predicting the behavior of two-phase CO₂ in supersonic flows, including the presence of shock waves, is essential for optimizing the design and the efficiency of these components.

The separation of CO₂ from gas mixtures has become a critical concern in addressing global environmental issues and pursuing sustainable energy solutions. Carbon capture and storage systems enable the capture of CO₂ emissions generated by power plants or industrial facilities and their transportation to a storage site for safe and secure containment. Supersonic separators efficiently separate CO₂ from other gases in a supersonic expansion process and find extensive use, for instance, in the oil and gas industry [17] and in carbon capture and storage [18].

Supersonic ejectors are emerging as key elements to enhance the performance of CO₂ transcritical refrigeration cycles [19]. These cycles offer environmental advantages over chlorofluorocarbon and hydrochlorofluorocarbon refrigerants. This is due to CO₂ having a zero ozone depletion potential and low global warming potential. Ejectors can increase the cooling capacity of the system by effectively expanding the high-pressure, high-temperature supercritical CO₂ from the gas cooler to a lower pressure.

The formation of liquid droplets and their subsequent evaporation can significantly impact the performance of sCO₂ centrifugal compressors utilized in power cycles, particularly when operating near the saturation curve and the critical point. Any reduction in compressor performance has negative effects on the overall cycle efficiency, leading to increased fuel consumption and emissions.

Shock waves are characterized by sudden and significant changes in pressure, temperature, and velocity. They exhibit complex behaviors when propagating through multiphase flows, which can lead to various phenomena when interacting with CO₂ mixture, including phase-change, real gas effects, and complex wave structures [20]. Understanding and predicting these effects are essential in the aforementioned applications.

The rapid expansion of a flow near the saturation curve leads to thermodynamic nonequilibrium condensation, where the vapor phase becomes supersaturated, resulting in the formation of droplets that do not share the same temperature as the surrounding vapor phase. Nonequilibrium condensation models are required to simulate the nucleation and growth of the droplets in these flows [21–26]. In some cases, homogeneous equilibrium models can also provide a suitable solution, especially at higher pressures near the critical point, where the nonequilibrium effects are reduced [27,28].

In this study, we employ a nonequilibrium condensation and evaporation model based on the kinetic theory, also known as Hertz–Knudsen models [29], within the framework of computational fluid dynamics (CFD). Our aim is to explore the behavior of two-phase shock waves in CO₂ flows under various geometries and thermodynamic conditions near the critical point, capturing nonequilibrium effects associated with phase-change.

Additionally, we will present the key results obtained through our numerical simulations and draw conclusions that contribute to our understanding of two-phase shock waves in CO₂ flows and the boundary conditions that are relevant to the problem.

The research questions driving this study concern the shock-wave/boundary layer interaction and the abrupt evaporation of liquid droplets across shock-waves in two-phase CO₂ flows, especially within the thermodynamic region near the critical point. The objective of this work is to provide valuable insights into the numerical modeling of two-phase shock waves in CO₂ and their potential for improving the performance of components operating with supersonic flows. The findings of this research have the potential to advance the design and operation of these systems, contributing to more efficient and sustainable industrial processes.

Simulations of supersonic nozzle geometries involving two-phase CO₂ were conducted using Ansys Fluent 2022 R2. The simulations employed the multispecies user-defined real gas model integrated into the software. The thermodynamic and transport properties of CO₂ were determined using the Span and Wagner equation of state [30], accessed through Refprop [31]. These properties were defined within the solver using external look-up tables generated with an in-house Python code linked to Refprop. Pressure and specific enthalpy served as independent variables for the vapor phase, while only pressure was considered for the dispersed liquid phase, under the assumption of saturated conditions. Vapor properties were computed up to the spinodal limit, as depicted in Fig. 1 for density. The region between the saturation curve and the spinodal limit denotes the metastable region, which is crucial for investigating and simulating rapid expansions leading to homogeneous nonequilibrium condensation of supersaturated vapor flows.

Sections 2.1 and 2.2 will elaborate on the numerical framework and the phase change models employed to simulate nonequilibrium condensation resulting from rapid expansion and subsequent evaporation across the shock wave.

In this section, the nonequilibrium condensation model, which was previously validated in Ref. [34], is described. Additionally, the models for calculating mass transfer due to evaporation and denucleation of liquid droplets that occur across the shock wave are explained.

Two nozzle geometries were investigated in this study. The first geometry, introduced by Berana et al. [37], is a three-dimensional (3D) nozzle with a rectangular cross section, having a throat area of 0.40 mm². The convergent section was shortened compared to the real experimental geometry, ensuring that the Mach number at the inlet is very low. This resulted in a simulated total length of 13.5 mm. The divergent section profile has a constant slope of 0.48 deg with a flow cross section area at the outlet of 0.52 mm². The experiments conducted for this nozzle were utilized to validate the numerical model. The second nozzle, presented by Lettieri et al. [38], also has a rectangular cross section. However, in this study, the nozzle profile of Lettieri was used to create and investigate a two-dimensional (2D) axisymmetric geometry, commonly used in the design of supersonic components like supersonic separators. The simulated nozzle has a total length of 98.4 mm and a cross section area of 1.87 mm² at the throat and 3.16 mm² at the outlet. Illustrations of the geometries and grids can be found in Fig. 2, including details of the near-wall region and the cross section of the 3D geometry. It was ensured that the y⁺ at the wall remains below 30 for both grids. Attempts to further reduce the mesh size closer to the wall, aiming to resolve the viscous sublayer, posed convergence challenges due to the problem's complexity. A grid independence study was conducted using three grids with different levels of refinement in all spatial directions. In Table 1, the average liquid mass fraction at the outlet for Berana's nozzle are compared, and in Table 2, the shock location on the nozzle axis for Lettieri's nozzle are examined. The results indicate that grid independence is achieved with the middle-sized grids, with an absolute difference of 0.01% for the calculated outlet liquid mass fraction in Berana's nozzle and 0.1 mm for the calculated shock location in Lettieri's nozzle compared to the finest grid. The middle-sized grids are subsequently used in numerical simulations.

The results of the numerical simulations are presented in this section. Initially, the numerical and physical models were validated against experimental data using the 3D nozzle geometry from Ref. [37]. Following the validation, the same modeling approach was applied to investigate a 2D axisymmetric nozzle, adopting the profile from Ref. [38]. For the latter, various inlet thermodynamic conditions near the critical point and back pressures were explored.

The numerical and condensation models used in this study were previously validated in Ref. [34]. Furthermore, an additional experimental case is included in this paper to validate the evaporation model used for simulating phase-change across the shock wave. The experiments conducted by Berana et al. [37] involved a supersonic nozzle with supercritical CO₂, expanding the flow beyond the saturation curve and applying back pressures ranging between 36 and 61 bar to observe shock waves in the divergent section of the nozzle.

Experimental measurements were performed by measuring static pressures in the converging section and temperatures in the diverging section, using pressure gauges and thermocouples mounted along the centerline of the wall. Subsequently, saturated pressures in the diverging section were calculated based on the measured temperature. The simulated cases have a total temperature of 45 °C and a total pressure of 90 bar at the inlet boundary, with a reduced entropy of 1.13, and four different back pressures were set at the outlet boundary.

The validation results are presented in Fig. 3, showing the pressure distribution plotted along the streamwise direction. Here, the nozzle throat is located at x = 0 m. The pressure profiles resulting from the simulations were derived using the static pressure in the convergent section and the saturated pressure calculated from the temperature at the centerline of the nozzle wall for the divergent section. The pressure profile is well captured upstream of the throat but underestimated downstream of the throat, where condensation occurs. It is worth noting that wall temperature is highly sensitive to the turbulent boundary layer, which can be approximated as isobaric but not isothermal. This complexity makes the calculation of the saturated pressure more challenging and uncertain. However, the shock locations computed numerically align with the experiments, shifting upstream as the back pressure increases, even though the measured pressure rise is located upstream and not so abrupt.

To explore the effects of geometry type and inlet thermodynamic conditions, additional simulations were conducted. These simulations include a 2D geometry with the same nozzle profile and identical inlet boundary conditions as the validation cases. Both 3D and 2D geometries were also simulated with inlet total conditions set to 90 bar and 100 °C to prevent condensation throughout the expansion process, ensuring that the flow remained supercritical or superheated. The results, presented in Fig. 4, underscore the need for 3D simulations when simulating shock waves in both single-phase and two-phase flows in this nozzle. This is due to the interaction of the shock with the wall, which leads to the formation of a recirculation zone, inducing significant 3D effects, particularly in the corner of the cross section. Figures 5 and 6 illustrate this phenomenon through contour plots of velocity magnitude and liquid mass fraction for two cross section in the diverging section: one upstream of the shock and one downstream of the shock. The recirculation zone is clearly identifiable by the low velocities in the corner, resulting in higher temperatures and evaporation rates in this region.

The validated numerical model was employed to simulate two-phase flows in a 2D axisymmetric nozzle, adopting the profile from Ref. [38]. To investigate the impact of stagnation thermodynamic conditions at the inlet on the behavior of two-phase flow both upstream and downstream of the shock wave, four thermodynamic inlet points were selected near the critical point. These points are summarized in Table 3 in terms of temperature, pressure, and entropy. Additionally, they are depicted on a nondimensional temperature-specific entropy diagram in Fig. 7, superposed to the compressibility factor contours. The compressibility factor provides a quantitative measure of the nonideality of the gas in this region. The entropy and compressibility factor are calculated using Refprop, which solves the Span and Wagner equation of state [30] up to the spinodal limit. The simulations were conducted using three different back pressures (⁠$pout$⁠) set at 0.5, 0.6, and 0.7 times the inlet total pressure of each case, and the outlet points are also illustrated in Fig. 7 under the assumption of an isentropic process. The chosen back pressures enable the study of the generated shock waves at various Mach numbers, thereby exploring the effects of the shock intensity on the two-phase flow.

The pressure distributions calculated along the nozzle wall for each simulation are presented in Fig. 8. In the converging section, the pressures remain unaffected by the different thermodynamic conditions because the droplet nucleation occurs downstream of the throat. This is evident from the pressure rise due to the latent heat release to the vapor phase. However, the inlet conditions have a significant impact on the pressures in the diverging section upstream of the shock. When the inlet conditions are shifted farther from the critical point, the onset of condensation is moved downstream along the nozzle. This occurs because the so-called Wilson point, representing the thermodynamic conditions at which condensation takes place, is located deeper within the metastable region at higher specific entropies. Additionally, the saturation curve is crossed farther downstream in the nozzle, contributing to this shift. As the back pressure increases, the shock wave is consistently shifted upstream across all cases. Furthermore, the closer the operating conditions are to the critical point, the farther downstream the shock wave occurs. This trend is particularly evident for $pout/p0$ = 0.5 and $pout/p0$ = 0.6, although it is less pronounced for $pout/p0$ = 0.7. In this latter case, the shock waves generated are much closer to each other, with cases 1 and 2 almost superposed. It is important to note that, especially for case 1, the pressures across the shock exceed the critical pressure over a very small distance. This can lead to numerical anomalies in the physical model developed for phase-change occurring at subcritical pressures, potentially resulting in inaccuracies in the calculation of mass and heat transfer across the shock wave.

The liquid mass fraction and droplet density distributions along the nozzle axis for $pout/p0$ = 0.6 are shown in Fig. 9. In accordance with the pressure behavior, the nucleation and mass transfer due to condensation occur closer to the nozzle throat for cases with thermodynamic conditions approaching the critical point. Additionally, the liquid mass fraction and the number of droplets decrease moving from case 1 to case 4. In correspondence with shock waves, the resulting pressure and temperature rise drive the vapor into a superheated state, leading to the collapse of the droplets and subsequent mass transfer from the liquid to the vapor phase, that is evidenced by the sudden drop of the curves along the axial direction.

Figures 10 and 11 display the contours of liquid mass fraction and vapor Mach number for test cases 1 and 4 with the three back pressures. In situations where the flow contains a high wetness fraction, the liquid phase does not completely evaporate across the shock, and it ends up evaporating downstream of the shock in the compressive process within the divergent section. This phenomenon is evident in case 1 for $pout/p0$ = 0.5, where the liquid mass fraction upstream the shock reaches 22%, and the flow still retains a high wetness fraction even downstream of the shock. Notably, the Mach numbers are higher in case 4, given that the latent heat release due to condensation is smaller and occurs farther downstream in the supersonic flow, reaching a maximum value close to 1.55 for the lowest back pressure.

The contours of the mixture turbulent kinetic energy (TKE) and the pressure isolines in the region across the shock are shown in Fig. 12 for $pout/p0$ = 0.5. TKE provides a measure of the shock wave's intensity and its interaction with the nozzle wall, quantitatively estimating the turbulence generated in the recirculation zone. The TKE values calculated for case 4 are slightly higher than those for case 1, which is directly linked to the fact that the Mach numbers involved are larger. Additionally, Fig. 13 presents the mixture velocity contours and pressure isolines for the same fluid region and cases. Even though the calculated chocked mass flowrate is 1.5 times larger in case 1 than in case 4, the velocity magnitudes are higher in case 4 because the flow density is significantly lower, approximately 0.6 times the density in case 1. The isobars across the shock exhibit similar behavior in both cases, even though the pressure gradients at the shock front are slightly larger in case 1.

In this study, a numerical model was developed that presents an advancement in the field of thermodynamic nonequilibrium condensation and evaporation simulations in real gases, particularly in the context of simulating two-phase shock waves in supersonic nozzle geometries. The model includes the nucleation of liquid droplets, which is based on a corrected version of the classical nucleation theory, and the mass transfer rate due to phase changes, modeled using a modified Hertz–Knudsen equation.

Validation of the model against experimental data demonstrates its reliability in capturing the fundamental physics involved in shock wave phenomena combined with two-phase flows. A key insight from this research is the importance of conducting 3D computational fluid dynamics simulations when simulating shock waves in geometries with rectangular cross section. The interactions of the shock with the walls generate 3D flow effects that cannot be captured with 2D geometries.

Furthermore, the application of the developed numerical model to a 2D axisymmetric geometry, considering varying thermodynamic inlet conditions and nozzle back pressures, highlights the influence of these parameters on the location and intensity of the shock waves, as well as the behavior of immediate evaporation and denucleation of liquid droplets.

In conclusion, this paper contributes to the field of two-phase shock wave modeling by providing a robust numerical tool that can aid in the design, analysis, and optimization of systems where nonequilibrium condensation and evaporation play a significant role. Additionally, future work could explore the impact of turbulence models on the intensity and position of two-phase shock waves.

The authors gratefully acknowledge the financial support of the Doctoral School of Lappeenranta-Lahti University of Technology LUT.

• Lappeenranta University of Technology (Funder ID: 10.13039/501100004105).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

h =

specific enthalpy (J/kg)

J =

nucleation rate (1/(m³ s))

k =

thermal conductivity (W/(m K))

K =

turbulence kinetic energy (J/kg)

k_b =

Boltzmann constant (J/K)

Kn =

Knudsen number

m =

molecular mass (kg)

mfp =

mean free path (m)

p =

pressure (Pa)

Pr =

Prandtl number

Q_conv =

convective heat transfer (J/(m³ s))

R =

gas constant (J/(kg K))

$r¯$ =

mean droplet radius (m)

$r˙$ =

droplet growth rate (m/s)

s =

specific entropy (J/(kg K))

S =

supersaturation ratio

S_q =

source term in a transport equation for scalar q

t =

time (s)

T =

temperature (K)

u =

velocity vector (m/s)

u,v,w =

velocity components (m/s)

x,y,z =

Cartesian coordinates (m)

z =

compressibility factor

Greek Symbols

$α$ =

volume fraction

$η$ =

droplet density (1/m³)

$μ$ =

dynamic viscosity (Pa·s)

$μT$ =

dynamic eddy viscosity (Pa·s)

$ρ$ =

density (kg/m³)

$σ$ =

planar surface tension (N/m)

$τ¯¯$ =

stress tensor (N/m²)

$τ¯¯T$ =

Reynolds stress tensor (N/m²)

Subscripts

c =

critical point

denuc =

denucleation

eff =

effective

fg =

latent heat

int =

interfacial

l =

liquid

m =

mixture

nuc =

nucleation

s =

saturation

v =

vapor

vs =

saturated vapor

0 =

total quantity

Acronyms

CFD =

computational fluid dynamics

CNT =

classical nucleation theory

sCO₂ =

supercritical carbon dioxide

TKE =

turbulent kinetic energy