1. Introduction
High-speed air-breathing vehicles flying faster than Mach 6 are typically equipped with supersonic-combustion ramjet (scramjet) propulsion systems; see figure 1. If the high-speed flow in the inlet of the combustor of a scramjet propulsion system slowed down to subsonic speed, then temperatures, pressures and densities would be encountered that are too high for an efficient combustion (Roberts Reference Roberts2008), and materials might melt or lose integrity. Therefore, for facilitating efficient combustion and preventing material damage, the inlet flow has to slow down from a high supersonic speed to a suitable, but lower, supersonic speed. In general, this will hold for vehicles that can achieve speeds 5–12 times the speed of sound.
In a scramjet, compressibility effects delay shear layer mixing (Landsberg et al. Reference Landsberg, Wheatley and Veeraragavan2016). To compensate for this effect, long combustors are required to achieve the longer flow residence time required for higher combustion efficiency. However, a longer combustor increases the internal skin friction of the combustor, therewith increasing the overall drag of the engine (Tanimizu et al. Reference Tanimizu, Mee, Stalker and Jacobs2013). A longer combustor also increases the mass of the vehicle. Thus achieving faster mixing and more efficient combustion is of utmost importance for scramjets.
To mix fuel with the cross-flow, several approaches have been proposed (Seiner et al. Reference Seiner, Dash and Kenzakowski2001; Urzay Reference Urzay2018). An often-researched design solution is injecting the gaseous fuel via an under-expanded sonic jet from the wall into the supersonic cross-flow (Mahesh Reference Mahesh2013; Huang Reference Huang2016; Sun et al. Reference Sun, Liu and Hu2019). When such a jet is injected in the cross-flow, a bow shock arises upstream of the penetrating jet. The bow shock decelerates the cross-flow (see figure 2 a). Further downstream, the large-scale coherent structures in the jet upper shear layer (Ben-Yakar et al. Reference Ben-Yakar, Mungal and Hanson2006) and the counter-rotating streamwise vortex pair of the jet plume (VanLerberghe et al. Reference VanLerberghe, Santiago, Dutton and Lucht2000) are instrumental in the mixing of fuel and air. Subsequently, the air–fuel mixture is ignited (Gamba & Mungal Reference Gamba and Mungal2015).
An important parameter of the flow is the jet-to-cross-flow momentum flux ratio
$J$
, given by
\begin{align} J\equiv \frac {\rho _{\kern-1pt j} U_{\kern-1pt j}^2}{\rho _{c} U_{c}^2}=\frac {\gamma _{\kern-1pt j} p_{\kern-1pt j} {M}_{\kern-1pt j}^2}{\gamma _c p_c {M}_c^2}, \end{align}
where a large
$J$
results in a strong jet. Here,
$U$
is the velocity of the flow,
$\rho$
is the density of the fluid,
${M}$
is the Mach number,
$\gamma$
is the ratio of specific heats, and
$p$
is the static pressure. The subscripts
$c$
and
$j$
denote the cross-flow and the jet, respectively. An increase of the vertical spreading of the jet plume is desired for optimal mixing, but increasing
$J$
leads to higher total-pressure losses due to the stronger bow shocks (Gruber et al. Reference Gruber, Nejad, Chen and Dutton1995; Smink et al. Reference Smink, Hoeijmakers and Venner2022a
). Therefore, other options to enhance penetration have been proposed, such as fin-guided injection (Aguilera & Yu Reference Aguilera and Yu2015), extruded jet orifices (Seraj et al. Reference Seraj, Hosseinnejad, Rostamiyan and Fallah2023), coaxial air–fuel injection (Zhang et al. Reference Zhang, Rana, Moradi and Li2021a
) and the use of multiple jets. Multiple laterally arranged jets (e.g. Sebastian & Schreyer Reference Sebastian and Schreyer2022) made the spanwise penetration of the fuel in the cross-flow more uniform, but did not lead to enhanced penetration depth (Maikap & Rajagopal Reference Maikap and Rajagopal2024). A second emanating jet from the opposite wall showed enhanced mixing efficiency by benefiting from the interacting bow shocks (Liang et al. Reference Liang, Sun and Xiong2022). When jets are placed in a tandem configuration (see e.g. Lee Reference Lee2006a
; Pudsey & Boyce Reference Pudsey and Boyce2010; Zhao et al. Reference Zhao, Bian, Li and Ye2019), the bow shocks also enhance the penetration depth. Figure 2(b) shows the structure of the flow in the case of injection of two – in tandem arrangement – under-expanded sonic jets (
${M}_{\kern-1pt j}=1$
) perpendicular to the supersonic cross-flow. The spacing between the centres of the jets (with orifice diameters
$D_1$
and
$D_2$
) is
$\Delta x$
, which is non-dimensionalised using the total diameter
$D_t$
(an equivalent diameter taking into account the injection from two jets, where it is assumed that both jets have the same
$J$
) as the dimensionless distance:
\begin{align} S\equiv \frac {\Delta x}{D_t} = \frac {\Delta x}{\sqrt {D_1^2+D_2^2}}. \end{align}
In front of each of the two jets, a three-dimensional bow shock forms, which decreases the momentum of the cross-flow, and allows the most downstream jet to penetrate deeper into the cross-flow (Lee Reference Lee2006a ; Landsberg et al. Reference Landsberg, Wheatley and Veeraragavan2016).
Supersonic combustion ramjet. By decelerating the flow from high supersonic speed to some lower supersonic speed, the temperature of air entering the combustion chamber will remain below the ignition temperature of the fuel (Smink et al. Reference Smink, Hoeijmakers and Venner2022a ).
Dual-jet injection into a supersonic cross-flow. (a) Example schlieren image of sonic jet injection in
${M}_c=1.6$
air stream. Here,
$D_1=1$
mm,
$D_2=2$
mm,
$S=7.17$
and
$J=3.8$
. (b) Sketch of some flow features, adapted from Lee (Reference Lee2006a
). The plane of visualisation shows the focal plane for experiments, and the plane on which the
$(x,z)$
results in the numerical simulation are visualised.
Several studies have investigated dual-jet injection into a supersonic cross-flow, such as Lee (Reference Lee1994), Lee (Reference Lee2006a
,Reference Lee
b
), Landsberg et al. (Reference Landsberg, Wheatley and Veeraragavan2016), Liu et al. (Reference Liu, Baccarella, McGann, Lee and Do2018) and de Maag et al. (Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2020). Lee (Reference Lee2006a
,Reference Lee
b
) investigated the mixing and combustion characteristics of transverse injection of two equal-sized sonic under-expanded jets into a supersonic cross-flow of
${M}_c=2$
,
$J\in [1.0,2.0]$
. The use of dual jets resulted in a higher mixing rate and a deeper penetration compared to an equivalent single jet. Furthermore, an increase of
$J$
resulted in deeper penetration, but slower mixing. For optimal mixing, for each
$J$
there is an optimal non-dimensional spacing
$S$
between the jets, which increases with increasing
$J$
. Although in the combustor the total-pressure losses increase due to dual-jet injection, the energy losses are subordinate to the mixing enhancement. Landsberg et al. (Reference Landsberg, Wheatley and Veeraragavan2016) studied injection of two sonic under-expanded hydrogen jets at a 45 degrees inclination into a supersonic cross-flow of
${M}_c \in [2.6,5.3]$
, with
$J=1$
, by means of numerical simulations using the Reynolds-averaged Navier–Stokes (RANS) equations. The upstream jet (
$D_1 = ( {1}/{2})D_2$
) served to create a shielding effect for the downstream main jet. The shielding effect improved the mixing characteristics of the jet (i.e. the penetration depth and its spanwise spreading) compared to single-jet injection. This improvement depends on the distance between the jet orifices, while the optimum distance depends on the Mach number of the cross-flow. If the distance between the jets becomes too large, then they act as spatially separated single jets (Pudsey et al. Reference Pudsey, Wheatley and Boyce2015). A relation for the penetration depth of injection of air as a function of
$J$
and
$S$
(
$D_1 = ( {1}/{2})D_2$
) for
${M}_c=1.6$
has been established in an experimental study by Smink et al. (Reference Smink, Hoeijmakers and Venner2022a
,Reference Smink, de Maag, Lerink, Giskes, Hoeijmakers, Venner, Segerink and Offerhaus
b
).
The bow shock interacts with the flow in the subsonic part of the incoming wall boundary layer, which results in separation of the boundary layer upstream of the bow shock. Outside the boundary layer, the resulting increase in thickness of the boundary layer causes the formation of the (oblique) separation shock that intersects the bow shock. The recirculating flow in the separated flow region forms a horseshoe vortex wrapping as a scarf around the bow shock (Lee Reference Lee2006a ) (see also figure 2 b).
Inside the first portion of the plume of the jet, the interaction of the under-expanded sonic jet with the surrounding flow results in a barrel shock. The jet flow inside the barrel shock accelerates to supersonic speeds. The barrel shock is closed off by a Mach disc that decelerates the jet flow. Furthermore, the jet is forced to tilt from a normal direction at the exit of the jet to a direction more parallel to the supersonic cross-flow. Due to the interaction with the main stream, the edge of the tilted jet deforms and rolls up in a counter-rotating pair of longitudinal vortices; see figure 2(b) (Lee Reference Lee2006a ). The mixing of the fuel injected by the jets and the cross-flow fluid depends strongly on these flow structures (VanLerberghe et al. Reference VanLerberghe, Santiago, Dutton and Lucht2000). The resolution of the highly time-dependent flow features in the jet plume (Ben-Yakar et al. Reference Ben-Yakar, Mungal and Hanson2006; Zhang et al. Reference Zhang, McCreton, Awasthi, Wills, Moreau and Doolan2021b ) requires time-resolved results with a time scale down to the order of microseconds. However, this is experimentally challenging since it requires advanced recording techniques. In the present study, experimental work is complemented by numerical studies, providing deeper understanding of the time-dependent flow features and the three-dimensional aspects of dual-jet injection into a supersonic cross-flow. How large are the variances of the location of the jet shear layer in time, and is its behaviour captured in the experiments? In addition, the validity of the experimental results will be assessed based on the fundamental flow equations in the modelling. The time-averaged results are important for characterisation of the efficiency of jet mixing, needed e.g. for scramjet designs. However, do we need time-resolved solutions for assessing the time-averaged behaviour, or is an averaged method such as RANS adequate?
Therefore, the present study comprises experimental work, using schlieren flow visualisation in the test section of the supersonic wind tunnel at the University of Twente, complemented with numerical simulations of the flow inside the test section of the wind tunnel, using the open-source code SU2 (Palacios et al. Reference Palacios, Colonno, Aranake, Campos, Copeland, Economon, Lonkar, Lukaczyk, Taylor and Alonso2013; Economon et al. Reference Economon, Palacios, Copeland, Lukaczyk and Alonso2016). The RANS equations and (time-resolved) Delayed Detached Eddy Simulations (DDES), a hybrid unsteady RANS-Large-Eddy Simulation (LES) method, are used for the numerical simulations. The goal is to study dual injection of a tandem of under-expanded sonic air jets (
${M}_{\kern-1pt j} = 1.0$
) into a supersonic cross-flow of Mach 1.6, to compare results of RANS and DDES with experimental results, and to provide deeper understanding of the flow phenomena, which was not possible with experiments alone. The downstream jet, diameter
$D_2$
, is the main jet, where
$D_2 = 2D_1$
. The jet-to-cross-flow momentum flux ratio
$J$
in (1.1) is varied, and for given
$J$
, an optimal dimensionless distance
$S$
(1.2) between the two jets is focused on. The injection of a single jet, investigated in Giskes et al. (Reference Giskes, Verschoof, Segerink and Venner2016) and de Maag et al. (Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2020), is taken as reference. Numerical results are analysed in detail, and wind-tunnel schlieren images are compared with numerical equivalents (Luthman et al. Reference Luthman, Cymbalist, Lang, Candler and Dimotakis2019), or rather an approximation, coined ‘numerical schlieren’ images. The latter is the distribution of a derivative, i.e. in the
$x$
- or
$z$
-direction of the density field in an
$(x,z)$
-plane.
The method section (§ 2) describes the experimental set-up (§ 2.1) and the numerical method (§ 2.2). The results section (§ 3) describes results of the the wind-tunnel experiments (§ 3.1), time-averaged numerical results (§ 3.2), the comparison between numerical and experimental results (§ 3.3), and time-resolved numerical results (§ 3.4). Conclusions are presented in § 4.
2. Method
2.1. Experimental set-up
In the present study, flow structures are visualised employing the schlieren set-up (Settles Reference Settles2001) installed in the supersonic wind tunnel at the University of Twente; see figure 3.
Set-up of the University of Twente continuous-flow supersonic wind tunnel with schlieren visualisation system. The flow is from right to left. The atmospheric inlet, with silencer, is not visible in the picture. The exit of the channel at the left is connected to a (0.1 atm) vacuum pump system.
This supersonic flow facility is capable of achieving flows up to Mach 1.7 continuously. In the present study, the flow in the test section (located between lens 1 and lens 2 of the schlieren system) is set to Mach 1.6. The total pressure and total temperature of the cross-flow are equal to ambient conditions (
$p_{0,c} = 101$
kPa,
$T_{0,c} = 292$
K). The test section measures
$w\times h = 45$
mm
$\times$
48 mm perpendicular to the upstream flow direction, and has length
$\sim 100$
mm. The upstream and downstream sonic jets (
${M}_{\kern-1pt j}=1$
) have diameters
$D_1=1$
mm and
$D_2 = 2$
mm, respectively. Therefore, the equivalent total diameter is
$D_t \approx 2.23$
mm.
A schlieren system visualises density gradients in a transparent medium (Settles Reference Settles2001). The light rays from the light source pass through lens 1 (focal distance
$f=0.5$
m), which results in a parallel bundle entering the test section through a high-quality-glass window. In the case of undisturbed constant-density flow, the bundle leaves the test section through the glass window on the opposite site. Then lens 2 (focal distance
$f=1.0$
m) focuses the light into a focal point. A variation in the density field inside the test section refracts the light, with the direction of the light rays depending on the local density gradient. This light leaves the test section through the glass window on the opposite side of the test section, and subsequently passes through lens 2. However, now not all refracted light passes through the focal point; a part of the refracted light meets a barrier, the ‘knife’, located in the same plane as the focal point (see figure 3). This results in a dark region in the image. Another part of the refracted light will not be blocked by the knife, resulting in a lighter region in the image. These darker and lighter regions therewith visualise features that caused the refraction, i.e. variations in the density field. A horizontal knife edge is used in the experiments, resulting in the visualisation of the derivative of the density in the vertical direction, i.e. the
$z$
-direction.
A Phantom V611 camera equipped with a Nikon TC-20 E AF-S Teleconverter II plus Nikkor AF-S 200 mm f/2.0 ED VR lens captures the schlieren images. For obtaining high-quality images, a vertical-cavity-surface-emitting laser (VCSEL) from Tyson Technology is used (de Maag et al. Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2020). The emitted wavelength is
$\lambda =808$
nm, with power
$P = 4.5$
W emitted from a laser surface of 1.19 mm
$\times$
1.19 mm (de Maag et al. Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2019). The camera frame rate is limited to 1000 fps. This limitation would result in motion blur due to rapidly propagating flow features in combination with a long exposure time. Therefore, the VCSEL is pulsed with frequency 1 kHz and pulse width 100 ns. This is easily captured by the camera, whose shutter is open for a long period (990
$\unicode{x03BC}$
s). Effectively, the exposure time is reduced to 100 ns, preventing motion blur. For more details on the VCSEL, see de Maag et al. (Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2019, Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2020).
The present investigation employs a largely automated post-processing procedure developed by Smink et al. (Reference Smink, Hoeijmakers and Venner2022a
) that, for each
$(S,J)$
combination, can handle a relatively large number of schlieren images fast and reliably. The first step of this procedure, which is implemented in
$\textrm {MATLAB}^{\circledR }$
, is to isolate the jet plume by determining a window of interest, providing a binary mask that excludes all flow phenomena that can cause noise in the images, such as bow shocks, reflected shocks and boundary layers. This window is applied to all images of the same experiment. Second, for every image, the flow features of the jet, which are visible due to local density differences, are captured. This is achieved by first increasing the contrast of the image (function localcontrast with edge threshold 0.3, and enhancement 1; function imadjust; function adapthisteq with contrast enhancement limit
$0.4\pm 0.05$
,
$8\times 16$
tiles) and thereafter applying binarisation (function imbinarize with locally adaptive image threshold setting and sensitivity
$0.36\pm 0.04$
). The noisy pixels are filtered from the image using the function bwmorph, resulting in an image with a collection of coherent flow features. These data are collected from thirty consecutive images, 1 ms apart in time, for each
$(S,J)$
combination. The end result of the procedure, which manipulates pixel-wise light-intensity maps in the region of the jet plume, is the iso-incidence plot. Here, per pixel of each image, the number of times a feature occurs is determined, which is equivalent to a time-averaged quantity. More details can be found in Smink et al. (Reference Smink, Hoeijmakers and Venner2022a
).
2.2. Numerical simulations – wind-tunnel test section (RANS and DDES)
The open-source code SU2 (Palacios et al. Reference Palacios, Colonno, Aranake, Campos, Copeland, Economon, Lonkar, Lukaczyk, Taylor and Alonso2013; Economon et al. Reference Economon, Palacios, Copeland, Lukaczyk and Alonso2016) has been used for the numerical simulations, at three
$(S,J)$
combinations (
$S=7.17$
and
$J\in \{2.8, 3.8, 4.8\}$
), making use of both a RANS simulation and DDES. The DDES are a hybrid LES/unsteady RANS (URANS) model of turbulent flow, developed by Spalart et al. (Reference Spalart, Jou, Strelets and Allmaras1997, Reference Spalart, Deck, Shur, Squires, Strelets and Travin2006). The turbulence model used for the RANS simulations is the standard Spalart–Allmaras (SA) model (with transport variable
$\tilde {\nu }$
, which is related to the eddy viscosity). This model is chosen because of the compatibility with the DDES. The DDES version is a further development of the original DES model, which reduces the dependency of the numerical solution on the computational grid. The model uses a URANS formulation near the wall, but an LES formulation away from the wall. By using a URANS formulation near the wall, the huge grid requirements of fully wall-resolved LES are circumvented, while still capturing the transient behaviour of the turbulence in the solution. This is accomplished by modifying, away from the wall, the original URANS turbulence model, i.e. the SA model. For more details, see Spalart et al. (Reference Spalart, Jou, Strelets and Allmaras1997, Reference Spalart, Deck, Shur, Squires, Strelets and Travin2006).
Computational domain of (a) the entire wind tunnel and (b) the test section. (a) Geometry of the entire wind tunnel, which provides the inlet condition for the test section. The shown plane with the Mach number distribution is the location of the inflow boundary of the computational domain of the test section. (b) Grid of the midplane of the test section of the wind-tunnel channel. The size of the test section is
$(x,y,z) = (81,45,48)\ \text{mm}=(36.3D_t, 20.2D_t, 21.5D_t)$
. The
$y$
-coordinate points inwards. The refined grid is twice as fine in all dimensions. Here,
$S=7.17$
.
For obtaining detailed results of the numerical simulations, a large number of grid points are still needed, especially for the DDES calculations. Therefore, it is desirable to limit the computational domain as much as possible, here to only the test section. For the then required inflow condition, an exploratory full-channel RANS simulation has been performed, confirming the
${M}_c=1.6$
inflow condition (see figure 4
a). From the solution plane in the simulation of the entire wind tunnel (figure 4
a), the part corresponding to the start of the test section is extracted, and this solution is used as inflow condition for the simulations on the grid of the computational domain limited to the test section (figure 4
b). The justification for this approach is that in the wind tunnel, the core flow is supersonic and the boundary layers are nicely attached. The boundary layer thickness is 3 mm. A comparison between the solution for the empty test section (figure 4
b) using this approach and the solution for the entire wind tunnel (figure 4
a) showed negligible differences. In addition, for the incoming cross-flow, the total pressure and total temperature are equal to ambient conditions (
$p_{0,c} = 101$
kPa,
$T_{0,c} = 292$
K), identical to the experimental conditions (§ 2.1). As for the experiments, the upstream and downstream sonic jets have diameters
$D_1=1$
mm and
$D_2 = 2$
mm, respectively. The flow at the exit is still supersonic in the core flow, hence no boundary conditions are needed there. However, near the walls, the flow is subsonic and the static pressure is described in those regions to keep the solution stable. This pressure is taken equal to 90 % of the static pressure at the inlet of the test section. Variation of this value showed that the solution upstream of the outlet is insensitive to it.
The computational grid used for the simulation of the flow inside the test section is a multiblock structured grid composed of hexahedra. This results in better grid quality than unstructured grids. The four wind-tunnel walls are specified as no-slip walls. For the original grid, there are 54 497 635 points and 54 050 816 volume elements; for the refined grid, there are 434 191 045 points and 432 406 528 volume elements. So the refined grid has a uniform refinement of factor 2 in each of the three directions.
Schlieren images are typically used for the investigation of flows in wind tunnels (as well as in flight experiments) (Settles & Hargather Reference Settles and Hargather2017). In experiments, schlieren images of the flow field in the wind-tunnel section are created optically due to the non-uniform density field (see § 2.1). Therefore, the schlieren images of compressible flow visualise shock waves, expansions, contact surfaces, jets and boundary layers. In the present study, the
$(x,z)$
-plane is chosen as the vertical plane, parallel to the main flow direction, which passes through the middle of the test section, designated as midplane. In this respect, the distribution of
$\partial \rho /\partial z$
is associated with the horizontal position – parallel to the
$x$
-axis – of the knife edge of the wind-tunnel schlieren system (see § 2.1). The numerical schlieren images reveal, as the wind-tunnel schlieren images, most of the topology of the compressible flow inside the wind-tunnel test section. In addition, the numerical solution provides detailed flow-field distributions of all three velocity components, Mach number, total pressure and derived quantities such as vorticity.
2.2.1. Wind-tunnel test section simulation – RANS
Before DDES are performed, a (steady) solution based on the RANS equations, in combination with the SA turbulence model, is obtained. The chosen spatial discretisation scheme is a second-order vertex-based upwind finite-volume method in which the gradients of the flow variables are computed using a weighted least squares algorithm in combination with the Wang (Reference Wang2000) formulation of the Venkatakrishnan limiter and Roe’s approximate Riemann solver (Roe Reference Roe1981) for the advective fluxes. For the viscous terms, a central discretisation scheme has been employed. The iterative algorithm is a quasi-Newton method in which the linear system of algebraic equations, a result of an approximate linearisation, is solved by a flexible generalised minimal residual algorithm. In the present application, this iterative algorithm typically reduces the residuals three orders of magnitude before oscillations set in, indicating that the flow shows (some) transient effects, even within the RANS approximation. However, inspection of the solutions at multiple points in the iteration process shows that the transients in the larger flow structures are negligible. The obtained results of this simulation serve two goals: (i) to be used as an initial solution for the DDES (see § 2.2.2), and (ii) to compare the time-averaged DDES results with the RANS solution. The flow features occurring because of jet injection into a supersonic cross-flow are of an unsteady nature, therefore it will be investigated whether the RANS solutions are similar to the time-average of the time-resolved DDES solutions.
2.2.2. Wind-tunnel test section simulation – DDES
The time-dependent DDES on the test-section grid were initialised using a (steady) RANS solution (see § 2.2.1), followed by a restart using DDES, with time step 1
$\unicode{x03BC}$
s on the original grid, and 0.5
$\unicode{x03BC}$
s on the refined grid. The time integration scheme used is the second-order implicit backward Euler scheme (Jameson Reference Jameson1991). The spatial discretisation is the same as for the RANS simulations, described in § 2.2.1. To wash out the transients caused by the restart from the RANS solution, 4000 time steps of 1
$\unicode{x03BC}$
s on the original grid and 8000 time steps of 0.5
$\unicode{x03BC}$
s on the refined grid were performed, corresponding to 25 test section convective flow-throughs. Thereafter, the actual data were collected for 2 ms in real time, which corresponds to 2000 and 4000 time steps on the original and refined grids, respectively. The choice of 2 ms is justified in § 2.2.3. Based on the free-stream velocity, 2 ms corresponds to approximately 13 test section convective flow-throughs.
For every time step, the system of nonlinear equations is solved using the dual time-stepping procedure in combination with an implicit Euler scheme, solving the linearised equations with two steps of the generalised minimal residual method (Jameson Reference Jameson1991). A total of 20 inner iterations were sufficient to converge the solution of the nonlinear equations for every time step 2–3 orders of magnitude. The simulation was run on the Snellius machine of the Dutch Supercomputer Centre. For the simulation on the original grid, 384 cores (2 nodes) were used, and for the simulation on the refined grid, 3072 cores (16 nodes) were employed. The size of the resulting data set for the time-resolved series of 2
$\unicode{x03BC}$
s is 7 TB for the original grid, and 112 TB for the refined grid (refinement of factor 2 in both space and time). For the visualisation of the flow solution, the open-source post-processing tool ParaView (Ayachit Reference Ayachit2015) has been used.
2.2.3. Convergence of time-averaged DDES
In this subsubsection, it is investigated whether the grid convergence both in space and time of the DDES results is sufficient. Therefore, the results for time-averaged solutions and for different grid refinements are compared.
Comparison of the DDES results (time step 1
$\unicode{x03BC}$
s on the original grid, time step 0.5
$\unicode{x03BC}$
s on the refined grid) in the
$(x,z)$
midplane, time-averaged over 1 and 2 ms. The first and second rows compare time averaging on the refined grid for 1 and 2 ms. The second and third rows compare the time-averaged results on the refined and original grid. (a) Numerical schlieren images (
$\partial \rho /\partial z$
). (b) Vorticity magnitude (
$|\boldsymbol{\omega }|$
) in s
$^{-1}$
. Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
Figure 5. Long description
Panel A: Three numerical schlieren images are displayed. The first row shows time-averaged results over 1 millisecond on a refined grid. The second row shows time-averaged results over 2 milliseconds on a refined grid. The third row shows time-averaged results over 2 milliseconds on the original grid. The density gradient is represented with a grayscale color bar ranging from -100 to 100. Panel B: Three vorticity magnitude images are displayed. The first row shows time-averaged results over 1 millisecond on a refined grid. The second row shows time-averaged results over 2 milliseconds on a refined grid. The third row shows time-averaged results over 2 milliseconds on the original grid. The vorticity magnitude is represented with a color bar ranging from 0 to 5 times 10 to the power of 4.
Comparing the results for time averaging over 1 and 2 ms (first and second rows in figure 5, respectively), it becomes clear that the results are almost identical. A slight difference is observed in the vorticity magnitude, more pronounced due to the chosen small range for the scale. But all in all, the differences between time averaging for 1 and 2 ms are almost negligible, which confirms that 2 ms is a sufficient time span for time averaging of DDES the results.
When comparing the results for the refined and the original grid (second and third rows in figure 5, respectively), it appears that the time-averaged DDES results are already well represented on the original grid. Of course, the capturing of the shock waves is spread over a somewhat larger distance in the original grid compared to that in the refined grid, but compared to the original grid, there is not much more detail visible in the refined grid. Furthermore, the different flow features are almost the same when overlaying the results. Therefore, it is justified to use the original grid for the calculations, especially for obtaining the averaged flow field, as it provides results of adequate accuracy, at only a fraction (less than 10 %) of the computational cost for a simulation using the refined grid.
3. Results
3.1. Results – experimental study
The schlieren visualisation of the injection of dual sonic jets has been carried out for a number of conditions. Figure 6 presents, for
$J = 3.8$
, one sample schlieren image for each value of
$S$
that has been investigated. These results show that for smaller values of
$S$
(
$S\lt 4.732J^{0.288}$
,
$\sim 7$
for
$J=3.8$
), the shielding effect due to the upstream jet increases with
$S$
for which the penetration height reaches a maximum. For still higher values of
$S$
, the shielding effect decreases with increasing
$S$
.
Schlieren images for
$J=3.8$
, for a range of values of
$S$
. For every
$(S,J=3.8)$
combination, one of the 30 images available is shown. Here,
${M}_c=1.6$
,
$D_1 = 1$
mm,
$D_2= 2$
mm. Red square: schlieren image for
$(S=7.17,J=3.8)$
, featuring in figure 2 and most other figures in § 3. The vertical dashed white line indicates the location of the main orifice.
The upper boundary of the jet plume, with non-zero values of the incidence (see § 2.1), represents the time-averaged location of the upper edge of the plume, which in the present study is used as measure for the local penetration depth of the jet. A parameter study on single-jet injection has been carried out by Papamoschou & Hubbard (Reference Papamoschou and Hubbard1993), and the corresponding empirical relations for the penetration of a jet into a supersonic cross-flow as a function of especially
$J$
and
${M}_c$
have been studied extensively by e.g. Portz & Segal (Reference Portz and Segal2006) and Fries et al. (Reference Fries, Ranjan and Menon2021). For dual-jet injection, empirical relations are not readily available. However, in previous work by Smink et al. (Reference Smink, Hoeijmakers and Venner2022a
,Reference Smink, de Maag, Lerink, Giskes, Hoeijmakers, Venner, Segerink and Offerhaus
b
), the penetration depth has been found to satisfy an empirical relation that fits the data from all
$(S,J)$
combinations considered at
${M}_c=1.6$
:
\begin{align} \frac {z}{D_t} \left (\frac {x}{D_t} ;S,J\right ) & =0.432J^{0.461} \left[2+\textrm{e}^{-\tfrac {1}{2} \left [\tfrac {S-4.732J^{0.288}}{6.103J^{-0.521} }\right ]^2}\right] \nonumber \\ & \quad \times \left [\left (\frac {x}{D_t}+0.68\right )^{0.333}+0.80\right ] \!, \end{align}
with
$z/D_t$
the (time-averaged) vertical penetration depth as a function of the horizontal distance
$x/D_t$
from the centre of the downstream jet orifice, for given
$J$
(see (1.1)) and
$S$
(see (1.2)). The relation is based on data for
${M}_c=1.6$
,
$2.8\lt J\lt 4.8$
and
$0\lt S\lt 9.87$
. This relation describes the data satisfactorily with a coefficient of determination
$R^2=0.73$
(satisfactory
$0.65\leqslant R^2\leqslant 0.8$
, good
$R^2\gt 0.8$
), where
$R^2$
is defined as
\begin{align} R^2\equiv 1- \frac {\sum ^N_{i=1}(f(x_i)-z_i)^2}{\sum ^N_{i=1} (\bar {z}-z_i)^2},\quad \text{with} \ \bar {z}=\frac {1}{N}\sum ^N_{i=1}z_i. \end{align}
Here,
$(x_i,z_i),\ i=1,2, \ldots, N$
, are the coordinates of
$N$
data points,
$\bar {z}$
is the mean value of all values of
$z_i$
, and
$f(x_i)$
is the function value of the fit evaluated at
$x_i$
.
An overall measure of the penetration depth of the upper jet shear layer of the plume is the spatial average of (3.1) from the downstream edge of the main orifice at
$x\approx 0.5D_t$
to
$x=15D_t$
downstream. The upper value
$15D_t$
is chosen such that the shock reflections from the upper and side walls do not have an influence on this measure. The average penetration depth yields
\begin{align} \frac {z_{avg}}{D_t} (S,J) & \equiv \frac {1}{14.5} \int _{0.5}^{15} \frac {z}{D_t} \left (\frac {x}{D_t} ;S,J\right )\, {\rm d}\left (\frac {x}{D_t}\right ) \nonumber \\ & = 1.194J^{0.461}\left [2+\textrm{e}^{-\tfrac {1}{2} \left [\tfrac {S-4.732J^{0.288}}{6.103J^{-0.521}} \right ]^2} \right ]\!. \end{align}
Because of the shielding effect of the jet from the upstream orifice, for the dual-jet injection, the penetration depth of the jet plume into the cross-flow exceeds the penetration depth of a single-jet configuration, at the same
$J$
. However, for larger distances
$S=\Delta x/D_t$
, the shielding effect decreases and the penetration depth of the main jet plume decreases towards that of a single jet with diameter
$D_2$
. Similarly, for
$S\rightarrow 0$
, the value converges to that of a single jet with diameter
$D_t$
. Therefore, there is an optimum value of
$S$
for which the penetration will be maximal (Landsberg et al. Reference Landsberg, Wheatley and Veeraragavan2016; Liu et al. Reference Liu, Baccarella, McGann, Lee and Do2018), equal to
$S_{opt}$
, which increases with increasing
$J$
(Lee Reference Lee2006a
). It follows from (3.3) that
\begin{align} S_{opt}(J) = 4.732 J^{0.288},\quad \text{ for which } \frac {z_{avg,opt}}{D_t}(J) = 3.582 J^{0.461}. \end{align}
The
$(S,J)$
parameter plane of the experimental study. Black open circles: de Maag et al. (Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2019, Reference de Maag, Hoeijmakers, Venner, Segerink and Offerhaus2020). Red open circles: Smink et al. (Reference Smink, Hoeijmakers and Venner2022a
,Reference Smink, de Maag, Lerink, Giskes, Hoeijmakers, Venner, Segerink and Offerhaus
b
). Big red circles: selected conditions for RANS numerical simulations and DDES. Solid greyscale curves: iso-contours of
$z_{avg}/D_t$
from (3.3), with increasing brightness, for
$z_{avg}/D_t = 2.5,3.0,3.5,\ldots , 8.0$
. Dashed black line:
$S_{opt}(J)$
from (3.4).
Figure 7 provides an overview of the
$(S,J)$
combinations for which schlieren images have been obtained. The data used to develop the empirical scaling relation for
$z/D_t$
, given in (3.1), have been extracted from the schlieren images based on the post-processing procedure described in § 2.1, which utilises iso-incidence plots. The resulting empirical scaling relation (3.3) for the overall measure
$z_{avg}/D_t$
of the penetration depth has been plotted in figure 7 in terms of iso-contours in the
$(S,J)$
parameter plane, for values of
$z_{avg}/D_t$
between 2.5 and 8, with steps of 0.5. Included in figure 7 is the dashed curve for the optimal value of
$S$
equal to
$S_{opt} (J)$
given in (3.4). Figure 7 shows that the region in the
$(S,J)$
-plane for which the
$(S,J)$
combination is near-optimal becomes narrower with increasing strength
$J$
of the jets. This indicates that the choice of the parameters for the near-optimal
$(S,J)$
combination becomes more critical with increasing
$J$
.
3.2. Time-averaged results – numerical simulations (RANS and DDES)
In this subsection, RANS results are compared with time-averaged DDES results for
$J=3.8$
in § 3.2.1. Thereafter, for different cutting planes, the time-averaged DDES results are discussed for
$J = 2.8$
,
$3.8$
and
$4.8$
in § 3.2.2.
3.2.1. Comparison of RANS and time-averaged DDES results
To assess the need of time-resolved solutions for the representation of time-averaged results, the results for RANS are compared with the time-averaged DDES results. Figure 8 presents, for the
$(x,z)$
midplane, the (steady) RANS solution (first row) and the time-resolved DDES solution (second row), time-averaged over the final 2 ms in real time (approximately 13 convective flow-throughs) carried out in the numerical simulation. The vorticity magnitude and total-pressure distribution are the most striking representations for showing the differences between RANS and DDES. The RANS solutions show behaviour similar to the time-averaged DDES results, except for the jet plume. Here, the time-averaged DDES solution shows mixing of the jet plume with the cross-flow, whereas in the RANS solutions, the jet plume remains separated from the cross-flow through a layer with a large density gradient, resulting in significant total-pressure losses (see figure 8
b) and sharp vorticity shear layer (see figure 8
a). This stresses that although the RANS solution shows most features satisfactorily, the DDES solution describes the mixing of the jet with the cross-flow more accurately, as these results correspond with the experiments (see § 3.3). The RANS simulations predict a deeper penetration into the cross-flow, which could be due to the larger diffusion as a result of higher eddy viscosity. Therefore, for investigation of the jet penetration and mixing with the cross-flow, it is important to have time-resolved solutions such as DDES.
Comparison of the time-averaged (over 2000 time steps of 1
$\unicode{x03BC}$
s) DDES solution and (steady) RANS solution in the
$(x,z)$
midplane in the test section, for the original grid. (a) Vorticity magnitude (
$|\boldsymbol{\omega }|$
) in s
$^{-1}$
. (b) Total pressure
$p_0$
in Pa. Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
Figure 8. Long description
Panel A: Two heat maps showing vorticity magnitude in the midplane of a test section. The top heat map represents the RANS solution, while the bottom heat map represents the time-averaged DDES solution. The color scale ranges from 10^4 to 5*10^4, with darker colors indicating lower vorticity and lighter colors indicating higher vorticity. The heat maps show the distribution of vorticity around two jets, with noticeable differences in the patterns between the RANS and DDES solutions. Panel B: Two heat maps showing total pressure in the midplane of a test section. The top heat map represents the RANS solution, while the bottom heat map represents the time-averaged DDES solution. The color scale ranges from 2*10^4 to 1.4*10^5 Pa, with darker colors indicating lower pressure and lighter colors indicating higher pressure. The heat maps show the distribution of total pressure around two jets, with noticeable differences in the patterns between the RANS and DDES solutions.
This observation is confirmed by figure 9, which presents results in the
$(y,z)$
exit plane of the test section, for again the RANS solution, together with the DDES solution time-averaged over 2 ms. Figure 9 shows (i) a colour plot of the distribution of the
$x$
-component of vorticity (
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
) and (ii) a colour plot of the distribution of the total pressure. Figure 9(a) shows that though the time-resolved distribution of
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
is not left–right symmetric at all (as will be discussed in § 3.4.2), the time-averaged distribution of
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
in figure 9(a) is nearly perfect left–right symmetric (with positive
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
on the left, and negative
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
on the right). The distribution of the total pressure, presented in figure 9(b), also nicely shows the cross-section of the jet plume as a mushroom-like structure, corresponding to the vortex dipole structure in figure 9(a). This time-averaged structure, and the secondary vortical structures at the bottom wall, are stationary (showing primary and secondary separation of the boundary layer). These are caused by the vortical-flow regions being close to the lower wall, in effect causing the presence of a virtual, mirror-imaged, vortical structure below the wall, such that the zero-normal-velocity boundary condition is satisfied at the wall.
The RANS solution in figure 9 shows a similar mushroom-like vortex dipole structure with positive
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
on the left and negative
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
on the right, but somewhat further away from the bottom wall. This was also observed upon inspection of figure 8, where the penetration depth of the jet was larger for RANS than for time-averaged DDES. As a result of the larger distance from the wall for RANS, the boundary-layer separation induced by the vortex dipole induced flow is weaker than for time-averaged DDES, where RANS only shows primary separation. The comparison of RANS and time-averaged DDES results in figure 9 stresses the importance of time-resolved simulations in order to obtain accurate time-averaged results.
Comparison of the time-averaged DDES (over 2000 time steps of 1
$\unicode{x03BC}$
s) and (steady) RANS flow solution in the
$(y,z)$
exit plane, both for the original grid. (a) Vorticity
$x$
-component (
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
) in s
$^{-1}$
; red indicates positive (anticlockwise), and blue indicates negative (clockwise). (b) Total pressure
$p_0$
in Pa. Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
3.2.2. Time-averaged DDES results for varying
$J$
Figure 10 presents the time-averaged DDES results for the test section (
$(x,z)$
-plane). Four types of plots are shown: numerical schlieren image, colour plot of the vorticity magnitude distribution (
$|\boldsymbol{\omega }|$
), colour plot with the distribution of local Mach number, and colour plot of the total-pressure distribution
$p_0$
. The shock wave pattern, including the
$\lambda$
-shock at the foot of each of the two primary bow shocks, the shock reflections at the walls, and the barrel shocks, are grossly stationary. Figure 10 shows that the unsteady internal dynamics of the jet plumes of the dual-jet injection is averaged out within the envelope of the jet plumes, when compared with the transient results in § 3.4.
Comparing the time-averaged DDES results for
$J=3.8$
(second row in figure 10) with those for
$J=2.8$
(first row in figure 10) and
$J=4.8$
(third row in figure 10), some differences become apparent. As expected, an increase in
$J$
increases the penetration depth of the (main) jet into the cross-flow. Consequently, the bow shocks in front of the jets are steeper and will therefore be stronger and lead to higher total-pressure losses (figure 10
d). The vorticity magnitude within the jet plume is higher for lower
$J$
because of the jet being more confined.
Comparison of time-averaged (over 2000 time steps of 1
$\unicode{x03BC}$
s) DDES flow solution in the
$(x,z)$
midplane in the test section for the original grid for
$J = 2.8$
,
$J=3.8$
and
$ J= 4.8$
. (a) Numerical schlieren images (
$\partial \rho /\partial z$
). (b) Vorticity magnitude (
$|\boldsymbol{\omega }|$
) in s
$^{-1}$
, (c) Local Mach number
${M}$
. (d) Total pressure
$p_0$
in Pa. Here,
${M}_c = 1.6$
and
$S=7.17$
.
Comparison of time-averaged (over 2000 time steps of 1
$\unicode{x03BC}$
s) DDES flow solution in the
$(y,z)$
exit plane for the original grid for
$J = 2.8$
,
$J=3.8$
and
$ J= 4.8$
. (a) Vorticity
$x$
-component (
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
) in s
$^{-1}$
; red indicates positive (anticlockwise), and blue indicates negative (clockwise). (b) Total pressure
$p_0$
in Pa. Here,
${M}_c = 1.6$
and
$S=7.17$
.
Figure 11 presents the time-averaged DDES results for the
$(y,z)$
exit plane for varying
$J$
. Comparing the results for different
$J$
, the observation from figure 8 that the jet plume penetrates further into the cross-flow with increasing
$J$
is confirmed. In the exit plane, it also becomes clear that the total-pressure losses within the mushroom-shaped core of the jet plume become higher for increasing
$J$
. At the same time, for larger
$J$
, the mushroom vortex pair becomes stronger, resulting in a more pronounced left–right symmetry observed for increasing
$J$
. The vorticity magnitude in the jet plume for
$J=2.8$
(see figure 10
b) is also significantly larger than for higher values of
$J$
. This might be the reason for the less-symmetric shape of the jet plume for
$J=2.8$
compared to the shape for larger values of
$J$
.
The DDES flow solution in the test section. Time-averaged (over 4000 time steps of 0.5
$\unicode{x03BC}$
s) flow solution on the refined grid, numerical schlieren image (
$\partial \rho /\partial z$
) in the
$(x,y)$
cutting plane, for (a)
$z/D_t =0.90$
, (b)
$z/D_t =3.14$
, (c)
$z/D_t =6.73$
, (d)
$z/D_t =10.76$
. Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
Figure 12. Long description
Panel A: A grayscale image showing the flow pattern with a density gradient scale at the bottom. The image depicts two jets injecting into a supersonic cross-flow, creating distinct flow structures. Panel B: Another grayscale image with a similar density gradient scale, showing a different flow pattern with the jets injecting at an angle. Panel C: A grayscale image with a density gradient scale, illustrating a more spread-out flow pattern. Panel D: The final grayscale image with a density gradient scale, showing a different configuration of the dual-jet injection into the supersonic cross-flow.
Figure 12 presents, for four
$(x,y)$
cutting planes (
$z/D_t \in \{0.90,3.14,6.73,10.76\}$
), the numerical schlieren image of the DDES solution, time-averaged over the final 2 ms of the simulation. The planes cut through the bow shocks and their reflections on the walls. The white shock waves in figures 12(c,d) are the reflections from the upper wall. The shock waves reflect at the wall via a Mach reflection, and intersect the jet again. Here, a non-trivial reflection at the jet appears, visible in the subsequent cutting planes. The jet penetrates deeper than
$z/D_t=6.73$
, as is clear from figure 12(c).
3.3. Comparison of numerical and experimental results
Iso-incidence plots obtained from (a,b) experimental and (c–e) (numerical) schlieren images (
$\partial \rho /\partial z$
) generated from the time-resolved DDES solution. Each plot is based on 30 snapshots. The value of
$J$
is varied: (a,c)
$J=2.8$
, (b,d)
$J=3.8$
, (e)
$J=4.8$
. The scale gives the iso-incidence values. Black solid curves indicate the similarity relation for the time-averaged location of the upper edge of the jet shear layer in (3.1). Coefficients of determination are (a)
$R^2= 0.70$
, (b)
$R^2= 0.73$
, (c)
$R^2= 0.72$
, (d)
$R^2= 0.81$
, (e)
$R^2= 0.79$
. Here,
${M}_c = 1.6$
and
$S = 7.17$
.
In the present study, the post-processing procedure for experimental schlieren images (see § 3.1) has also been applied to the numerical schlieren images. The obtained iso-incidence plot from both the experimental and numerical results (all for the original grid) are presented in figure 13. Figures 13(a,b) are iso-incidence plots that have been generated from the experimental data using 30 consecutive images with 1 ms time spacing. Figures 13(c–e) are iso-incidence plots that have been generated from numerical schlieren images obtained from the time-resolved DDES flow solution, using 30 consecutive images with 60
$\unicode{x03BC}$
s time spacing as input. The time step 60
$\unicode{x03BC}$
s is chosen such that the schlieren images are not correlated in time (convective time scale
$D_2/U_c \approx 5$
$\unicode{x03BC}$
s), to follow an approach similar to that for the experimental results. In the iso-incidence plot, at a specific location in the jet plume, the colour indicates the incidence – an integer between 0 and 30 – i.e. the number of times a feature is present at a certain location. When no features have been present in the snapshots (incidence of 0), the background schlieren snapshot is visible. The iso-incidence plot reveals either largely stationary or frequently occurring flow structures. For instance, in figure 13, the region just above the barrel shock indicates high values of the incidence (dark red colour), because the Mach disc shows up at the same location in almost all of the 30 images processed, and is therefore largely stationary. On the other hand, a blue or purple colour indicates that flow features were not frequently encountered at that location. The lower value of the incidence in the plume upper-edge region shows that there is a large spatial variance in the jet shear layer, which multiple snapshots are needed to capture. This stresses the importance of the use of multiple snapshots for capturing the time-averaged location of the upper shear layer of the jet plume. Also, note the empirical similarity relation (3.1) for the time-averaged location of the upper shear layer of the jet plume, derived from wind-tunnel schlieren images in all iso-incidence plots.
Figures 13(a,b) are iso-incidence plots obtained from the experiments. In these plots, the Mach disc and the area around the recompression shock show up as dark red regions of high incidence. Further, the absence of clearly evenly coloured areas suggests that the flow features within the jet plume are not clearly correlated and originate from unsteady flow behaviour. This is even clearer in figures 13(c–e), which are iso-incidence plots obtained from the results of the DDES. Here, only the Mach disc and barrel shock are clearly identified flow features, whereas the rest of the jet plume suggests unsteady flow behaviour.
Comparing for
$J=2.8$
the experimental iso-incidence plot (figure 13
a) with its numerical counterpart (figure 13
c), it is observed that the same flow features, such as the bow shocks arising in front of the jets, the barrel shock and the Mach disc, are at the same position. The same is observed when comparing the experimental (figure 13
b) and numerical (figure 13
d) iso-incidence plots for
$J=3.8$
. Figure 13(e) shows an iso-incidence plot of the plume of the main jet for
$(S,J)=(7.17,4.8)$
. Here, it is clear that comparing with the plots for
$J=2.8$
(figures 13
a,c) and
$J=3.8$
(figures 13
b,d), the jet penetrates further into the cross-flow. The bow shocks arising in front of the jets have a larger part that is a normal shock, decelerating the cross-flow to a larger extent near the jet orifices, and therefore allowing the jets to penetrate deeper into the cross-flow. This, however, is accompanied by an increase in total-pressure losses induced by the shock waves.
The position of the jet upper shear layer, as calculated from the empirical similarity relation (3.1), has been added to figure 13. This empirical relation shows good agreement with the position of the upper side of the jet plume, validating (3.1) for these conditions of
$J=2.8$
and
$J=3.8$
, but also for a condition that was not investigated experimentally (
$J=4.8$
). The experimental (figures 13
a,b) and numerical (figures 13
c,d) results show a behaviour that is very similar and well-described by the empirical similarity relation.
To conclude, the experimental and numerical results show good agreement with respect to the obtained flow field, and (3.1) describes the jet penetration accurately for different
$(S,J)$
conditions.
3.4. Time-resolved results – numerical simulation
To obtain a more detailed insight into the penetration of the jets into the cross-flow, and the interaction of the jets with the cross-flow, this subsection discusses the time-resolved numerical results. The results for
$(S,J)=(7.17,3.8)$
will be discussed in detail, i.e. where the results on the refined grid have been obtained. To show the flow behaviour, flow solutions at five discrete times are shown, with equal time spacing 15
$\unicode{x03BC}$
s. This time spacing is chosen to be of the order of the convective time scale of the propagation of the flow features in the jet plume (
$D_2/U_c \approx 5$
$\unicode{x03BC}$
s). This allows us to follow the flow structures over time, which would not be possible for larger time spacings. Simultaneously, the time spacing is sufficiently large such that relevant flow evolution is visible. The discussion focuses on the results in the midplane (§ 3.4.1), the exit plane (§ 3.4.2) and the horizontal plane through the jet plume (§ 3.4.3).
3.4.1. Results in the
$(x,z)$
-plane
The DDES flow solution in the test section in the
$(x,z)$
midplane on the refined grid at five successive points in time: (a) numerical schlieren image (
$\partial \rho /\partial z$
); (b) vorticity magnitude (
$|\boldsymbol{\omega }|$
) in s
$^{-1}$
; (c) local Mach number; (d) total pressure
$p_0$
in Pa. Images are part of supplementary movies 2–5 are available at https://doi.org/10.1017/jfm.2026.11982. Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
The results of the DDES time-resolved numerical simulation of the flow in the wind-tunnel test section show an unsteady flow pattern containing features ranging over a large number of scales. For five discrete times, figure 14 presents the flow solution as a numerical schlieren image, and as colour plots of the magnitude of the vorticity, of the local Mach number and of the total pressure.
The numerical schlieren images are quite similar to the corresponding wind-tunnel schlieren images shown in figure 6 (with the red dashed square around the corresponding condition). In the snapshots in figure 6 and in figure 14, as well as in supplementary movies 1, 2, 9, 10 and 11, similar shock patterns are evident. Specifically, the two bow shocks associated with the two sonic jets emanating from the bottom wall of the test section, as also sketched in figure 2(b), are evident. The same holds for the reflections of these shocks at the top wall. However, in figure 14, considered as the numerical counterpart of figure 6, many more details of the interaction of the two jets with the cross-flow are evident, specifically due to the refined time resolution. Also revealing is the combination with the pictures of the distribution of the magnitude of the vorticity
$|\boldsymbol{\omega }|$
and of the local Mach number. Moving from upstream to downstream, one meets the following (see figure 15 for the locations of the features encountered).
Enumerated locations of the flow features in dual-jet injection into a supersonic cross-flow in the midplane (refined grid), corresponding to the main text: (a,c) numerical schlieren image (
$\partial \rho /\partial z$
); (b) distribution vorticity magnitude; (d) distribution local Mach number. The white dashed box in (c,d) indicates the region of the zoom of (a,b). Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
-
(i) The flow upstream of the first bow shock appears to be steady. Specifically, also in the bottom-wall boundary layer approaching the jet from the upstream orifice, the solution of DDES, i.e. the URANS equations solved in that part of the computational domain, is steady.
-
(ii) The separation shock is upstream of the first bow shock, due to the interaction of the subsonic part of the incoming thickening boundary layer and the pressure rise due to the bow shock. The separation shock and the bow shock are clearly part of a
$\lambda$
-shock that forms on top of the incoming boundary layer. -
(iii) Embedded within the subsonic-flow region on the downstream side of the normal part of the bow shock, one observes: (a) the tilted barrel shock with high-supersonic flow inside; (b) closed off by the Mach disc; and (c) the region with lower supersonic flow between the barrel shock and the upstream facing part of the mixing shear layer of the jet from the upstream orifice, which forms the contact surface (layer) between the fluid of the jet and that of the cross-flow. This results in a weak secondary bow shock that develops just downstream of the most upstream bow shock. The foot of the secondary bow shock is attached to the right-hand side of the Mach disc.
-
(iv) Downstream of the location of the Mach disc, the jet plume from the upstream jet is supersonic, with a narrow subsonic-flow region at the centre (a) originating at the Mach disc. The flow inside this jet turns unsteady (b), which is evident from successive snapshots of schlieren images in figure 14(a) and from supplementary movie 2. Actually, in the supplementary movies, it is observed that the trigger for unsteadiness appears to be the subsonic (cross-)flow around the Mach barrel becoming unsteady, similar to the subsonic flow around a circle-cylinder perpendicular to an oncoming flow.
-
(v) Further downstream, the jet (a) from the first bottom-wall orifice interacts with the larger, main, jet emanating, at sonic speed, from the second bottom-wall orifice (b). This main jet, tilted by the cross-flow, also features a high-supersonic flow inside its barrel shock (c), which is enclosed by a lower supersonic flow between the barrel shock and the mixing shear layer (d), forming the interface between the fluid from the second jet and the surrounding fluid. The barrel shock is closed-off by the Mach disc (e). In the second, main, jet plume, the flow is likewise the flow in the upstream jet plume, supersonic with a core (f) of subsonic flow of jet fluid that passed through the Mach disc.
-
(vi) The jet from the upstream orifice causes the bow shock in front of the downstream jet to lift off and form a
$\lambda$
-shock on top of the jet plume from the first jet. -
(vii) Downstream of the Mach disc of the second barrel shock, the plume of the jet from this orifice develops as an unsteady supersonic flow with patches of subsonic flow (a). As is clear from the numerical schlieren images, but especially from the plot with the distribution of the magnitude
$|\boldsymbol{\omega }|$
of the vorticity within the upper part of the plume, coherent vortical structures form in the jet plume. These vortical structures (b) travel with the flow, as observed in the time series in figure 14. -
(viii) In the emanating main jet, on the downstream side of its barrel shock (a), the fluid from the jet enters the region between the barrel shock and the mixing shear layer where jet fluid mixes with cross-flow fluid. Here, the flow, upwards along the barrel shock, is supersonic until at the position of the Mach disc. Here, a normal shock (b) slows down the fluid to subsonic speed in order to match up with the flow that passed through the Mach disc. Furthermore, the flow on the downstream side of the mixing layer, slightly supersonic, meets a normal shock (c) that slows down the flow to subsonic speed in order to match with the flow passing through normal shock (b).
-
(ix) In the numerical schlieren images in figure 14(a), immediately downstream of the second bow shock, more or less horizontal waves are distinguishable that travel to the top wall. These waves are caused by the unsteadiness in the first part of the upper shear layer of the main jet plume, i.e. the region where the jet plume from the first orifice collides with the upper shear layer of the main jet.
-
(x) The two primary bow shocks reflect at the upper wall of the test section through a Mach reflection. These show the same behaviour as in the experimental results of the flow (figure 6).
-
(xi) In the numerical schlieren images as well as in the plot with the distribution of the local Mach number, another pair of shocks is discernible, which reflect at the upper side of the main jet plume. These incoming shocks have been identified as reflections of the primary bow shocks at the side walls; see § 3.4.3.
-
(xii) The colour pictures of the vorticity magnitude
$|\boldsymbol{\omega }|$
show that vorticity is generated at the curved parts of the primary bow shocks, at the Mach discs, within the boundary layer along the bottom wall, as well as in the mixing shear layers, where the fluid from the jet hits the cross-flow fluid.
3.4.2. Results in the
$(y,z)$
-plane
Figure 16 presents, for five successive points in time, the time-resolved DDES solution in terms of the distribution of the
$x$
-component of the vorticity
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
(figure 16
a) and the distribution of the total pressure (figure 16
b), both at the exit plane of the computational domain. In figure 16(a), two types of vortical-flow regions are distinguished:
-
(i) vortical flow structures in the field, i.e. in the jet plume, that form during the deflection of the plume by the cross-flow;
-
(ii) vortical structures in the bottom-wall boundary layer with main components of
$\boldsymbol{\omega }$
in the
$y$
- and
$z$
-directions.
The DDES flow solution in the test section at the
$(y,z)$
exit plane at five successive points in time. (a) Vorticity
$x$
-component (
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
) in s
$^{-1}$
; red indicates positive (anticlockwise), and blue indicates negative (clockwise). (b) Total pressure
$p_0$
in Pa. Images are part of supplementary movies 6 and 7. Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
Concerning (i), in the field, the mixing-promoting dynamics of the flow is evident. Vortical structures with positive (red, anticlockwise) and negative (blue, clockwise)
$x$
-components of the vorticity (
$\boldsymbol{\omega }\boldsymbol{\cdot } \boldsymbol{e}_x$
) interact with each other. In addition, the whole collection of vortical structures appears to sway from left to right.
When the vortical structures are in close proximity, there is a strong mutual interaction: two close vortical structures of the same sign rotate around their centre of vorticity, while two close vortical structures of the same strength but of opposite sign translate in an upward direction, perpendicular to the line connecting their centres.
The collection of vortical structures seen in the exit plane can be represented by a pair of counter-rotating vortical structures, a vortex dipole. This corresponds to the mushroom-like cross-section of the jet plume, as indicated in figure 2(b) and as observed in the time-averaged DDES solution presented in § 3.2.
Concerning (ii), in the boundary layer along the bottom wall, the net effect of the presence of the field vortical structures results in an adverse pressure gradient in the direction of the inward-directed flow induced by the field vorticity at the bottom wall. The latter induced flow is directed from the side walls to the midplane. The adverse pressure gradient causes the boundary layer to separate, resulting in additional vortical structures of opposite sign, which interact mutually as well as with the vortical structures in the field.
Included in figure 16(b) are pictures of the distribution of the total pressure within the exit plane. The results confirm that the vortical structures and boundary layers are associated with losses in total pressure. Therefore, regions with total-pressure losses are often used to identify (strong) vortical structures in the form of free or wall-bounded shear layers.
3.4.3. Results in the
$(x,y)$
-plane
Figure 17 presents, for the same five successive points in time, for the time-resolved DDES solution, numerical schlieren snapshots on the plane
$z=3.14D_t$
, a cutting plane that is parallel to the bottom wall at a distance just above the two barrel shocks. The plane cuts through the bow shocks and their reflections at the side walls, and cuts lengthwise through the main jet plume. Most of the flow features identified in the midplane in figure 14 can be identified in figure 17, which provides useful insight into the structure of the flow field in the lateral direction, and therewith into the three-dimensional structure of the flow field.
The DDES flow solution in the test section in the
$(x,y)$
cutting plane at five successive points in time. Numerical schlieren images (
$\partial \rho /\partial z$
) on the cutting plane
$z/D_t=3.14$
. Images are part of supplementary movie 8. Here,
${M}_c = 1.6$
,
$S=7.17$
and
$J=3.8$
.
-
(i) The first, upstream, primary bow shock, as well as the second primary bow shock, is well-resolved. The lateral extent of the weaker secondary bow shock that forms just downstream of the first primary bow shock is clearly visualised.
-
(ii) The horizontal cutting plane intersects the second bow shock just below the merging point of the two legs of its
$\lambda$
-shock, evidenced by the blackish area immediately upstream of the second primary bow shock. -
(iii) At the side walls, the two primary bow shocks feature a Mach reflection with a Mach stem. The reflected shocks hit and interact with the unsteady flow inside the main jet plume. Subsequently, these shocks reflect and continue their path, back towards the side walls, where they reflect once more.
-
(iv) The flow topology is mostly stationary, specifically with respect to the shocks, except for the main jet plume, where the mixing takes place.
4. Conclusions
In this study, dual, tandem, under-expanded sonic jet injection, normal to the supersonic cross-flow of Mach 1.6, was studied. Experimental results of schlieren images in the midplane were complemented with results of numerical simulations, which offered time-averaged and time-dependent results in additional cutting planes. The numerical results matched the experimental results.
The schlieren images from the experimental study have been utilised to investigate the penetration depth of the jets into the cross-flow. This resulted in an empirical scaling relation of the time-averaged location of the upper shear layer of the jet, i.e. the penetration depth of the dual-jet injection. For
$J$
-dependent intermediate values of
$S$
, i.e. the dimensionless distance between the two orifices, the upstream jet imposes a shielding effect on the downstream jet. Compared to single-jet injection, this results in an enhanced penetration depth of the upper shear layer of the main jet plume.
The results of the numerical study of the flow in the wind-tunnel test section illustrate that the computational results of (steady) RANS and the time-averaged ones of time-resolved DDES match the results derived from the experiments, except for the jet plume mixing. For the jet plume mixing, the time-averaged DDES show better correspondence with the experiments than RANS. This stresses the importance of the use of time-resolved simulations for jet injection into a supersonic cross-flow. The numerical results, specifically the time-resolved ones from DDES, are clearly very powerful in enabling a thorough analysis of the complex flow field of dual-jet injection in a supersonic cross-flow. The numerical results reveal many details of the way the jet plume from the main orifice is affected by the shielding flow field induced by the jet emanating from the upstream orifice. In addition, the post-processing procedure developed for extracting iso-incidence plots from wind-tunnel schlieren images, from which the penetration depth of the main jet plume has been determined, has also found fruitful application in processing numerical schlieren images.
For future work, utilising an improved wind-tunnel schlieren system, and carrying out wind-tunnel experiments at further conditions in the
$(S,J)$
parameter plane and for other Mach numbers
${M}_c$
, will enhance the empirical similarity relation for the penetration of the jet into the cross-flow. Extending the similarity relation also to include numerical simulations with DDES and wall-modelled LES for additional conditions in the
$(S,J)$
parameter plane will reveal the flow structures responsible for the enhanced penetration for certain values of
$S$
. For obtaining enhanced jet penetration without excessive total-pressure losses, a measure for the efficiency of the mixing in relation to the loss in total pressure in the test section should be determined.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11982.
Acknowledgements
The authors acknowledge W. Lette, S. Wanrooij, H. Stobbe, E.M.P. Leusink, L. Tiemersma and F.B. Segerink for their contributions in the wind-tunnel experiments, and M.H.G. Wilens for his contribution in the numerical part of the study. TNO is thanked for their cooperation in the realisation of the present publication. The reviewers are acknowledged for their careful reading and constructive feedback.
Funding
This work used the Dutch national e-infrastructure with the support of the SURF Cooperative using grant no. EINF-14102.
Declaration of interests
The authors report no conflict of interest.
Mc=1.6
D1=1
D2=2
S=7.17
J=3.8
(x,z)
(x,y,z)=(81,45,48) mm=(36.3Dt,20.2Dt,21.5Dt)
y
S=7.17
μ
μ
(x,z)
∂ρ/∂z
|ω|
−1
Mc=1.6
S=7.17
J=3.8
J=3.8
S
(S,J=3.8)
Mc=1.6
D1=1
D2=2
(S=7.17,J=3.8)
(S,J)
zavg/Dt
zavg/Dt=2.5,3.0,3.5,…,8.0
Sopt(J)
μ
(x,z)
|ω|
−1
p0
Mc=1.6
S=7.17
J=3.8
μ
(y,z)
x
ω⋅ex
−1
p0
Mc=1.6
S=7.17
J=3.8
μ
(x,z)
J=2.8
J=3.8
J=4.8
∂ρ/∂z
|ω|
−1
M
p0
Mc=1.6
S=7.17
μ
(y,z)
J=2.8
J=3.8
J=4.8
x
ω⋅ex
−1
p0
Mc=1.6
S=7.17
μ
∂ρ/∂z
(x,y)
z/Dt=0.90
z/Dt=3.14
z/Dt=6.73
z/Dt=10.76
Mc=1.6
S=7.17
J=3.8
∂ρ/∂z
J
J=2.8
J=3.8
J=4.8
R2=0.70
R2=0.73
R2=0.72
R2=0.81
R2=0.79
Mc=1.6
S=7.17
(x,z)
∂ρ/∂z
|ω|
−1
p0
Mc=1.6
S=7.17
J=3.8
∂ρ/∂z
Mc=1.6
S=7.17
J=3.8
(y,z)
x
ω⋅ex
−1
p0
Mc=1.6
S=7.17
J=3.8
(x,y)
∂ρ/∂z
z/Dt=3.14
Mc=1.6
S=7.17
J=3.8
Open access
