3.1.1. Shock Tube Problem.
Shock tube is a well-known benchmark problem (also known as Sod shock tube problem [
52]) for the validation of compressible codes. It is an attractive test case for model validation due to the simple one-dimensional geometry of computational area and the possibility to simulate a wide range of flow physics, including shock and rarefaction waves, contact discontinuities, and complex wave interactions. The available analytical and computational solutions in the literature allow for the assessment of the code's accuracy and reliability. The shock tube problem for two cases have been solved: (1) a shock tube with an inert multicomponent gas mixture, and (2) a shock tube with a chemically reactive multicomponent gas mixture. In these validations, the initial discontinuity of the gas in the shock tube is decayed into different gas configurations such as shock waves, contact discontinuous, and rarefaction waves spreading into different directions (See
Figure 1).
The first validation case is a multicomponent inert shock tube simulation without the chemical source term in order to avoid shock-induced autoignition of the gas mixture. A homogeneous premixed gas mixture of 20%H
2/10%O
2/70%Ar (in mole%) is used in the simulation to verify the numerical accuracy of the model implementations for thermal and mass transport of multicomponent mixture. The shock tube is 0.10 m in length, and it is discretized with 400 uniform elements. The time step is equal to
. In the middle of the tube is a fixed diaphragm that separates the gas between the left side (
and
) and the right side (
and
). At time
the diaphragm is eliminated, and the shock wave starts to propagate from right to left and expansion wave moves in the opposite direction. The initial configuration of the gas mixture generates the propagation of shock and discontinuous waves to the left (low-pressure) and a rarefaction (expansion) wave to the right side (high-pressure) of the tube.
Figure 2 shows the current simulation results (solid line) compared with the simulation results of Huang et al. [
53]. The modeling results have good agreement. It is noteworthy that the
compressibleCentralReactingFoam solver's accuracy is 2nd order for both time and spatial convective terms.
A reactive multicomponent gas mixture in the shock tube is also modeled using the conditions of Martínez-Ferrer et al. [
34] for comparison. The molar composition of the gas mixture is 20% H
2/10% O
2/70% Ar. The length of the tube is equal to
with the initial conditions corresponding to the left and right sides as follows:
and
, where
are measured in (
), (
), and (
), respectively. The length of the tube is discretized by 400 uniform elements with 0.01
timestep. The simulation used the reaction mechanism of Conaire et al. [
54] for hydrogen which consists of nine (9) chemical species (
H2, O2, H, O, OH, HO2, H2O2, H2O, ), and 18 elementary reactions. A wall boundary condition is implemented on the left side of the tube, and on the right side, a non-reflected boundary condition is implemented. The initial configuration of non-stable discontinuity is decayed into a shock propagated to the tube's left side. The solutions of spreading discontinuities in the shock tube for velocity, mass fraction of
, and temperature:
,
, and
are shown in
Figure 3 for 230
simulation time. Comparisons of current solution (solid lines) in
Figure 3 with computational data of Martínez-Ferrer et al. [
34] (symbols) show good agreement. The visible discrepancies are seen in the vicinities of local maximums around x = 0.06 m at 230
that is explained of using by Martínez-Ferrer et al. [
34] a seventh-order accurate Weighted Essentially Non-Oscillatory (WENO) scheme to discretize the non-linear advective terms.
The decay process of the initial multicomponent gas mixture, and how the shock wave propagate in the tube was investigated in detail by comparing the profiles of inert and reactive mixture of H
2/O
2/Ar.
Figure 4 shows the simulation results for the temperature profiles of inert and reactive mixture at various simulation times.
The shock starts to propagate from its initial position in the middle of the tube from right to left toward the wall (
i.e., x = 0) with a velocity
= 8100 m/s and reaches the wall around 75
. After reflection from the wall, the shock moves in the opposite direction from left to right. The reflected shock moves with a velocity of 4500 m/s that is almost half of
as the gas flows right to left with a velocity of 490 m/s. Due to the interaction of the shock with the wall, the temperature of the reflected shock is increased from 750 K to 1200 K for both inert and reactive cases. As seen in
Figure 4b, the shock-induced autoignition starts to occur (
i.e., a slight increase in temperature at x = 0) at 120
for the reactive case as the simulation time is longer than the chemical induction time to generate combustion radicals. Due to the combustion heat release, the gas temperature increases to 2100 K as the simulation time is increased to 150
. The gas temperature approaches equilibrium state as the simulation time is further increased to 220
.
3.1.2. Simulation of Ladenburg Jet Problem.
Simulation of a well-known experiment of underexpanded jet problem by Ladenburg et al. [
55] is conducted to verify code accuracy for the shock formation of supersonic jets. Two-dimensional axisymmetric simulation of the problem using
rhoCentralFoam was reported previously [
41,
56]. In the current work, 3D simulation of the Ladenburg experiment is performed with the
compressibleCentralReactingFoam solver. As the maximum temperature of the experiment is no more than 300 K, the chemical reactions will be at frozen state and hence, the gas mixture is considered be inert. The computational 3D mesh (see
Figure 5) has 5 mm inlet radius, 10 mm free surface radius, and 30 mm length. The following inlet conditions for pressure, axial velocity and temperature are used to ensure a sonic condition with Mach number of one (1) at the inlet surface:
,
and
. The conditions at the free stream surface are:
,
and
. Non-reflected boundary conditions are used at the outlet. The simulations are performed using a mesh with total hexahedral elements equal to 1.72 million. The sketch showing wave structures and the simulation flow of Mach number contours is shown in
Figure 6a. The pressure at the inlet of the computational region is the pressure of the reservoir, which is greater than the ambient pressure; thus, the flow is defined as underexpanded. When air accelerates from high inlet pressure to the atmosphere, the gas flow declines due to its expansion and acceleration. The gas jet expands to the atmospheric pressure through an expansion fan (
Figure 6a). It is seen in
Figure 6b in the free jet boundary above which the gas flow is at rest (i.e., M=0). In
Figure 6a, numbers 1-5 in the circles indicate different gas structures, such as expansion and shock waves. The expansion waves (denoted by 1) are reflected from the free boundary as compression waves. The accumulated compression waves at a point then form an oblique shock wave, shown by the solid line (denoted by 2). The formed oblique shock waves are again reflected off the center line of the flow. This reflection causes the formation of a Mach reflection before the intersection point with the centerline, producing a triple shock point where the compression waves, oblique shock, and centerline coincide. The incoming oblique shock wave (denoted by 2) is the first Mach reflection. The second shock (denoted by 3) is a Mach reflection, perpendicular to the centerline termed a normal (Mach disk) shock. The third Mach reflection (denoted by 4) is another oblique shock wave that is reflected off the constant pressure-free air-jet boundary as an expansion fan (denoted by 5). The escalation of the pressure and temperature occurs as the flow passes through the normal shock or Mach disk. So, the Mach disk is the high-pressure region in an underexpanded air jet that is formed through shock waves and expansion waves created by the vast difference between the inlet pressure and the ambient pressure.
Figure 7 shows the comparison of the current simulation with experimental data of Ladenburg et al. [
55] for density iso-lines. The simulation results have a good agreement with the experimental data, especially the location and height of the Mach disk. The experimental results show that a weak shock is produced due to the expansion of air from the nozzle orifice (close to density
) and then extends toward the nozzle axis. Thus, it creates a Mach disk feature and a triple point (i.e., the intersection of Mach disk, barrel shock and reflected shock). The experimental data in
Figure 7 (i.e., lower panel) shows that the position of the triple point is at 1.7 mm radial distance and 13.3 mm axial distance. The computed radial and axial distance of the triple point are 1.65 mm and 13.5 mm, respectively. Thus, the computational results yield less than 3% error compared to the experimental data. In summary, the computational results are in good agreement with the experimental data as well as the physics of the shock formation, expansion wave behavior, and their interactions with jet boundary and slip lines.
3.1.3. Simulation of DLR Scramjet Combustor
The combustion experiment of the Waidmann et al. [
57] is selected to validate the
compressibleCentralReactingFoam solver to model the turbulence combustion interactions in a scramjet combustor. The experimental data were obtained for hydrogen combustion using the Scramjet test facility at the German Aerospace Center (DLR) at Mach number two. The experimental system was modeled using 2D and 3D computational studies by Oevermann [
58] and Haung et al. [
35], respectively. Oevermann [
58] used RANS (i.e., k-ε turbulence model) with laminar flamelet whereas Huang et al. [
35] used LES with PaSR sub-grid scale combustion model to simulate the Scramjet. Both computational studies produced reasonable agreement with the experimental data. In the current work, 2D LES simulation is conducted with PaSR sub-grid scale combustion model to validate
compressibleCentralReactingFoam solver. Although the two-dimensional assumption may introduce some uncertainty when directly comparing it to the experiment, the 2D computational study will help to verify the accuracy of the model implementation in
compressibleCentralReactingFoam qualitatively.
Figure 8 shows the geometry of DLR test rig. The combustor has a height of
and width of
at air inflow. The divergence angle of the upper wall is
to compensate for the expansion of the boundary layer. The following boundary conditions for air are established to ensure the Mach number
:
,
,
with the mass rate of airflow
. The injector of fuel (
) is a slot in the center of the wedge-shaped strut with a height
that gives the surface area of injection
. The following boundary conditions are used to establish fuel injector Mach number
:
,
,
with the hydrogen mass flow
. The no-slip boundary conditions are implemented on the solid walls. The computational region is discretized on 100,000 hexahedrons (See
Figure 9). A detailed kinetic mechanism for hydrogen combustion with 9 species and 18 elementary reactions [
59] is used.
Initially, the cold flow is simulated.
Figure 10 shows a Schlieren photograph of the channel flow and computational data of density gradient contours corresponding to the experiment. Comparisons of these two pictures show that the current compressible solver is able to accurately capture the aerodynamics features such as shock waves, expansion waves, interaction, and reflection of these structures from the walls and jet-shear layer. The computational results shown in
Figure 10b also identify various flow characteristics with labels 1 to 7. Label (1) shows the leading shocks formed on the tip of the edge (see
Figure 10). Label (2) refers to the expansion waves due to the supersonic flow past the corners of the edge. Label (3) is the reflected shock of (1) from the upper wall whereas label (4) refers to the shocks formed due to collision of gas flow passing through the corner of the edge. Label (5) is an oblique shock formed from merging of the reflected shock and shocks in (4). Lebel (6) is a reflected from the upper wall oblique shock whereas label (7) is the jet shear-layer formed from the Kelvin-Helmholtz instability. The oblique shocks that are reflected from the upper and lower walls interact with the jet shear layer to get reflected again without penetration as shown
Figure 10b.
Figure 11a shows the normalized contours of density gradient
|.
Figure 11b shows the Mach numbers of the cold flow. Shocks are generated from the tip of the wedge and are further reflected from the lower and upper walls of the combustor. Also, the expansion wave on the upper wall at the starting divergence is seen. Gas flow passing the corners of the wedge are subjected to expansion by forming two expansion waves. The jet vortex structure arises from the Kelvin-Helmholtz instability of the jet-shear layer [
14]. It is noteworthy that the simulation with the
turbulent model did not generate shear-layer instabilities as reported by Oevermann [
58]. This finding confirms that
models are limited when predicting such phenomena and
models are more appropriate in scramjet combustor simulations.
Figure 12 compares the cold flow simulation results for pressure distributions with the experimental measurements of Waidmann et al. [
57].
Figure 12a show the distribution of pressure on the lower (solid line) and upper (dashed line) walls along with experimental data obtained on the lower wall (symbols).
Figure 12b compares the pressure along the center line with the experimental data. The model predictions are able to capture the pressure profiles observed experimentally.
Simulations were performed for reactive flow with detailed hydrogen kinetic mechanism for experimental conditions of Waidmann et al. [
57]. Combustion was initiated with the help of an igniter during the test as no autoignition was observed experimentally. Our initial simulation also showed that combustion did not occur without an external ignition source. Hence, ignition zones at the two corners of the strut with temperature
are implemented.
Figure 13 compares Schlieren's photograph of the reactive flow with the computational data for the density contours obtained in the current work for the experimental condition of Waidmann et al
. [
57]. The comparison shows similar gas dynamic features such as shock waves, expansion waves, interaction, and reflection of the flow structures from the walls and jet-shear layer. Significant gas dynamic differences can be seen between cold flow (
Figure 10) and hot flow (
Figure 13). For example, expansion waves (i.e., labeled 2 in
Figure 10) at the corners of the wedge in cold flow are replaced with the shocks generated due to the combustion in
Figure 13. It should be noted that a considerable subsonic region is developed behind the hydrogen injection for reactive flow.
A discussion on the influence of different turbulent models on supersonic combustion is presented in this section. Scramjet simulations for the DLR Scramjet experiment [
57] are performed with the following turbulent models described in Subsection 2.1: SMG, LDkEqn, and WALE.
Figure 14 compares the experimental data for temperature at three cross sections along the combustor with the simulation results obtained with three turbulent models. The figure also shows the location of the temperature measurements in the DLR Scramjet combustor. The simulation results close to the strut (
i.e.,
x = 62 mm) show that LDkEqn and WALE models have similar profiles. Further downstream (
i.e.,
x = 216 mm), all three models have similar profiles. Overall, WALE model has slightly better agreement with the experimental data.
Figure 15 compares the current WALE model simulation results with other DLR Scramjet simulation data reported in the literature [
12,
30,
58,
60]. Modeling result of Berglund and Fureby, [
60] closely follow the experimental profile except at
x = 62 mm. Overall, the present model has better agreements with the experimental data at all locations, but it does show relatively higher diffusivity at downstream locations (
i.e.,
x = 108 mm and 216 mm).
Based on the analysis of different turbulent models, the WALE model is selected as the most suitable for Scramjet simulation. The effect of turbulence-combustion coupling on gas dynamics and flame structure is discussed in this section using the DLR Scramjet simulation results obtained with WALE model.
Figure 16 shows the Mach number (Ma), rate of heat release per unit volume (Qdot), instantaneous temperature (T), and normalized contour lines of densities’ gradient (
|) obtained from the simulation. The shock structure, shown in
Figure 16a,d, is characterized by strong interaction with the shear layers characterized by large density gradients. Behind the strut, a central recirculation bubble is generated with a subsonic region. The formation of the recirculation zone is very important for flame holding and, hence, stabilizing the flame in the scramjet combustor. The rate of heat release contours in
Figure 16b show that combustion is intensive within the thin shear layers. The combustion is initiated within the shear layers separated by the recirculation zone close to the edges of the strut. This observation is supported by the temperature data presented in
Figure 14a.
The direction of diffusive spreading of fuel and oxidizer helps to analyze turbulence-chemistry interactions. It can be estimated by the interaction between the gradients of fuel (
) and oxygen (
) mass fractions. Takeno and co-workers [
61] suggested the use of a flame index (known as the Takeno Flame Index - TFI) based on scalar productions of the gradients of the fuel and oxidizer mass fraction, which help to distinguish between the premixed and non-premixed zones [
61]. The TFI is defined as
. Positive and negative values of TFI indicate non-premixed or premixed modes, respectively [
61].
Figure 17 shows the TFI calculated for the DRL Scramjet conditions. Pockets of non-premixed zones shown in red are surrounded by large premixed zone shown in dark blue.
Figure 18 shows the species mass fraction contours for
,
,
and
. It is observed [
9] that the OH field is far from uniform in the mixture due to turbulence, which would have been expected from fast chemistry arguments. High OH regions are found where the mixture fraction iso-surfaces are highly convoluted, and low values are in the areas where the mixture fraction iso-surfaces are stretched and not wrinkled. The mass fraction of OH is relatively small at the beginning, while it is most prominent in the downstream regions, especially in the transition zone. The OH radical represents the existence of a flame front which correlates with heat release rate. Field of H2O correlates with the field of instantaneous temperature while the field of HO2 correlates with OH and Qdot on the shear layer.
Figure 19 shows mixture fraction (
) and scalar dissipation rate (
. The mixture fraction measures the local fuel/oxidizer ratio, and calculated as
where
is equivalence ratio. The scalar dissipation rate is calculated as
, where
is the mixture diffusion coefficient. The scalar dissipation rate is essentially the rate of mixing between fuel and oxidizer [
9] and
represents the local mixing rate at the molecular level [
62]
. The scalar dissipation rate
is critical for modeling the effect of turbulence on reaction rates as it has significant influence on non-premixed combustion. It often provides the connection between the molecular mixing and the combustion. In turbulent flows, the scalar dissipation is seen as a scalar energy dissipation. Its role is to destroy (dissipate) scalar variance (scalar energy) analogous to the dissipation of the turbulence. Unlike the kinetic energy dissipation, most of the scalar dissipation occurs at the finest scales. The scalar dissipation has higher values near the strut and decreases further downstream as shown in
Figure 19b.
A time-scale analysis based on the estimation of the local Damköhler Number (Da) is also calculated to examine the reactive zones in more detail. The Da number, quantified as the ratio between the turbulent mixing time scale
and the chemical time scale
:
[
63,
64]. The Borghi diagram [
63] shown in
Figure 20a helps to recognize different regimes of turbulent combustion.
Figure 20b,c show the Damköhler number (Da) on a logarithmic scale and the fraction of fine structure (
), respectively. The bright white lines indicate the locations where
. The Damköhler number is relatively high (
) in the core of the jet flow. This indicates that the time scale of chemical reaction is smaller than the turbulent mixing time scale which means combustion regime with fast chemistry. Thus, the fraction of reactive cells (
) (see Section 0 for details) is mostly small in this zone. In contrast, in the shear layer where the rate of the heat release (see
Figure 16b) and the scalar dissipation rate (see
Figure 19b) have the highest values, the flame is controlled by finite rate chemistry with
and
In the core flow, it is also seen small pockets where
indicating the combustion zone is in the “island formation” regime as shown in
Figure 20a. These zones are characterized by the combustion regime with a finite rate chemistry with a rather high fraction of reactive cells (
) and are called partially-stirred reactor.