Role of Inlet Boundary Conditions on Fuel-Air Mixing at Supercritical Conditions

Fuel injection into the combustor can be described by two extreme phenomena. At low ambient pressures, typically below the thermodynamic critical condition of the fuel, a classical situation exists where a well-defined interface separates the injected liquid fuel from the surrounding gaseous fluid due to the presence of the surface tension. Under these conditions, surface tension forces form a discontinuous interface that promotes primary atomization, secondary breakup, and the resultant spray phenomenon that has been well characterized and is widely assumed. As the ambient pressure exceeds the critical pressure of the injected fuel, a distinct gas–liquid interface might not exist depending upon the local fluid temperature. Instead, the injected fuel undergoes a continuous change in density as it mixes with surrounding gas if the local mixture temperature exceeds the local mixture critical temperature. In such a case, the surface tension is diminished, and the fuel atomization is suppressed. Instead, diffusion dominated mixing processes without fuel atomization occur. Fuel injected at such conditions has been observed to evolve in the presence of large (but continuous) thermophysical gradients, with no drops present in the flow field. However, the transition of classical spray atomization processes to single-phase continuous dense-fluid mixing dynamics with diminished surface tension forces has not been fully understood.

Advanced engine design and alternative fuels present the possibility of fuel injection at purely supercritical conditions in diesel-like applications and gas turbines, as is already the case for rocket engines. Oefelein [1] identified seven topics in need of basic research to develop advanced predictive capabilities in high pressure, multiphase reactive flows. Among these, the topic of “Supercritical fluid mixing and combustion” includes “Real-fluid nonideal thermodynamics,” “State equation correlations and closure,” and “Pressure dilatation and compressibility effects” as important subtopics [1]. Many simulations of supercritical fuel injection have been presented in the literature, with the majority of them assuming one of the cubic equations of state and the associated mixing law [2–9]. Proper treatment of the inlet boundary conditions is an important aspect of modeling that has not been investigated in past studies and requires further attention.

For traditional liquid fuel injection into a gaseous ambient, the two phases are treated separately in the Lagrangian–Eulerian approach. Therein, liquid droplets are treated as a discrete phase, and their trajectories are tracked in the Lagragian reference frame, while the surrounding gas is treated as continuum, and is governed by the conservation laws in the Eulerian reference frame. Drag force and heat flux source terms associated with each reference frame account for momentum and energy transfer between the discrete and continuum phases [10–12]. The Lagrangian–Eulerian formulation is adequate for mixing with two distinct phases, where liquid fuel breakup into evaporating droplets is driven by surface tension, and the vapor phase is nearly an ideal gas because of the high ambient temperature compared with the temperature of the injected liquid.

For supercritical injection, a purely Eulerian approach can be used, where all the associated fluids are treated as a continuum, and all of the transport phenomena are described by the conservation laws in the Eulerian reference frame, where each local mixture has homogeneous properties calculated by the equation of state (EOS) and fluid mixing law. This approach is effective and well justified when the conditions throughout the domain remain above the critical pressure and/or temperature of the fluid mixture throughout the mixing process, such that the surface tension forces no longer form a discontinuous phase transition. In such a framework, no momentum and temperature coupling between different working fluids must be considered, and the species mixing is controlled by convection, as well as molecular and turbulent diffusion.

For Eulerian simulation of fuel at supercritical conditions, a real-gas equation of state must be used to account for the compressibility effects that occur inside and directly after the injector. The most popular real-gas equations of state for computational fluid dynamics (CFD) are the cubic equations of state, which comprise a family of equations based on specific properties for gases at their respective critical points. Cubic equations of state assume a uniform behavior among all fluids, accounting for differences between fluids using two to five constant property values. Because of their simplifying assumptions, some degree of error is inherent in any cubic equation of state formulation, and that error is compounded with the mixing of multiple fluids. With mixing, the cubic equations require that the mixture be treated as a single pseudo-fluid, with critical properties calculated as a function of the mixture fraction and of the critical properties of the respective fluids. Despite inherent limitations to accuracy, cubic equations of state are popular in CFD because they are relatively efficient to compute and are usually sufficiently accurate, particularly if no species mixing occurs in the domain. More accurate equations of state, such as the Benedict-Webb-Rubin equation and models based on Helmholtz energy, can account for nonidealities introduced by mixing of fluids, can work together with chemical kinetics mechanisms for combustion to improve accuracy, and can be tuned to fit the experimental data [13]. These formulations have been used in CFD simulations of fuel injection and other simulation applications, but they increase the required computational time dramatically [1,14–17]. Due to the increased computational cost, particularly as more molecular species are added, fuel injection simulations that include combustion models using many species and multiphase simulations usually make complex EOS formulations infeasible, so simpler cubic EOS formulations are more popular [18–20].

Boundary conditions for fuel injection are critical to accurate simulation. Even so, the internal geometry of the injector itself is often omitted in favor of assumed inlet (or nozzle exit) boundary conditions of velocity or mass flux, temperature, and turbulence levels. However, the effect of different boundary conditions on CFD predictions has not been well documented in the literature. The primary objective of this study is to implement comprehensive real-gas and fluid mixing models to predict fuel-air mixing at purely supercritical conditions. Thus, the fuel at supercritical pressure and temperature is injected into an ambient at supercritical pressure and temperature, to ensure no surface forces are present and the mixing process is determined only by convection and diffusion calculated using the Eulerian approach. First, results are compared with CFD predictions of Qiu and Reitz [9] to verify the present implementation of real-gas and fluid mixing models. Next, experimental data presented by Roy et al. [21] are used to investigate the effect of the boundary conditions and to show how prediction could be improved by the knowledge of the correct inlet boundary conditions.

The governing equations are solved in an axisymmetric computational domain to match with Qiu and Reitz [9]. Figure 1 displays the domain, the boundary surfaces, and the mesh. The domain has a radius of 40 mm and length of 120 mm. The fluid is injected at the center through a 1 mm radius nozzle, or in the second case in this work, a 2 mm radius pseudo-nozzle. A constant pressure is specified at the outlet boundary, and all other surfaces of the chamber are insulated with no-slip walls. The injector radius is represented by 45 uniform grids, and 200 divisions span the remaining chamber radius, resulting in a total of 245 grids in the radial direction. The axial direction is represented by 200 divisions, resulting in a total of 98,000 grids as compared with a total of 2100 grids used by Qiu and Reitz, which they found to provide sufficiently grid independent results [9].

Qiu and Reitz [9] sought to replicate the experimental results of Roy et al. [21], who used planar laser induced fluorescence (PLIF) in an optically accessible pressure chamber to study supercritical mixing of fluoroketone (C₆F₁₂O) injected into a quiescent nitrogen environment. In the selected case from Roy et al. [21], the reduced temperature (temperature normalized by the critical temperature of fluoroketone) at the injector exit and in the chamber is 1.08 and the reduced pressure (pressure normalized by the critical pressure of fluoroketone) is 1.5 at the injector exit, and 1.41 in the chamber. The thermodynamic and flow properties of fluoroketone are not available in the literature. Thus, Qiu and Reitz [9] used perfluorohexane (C₆F₁₄) as surrogate fluid, noting that previous work had identified the critical point, latent heat, and surface tension at normal conditions to be similar for perfluorohexane and fluoroketone [23].

In this study, perfluorohexane is used as the fuel surrogate, as it was for Qiu and Reitz [9]. The critical temperature, critical pressure, acentric factor, and molecular weight of perfluorohexane, required for the Peng–Robinson equation of state, were set to the values reported by Qiu and Reitz [9]. However, the critical specific volume was not reported, and so was obtained from Refprop 10.0 [24]. Likewise, the transport properties and the specific heat of perfluorohexane were obtained from Refprop for a range of temperatures at the relevant pressure, and curves were fitted to that data [24]. A similar procedure was also performed for N₂, using the same Peng–Robinson parameters as reported by Qiu and Reitz [9], although transport properties and specific heat were set to default correlations instead of Refprop-derived correlations [9,24]. Mixing rules for transport properties and specific heat were set to the default kinetic theory-based mixing laws, and the rules for determining the Peng–Robinson pseudo-fluid properties for a given mixture fraction were set to the one-fluid van der Waals mixing law.

Figure 2 shows the density contours for perfluorohexane with respect to pressure and temperature. Each successive line represents an increase of 50 kg/m³, but some labels are skipped for clarity. The saturation curve is shown starting from the lower left and extending up into a high concentration of lines, where to the right, perfluorohexane is a vapor, and to the left, it is a liquid. The end of the saturation line represents the critical point, where the density is highly sensitive to the changes in the pressure and temperature, as expected. Similar curves showing sharp gradients in the region adjacent to saturation line and critical point were obtained for other thermophysical properties including specific heat, thermal conductivity, and viscosity.

Qiu and Reitz assumed top hat profiles for velocity (14.47 m/s) and temperature at the injector exit. Using critical properties of the perfluorohexane based on the reduced temperature of 1.08 and reduced pressure of 1.41, the chamber temperature and pressure were 485.5 K and 2.468 MPa, respectively, and the fuel temperature was 485.5 K [9,21]. Figure 3 depicts that the predictions by Qiu and Reitz [9] had large deviations from experimental measurements of density. Qiu and Reitz noted that some shortcomings of their formulation include the RANS model’s inadequacy for capturing wavy structures, and the smaller spreading angle relative to the experiment [9]. However, a detailed explanation of these large deviations remains unknown.

In this work, the difficulty in accurately replicating experimental data on supercritical mixing is examined using a two-step approach. First, real-gas and fluid mixing models implemented in this study are used to replicate previous CFD results and to ensure that the present models are consistent with those of Qiu and Reitz [9]. Next, the developed models are used to assess the applicability of different boundary conditions in reducing discrepancies with experimental data. Finally, it is hypothesized that the experimental data itself is subject to significant measurement uncertainties.

Simulations were performed for pressures and temperatures at the injector exit and in the chamber taken from Qiu and Reitz [9]. First, for the verification of the CFD models, a top hat velocity profile was used at the injector exit. Qiu and Reitz [9] report a Reynolds number of 8.1 × 10⁵, which would indicate a turbulent intensity of 2–3%. In this study, the turbulence intensity at the injector exit was fixed at 2.5%, and the length scale was assumed to be the diameter at the injector exit. Each simulation was run for an injection duration of 10 ms, matching Qiu and Reitz [9]. Mass flowrate was calculated as 0.01761 kg/s, which is within 1.5% of the reported value of 0.01735 kg/s from Roy et al. [21].

Figure 4 presents the computed density contours and reveals the similarity of these results with those by Qiu and Reitz [9] shown in Fig. 3(b). Although the present results match those by Qiu and Reitz [9], neither replicates the experimental results shown in Fig. 3(a). Turbulence intensity could be one source of error, as hypothesized by Qiu and Reitz [9]. In this study, low and high turbulence intensity levels were investigated, and grid effects were examined through simulation using a coarse grid with equivalent resolution to Qiu and Reitz [9]. For a quantitative comparison, radial profiles of density are compared in Fig. 5 at two axial locations x = 15 mm and x = 35 mm for five cases. The first two cases pertain to this study with turbulence intensity of 2.5% and 0%, respectively. The third case is similar to the 2.5% turbulence intensity case, but with the coarse grid. The fourth case is the result from Qiu and Reitz [9], and the fifth case is the experimental data from Roy et al., as presented by Qiu and Reitz [9,21]. These two axial locations were selected because detailed information was provided by Qiu and Reitz [9].

Figure 5 shows that without turbulence at the inlet, the density is significantly overpredicted at the centerline, particularly at the selected downstream locations. The choice of 2.5% for turbulent intensity is seen to be reasonable, although the centerline density is slightly overpredicted at x = 15 mm, and slightly underpredicted at x = 35 mm. For higher turbulence intensity at the inlet, the density is expected to drop too quickly along the centerline. Also shown is the marginal difference in density at the selected locations between the full resolution and coarse resolution cases for 2.5% turbulence intensity.

For the remainder of this work, the 2.5% turbulent intensity case is referred to as the top hat profile case. At x = 15 mm, the density from Qiu and Reitz [9] is significantly higher at the center than either of the top hat profile cases from this work, but at x = 35 mm, the center point density from Qiu and Reitz [9] is lower than the case from this work that has no turbulence at the injector exit [9]. From this information, we infer that Qiu and Reitz [9] assumed a turbulence level at the inlet. The higher density at the inlet for Qiu and Reitz [9] is attributed to the differences in the method to calculate the critical specific volume: an unknown source for Qiu and Reitz [9], and Refprop in this work [24]. Nevertheless, the general trends of the present results match well with the previous work, and thus, they can be used for parametric analysis presented in Sec. 3.2.

With the computational framework and domain validated as such, a second case is introduced to improve upon the results relative to the experimental data. The experimental data given in Roy et al. [21] apparently do not start at the injector exit, and instead start somewhere downstream, presumably due to optical limitations of the experimental setup. For this reason, direct comparison with CFD is difficult as the exact axial location for the inlet boundary is unknown. The PLIF method used by Roy et al. [21] is also subject to uncertainties, as evidenced by the fact that, near the nozzle exit, radial profiles of the reported density are nonmonotonic, with values at the center being lower than those in the immediate surroundings, which is clearly not physical.

To replicate the experimental results as closely as possible, the experimental density data from Roy et al. [21] is used to determine the mass fraction of fuel as a function of radial location at the earliest reported axial location, and a velocity profile is determined by assuming that the velocity and mass fraction profiles are similar. The radius at the inlet boundary (downstream of the injector exit as explained above) is adjusted to encompass the width of the profile and supply the known mass flowrate. The mass fraction at the center point is assumed to be 1.0 to rectify the experimentally presented nonphysical density at center being less than that of the immediate surroundings. Away from the center, after the point where the experimental density monotonically decreases in the radial direction, the density is matched with the experimental data by assuming a mass fraction profile at the temperature and pressure given in the boundary condition, and this profile is connected to the value of 1.0 at the center by a simple spline.

Note that the experimental data are asymmetrical, and the center of the jet is not at the centerline of the presented graph (Fig. 3). Thus, the density profile at the inlet is obtained by locating the true centerline and averaging the density profiles on either side. Because the inlet in this simulation is located downstream of the physical injector exit, a constant total inlet mass flowrate is used instead of a velocity profile, and the profile shape of the mass flux is determined by assuming that the axial velocity and mass fraction profiles are similar. This approach results in a mass flux profile equivalent to the profile obtained from the product of experimentally measured density and imposed mass fraction at the inlet, and also fixes the velocity profile to be equivalent in shape to the mass fraction profile of the injected fluid. Figure 6 compares the corrected profile to the top hat profile at the domain inlet. The boundary condition with gradually changing mass fraction profile derived from the experimental data creates a well-developed density profile at the inlet. As before, the turbulence intensity at the inlet is set to 2.5%, with length scale equal to the original injector diameter. With these modified boundary conditions, the results in Fig. 7 are obtained.

In Fig. 7, the spreading angle of the jet is not apparently significantly different, but the overall profile matches the experimental results in Fig. 3(a) much better at almost every point than the previous test cases [9,21]. Qiu and Reitz noted that the spreading angle in their simulation appeared to be too low, but the fact that the spreading angle is not significantly increased in the corrected inlet boundary profile case, while the density contour profile matches more closely in overall shape, indicates that the problem is not so much underprediction of the spreading angle, as Qiu and Reitz theorized, as it is imposition of a boundary condition inconsistent with the experimental data. For a quantitative comparison, Fig. 8 shows the radial profiles of density for the top hat profile case with 2.5% turbulence, the experimentally corrected inlet profile case (as described in this section), the experimental data [21], and the top hat case from Qiu and Reitz [9] at x = 15 mm and x = 35 mm.

The radial profiles of density shown in Fig. 8 at x = 15 mm demonstrate the advantage in accuracy and profile shape of obtaining the boundary condition as directly as possible from the experimental data. The shape of the density profile in both top hat cases is too narrow, overpredicting density significantly at the center, and then underpredicting density after approximately 0.8 mm radially. In contrast, the profile predicted by corrected boundary conditions closely matches the downstream density at the center. Although density for most of the remaining mixing layer is still underpredicted, the prediction is much closer to the experimental result than either top hat case.

At x = 35 mm, all CFD cases underpredict density, except for the top hat case of Qiu and Reitz, which significantly overpredicts density near the center [9]. While none of the CFD results match the experimental data well, the case with the corrected inlet profile of this section provides a significantly closer match with experimental data. The density at the center for both axial locations shown has a much better match with the experiment, indicating that the rate at which density decreases along the centerline is consistent, unlike the two top hat profile cases.

Figure 8 shows that the top hat profile case from this effort overpredicts density at the center at x = 15 mm and underpredicts it at x = 35 mm, indicating that the predicted rate of decrease in density is too high. This discrepancy could be reduced by decreasing the turbulent intensity at the injector exit, but that would also make the overall agreement worse. Further, the turbulent intensity is already less than the predicted intensity based on the experimental Reynolds number. For all cases, this discrepancy in the rate at which the center point density decreases, is also clear in the contour plots in Figs. 3, 4, and 7.

In Figs. 3 and 7, respectively, center point density in the experimental data and the corrected profile cases both apparently decrease more quickly at first, up to around x = 30 mm axially, while the dense core in both of the two top hat cases in Figs. 3 and 4, respectively, is still at a higher density up to that point. Further downstream in these contour plots, the center point density in the top hat cases decreases significantly and quicker than the experimental data, whereas the corrected profile case matches relatively well in the whole domain. If the effective impact of using a corrected profile boundary condition is simply adjusting the inlet boundary to a downstream location to match the experimentally reported field of view, we should expect the top hat case to match the inlet of the corrected case at some near-field location, which clearly is not the case. These results indicate that in the near-field of the injector either the model does not capture the key physics or the measurement uncertainty is too high. Although this discrepancy is unresolved, the top hat and corrected profile cases are more similar at x = 35 mm than at x = 15 mm, and generally agree better at downstream locations.

Velocity also plays an important role in determining the density field. Figure 9 shows radial profiles of axial velocity at x = 15 mm and x = 35 mm. The axial velocity in the center region for the top hat profile case is nearly twice that for the corrected profile case. In spite these large differences, the density in the center region is similar for the two cases since the actual mass flow in the center is a small fraction of the total mass flow. Figure 10 shows axial profiles of density and axial velocity along the centerline for the top hat and corrected profile cases. The sharp drop in density and velocity at the downstream edge indicates that the jet has reached x = 60 mm for the corrected profile case and x = 80 mm for the top hat profile case, which is expected because the top hat profile case pertains to higher axial velocities.

Good agreement between density from experiments and the corrected profile case of the present study indicates that the present implementation of the inlet boundary conditions is a significant step forward compared with the top hat profile case. The results highlight the need to choose the correct variable to match when determining the inlet boundary conditions from the experimental data. A conserved variable such as the mass flowrate (as in the present study) might be more accurate to match with the experimental data compared with the derived variable such as velocity, when it is based on an assumed density and mass flowrate.

Supercritical fuel injection is an increasingly relevant field of study and presents new challenges to simulation efforts. This work implemented real gas models and introduced two simple CFD cases as a starting point to develop a robust computational framework and demonstrated the importance of matching experimental data in specifying the inlet boundary conditions. The test cases were designed to replicate experimental data for supercritical fluid injection by Roy et al. [21]. The first case was similar to a previous effort to replicate the experimental data by Qiu and Reitz and showed reasonable agreement with the previous CFD results [9]. The second case improved upon the boundary conditions by adjusting the mass fraction profile at the inlet to match the experimentally obtained density. Predictions from this second case matched the experimental density profile more closely than either of the top hat cases and indicate the need to infer boundary conditions directly from experimental data. For supercritical injection, the strong dependency between temperature, density, mixture fraction, and velocity in Eulerian reference frame requires careful consideration of boundary conditions, and conventional methods are not appropriate. Future work will include a study of propane fuel injection and the effects of decoupling the fuel injector from its internal geometry at the nozzle exit, as well as further experimental validation.

This research was supported in part by the Office of Energy Efficiency and Renewable Energy Grant DE-EE0007301.

There are no conflicts of interest.