2. Materials and Methods
The considered geometry consists of an axisymmetric convergent–divergent nozzle, which is illustrated in
Figure 2. Twelve cylindrical pipes, inclined
to the nozzle centre-axis, are arranged equidistantly on its circumference in the divergent section. The design Mach number at the nozzle exit is
, and the area ratio is
, which is defined as the nozzle exit area
to the throat cross-sectional area
. A total pressure source,
, drives the nozzle flow, which is four times the ambient pressure
. The employed total temperature
at the domain inlet is set slightly higher than the ambient temperature
to mimic the operating conditions of the experimental measurements at the University of Cincinnati [
21] and allow thereby direct comparison of the results. Downstream of the nozzle exit, the stream expands into a quiescent ambient environment. The applied boundary conditions are summarised in
Table 1.
Three injection locations
, i.e.,
,
, and
upstream from the nozzle exit, are considered to investigate the sensitivity of the shock pattern on the injection location found beneficial by Semlitsch & Mihăescu [
16]. The injector diameter is
of
mm. The injector total temperature is set to the ambient temperature, and the flow is aligned with the pipe axis. The injector total pressure is an investigation parameter referred to as Injection Pressure Ratio (IPR), relating the applied total pressure at the injector inlet to the ambient pressure.
2.1. Numerical Simulation Procedure
The three-dimensional nozzle flow is simulated by solving the Navier–Stokes equations for compressible fluids numerically, which can be written as,
where
t is the time,
the spatial coordinates,
is the fluid density,
u the velocity,
p the static pressure defined by the ideal gas law,
, and
the viscous stress tensor. For a Newtonian media, the viscous stress tensor can be written as,
where
is the Kronecker delta,
is the dynamic viscosity, and
is the strain rate tensor. The temperature dependence of the dynamic viscosity,
, is modelled using Sutherland’s formula,
where
is the reference dynamic viscosity,
is the reference temperature, and
is the Sutherland temperature. The values set for the constants are;
kg/s m,
273.15
, and
110.4
. The strain rate tensor can be defined as,
The conservation of mass and energy is guaranteed, by solving the conservation equations
and
respectively, where
is the total internal energy and
is the heat flux. The total internal energy can be linked to the other primary variables by,
The heat flux is calculated by employing Fourier’s law,
where
K is the heat conductivity.
An explicit low-storage four-stage Runge–Kutta scheme of second-order accuracy was utilised for temporal integration. A constant time-step
was set for time advancement to satisfy the Courant–Friedrichs–Lewy condition (see
Table 2). The inviscid flux is computed for finite volumes using a second-order central difference scheme, where a Jameson-type artificial dissipation [
22] is added to avoid spurious numerical oscillations near sharp gradients or discontinuities, such as shock waves. The added dissipation is based on a blend of second and fourth-order differences triggered by a pressure difference sensor. The viscous stresses
in the governing equations can be expressed by a vector Laplacian of the flow-field. The identity of the Laplacian equation is used to calculate the second-order gradient terms. The remaining terms for a fully viscous approach are obtained using a Green–Gauss formulation. Thereby, a second-order accurate discretisation for the viscous terms is achieved. A more detailed description of the utilised solver,
edge, has been provided by Eliasson [
23].
The smallest flow scales, i.e., shock wave thickness or smallest turbulent length scales, are lower than a mesh resolution with efficient use of computational resources could handle. Thus, the governing equations are spatially filtered, and the Favre-averaging of the equation set is performed. The arising subgrid-scale terms account for unresolved physics. The nature of the smallest flow scales is to dissipate the kinetic energy into heat at the molecular level. This effect needs to be endorsed employing some model since insufficient dissipation would promote an unphysical energy build up [
24]. Smagorinsky [
25] introduced an artificial viscosity model, similar to the eddy viscosity concept, acting as a subgrid-scale model to stabilise the numerical approach. Calibration coefficients are required, which can affect the flow-field substantially, e.g., shown by Uzun et al. [
26], and need to be obtained empirically. Despite the recent development of novel subgrid-scale models [
27,
28,
29], the compressible nature of nozzle flow challenges the choice of an adequate explicit subgrid-scale model. At a distance from solid walls, the dissipation is dominated by the kinetic energy transfer from the large to the small flow scales. This behaviour is statistically universal at sufficiently high Reynolds number accordingly to Kolmogorov’s hypothesis, which is promising for modelling. In the present approach, the subgrid scales are represented implicitly by the intrinsic numerical dissipation of the solver. The approach is also known as monotonically integrated Large Eddy Simulation or implicit LES. For a non-oscillatory finite volume discretisation of at least second-order accuracy, the numerical truncation error can be interpreted as a Clark-type subgrid-scale model for the momentum equation and as a Smagorinsky-type subgrid-scale model for the energy conservation Equation [
30]. Hence, the kinetic energy fluctuations are absolutely decreasing with this approach and no spurious energy build-up at the high frequencies occurs, which is proven by the spectra shown in
Figure 3a. The entire discretisation resolution of the mesh grid is exploited, and no information is truncated by an explicit filter [
31,
32]. Further, no calibration constants are required, which might be difficult to obtain for internal supersonic nozzle flows.
2.2. Meshes, Verification and Validation
The numerical grids consist of 259 hexahedral blocks arranged in a centred O-structure, which is shown in
Figure 4. The discretisations of all injectors are structured with an individual O-grid to enable consistent near-wall refinement. The mesh in the investigation section is retained to be as uniformly fine spaced as possible. The inlet and outlet sections assist to buffer undesired reflections by gradual cell size increments the boundary conditions. The buffer region downstream of the nozzle exit covers an expansion zone into stagnant ambient conditions.
A mesh resolution study was performed using three grid levels while conserving the block structure. The specific properties of the meshes are listed in
Table 2.
Figure 3a illustrates the velocity spectra obtained from static pressure signal recorded in the probe
, which is located in the wake of the injection as indicated in
Figure 4b. The spectra reveal the same broadband content at low frequencies (consistent with the reported results by Kawai & Lele [
12], and Genin & Menon [
11]), while the captured turbulent cascade is prolonged with the mesh refinement.
The methodology suggested by Celik et al. [
33] has been chosen to analyse the effect of grid refinement quantitatively, which is based on the Richardson extrapolation. The numerical solution obtained on three grid levels is extrapolated to an infinite fine grid, and the resulting error is assessed. For an injection pressure ratio of
at the intermediate jet injection location (i.e., on the divergent slope of the nozzle at
upstream from the nozzle exit), the solutions of the estimated static pressure on the nozzle centreline with three grids and the Richardson extrapolation are shown in
Figure 3b. The solutions show nearly overlapping contours in the smooth flow-field regions and a shift (in the order of a mesh cell) of the shock wave in the upstream direction with increased mesh resolution.
Table 3 shows that the mean relative error between intermediate and fine grid is small, although the maximal relative error is high. Shock waves represent steep discontinuous gradients, which induce even for a small shift a large maximal relative error. The relative errors between the numerical solution and the extrapolation are even lower (as shown in
Table 3), which suggests that the grids are properly designed to capture shock pattern changes with fluidic injection accurately.
The experimental visualisation of confined supersonic flow is challenging. Instead, the numerical simulation data are compared to measurements of the external shock pattern performed at the University of Cincinnati [
34].
Figure 5 reveals that the numerical approach captures (on a coarser mesh resolution than the intermediate mesh) accurately the shock wave angles and locations, which are the focus of the present investigation. More quantitative validation with particle image velocimetry data has been documented by Semlitsch et al. [
4] for external fluidic injection, where a good agreement regarding the streamwise velocity magnitudes and shock wave angles and locations are reported.
3. Results
Figure 6 illustrates the shock pattern transformation in the divergent nozzle section for different injection pressures and locations. We define the shock pattern occurring without injection as the baseline, where an oblique shock anchors slightly downstream of the nozzle throat and merges in the core in the form of a Mach-disk. There, shock wave reflections and slip-lines arise. This shock structure can be split into two separate employing fluidics. Bow shocks establish upstream of the injection governing the shock structure. A second shock structure establishes on top of the injected flow. The upstream shock structure steepens with increasing injection pressure, and the downstream shock structure shifts the nozzle exit.
The oblique shock angle,
, is governed by the pressure ratio across it, the upstream Mach number,
, and the ratio of specific heats,
,
where the indices 1 and 2 refer to the quantities to the upstream and downstream locations, respectively.
Figure 6 shows that the static pressure distributions (without injection) differ on the nozzle circumference for the individual injection locations. As a consequence, the crossflow momentum is there different, which pushes against the injected stream. The force balance of the injected flow momentum blocking the crossflow induces locally high static pressure regions at the injector orifice (for all
). The mass flow rate and the effective velocity ratios,
, develop differently for individual injector configurations with equal IPRs accordingly to the static pressure at the injector orifice. The pressure ratio over the bow shock increases with increasing IPR. This also amplifies the static pressure downstream of the bow shock continuation, and the angle of the first upstream shock structure steepens.
The flow accelerates in the divergent nozzle section until hitting the first upstream shock. Hence, the further downstream the first upstream shock is situated, the higher the upstream Mach number, , and the lower the upstream static pressure, , of the shock. The pressure downstream of the first upstream shock, , is governed by the injection. The pressure ratio is too high for an oblique shock, and a Mach disk establishes for the downstream injection location. On the contrary, the flow did not accelerate enough to reach low enough pressures to support a normal shock wave for the upstream injection location.
The fluidic injection contracts the core flow. Thereby, the penetration depth governs the location of the second downstream shock structure. For low injection pressure ratios, the injected jets remain attached to the nozzle walls. The kinetic energy fluctuations indicate that once the injection separates from the nozzle walls, the penetration depth into the crossflow scales with the effective velocity ratio, , until the injectors choke (). The radial penetration ceases beyond choke.
The evolution of the injection shock pattern with increased
is sketched in
Figure 7. The injection remains subsonic for
below
. A leeward shock arises for
separating the injected jet and the flow recirculation downstream of it. Mach disks appear spontaneously on top of the leeward shock. The crossflow shears the windward side of the injection, and windward shocks do not emerge until
reaches
. These windward shocks remain weak and exhibit intense motion. Genin & Menon [
11], and Kawai & Lele [
12] simulated higher effective velocity ratios and observed that this motion of the shocks in the injection causes vortex breakdown or bursting.
Figure 7 shows the counterrotating vortex pairs in terms of streamwise vorticity iso-surfaces, which reveal a lower penetration depth in radial direction for injection with shocks (
) than injection with sporadic shocks (
). The counterrotating vortex pair results relatively weak for
, but an additional counterrotating vortex pair establishes in between the injection (highlighted in white in
Figure 7). In proximity to the nozzle wall, the nozzle flow is diverted and partially forced in between the injection. Thereby, the incoming stream is contracted and pushed partially away from the walls. Such an additional counterrotating vortex pair can also be observed for other injector operating conditions but evolves less clearly due to the enhanced unsteadiness.
Several compression waves form the second downstream shock structure (see
Figure 6). Illustrating the time-averaged pressure fluctuations as transparent iso-surface on top of the counterrotating vortex pair in
Figure 8 demonstrates that the compression waves anchor on the injected stream. Especially for higher injection pressure ratios, these compression waves display strong motion as the kinetic energy fluctuations show in
Figure 6, which is related to the injection unsteadiness. For intermediate injection pressure ratios (3.0–4.4) at the upstream and intermediate injection location, compression waves cannot always focus as a second downstream shock pattern (the transient behaviour of the compression waves can be observed in the animation provided as
Supplementary Material.). It is noteworthy that the first upstream shock structure is not reflected (at all times), and therefore, the exhaust remains without persistent shock waves.
Fluidics can be potentially used to compensate for the under-expansion of an exhausting jet and match the static with the ambient pressure at the nozzle exit. The profiles of the time-averaged static pressures shown in
Figure 9 exhibit significant radial variations at the nozzle exit, which can be attributed to the shock waves and injected wakes penetrating this plane. Without fluidic injection, the static pressure is high in the interior flow while being significantly lower at the outer circumference.
Figure 9 demonstrates that injection augments the static pressure successively at the circumference, where the injection location affects its impact. Significant differences between the wake and in between wake profiles can only be noted for the most downstream injection location. Hence, the injected flow is not the principal contributor to the static pressure rise at the circumference. Moreover, the downstream shift of the second, downstream shock pattern towards the nozzle exit is the cause. The best pressure match (uniform around the circumference) is achieved for injection pressure ratios of approximately
for the intermediate injection location.
Shock, mixing and blockage losses reduce the nozzle efficiency to generate thrust. The mean and variance of the thrust (samples over simulation time) are illustrated in
Figure 10. Even low injection amounts increase the thrust, while a little additional benefit is achieved with higher injection pressures. This is not unexpected because the nozzle is operating an underexpanded exhaust, as shown in
Figure 5.
Figure 6 shows that the Mach-disk is most extensive for the downstream injection configuration. Although the most considerable shock associated losses may be expected for this configuration (neglecting bow shock losses), the thrust is the highest, as revealed by
Figure 10. (the cross-sectional area of the Mach-disk remains small compared to the total cross-section.) The high thrust results are due to the high flow velocities in the core flow and beneficial static pressures at the circumference for this configuration. Thus, shifting the shock pattern is more beneficial for thrust optimisation than shock strength reduction.
The thrust variance can be reduced using fluidics with low injection pressure ratios, i.e., . The injection introduces unsteadiness in the divergent nozzle section for higher injection pressures. The generated thrust unsteadiness arises from induced shock pattern motion rather than the flow unsteadiness in the injection wakes.
4. Conclusions
The ability to transform shock patterns in convergent–divergent nozzles with fluidic injection has been investigated using Large Eddy Simulations. The numerical approach has been validated with Shadowgraph visualisations, and the uncertainty due to the grid resolution has been assessed. We showed that fluidics could transform a single shock pattern manifesting in a converged-divergent nozzle into two weaker shock structures. The bow shocks (upstream of the injection) form the first shock structure, and a second shock structure establishes on top of the penetrating injection. The first upstream shock pattern steepens with increased injection amount, while the second shock pattern location shifts downstream. The sensitivity of the shock pattern to injection pressure changes is linked to the injection location. The surface static pressure is lowest towards the nozzle throat. Consequentially, higher mass flow rates arise for injection close to the nozzle throat, which causes deeper injection penetration. Once the injected jets choke, the shock waves compensate partially for the pressure difference.
It was shown that even low injection pressures augment the thrust, compensating for the underexpanded nozzle operating condition. The most favourable injection configuration was found when a (second downstream) shock pattern was shifted towards the nozzle exit plane. Although the shock associated losses are high for this configuration, the high core-flow velocities in combination with high static pressures at the nozzle circumference generate the highest thrust.
The closely spaced fluidic injection can provoke additional counter-rotating vortex pairs in between injectors. The generated vorticity forms the source of the additional established counter-rotating vortex pairs within the injection ports for some of the injection pressure ratios considered. The formation of additional counter-rotating vortex pairs was linked with compression waves forming in between the upstream shock wave and the consecutive shock wave reflection.
We have shown that the shock pattern could be transformed into a transient compression wave pattern for some injection configurations. The compression waves exhibited intense motion and did focus only at times to a shock pattern. The fluidic injection can also be used to match the static pressure at the nozzle exit plane for over-expanded exhausts to evade shock formation at the nozzle exit. We find that the injection pressure is not of primary importance to achieve pressure distributions favouring ideal expanded operating conditions at the nozzle exit. The shift of the shock pattern inside of the nozzle is the critical parameter. Hence, shock associated noise can be attenuated most efficiently using internal fluidic injection, as demonstrated by Morris et al. [
17] and Cuppoletti et al. [
21].