1. Introduction
At a certain mixture-dependent pressure and temperature, a given thermodynamic system reaches a singularity called the critical point (Mayer et al. Reference Mayer, Schik, Vieille, Chaveau, Gökalp, Talley and Woodward1998; Bellan Reference Bellan2000). Here, the enthalpy of vaporisation and interfacial tension both vanish, with the result that two-phase boundaries thicken to the continuum limit (Prausnitz, Lichtenthaler & de Azevedo Reference Prausnitz, Lichtenthaler and de Azevedo1998; Dahms & Oefelein Reference Dahms and Oefelein2013). Hence beyond the critical point, a two-phase solution for a given system is impossible to achieve.
Though the field achieves some simplicity in the phases, the system remains thermodynamically complex. For instance, strong nonlinear variations in the thermodynamic and transport properties are found at thermodynamic states near the critical point, and the system exhibits significant susceptibility to perturbations in the state as a result. A visualisation of this behaviour for a system composed of pure carbon dioxide is provided in figure 1.
Selected thermodynamic property variations as functions of pressure and temperature for carbon dioxide. The isobar identified by the solid red line in each graph clearly captures the order-of-magnitude change in these properties exhibited near the critical point. Variations of similar magnitudes are found in all thermodynamic and transport properties. In each graph, broken black lines represent relaxations of these variations at different conditions, while ideal gas behaviour is identified by the solid black line in each case.
The large variations shown in figure 1 occur when the thermodynamic state transitions across the so-called pseudo-boiling line (sometimes referred to as the Widom line), which is an extension of the liquid–vapour saturation line. It is defined based on the locus of points where some thermodynamic property achieves a characteristic value (e.g. loci of maxima in isobaric heat capacity). As this transition occurs, the fluid switches between taking on ‘gas-like’ and ‘liquid-like’ characteristics. This behaviour is unique to supercritical fluids.
The strong variations in thermodynamic and transport properties as the Widom line is crossed couples with, and influences, fluid dynamic motions. Hence turbulent flows of supercritical fluids contain not only the complexities associated with classical turbulence, but also those arising from supercritical nonlinear thermophysics. For instance, supercritical turbulent flows at thermodynamic states near the critical point are compressible (see figure 1 c), and contain variations in momentum and thermal diffusivities that are absent in ideal gas flows. Figure 2 shows this behaviour graphically.
Variation in selected transport properties for CO
$_2$
: (a) dynamic viscosity; (b) Prandtl number.
These considerations make analysis of supercritical turbulence challenging. However, these flows are of direct engineering relevance to propulsion and power platforms such as cryogenic rocket injection systems, as well as diesel engines. Much research and development has already occurred in this area, using both experimental and computational techniques (see Oefelein & Yang (Reference Oefelein and Yang1998),Yang (Reference Yang2000), Mayer et al. (Reference Mayer, Schik, Schäffler and Tamura2000) and Branam & Mayer (Reference Branam and Mayer2003) for pertinent cases). Recently, ground power systems have been developed that are able to successfully exploit the supercritical nonlinear thermophysics shown in figures 1 and 2 to achieve higher cycle efficiencies than traditional ground power systems. For instance, oxy-combustion has historically been challenging to implement in practical systems given the high resultant turbine inlet temperatures that occur due to the absence of the diluent N
$_2$
species in the flue gases. Injection of supercritical CO
$_2$
into the combustion chamber itself allows the turbine inlet temperature to be successfully modulated to minimise thermal fatigue on the turbine blades. This is possible due to the peak in isobaric heat capacity found near the critical point of CO
$_2$
, which is shown in figure 1(b). Additional benefits include easier carbon capture and hence lower emissions overall. This example provides insight regarding the practical relevance and promise of supercritical fluids to both current and next-generation engineering devices. Fundamental understanding regarding the dominant driving mechanisms of turbulent flows of such fluids is still required, however, and this fact provides the motivation for the present research.
Early experimental research in supercritical fluid dynamics was focused primarily on establishing the absence of two phases at supercritical conditions, and showing that supercritical shear layers and jets behave similarly to their variable density gaseous counterparts (Newman & Brzustowski Reference Newman and Brzustowski1971; Mayer & Tamura Reference Mayer and Tamura1996; Mayer et al. Reference Mayer, Schik, Vieille, Chaveau, Gökalp, Talley and Woodward1998). Later research employed Raman spectroscopy to examine canonical configurations (Oschwald & Schik Reference Oschwald and Schik1999; Chehroudi, Talley & Coy Reference Chehroudi, Talley and Coy2002). In particular, the research of Oschwald & Schik (Reference Oschwald and Schik1999) and Chehroudi et al. (Reference Chehroudi, Talley and Coy2002) used a single-component fluid to focus on mixing characteristics in the absence of complexities such as multiple critical points associated with multiscalar mixing and chemical reactions. This is the approach preferred in this study, since the problem remains tractable yet appreciably complex through this choice.
The aforementioned experimental research has provided valuable insights regarding jet divergence angle as a function of the injector configuration, chamber density ratio and chamber pressure. A particularly insightful contribution is that of Chehroudi et al. (Reference Chehroudi, Talley and Coy2002), which includes quantitative verification that supercritical mixing is akin to variable density gaseous mixing. The verification was done by comparing supercritical jet and shear layer growth rates with the classical gas-based theory of Dimotakis (Reference Dimotakis1986), that of Papamoschou & Roshko (Reference Papamoschou and Roshko1988), and the seminal incompressible variable-density shear layer data of Brown & Roshko (Reference Brown and Roshko1974). The variable-density shear layer experiments of both Brown & Roshko (Reference Brown and Roshko1974) and Papamoschou & Roshko (Reference Papamoschou and Roshko1988) in particular provided valuable insight regarding the weak dependence of shear layer growth rates on the density ratio between the two streams. Compressibility effects were found to more appreciably influence shear layer growth than both density and velocity ratios (Papamoschou & Roshko Reference Papamoschou and Roshko1988). This concept was further refined by Dimotakis (Reference Dimotakis1989), whereby shear layer mixing is characterised in terms of an entrainment ratio as opposed to velocity and density ratios.
Despite these contributions and the confirmation that supercritical shear layers are found to behave quantitatively as variable density layers, there is a lack of data particularly at near-critical conditions. At these conditions, the mixing will occur between thermodynamic states that bracket the Widom line, which exposes the nonlinear behaviour seen in figures 1 and 2. Given that such conditions are often encountered in supercritical devices of engineering interest, this research effort aims to provide a quantification of the mixing dynamics at precisely these conditions.
Computational analyses of supercritical fluid dynamics have occurred over approximately the past two decades. This work has used the direct numerical simulation (DNS) method to investigate fundamental questions, and assess and extend existing turbulence closures to enable accurate large-eddy simulations (LES) of supercritical fluid flows. Multiple DNS calculations of a temporally evolving shear layer were conducted by Bellan and colleagues (see e.g. Miller, Harstad & Bellan Reference Miller, Harstad and Bellan2002; Okong’o & Bellan Reference Okong’o and Bellan2002; Okong’o et al. Reference Okong’o, Harstad and Bellan2002; Bellan Reference Bellan2006). These studies uncovered the existence of regions of high density gradient magnitude (HDGM) within the mixing region. These regions are initially introduced to the field by the initialisation of the temporal mixing layer, and are a manifestation of the anisotropic nature of the shear layer itself. As the layer develops, these HDGM regions become distorted and are exacerbated by preferential species diffusion through the layer. This finding is also corroborated by Zong et al. (Reference Zong, Meng, Hsieh and Yang2004). The HDGM regions are found to contribute significantly to the development of turbulence anisotropy within the field, and are a distinctly unique phenomenon to supercritical mixing layers. The regions of HDGM act similarly to local material boundaries in a coordinate system defined by the local interface normal and tangent vectors. This dampens turbulent fluctuations normal to the HDGM region, and the turbulent energy in this component is instead preferentially distributed to the remaining two. This in turn causes the HDGM regions to act in a positive-feedback loop, whereby they amplify existing field anisotropy. These regions are also known to form in single-species supercritical shear layers when there are differences in the stream temperatures. Bellan (Reference Bellan2006) remarks that this renders the study of isotropic turbulence for supercritical combustion applications largely immaterial, and motivates the use of a spatially evolving shear layer in this research.
The DNS databases allow careful examination of LES subfilter models, and significant attention has been devoted to this effort by the computational research community in recent years. In particular, the research has been focused on answering an open question regarding the validity of certain LES modelling assumptions under supercritical conditions. The LES modelling assumptions are those that imply that the filtered values of secondary quantities (e.g. mass diffusion) can be calculated equivalently using individually filtered primary quantities. These assumptions break down for nonlinear functions such as those shown graphically in figures 1 and 2. The assumptions are found to be justifiably negligible under low-pressure ideal gas settings; however, they break down under supercritical conditions (Selle et al. Reference Selle, Okong’o, Bellan and Harstad2007; Ribert et al. Reference Ribert, Domingo, Petit, Vallée, Blaisot and Bellan2020). The invalidity of these assumptions in high-pressure settings necessitates the use of additional subfilter models to close the filtered system of equations that govern an LES flow field. In particular, terms related to subfilter pressure gradients, heat flux divergences and subfilter density were found to require treatment (Selle et al. Reference Selle, Okong’o, Bellan and Harstad2007; Ovais, Kemenov & Miller Reference Ovais, Kemenov and Miller2021; Unnikrishnan et al. Reference Unnikrishnan, Huo, Wang and Yang2021, Reference Unnikrishnan, Oefelein and Yang2022). Various types of closure have been proposed to treat these terms, ranging from Taylor series expansions (Selle et al. Reference Selle, Okong’o, Bellan and Harstad2007) to presumed
$\beta$
probability density function distributions (Lapenna & Creta Reference Lapenna and Creta2017).
This research aims to address some of the existing knowledge gaps regarding supercritical fluid mixing in a joint fashion. An experimental platform is constructed for this purpose, and results are gathered at a supercritical thermodynamic state near the critical point. This represents a unique and novel setting of engineering interest. It provides a way to characterise turbulent mixing both qualitatively and quantitatively at conditions not studied in past research. Experimental data are collected using state-of-the-art Raman spectroscopy. The experimental results are complemented by corresponding DNS, which provides further details regarding the dynamics of turbulence at these conditions.
2. Methodologies
2.1. Selection of operating conditions
Many experimental efforts documented in the literature have focused on either transcritical injection or supercritical injection at states far away from the critical point (e.g. Oschwald et al. Reference Oschwald, Smith, Branam, Hussong, Schik, Chehroudi and Talley2006; Mayer et al. Reference Mayer, Telaar, Branam, Schneider and Hussong2003). Experiments at such conditions were motivated either by cryogenic injection processes within rocket combustion chambers, or by supercritical combustion processes where the temperature is very high. As stated earlier, however, injection of supercritical CO
$_2$
and subsequent turbulent mixing between reactants is directly relevant to next-generation power systems, and understanding the processes at play is key to their design. With this application area in mind, and also to ensure novelty and interest in this research from a scientific viewpoint, attention is specifically focused on near-critical conditions. Doing so exposes the nonlinearities shown in figures 1 and 2. Given the shear layer topology of the set-up, selecting inlet streams at thermodynamic states that lie on either side of the Widom line ensures that as the streams mix, the Widom line is crossed. In this way, the nonlinear thermophysics is spatially confined and concentrated within the shear layer region.
Practical devices of course require multiple chemical species to mix and react; however, under such conditions, the mixture ceases to contain a single critical point and is influenced by additional processes such as species mass diffusion. This adds complexity to the detailed analysis of such flows. Thus to keep the problem focused while ensuring novelty and practical grounding in the study, a single-component flow of CO
$_2$
is considered.
The critical point of carbon dioxide is 7.38 MPa and 304.14 K. Thus the two fluid streams are injected at nominal pressure 8 MPa and nominal temperatures 308 K and 318 K. Figures 1 and 2 indicate that at this pressure, which is indicated by the red line in each graph, the selected temperatures bracket a region where nonlinear variations in properties are found. Shear is induced in the flow field by setting an upper/lower velocity ratio 2.3 : 1 between the two streams. Bulk velocities are set at 0.39
$\rm m \, s^{-1}$
for the upper stream, and 0.17
$\rm m \, s^{-1}$
for the lower stream. For a summary of the channel inlet boundary conditions applied in this work, the reader is referred to table 1. Note that the Reynolds number values in the table are reported using the channel half-height
$h$
as the length scale.
Channel inlet boundary conditions for the computations. The Reynolds number (
$Re$
) is based on the channel half-height
$h$
. Note that the upper and lower stream flow and boundary conditions for the experiment are inverted relative to the computations due to configurational considerations.
2.2. Experimental approach
2.2.1. Test bench
(a) Three-dimensional model of the experimental facility; colour is used to represent the inflow and outflow of CO
$_2$
. (b) Piping and instrumentation diagram.
Serially connected hardware components are used to ensure that the turbulent flow within the test section is fully developed and at the correct nominal operating conditions. A diagram of the experimental facility is provided in figure 3. Note that the accumulators shown in figure 3(a) are connected to the inlet and outlet of the compressor to dampen any pressure fluctuations. This allows the pressure within the test section to vary by no more than 0.83 % from its nominal value. The flow exiting the compressor and high-pressure accumulator is split into two streams. The mass flow rate of each stream is controlled using needle valves, and the temperature is set by the 10 kW heaters shown in figure 3. Settings are verified using both thermocouples and pressure transducers just upstream of the optical test section.
(a) Three-dimensional model of the test section. Note the presence of two (green) optical access panels. (b) Image of the test bench showing the inlet and optical sections. The inset shows the test section field of view as seen through the first (leftmost) optical access panel.
The test section itself is composed of two parts: an inlet section and an optical section, labelled accordingly in figure 4(a). The inlet section contains two stainless steel rectangular channels, each with a
$25.4\times12.7\,\text{mm}^2$
cross-section. The two channels are separated by a 1 mm thick plate. All solid walls are hydraulically smooth. Note that the pipe housing the inlet section is included for safety, given the extreme operating pressure used in this study. This pipe also houses the thermocouples and pressure transducers required for temperature and pressure control.
The optical section is 300 mm long, with outer diameter 127 mm and inner diameter 51 mm. A custom liner is used to change the cross-section shape from circular to square with side length 51 mm. Rectangular windows are fitted to permit optical access for the laser diagnostics. Four windows are included, with two in an upper–lower configuration, and two in a fore–aft configuration. The laser sheet used in the diagnostics passes through the upper and lower windows. Each of these windows has dimensions
$101.6\times6.35\,\text{mm}^2$
. The front and back windows are wider, having dimensions
$101.6\times25.4\,\text{mm}^2$
. The splitter plate separating the two fluid streams ends 51 mm from the entry of the viewable area. This leaves a
$51\times25.4\,\text{mm}^2$
viewable field where the shear layer dynamics may be observed.
2.2.2. Optical diagnostics
(a) Isometric and (b) front views of the spontaneous Raman scattering set-up in the current test section configuration.
Near the critical point, many thermodynamic and transport properties exhibit sharp variations. In particular, density represents a suitable property for optical measurement since field information may be extracted in a non-intrusive manner. Commonly used techniques include schlieren, interferometry and planar laser-induced fluorescence. However, these techniques are unable to provide accurate mixing quantification at the high operating pressure and associated critical point proximity, which is a key characteristic of this study. Instead, spontaneous Raman scattering is used to take advantage of the rotational and vibrational degrees of freedom of the carbon dioxide molecules. In addition, sCO
$_2$
exhibits higher than normal near-critical density relative to other commonly used gases. The relatively high scattering cross-section of CO
$_2$
results in a strong scattered signal, further underscoring the benefits of using Raman scattering for density measurements.
The optical diagnostic set-up for this study is visualised in figure 5. The technique measures the intensity of inelastically scattered light (see the inset of figure 5
a), from which density information can be extracted. The intensity of the Raman scattered signal (
$I_S$
) is related to various parameters according to
\begin{equation} I_{S} \propto \nu ^4I_{I}N\left ({\frac {\partial \alpha }{\partial Q}}\right )^2 \mbox {.} \end{equation}
In (2.1),
$\nu$
is the frequency of the exciting laser,
$I_I$
is the incident laser intensity,
$N$
is the number of scattering molecules,
$\alpha$
is the polarisability of the molecules, and
$Q$
is the vibrational amplitude (Larkin Reference Larkin2018). The angle of scattered light from the incident is
$90^\circ$
. The intensity of the Raman signal is proportional to the number density of the fluid under the condition that the flow is composed of a single chemical species. Equation (2.1) can be simplified down to
\begin{equation} I_{S} \propto C \sigma \rho I_{I} \mbox {,} \end{equation}
where
$\sigma$
is the scattering cross-section,
$\rho$
is the species density, and
$C$
is an experimental constant dependent on detector efficiency and geometry. Here,
$I_S$
and
$I_I$
retain their definitions. The Raman signal is calibrated using known density data at multiple thermodynamic states. This is done by collecting numerous flat-field signals (i.e. with both streams at identical thermodynamic and hydrodynamic states) in both isothermal and isochoric conditions, to assess whether the Raman signal intensity is proportional to the number density. For the isothermal tests, two tests were run, one at 35
$^\circ$
C, and one at 45
$^\circ$
C, with the pressure increasing from 7.6 to 8.4 MPa. The isochoric tests were done by varying both the temperature and pressure (to maintain a constant density and hence signal intensity). Two tests were run, one at target density 241
$\rm kg\, m^{-3}$
, and another at 417
$\rm kg\, m^{-3}$
. Further details and discussion regarding the calibration techniques and results are documented by Lim (Reference Lim2022, chapter 2).
The Raman signals are collected by directing a laser sheet through the upper and lower windows of the test section. An image of the scattered signal is then taken through the front window. The laser sheet and camera are visible in figure 5 and show how a signal is gathered. The laser sheet itself is formed using a Galilean beam expander, which uses one cylindrical diverging lens and two cylindrical converging lenses. The laser used for molecular excitation in this experiment has wavelength 532 nm and pulse speed 10 Hz. The energy within the incident light is 450 mJ. Scattered Raman signals are observed by blocking any extraneous incident light using a 532 nm notch filter and passing the Raman signal itself through a 570 nm band-pass filter. To amplify the scattered signal, an intensified charge-coupled device camera is used. The camera used (see figure 5) has full resolution of 5.5 megapixels and a field of view spanning
$2560\times2160$
pixels. The camera is operated with a
$4\times4$
pixel binning mode, to obtain a high signal-to-noise ratio. This leads to an effective camera resolution of
$640\times540$
pixels.
2.3. Computational approach
All of the computational work performed in this study is done using the RAPTOR code framework, which has been employed previously in similar settings. The theoretical framework includes the fully-coupled compressible conservation equations of mass, momentum, total energy and species for a general multicomponent reacting flow system. The solver is designed to handle high-Reynolds-number, high-pressure, supercritical fluid phenomena. Details related to the code have been well documented in the past (see e.g. Oefelein Reference Oefelein2006a ,Reference Oefelein b ).
The governing equations used to model the physical system are presented below in non-dimensional form using Cartesian tensor notation. The definitions of all dimensionless variables appearing in the equations are presented in table 2. In the equations,
$M$
represents a reference Mach number, and
$Re$
is a reference Reynolds number. The quantity
$Pr$
is the Prandtl number. Superscripts are used to distinguish quantities specific to a given chemical species
$\alpha$
, or for other notational purposes only. They are not to be confused with tensor components, which are denoted using subscripts. The Einstein summation convention holds for repeated tensor suffixes.
Definition of dimensionless variables used in the governing equations. Reference quantities (subscripted ‘ref’) are defined based on physical parameters of the system being modelled computationally. Dimensional quantites in the table are identified by a superscripted asterisk.
-
(i) Conservation of mass:
(2.3)
\begin{equation} \frac {\partial \rho }{\partial t} + \frac {\partial }{\partial x_i} \left ( \rho u_i \right ) = 0 \mbox {.} \end{equation}
-
(ii) Conservation of momentum:
(2.4)where
\begin{equation} \frac {\partial }{\partial t} \left ( \rho u_i \right ) + \frac {\partial }{\partial x_j} \left ( \rho u_i u_j + \frac {p}{M^2}\, \delta _{ij} \right ) = \frac {\partial }{\partial x_j} \sigma _{ij} \mbox {,} \end{equation}
$\delta _{ij}$
is the Kronecker delta tensor, and the viscous stress tensor
$\sigma _{ij}$
is given as(2.5)
\begin{equation} \sigma _{ij} = \frac {\mu }{Re} \left [ -\frac {2}{3}\, \frac {\partial u_k}{\partial x_k}\, \delta _{ij} + \left ( \frac {\partial u_i}{\partial x_j} + \frac {\partial u_j}{\partial x_i} \right ) \right ] \mbox {.} \end{equation}
-
(iii) Conservation of total energy:
(2.6)where
\begin{equation} \frac {\partial }{\partial t} \left ( \rho e^t \right ) + \frac {\partial }{\partial x_j} \left [\left ( \rho e^t + p \right ) u_j \right ] = \frac {\partial }{\partial x_j} \left [ q^e_j + M^2\left ( \sigma _{ij} u_i \right ) \right ] \mbox {,} \end{equation}
(2.7a)
\begin{align} e^t &= e + \frac {M^2}{2}\, u_i u_i \mbox {,} \end{align}
(2.7b)
\begin{align} e &= h^\alpha - \frac {p}{\rho } \mbox {,}\end{align}
(2.7c)
\begin{align} h^\alpha &= \left ( h_{f}^\circ \right )^\alpha + \int _{p^\circ }^{p} \int _{T^\circ }^{T} c_{p}^\alpha (T,p)\, \mathrm {d} T\, \mathrm {d} p \mbox {,}\end{align}
represent the total energy, internal energy and enthalpy of species
$\alpha$
, respectively. The quantity
$ ( h_f^\circ )^\alpha$
represents the enthalpy of formation of species
$\alpha$
at a given reference state. The energy flux from heat conduction is given as(2.8)
\begin{equation} q^e_j = \frac {\mu }{Pr}\,\frac {c_p}{Re}\,\frac {\partial T}{\partial x_j} \mbox {.} \end{equation}
Equations (2.3)–(2.6) require an equation of state for closure. Previous computational treatments of supercritical fluids tend to either employ a cubic equation of state such as the Peng–Robinson (PR) state equation (see e.g. Bellan Reference Bellan2006; Ovais et al. Reference Ovais, Kemenov and Miller2021; Unnikrishnan, Oefelein & Yang Reference Unnikrishnan, Oefelein and Yang2022) or the more sophisticated Benedict–Webb–Rubin (BWR) state equation (e.g. Dahms & Oefelein Reference Dahms and Oefelein2015). The BWR state equation provides consistently accurate results over a wider range of thermodynamic states than cubic state equations. However, it is computationally intensive. Cubic equations of state, on the other hand, can be less accurate, especially for mixtures at saturated conditions, but are computationally efficient. The most common of the cubic equations of state are the van der Waals, Redlich–Kwong (RK), Soave–Redlich–Kwong (SRK) and PR equations. These ‘two-constant’ equations provide theoretical and/or semi-empirical corrections to the ideal gas equation of state that account for intermolecular forces of attraction and volumetric effects. The van der Waals equation is the oldest of the two-constant cubic equations and was first presented in 1873. Research conducted by van der Waals provided the theoretical form of the corrective constants and also formed the basis of the corresponding states principle (Leland & Chappelear Reference Leland and Chappelear1968; Rowlinson & Watson Reference Rowlinson and Watson1969). The RK equation (presented by Redlich and Kwong in 1949) adopted a similar procedure to that followed by van der Waals to derive semi-empirical corrections to the orginal constants. The SRK and PR equations (presented by Soave in 1972, and by Peng and Robinson in 1976, respectively) provided further refinments by converting the corrected constants to empirical coefficients with a functional dependence on the reduced temperature and Pitzer’s acentric factor.
Experience has shown that that both the SRK and PR equations, when used in conjunction with the corresponding states principle, can give accurate results over the range of pressures, temperatures and mixture states of interest in this study. In particular, the operating pressure and stream temperatures create conditions that are sufficiently removed from the saturation curve. The SRK coefficients are fitted to vapour pressure data and are thus more suitable for conditions when the reduced temperature is less than 1. The PR coefficients, on the other hand, are more suitable for conditions when the reduced temperature is greater than 1. In the current study, the cubic PR equation of state is used (Peng & Robinson Reference Peng and Robinson1976). The historic success of this equation in similar studies, when coupled with the extended corresponding states methodology, provides some confidence regarding its use for the purposes of this research.
Having established an analytical representation for real gas PVT behaviour through the PR equation of state, dependent thermodynamic properties are obtained in two steps. First, properties are combined at a fixed temperature using the extended corresponding states methodology to obtain the state at a given reference pressure. A pressure correction is then applied using departure functions of the form given by Reid et al.(Reference Reid, Prausnitz and Poling1987, chapter 5). These functions are exact relations derived using the Maxwell relations (see e.g. Van Wylen & Sonntag Reference Van Wylen and Sonntag1986, chapter 10) and make analytic use of the real gas PVT path dependencies dictated by the equation of state. Transport properties are evaluated using the extended corresponding states model of Ely & Hanley (Reference Ely and Hanley1981a ,Reference Ely and Hanley b ) that provide departure functions for deviations from given reference states.
Only a brief discussion regarding the evaluation of thermophysical properties has been provided here in an effort to balance relevant detail with focus on the specifics of this research. For full details regarding the property evaluation models used in the numerical solver, see Oefelein (Reference Oefelein1997), which provides a complete discussion in this regard.
Sketch of the computational domain. The red boxed region demarcates the approximate location of the shear layer. Axis labels are provided in non-dimensional terms, with the thickness of the splitter plate used as the reference length.
2.3.1. Computational domain
The DNS calculation takes place on a domain sized as shown in figure 6, which is intended to mimic the optical section of the experimental facility, and has the same channel heights and splitter plate thickness. The two fluid streams are fed through the 12 mm tall inlet planes shown on the far left of the image in figure 6. The fluid injected through the inlet planes advects through channels with rectangular cross-sections for the 5 mm length of the splitter plate. The streams are brought together, and mixing commences downstream of the 1 mm thick splitter plate tip. The flow is allowed to evolve over the remaining 35 mm of the computational domain. In the interest of maintaining computational tractability while maintaining statistically decoupled spacing between boundaries, the cross-stream width of the domain is sized as 6 mm.
Plots of the (a,b) longitudinal, (c,d) cross-stream and (e,f) spanwise two-point, one-time spatial autocorrelation functions, for (a,c,e) the upper channel and (b,d,f) the lower channel. The results indicate that at separation distance 6 mm, the turbulence becomes decorrelated, and this is used as justification for the domain width in the computations.
Periodic boundary conditions are applied along the fore and aft walls of the computational domain to accommodate the difference in width between the experimental and computational set-ups. The resulting implicit assumption is that the channel core flow in the experiment is not affected by the presence of the fore and aft walls. The width of the computational domain is verified by performing initial calculations on a wider domain to obtain the two-point, one-time spatial autocorrelation function, and sizing the domain down based on the results obtained. The results shown in figure 7 indicate that the correlations decay to at least a weak level within approximately 2 mm for the upper channel turbulence, and 5 mm for the lower channel. Note that this measure is susceptible to convergence issues due to large-scale flow features, thus is not a typical method for domain sizing. However, with the computational intensity of DNS at these conditions prohibiting iterative domain design, the results in figure 7 at least provide some indication regarding domain size sufficiency.
2.3.2. Boundary conditions and grid characteristics
A summary of the inflow boundary conditions is presented in table 1. Turbulent inflow signals are generated using the ‘ensemble synthetic eddy method’ of Schau et al. (Reference Schau, Johnson, Muller and Oefelein2022) due to its computational efficiency and ability to reconstruct a correlated time-dependent inflow signal with realistic turbulent structures. These structures play a large role in governing the mixing layer dynamics once they are advected off the tip of the splitter plate. Hence capturing them accurately is essential to the reconstruction of the overall macroscopic mixing behaviour in the domain. The generation of the coherent inflow signal relies on prior knowledge of the mean velocity profile, Reynolds stress profiles, and integral scales (i.e. the energetic eddy length and time scale distributions). For these simulations, Reynolds stress profiles are obtained from the DNS data of Moser, Kim & Mansour (Reference Moser, Kim and Mansour1999), and the integral scales are obtained directly from Jarrin (Reference Jarrin2008). A turbulent signal is generated with a temporal periodicity of five flow-through times calculated using the upper channel bulk velocity. This ensures statistical decoupling of the velocity fluctuations. In all, 200 instantaneous velocity frames are synthesised to recreate a realistic inflow signal, with numerical interpolation between frames accommodating the somewhat smaller simulation time steps.
No-slip adiabatic boundary conditions are applied to all of the solid walls in the domain, namely, the upper and lower confining walls and the splitter plate. A subsonic outlet boundary condition is prescribed at the exit of the domain, and the fore and aft walls of the domain are assigned periodic boundary conditions.
Channel-specific resolution data for the DNS grid. The subscript
$w$
denotes data immediately adjacent to any solid walls, while
$c$
denotes channel core data. Here,
$N_{y^+ \leqslant 30}$
denotes the number of computational cells within the buffer region of the turbulent boundary layer adjacent to the walls. The
$+$
superscript is used to denote data reported in wall-normal units.
Global grid characteristics. Here, AR denotes aspect ratio, SF the skewness factor, and GR the cell-to-cell growth ratio. Data reported include cells within the buffer region.
(a) Kolmogorov length scale, and (b) axial, (c) cross-stream, (d) spanwise grid spacing normalised by the Kolmogorov length scale, as functions of height above the domain lower wall. Calculations consider only the solenoidal component of the total dissipation rate. They involve ensemble averaging in the spanwise direction using one instantaneous data frame with the mixing layer in a statistically stationary state. Data are sampled at an arbitrarily selected axial location 5 mm downstream of the splitter plate tip.
Plots showing (a) the smallest turbulent thermal length scale (
$\eta _\theta$
) and (b) axial, (c) cross-stream, (d) spanwise normalised grid spacing data as functions of height above the domain lower wall. Scales smaller than
$\eta _\theta$
are dominated by diffusive mechanisms. Calculations still assume classical Kolmogorov turbulence and involve ensemble averaging in the spanwise direction for one instantaneous data frame with the mixing layer in a statistically stationary state. Data are sampled at an axial location 5 mm downstream of the splitter plate tip.
The canonical nature of the computational domain allows the grid generation process to occur with primary focus on achieving adequate resolution and topological quality. Grid resolution especially is of paramount importance. With this in mind, the grid is initially fabricated assuming a one-to-one correspondence between the classical Kolmogorov length scale and the viscous length scale. Confidence in this method is based on the success of research from Kim, Moin & Moser (Reference Kim, Moin and Moser1987) and Lee & Moser (Reference Lee and Moser2015). A parametric characterisation of the DNS grid created in this manner is presented in tables 3 and 4. The grid is also designed with a buffer region near the exit plane to dampen vortices as they advect out of the domain. This is to alleviate issues related to solution destabilisation resulting from flow re-entry through the exit plane.
Once the flow field reaches statistical stationarity, it is sampled in order to verify the grid resolution. The results of one such field interrogation are presented in figures 8 and 9. Figure 8 shows the Kolmogorov length scale and grid spacing normalised by the Kolmogorov length scale in all three coordinate directions at axial location
$x = 10\, \rm mm$
. The data presented are obtained using the relations
\begin{align} \eta _k &= \left \langle \left (\frac {\nu ^3}{\varepsilon }\right )^{{1}/{4}} \right \rangle,\quad \mbox {where} \end{align}
\begin{align} \varepsilon &= 2\nu S_{ij} S_{ij}\quad \mbox {and} \end{align}
\begin{align} S_{ij} &= \frac {1}{2}\left ( \frac {\partial u_i}{\partial x_j} + \frac {\partial u_j}{\partial x_i}\right ) \mbox {.} \end{align}
Here, angular brackets denote ensemble averaging in the spanwise direction, and the results shown were obtained for one instant in time. Although classical Kolmogorov turbulence has been assumed, these results provide some indication regarding the adequacy of the grid resolution. In particular, the grid spacing in each of the three coordinate directions is adequate to resolve at the very least non-dilatational fluid motions within the multiscale field.
A similar calculation is performed to ascertain the resolution adequacy with regards to turbulent thermal fluctuations. The results of this analysis are shown in figure 9. The plots were produced using the methods found in the work of Sciacovelli & Bellan (Reference Sciacovelli and Bellan2019), which itself borrows results from the analysis of Goto & Kida (Reference Goto and Kida1999). Employing these findings, the smallest relevant thermal scale, denoted
$\eta _\theta$
, in the turbulent field can be estimated as
\begin{equation} \eta _\theta = \left \langle \eta _k\, Pr^{-{1}/{2}} \right \rangle \mbox {.} \end{equation}
In this equation, the quantity
$Pr$
is the local molecular Prandtl number, and angular brackets once again denote ensemble averaging in the spanwise direction. The expression implies that at Prandtl numbers greater than unity, the smallest active thermal scale ought to be smaller than the smallest dynamic scale. With the boundary conditions selected for this study, the Prandtl numbers are well over unity, hence stressing the importance of resolution in this sense. The results indeed show that the smallest thermal scales are smaller than the smallest dynamic scales. The non-dimensional grid spacing found in the field, e.g. shown here for axial location
$x=10$
mm (see figures 9
b–d), indicates that the resolution still largely lies within one order of magnitude of this scale, except very close to the upper wall of the domain, which is sufficiently far away from the shear layer region.
3. Results and discussion
In this section, both experimental and numerical results are presented. Specifically, the data are used in a complementary fashion. The experimental data are used to obtain detailed profiles and mixing behaviour far downstream of the splitter plate, while the numerical results provide information closer to the plate tip, where obtaining experimental measurements is difficult due to the short time and length scales involved.
3.1. Instantaneous and ensemble averaged field visualisations
Visualisation of the shear layer, obtained using shadowgraphy. The coloured boxes represent two fields of view (FoV) for which spontaneous Raman scattering signals are obtained. In the images, the flow travels from left to right. Flow conditions are the same as outlined in table 1 but are flipped relative to the configuration specified in the table, and this is accounted for in all quantitative results.
Figure 10 presents a view of the shear layer obtained through the direct shadowgraphy method. The data are instantaneous, and the image is created by stitching together instantaneous snapshots from three different fields of view. Also shown in the image are two boxes (in blue and red) demarcating the two fields of view used to obtain the Raman scattering results presented later in this paper.
Notably, the shear layer as visualised in figure 10 does not show the development of the coherent structures that are hallmarks of such configurations (Brown & Roshko Reference Brown and Roshko1974). This arises from the relatively low density ratio between the two streams, approximately 1.6 under the current conditions, which complicates imaging coherent structures using backlight techniques. At distance 25 mm from the splitter plate, marking the start of FoV 1, the initiation of coherent structures becomes faintly visible. Much finer details and the structures can be seen in the instantaneous Raman images presented in figure 11. They appear as angular yellow structures penetrating into the purple field above, and are perhaps most evident in the leftmost image in figure 11.
Instantaneous visualisation of the plane shear layer. Results are obtained via Raman scattering, and presented in terms of the density difference ratio, defined according to (3.1). Each image corresponds to a different instant in time, with no correlation between images. The arrows indicate the direction of the flow, as well as an approximate indication of the velocity magnitude in each stream. The development of the classical coherent structures, first reported by Brown & Roshko (Reference Brown and Roshko1974), can be seen here, especially in the leftmost image.
To obtain the images shown in figure 11, the spontaneous Raman scattering signal is post-processed. Full details can be found in Lim (Reference Lim2022, chapter 2); however, a brief description of the process is provided here. Spontaneous Raman scattering enables density estimation by correlating the signal intensity with the species density. This is done by applying (2.1) and (2.2). By calibrating the Raman signal using known density measurements under varying thermodynamic conditions, accurate density measurements can be deduced. To process the images, black and white field references are first collected to both evaluate the health of the CCD camera chips and also obtain a baseline signal-to-noise ratio for the current set-up. Next, flat-field images are obtained by imposing uniform conditions (pressure, temperature and mass flow rates) between the two streams. This should result in identical signals corresponding to the same mean density. However, due to imperfections in the circular laser beam (which are exacerbated through the conversion to a laser sheet), this is not the case. An adjustment is made using a reference sheet signal. A two-dimensional Gaussian filter is then applied to correct the resultant intensity profile to maintain the mean intensity of the flat-field image. For this specific image processing pipeline, a range of standard deviation values was tested for Gaussian smoothing, and standard deviation 3 was found to be the most suitable. This value effectively reduced noise in the signal while preserving the overall intensity profile of the one-dimensional signal. After this, the resulting image and the original flat-field signal are used together to provide a map of intensity corrections. Instantaneous Raman signals are corrected using this intensity correction map to provide corrected images. Finally, the effects of beam steering, which results from the laser passing through media with different densities, must be corrected. This is done using the frequency-space methods outlined in Mohaghar et al. (Reference Mohaghar, Carter, Pathikonda and Ranjan2019) and Lim (Reference Lim2022, chapter 4). The resulting images, which represent the density of the fluid, are presented in figure 11. The raw Raman signal intensity is converted to a non-dimensional density difference ratio, denoted
$\xi$
, which is calculated as
\begin{equation} \xi = \frac {\rho - \rho _1}{\rho _2 - \rho _1} \mbox {,} \end{equation}
where
$\rho$
is a locally sampled density gathered from the raw Raman signal, and
$\rho _1$
and
$\rho _2$
are the nominal densities of the upper and lower streams (in the experimental configuration), respectively. The density
$\rho$
is obtained from the density ratios of the isochoric Raman signal calibration exercise as
\begin{equation} \rho = \frac {I_s(x, y)}{I_{{ref}}(x, y)}\, \rho _{{ref}} \mbox {,} \end{equation}
where
$I_s$
is the (raw) sampled Raman signal intensity,
$I_{{ref}}$
is a flat-field Raman signal intensity collected at 8 MPa and 45
$^\circ$
C, and
$\rho _{{ref}}$
is the density value implied from a valid equation of state at 8 MPa and 45
$^\circ$
C.
The presence of the coherent structures within the plane shear layer is also confirmed by the computational results, presented in figure 12. Here, a translucent bounding box identifies the domain extents and is coloured to highlight the temperature variation within the domain. To identify vortical structures within the field, the Q-criterion (Jeong & Hussain Reference Jeong and Hussain1995) is used. Isosurfaces of the Q-criterion are coloured by the local isothermal compressibility of the fluid. The isothermal compressibility is one of two measures of thermodynamically induced compressibility, the other being the thermal expansion coefficient. Both show order-of-magnitude variations at the selected conditions (see e.g. figure 1 c), thus both show order-of-magnitude changes within the shear layer. The compressibility induced by variations in these properties permits an exchange between kinetic and internal energies, thus affecting the overall growth rate of the shear layer.
Three-dimensional view of the spatially evolving mixing layer obtained via DNS. The Q-criterion isosurfaces are visualised at
$Q=110\,000$
$\rm s^{-2}$
. The isosurfaces are coloured by the local isothermal compressibility of the fluid. The translucent domain bounding box is coloured by temperature to highlight its spatial variation.
In addition to instantaneous characteristics, mean field data are also required to provide information regarding shear layer growth rates and profiles through the field. To obtain these mean field visualisations, 250 instantaneous images are collected and ensemble averaged, with the fields shown in figure 13. Figure 13(a) represents the field of view outlined by the blue box in figure 10 (i.e. FoV 1), while figure 13(b) corresponds to the red box (i.e. FoV 2). The images in figure 13 show contours of the density difference ratio. The images represent mixing with a velocity ratio larger than unity, and the mixing layer shows a preference for growing upwards, due to the higher momentum of the lower stream, in line with Dimotakis (Reference Dimotakis1986). Common features of shear layers in general can be observed in the images in figure 13, such as the relaxation of gradients as both convective (i.e. entrainment) and diffusive processes act. The shear layer is seen to grow appreciably within the two fields of view.
Ensemble-averaged flow fields of the sCO
$_2$
mixing layer. Density difference ratios (see (3.1) for a definition) are plotted as functions of space for a two-dimensional slice through the mixing layer. Results are obtained via spontaneous Raman scattering.
Ensemble-averaged density field, obtained through spontaneous Raman scattering. Five temperature measurements at distinct locations are obtained using RTDs. These are represented by scatter points and are overlaid on the contour. Results are shown for FoV 1.
Figure 14 provides a similar visualisation, this time of the dimensional density field within FoV 1. Circular overlays are provided on the graph to show the variation in temperature, measured using resistance temperature detectors (RTDs). The graph in figure 14 in particular highlights the strong gradient in density within the shear layer that occurs with only a very mild change in temperature.
It has been established that supercritical shear layers at these conditions exhibit regions of HDGM (Bellan Reference Bellan2000; Miller et al. Reference Miller, Harstad and Bellan2002). These regions act as local material interfaces, thus inhibit entrainment on the large scales. This occurs instantaneously within the shear layer, but will also impact the mean growth rate of the layer itself. Hence with the HDGM regions opposing the equalising tendencies of convective entrainment and momentum diffusion, supercritical shear layers may grow more slowly than an equivalent configuration at atmospheric pressure. Though additional experiments and/or simulations are required to confirm this, quantitative results regarding shear layer growth are available from this effort. This and additional quantitative data are presented next.
3.2. Density profiles
Figure 15 includes two plots with data obtained via both Raman spectroscopy and the RTD devices in the test section. The data are taken from FoV2 at two streamwise locations. The data shown in figure 15(b) are taken approximately 15 mm downstream of those shown in figure 15(a). In each graph, the solid black line corresponds to the data obtained via Raman spectroscopy, and the green shaded region corresponds to one standard deviation from the measured value. The RTD temperature measurements are mapped to values of
$\xi$
to enable plotting both results on the same axis. The RTD measurements are shown as red circles, and error in their measurements is indicated by a horizontal red bar in each case.
Density difference ratio
$\xi$
as a function of height above the lower wall of the domain. The black solid line is obtained via direct density measurements from Raman scattering. The red scatter points represent calculated density values, here represented in terms of
$\xi$
. The shaded green region brackets one standard deviation of the Raman-obtained density profile.
The density profiles indicate a region of large density gradient (at least in the cross-stream direction) at height approximately
$y = 5$
mm. This sharp transition occurs due to the mixing between newly entrained fluid from the lighter, gas-like CO
$_2$
stream, and the fluid within the shear layer, which has already undergone a level of molecular mixing via diffusive mechanisms. Note that the data also show a second region where the density gradient is strong. This is at a height between
$-5\,\text{mm}$
and
$-10\,\text{mm}$
, and is perhaps most clearly seen in figure 15(a). The second gradient in
$\xi$
is due to the interaction of the heavier liquid-like stream and the fluid in the shear layer. The results in figure 15 show that this gradient relaxes more readily as the flow evolves downstream than the region between the shear layer and the gas-like CO
$_2$
.
The sharp gradients in density observed at approximately
$y = 5$
mm in both plots in figure 15 occur between the lighter gas-like fluid and the mixed fluid. This trend runs counter to the trends observed by Ovais et al. (Reference Ovais, Kemenov and Miller2021), whereby regions of HDGM form between regions of heavy fluid and mixed fluid. Further opposition is confirmed by contrasting both plots in figure 15. The relaxation in the gradient in
$\xi$
between
$-5\,\text{mm}$
and
$-10\,\text{mm}$
as the flow advects downstream implies a greater level of entrainment and subsequent mixing between the liquid-like CO
$_2$
and the mixed fluid than between the gas-like CO
$_2$
and mixed fluid. The HDGM results presented by Ovais et al. (Reference Ovais, Kemenov and Miller2021) (and also in Bellan Reference Bellan2000) are instantaneous, and arise from species diffusion, not thermal inhomogeneity. These authors theorise that the development of the HDGM regions in their multicomponent simulations could be a result of low rates of molecular diffusion (also theorised by Okong’o et al. Reference Okong’o, Harstad and Bellan2002). The configuration in this research does not contain multiple chemical species, so as a substitute, we consider the thermal diffusivity, which for the high-density (liquid-like) stream takes value
$3.65\times 10^{-9}\,\text{m}^2\,\text{s}^{-1}$
, while for the low-density (gas-like) stream, it takes on value
$3.15\times 10^{-8}\, \rm m^2\, s^-{^1}$
. Thus we ought to expect the gradient between the upper stream fluid and the mixed fluid to relax more readily than the lower stream fluid; however, surprisingly, this is not observed. Instead, it appears as though the region of HDGM forms on this side, inhibiting entrainment and hence mixing, which leads to the maintenance of the gradient observed at
$y = 5$
mm in figure 15.
3.3. Shear layer growth rate
Perhaps the most readily quantifiable and relevant characteristic of the plane shear layer is the growth rate. The level of mixing between the two fluid streams can be quantified using many methods, and a calculation of the so-called ‘mixed material’ is used here. The mixed material within the shear layer is computed experimentally via the ensemble-averaged Raman results as
\begin{equation} \delta = \frac {\iint K_m Y_{gl} Y_{ll} \, \mathrm {d}x\, \mathrm {d}y}{\int \mathrm {d}x}, \end{equation}
where
$Y_{gl}$
is the fraction of low-density (i.e. gas-like) CO
$_2$
, and
$Y_{ll}$
is the fraction of high-density (liquid-like) CO
$_2$
within a given volume. The quantity
$K_m$
is called a density difference ratio constant. This constant is selected such that the integrand in the numerator produces a value unity when the fluids are deemed fully mixed. This fully mixed criterion is determined from a mass flux balance using the initial thermodynamic states and conditions of the upper and lower streams. For the conditions selected in this study, the fluids are deemed to be fully mixed when
$Y_{ll}=0.75$
. This results in
$K_m$
having value 5.35. In this way, the integrand ranges from zero to one, with zero corresponding to unmixed material, and one corresponding to fully mixed material. The integrand thus physically represents how much mixing has resulted from the high-momentum stream penetrating into the low-momentum stream. When considering the area integral in the numerator of (3.3), it is important to note that only pixels where the integrand is larger than some threshold value (denoted
$\Theta$
) are to be considered. This is because the calculation naturally ought to exclude regions where the two streams are unmixed. Mathematically, this criterion can be represented as
\begin{equation} K_m Y_{gl} Y_{ll} \gt \Theta \mbox {.} \end{equation}
This threshold has taken on various values across the literature. For instance, Konrad (Reference Konrad1977) considers only a threshold
$\Theta = 0.80$
, while Mohaghar (Reference Mohaghar2019) and Mohaghar et al. (Reference Mohaghar2019) use
$\Theta = 0.86$
. This work uses a threshold
$\Theta = 0.70$
to match the visual thickness observations obtained via schlieren imaging (as defined by Brown & Roshko Reference Brown and Roshko1974). This resulted in a minimum mean square error between the visual thickness and the thickness obtained by considering the mixed material.
Mixed material as a function of streamwise location. The colours represent the fields of view highlighted in figure 10. In the graph legend, MM stands for mixed material (using
$\Theta = 0.70$
), and VT represents the visual thickness obtained by using a full-width half-maximum (FWHM) technique. The grey dashed line represents the growth rate found via the visual thickness, while the purple dotted line represents the growth rate calculated through the mixed material.
Figure 16 presents the results of the mixed material calculation obtained from the experimental results. The coloured scatter in the figure represents the two different fields of view highlighted in figure 10. The standard deviations are reported as well. Figure 16 shows a monotonically growing mixing layer, as expected. The lines of best fit (plotted for the visual thickness data as well as the mixed material data) have gradients approximately 0.20. Note that the image also includes growth rates calculated through the visual thickness metric, as well as the full-width half-maximum method, used in past research of such flows (see Oschwald & Schik Reference Oschwald and Schik1999; Chehroudi et al. Reference Chehroudi, Talley and Coy2002).
3.4. Turbulence kinetic energy spectrum
Additional insight into the fundamental nature of the turbulence dynamics in the field may be found using the computational results. However, prior to proceeding, it is important to ensure that the turbulent field contains adequate characteristics. A one-dimensional spectrum of the turbulent kinetic energy, sourced from the computational data, is presented in figure 17 with this in mind. The spectrum is computed along the spanwise (
$z$
) direction since the fore and aft walls have periodic boundary conditions applied to them, which permits application of the Fourier decomposition. The data were collected at an arbitrarily selected axial location 20 mm downstream of the splitter plate tip, and 12.5 mm above the domain lower wall (i.e. within the shear layer).
One-dimensional spectrum of the turbulent kinetic energy
$\text{E}(\kappa _3)$
taken from the DNS data sampled in the shear layer (
$x = 25$
mm,
$y = 12.5$
mm). Spectral decomposition is performed on a velocity signal taken along the homogeneous
$z$
direction. A reference
$-5/3$
line, which corresponds to classical inertial range Kolmogorov turbulence, is added as a point of comparison.
Schematic of the AIM. The envelope of realisability is identified by the red line, while the three-dimensional greyscale graphics represent the shape of the Reynolds stress ellipsoids along the right and left limbs of the AIM. An ellipsoid is also included for the isotropic turbulence limit, for reference.
The turbulence cascade is shown in the results presented in figure 17, with the large scales containing the majority of the turbulent kinetic energy and the energy per wavenumber dropping sharply at higher wavenumbers. This indicates that the highest resolved wavenumber indeed lies within the dissipation range. These trends, coupled with the results presented in figures 8 and 9, indicate that the simulation accurately captures all active dynamic and thermal scales. Note that the flattening of the spectrum at low wavenumbers is a natural result of aliasing since the spectrum as presented is one-dimensional.
3.5. Turbulence anisotropy analysis
The research of Miller et al. (Reference Miller, Harstad and Bellan2002), as well as Okong’o & Bellan (Reference Okong’o and Bellan2002) and Okong’o et al. (Reference Okong’o, Harstad and Bellan2002), identified and explained the significance of the HDGM regions that occur in supercritical shear layers. The regions are found to intensify any initial anisotropy in the density field, and are found in both multicomponent flows and thermally inhomogeneous flows. The latter scenario is directly relevant to the configuration selected here, so this subsection aims to examine the influence of the observed density stratification on the fluid dynamics.
The anisotropy invariance map (AIM) provides a way to analyse the data in this regard and uncover dominant trends. Though originally developed for incompressible flow analysis, it has been used successfully to study compressible flows in the past (see e.g. Maidi & Lesieur Reference Maidi and Lesieur2005; Kim, Elliott & Craig Dutton Reference Kim, Elliott and Craig Dutton2020). The plots locate the invariants of the Reynolds stress tensor at a given spatial location on the so-called invariant plane. A schematic is provided in figure 18 and highlights key points and regions of interest. The schematic also provides visualisations of the Reynolds stress ellipsoids that correspond to three regions and locations in the map.
Anisotropy invariant maps (Lumley triangles) obtained from the computational results: (a) randomly down-sampled DNS data; (b) points in mixing layer only. Invariants of the anisotropic part of the Reynolds stress tensor (as defined in (3.5)) are plotted as functions of ensemble-averaged density and coloured by height above the lower wall of the domain. Each scatter point represents the values of the respective invariants at a given spatial location in the domain. The realisability envelope is identified by black dashed lines.
Figure 19 shows a three-dimensional view of the AIM produced using the DNS data. The third dimension is used to highlight the variation in ensemble-averaged density, denoted
$\overline {\rho }$
. Each scatter point in the figure corresponds to data at a given cell within the computational domain. Since the invariant map necessarily strips the data of spatial context, each point is coloured by the cross-stream location where it was sampled. The maps are created by first reconstructing the Reynolds stress tensor
$\overline {\rho u^{\prime}_i u^{\prime}_j}$
at all points in the domain. The overline denotes ensemble averaging, while primes denote fluctuations from a Reynolds average. The anisotropic part of this tensor is then isolated. To provide a level of generality and intuitiveness to their interpretation, the values are non-dimensionalised using the turbulent kinetic energy at each corresponding location. Precisely, the formulation is
\begin{equation} b_{ij} = \frac {\overline {\rho u^{\prime}_i u^{\prime}_j} }{2 k} - \frac {1}{3}\,\delta _{ij} \mbox {,} \end{equation}
where
$b_{ij}$
is the anisotropic part of the Reynolds stress tensor,
$\delta _{ij}$
is the Kronecker delta tensor, and the turbulent kinetic energy
$k$
is defined as
\begin{equation} k \equiv \frac {1}{2}\, \overline {\rho u^{\prime}_k u^{\prime}_k} \mbox {.} \end{equation}
The invariants of the anisotropy tensor provide information that is, by definition, unaffected by coordinate axis transformation. As a result, they remain insulated from arbitrary selections in this regard, providing valuable insights into the field anisotropy. Two invariants may be defined as
\begin{align} \eta &\equiv \left ( \tfrac {1}{6} b_{ij} b_{ji} \right )^{{1}/{2}} \mbox {,} \end{align}
\begin{align} \epsilon &\equiv \left ( \tfrac {1}{6} b_{ij} b_{jk} b_{ki} \right )^{{1}/{3}} \mbox {,} \end{align}
where the Einstein summation convention is implied across repeated tensor suffixes. These definitions differ from the original ones of Lumley & Newman (Reference Lumley and Newman1977). They ensure that the map is bounded by two straight ‘legs’ instead of the curved boundary that results from the original formulation.
Two-dimensional view of the Lumley triangle plotted in figure 19. Data are isolated at
$\overline {\rho } = 324.6\, \rm kg\, m^{-3}$
, which is halfway between the upper and lower limits of the ensemble-averaged density. The realisability envelope is demarcated by black dashed lines. All data points found within the DNS field matching this density value are plotted on the map, and no down-sampling is performed.
Figures 19 and 20 show the Reynolds stress anisotropy for certain points within the DNS domain. A known drawback of the AIM is the difficulty associated with visualisation when many data points are included (Emory & Iaccarino Reference Emory and Iaccarino2014). To alleviate this issue, the data in figure 19(a) include DNS data after random down-sampling has been performed to maintain dominant trends, yet sparsen the data. Figure 19(b) includes data for points within the shear layer, with the points identified using a 95 % threshold based on the density. Random down-sampling is also performed here to sparsen the data. Figure 20 presents a two-dimensional visualisation of the AIM, including all points in the domain for which
$\overline {\rho } = 324.6$
kg m
$^{-3}$
. These data provide novel visualisations for at least supercritical flows at these conditions, if not supercritical flows in general.
A few dominant trends and subsequent implications emerge from the data. First, the data appear to show similarities with the AIMs obtained from ideal gas shear layer experiments (Pope Reference Pope2000, chapter 11). Points are clustered near the right and left limbs of the triangle, and this is perhaps seen most clearly in figure 20. Points near either limb correspond to either prolate or oblate Reynolds stress ellipsoids, with points near the right corresponding to data taken near the shear layer centre, and points near the left corresponding to data taken at the edge of the shear layer. No data are found at the origin of the invariant plane in figures 19 and 20, which serves to confirm the claim made in Bellan (Reference Bellan2000) that a study of supercritical isotropic turbulence is of limited relevance. Of course, the shear layer configuration itself introduces a level of anisotropy in the field. However, the presence of the HDGM regions within the field will serve to amplify the existing field anisotropy and push points farther away from the origin of the invariant plane. The combined action of the HDGM regions and the inherent anisotropy of the configuration explains this observation. In fact, it is likely that the data for such supercritical flows would show a tendency to move towards the anisotropic one- or two-component limits in the AIM. The data in figure 20 show an interesting tendency for some points within the shear layer to statistically approach the one component limit, which is a surprising result that may be due to attenuation of turbulent fluctuations due to the density stratification. To confirm or reject this hypothesis, DNS calculations or experimental measurements of an atmospheric pressure shear layer with identical thermal inhomogeneity and momentum ratio are required.
The anisotropy found in the turbulence has important implications for the mixing in the shear layer. Fundamentally, mixing between the two streams is a result of both advective and diffusive processes. Advective phenomena (e.g. coherent structure dynamics and growth) drive turbulent mixing, which acts to increase the effective surface area over which flow gradients develop and resultant diffusion processes engage. At these conditions, however, HDGM regions form (Bellan Reference Bellan2006) and inhibit turbulent fluctuations in directions normal to their local orientation. This attenuation of turbulent fluctuations in turn slows turbulent mixing relative to conditions where the HDGM regions are not present. Since turbulent mixing helps to increase the effective surface area available for molecular diffusion, this represents one way in which molecular-level scalar mixing is inhibited.
In addition to these fluid mechanical advective considerations, molecular mixing will also be inhibited by thermodynamic nonlinearities at the selected conditions and corresponding property variations. This creates a second way in which molecular mixing rates can be diminished relative to ideal gas conditions. At near-critical conditions, the dynamic viscosity
$\mu$
rises with increases in pressure at a given temperature. This is clearly seen in figure 2(a). The density also shows the same behaviour (see figure 1). However, the rise in density dominates the rise in dynamic viscosity, with the result that the kinematic viscosity (a measure of momentum diffusivity) decreases with pressure. This is shown in figure 21. From figure 2(b), it can be seen that the Prandtl number rises with increasing pressure at a given temperature. This quantity is defined as
$Pr \equiv \nu / \alpha$
, where in this subsection only,
$\alpha$
is used to denote the thermal diffusivity. Hence if the kinematic viscosity
$\nu$
decreases with pressure, but
$Pr$
rises with pressure, then
$\alpha$
must also decrease with pressure and exhibit increased sensitivity to these changes relative to
$\nu$
. Overall, this leads to the conclusion that the thermal diffusivity decreases with increasing pressure, hence inhibiting molecular mixing between thermally inhomogeneous streams at these conditions relative to those at lower pressures. This argument, which is independent of advective considerations, provides a way in which thermodynamics at near-critical conditions affects molecular mixing. Since the governing equation for the diffusion of thermal energy is identical to that for the diffusion of chemical species, the above arguments may be extended analogously to the ratio
$D/Sc$
(where
$D$
denotes mass diffusivity, and
$Sc$
denotes the Schmidt number) for a given chemical species. Note, however, that this depends strongly upon the particular chemical species under consideration, so the analogy may not always hold. Nevertheless, considerations such as that outlined above carry important implications for multi-species and thermally inhomogeneous shear mixing in the context of combustion physics.
Variation of the kinematic viscosity of CO
$_2$
with temperature and pressure.
3.6. Enstrophy dynamics
To further study the dynamics of turbulence within the shear layer, the mean square vorticity (i.e. enstrophy, defined as
$\overline {\omega _i \omega _i}$
) is analysed using the computational data set. This quantity is related to the rotational energy in the flow as well as the viscous dissipation function (Lele Reference Lele1994), highlighting its importance to the dynamics in the field.
The budget of
$\overline {\omega _i \omega _i}$
provides insight into the precise sources of the gains to and losses from the rotational energy in the flow field, since it is a square quantity. The coherent structures within the shear layer possess significant quantities of rotational energy, and ultimately, the mixing and entrainment of fluids from each stream can be linked to the rotational energy in the shear layer. The (Reynolds-averaged) enstrophy transport equation can be written as
\begin{equation} \frac {1}{2}\,\overline {\frac {\mathrm {D}\left (\omega _i \omega _i\right )}{\mathrm {D}t}} = \underbrace {\overline {\omega _i\omega _j\,\frac {\partial u_i}{\partial x_j}}}_{\text{I}} - \underbrace {\overline {\frac {\partial u_j}{\partial x_j}\,\omega _i\omega _i}}_{\text {II}} + \underbrace {\overline {\frac {1}{\rho ^2}\,\omega _i\epsilon _{ijk}\,\frac {\partial \rho }{\partial x_j}\,\frac {\partial p}{\partial x_k}}}_{\text {III}} + \underbrace {\overline {\omega _i\epsilon _{ijk}\,\frac {\partial }{\partial x_j}\left ( \frac {1}{\rho }\,\frac {\partial \tau _{kl}}{\partial x_l} \right )}}_{\text {IV}} \mbox {.} \end{equation}
The left-hand-side of (3.9) represents the Lagrangian derivative of the enstrophy. Term I represents vortex stretching, term II represents enstrophy generation/destruction due to fluid element dilatation, term III represents enstrophy production due to baroclinicity, and term IV represents viscous diffusion of enstrophy. The thermodynamic conditions selected for this research can activate certain mechanisms of enstrophy transport (specifically dilatation and baroclinicity), which do not need to be considered in the analysis of incompressible, constant-density flows. The relative importance of these pathways is yet to be precisely characterised for these flows across the literature, and this motivates the present analysis.
Enstrophy budget obtained from the DNS calculations. Results are included for three streamwise stations: (a)
$x=10\,\text{mm}$
, (b)
$x=20\,\text{mm}$
, (c)
$x=30\,\text{mm}$
. The data are ensemble averaged.
Figure 22 presents data showing the contributions of each term (I–IV) to the Lagrangian derivative of enstrophy. All data are obtained at spanwise location
$z = 3\, \rm mm$
, and for three streamwise stations, namely
$x = 10\, \rm mm$
,
$x = 20\, \rm mm$
and
$x = 30\, \rm mm$
. Data are presented as functions of the cross-stream coordinate
$y$
. All data are obtained by ensemble averaging. Although the averages for these higher-order quantities are not fully converged due to computational cost, the dominant trends are still clear, and it is these trends that we discuss here.
The data in figure 22 show both expected and interesting trends. First, it is clear to see that vortex stretching and viscous diffusion dominate the gains to and losses from the enstrophy at each streamwise location. This is an expected result, since vortex stretching is responsible for the energy cascade seen in figure 17, and viscous diffusion is the physical mechanism limiting the unbounded growth of enstrophy. Note also that the ensemble-averaged data at each axial station do not show significant contributions to the budget from either dilatation or baroclinicity, implying that these are not governing physical mechanisms in this set-up. This most likely originates from the lack of a strong mean pressure gradient in the field.
Overall, the results in figure 22 confirm that baroclinicity and dilatation are both not dominant mechanisms governing the vorticity dynamics in the mean, which has not been shown previously at the conditions of interest in this investigation. Note that alternative configurations that permit a dominant mean pressure gradient may produce different results due to the presence of a stronger baroclinic torque term.
4. Conclusions
This research uses both experimental and numerical data to study turbulent mixing in a supercritical CO
$_2$
shear layer. Specifically, near-critical boundary conditions are selected such that the mixing takes place in the presence of thermophysical nonlinearities. Experimental data are obtained using Raman spectroscopy on a bespoke test bench, while the DNS technique is used to obtain complementary numerical results.
Key takeaways are as follows. The shear layer growth rate is approximately 0.2, and further experiments at identical conditions but lower pressure are required to assess the sensitivity of this value to changes in the thermodynamic state. Interestingly, the experimental results indicate that density gradients associated with the mixing of the liquid-like fluid and mixed fluid relaxes more readily as the flow evolves downstream than when associated with mixing of the gas-like fluid and mixed fluid. This finding opposes that found in the literature. The computational results are able to identify a lack of turbulence isotropy in the field, using anisotropy invariance maps (AIMs). Some parts of the shear layer are found to contain turbulence that approaches the one-component limit, indicating statistical attenuation of certain turbulence components. This attenuation is theorised to originate from regions of high density gradient magnitude, which are known to exist in such layers. Enstrophy analysis indicates that baroclinic vorticity amplification is not a dominant phenomenon in this shear layer, despite the large density gradients.
Further data sets at identical flow conditions, but at lower (e.g. atmospheric) pressure, will complement the analysis presented here. Later research can particularly focus on a characterisation of the anisotropy associated with the small scales only, since the AIMs presented here are not able to provide differentiation between large- and small-scale anisotropy.
Acknowledgements
The authors gratefully acknowledge the support of the Department of Energy (DOE), National Energy Technology Laboratory (NETL) and the University Turbine Systems Research (UTSR) programme. This research was also supported in part through research cyberinfrastructure resources and services provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology, Atlanta, Georgia, USA.
Funding
This work was funded in part through the University Turbine Systems Research (UTSR) programme under project FE0031772. Part of this work used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC05-00OR22725.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
Data will be made available upon request.
Author contributions
D.P.: computational data set generation and data curation, writing (overall composition, original draft, review and editing). C.H.L.: experimental data set generation and data curation, writing (original draft, review and editing). A.M.S.: optical diagnostic acquisition, supervision (experimental work), writing (review). D.R.: experimental test bench acquisition, supervision (experimental work), writing (review). J.C.O.: funding acquisition, supervision (computational work), writing (review and editing).
Open access
