Next Article in Journal
NGD-Transformer: Navigation Geodesic Distance Positional Encoding with Self-Attention Pooling for Graph Transformer on 3D Triangle Mesh
Next Article in Special Issue
A New Symbolic Time Series Analysis Method Based on Time-to-Space Mapping, through a Symmetric Magnetic Field, Quantized by Prime Numbers
Previous Article in Journal
Diamond-α Hardy-Type Inequalities on Time Scales
Previous Article in Special Issue
Critical Dynamics in Stratospheric Potential Energy Variations Prior to Significant (M > 6.7) Earthquakes
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

OpenFOAMTM Simulation of the Shock Wave Reflection in Unsteady Flow

by
Lucas Monaldi
1,
Luis Gutiérrez Marcantoni
2 and
Sergio Elaskar
3,*
1
Instituto de Estudios Avanzados en Ingeniería y Tecnología (IDIT), Universidad Nacional de Córdoba and CONICET, Córdoba 5000, Argentina
2
Facultad de Ingeniería y Ciencias Básicas, Funfación Universitaria los Libertadores, Bogotá 111221, Colombia
3
Departamento de Aeronáutica, Instituto de Estudios Avanzados en Ingeniería y Tecnología (IDIT), FCEFyN, Universidad Nacional de Córdoba and CONICET, Córdoba 5000, Argentina
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(10), 2048; https://doi.org/10.3390/sym14102048
Submission received: 29 June 2022 / Revised: 13 September 2022 / Accepted: 22 September 2022 / Published: 1 October 2022
(This article belongs to the Special Issue Symmetry in Nonlinear Dynamics and Chaos)

Abstract

:
This work studies the impact of a shock wave traveling with non-constant velocity over straight surfaces, generating an unsteady and complex reflection process. Two types of shock waves generated by sudden energy released are studied: cylindrical and spherical. Several numerical tests were developed considering different distances between the shock wave origin and the reflecting surface. The Kurganov, Noelle, and Petrova (KNP) scheme implemented in the rhoCentralFoam solver of the OpenFOAMTM software is used to reproduce the different shock wave reflections and their transitions. The numerical simulations of the reflected angle, Mach number of the shock wave, and position of the triple point are compared with pseudo-steady theory numerical and experimental studies. The numerical results show good accuracy for the reflected angle and minor differences for the Mach number. However, the triple point position is more difficult to predict. The KNP scheme in the form used in this work demonstrates the ability to capture the phenomena involved in the unsteady reflections.

1. Introduction

A blast wave is a strong perturbation generated by the sudden release of a large amount of energy. Among the different phenomena associated with research on shock wave propagation and their interactions and effects [1,2,3,4,5], the reflection of spherical shock waves over straight surfaces bears technical and scientific applications [6,7,8,9]. The reflections of spherical or cylindrical shock waves over straight surfaces represent unsteady reflection processes, and there is no analytical theory to describe them (p. 297 in [10]).
Brode numerically studied the propagation of a blast wave from the detonation of a spherical charge of TNT [11,12]. The peak over-pressures as a function of the shock radius were asymptotically obtained for a point source and an isothermal sphere. Some years later, Dewey and their collaborators used a photogrammetrical experimental technique to study the interaction of two identical spherical shock waves. They compared this interaction with the reflection of one of the spherical shocks from the ground [13,14]. They found that a smooth surface induces a stronger Mach stem and a higher triple-point trajectory. Furthermore, in 1981 Takayama and Sekiguchi [15] investigated the interaction problem of a spherical blast wave with a flat plate. In addition to the work of Dewey et al., Colella et al. [16] numerically studied the flow field of strong blast wave interactions with a plate considering different values for the height of burst (HOB). They observed the transition from a regular reflection to a double Mach reflection in these interactions. More recently, the problem of an unsteady cylindrical blast wave interaction with a flat plate was numerically investigated by [17,18]. Both works aim to understand the blast wave propagation, reflection, and its transition phenomenon as well as the flow features induced by the blast wave using a high-resolution Euler/Navier–Stokes solver. Previously, Ref. [15] studied a similar problem where the blast wave was generated by a shock tube. The authors found that analytical results did not agree well with those obtained from the experiments.
The shock wave reflection presents a very complex phenomenon. Ref. [10] defines the transitions and reflections of the Old State-of-the-knowledge (p. 154) and the New State-of-the-knowledge (p. 178). However, in this paper, we focus on the work of Hu and Glass [19] who theoretically analyzed the interaction of a spherical blast wave (in pseudo-steady flow) with a planar surface in perfect air under standard ambient conditions. When the height of burst (HOB) value is greater than a given limit, there are only two types of shock-wave reflection: regular reflection (RR) and single Mach reflection (SMR). When the HOB value is less than this limit, there are four types of shock-wave reflection: RR, SMR, transitional-Mach reflection (TMR), and double-Mach reflection (DMR). The unsteady interactions of shock waves propagating in gases are extensively studied in the work of Gvozdeva et al. [20,21,22]. Shock wave interactions with concave and convex corners as well as with curved surfaces are described along with the transitions between RR and IR. Dixon-Hiester et al. [23] experimentally investigated the RR to MR transition region and its flow characteristic. More recently, Kleine [24] presented an experimental and numerical investigation into the unsteady process of blast wave reflection from straight smooth surfaces. They found significant discrepancies between the numerical and experimental results. However, they conclude that further improvement and enhancement in the resolution and experimental detection of the Mach stem will reduce this discrepancy to acceptable levels. Recently, a study by Ridoux et. al [25] presented a simplified model for simulating blast wave propagation in different geometries at a low computational cost. They provide several comparisons with experimental results, including a case describing the reflection over a straight surface.
We are particularly interested in the reflection processes of an unsteady shock wave over a straight surface. The classic book of Ben-Dor [10] extensively studies these phenomena in steady, pseudo-steady, and unsteady flows, presenting analytical, numerical, and experimental results.
The main objective of this work is the evaluation of the Kurganov, Noelle, and Petrova (KNP) scheme [26,27,28] in simulating the unsteady shock wave reflection on a straight solid wall and for capturing the different shock wave reflection transitions. For this purpose, the solver rhoCentralFoam of the OpenFOAMTM software is employed. It is worth noting that this scheme has shown to be a useful tool for simulating supersonic flow patterns [29,30,31,32]. This software is widely used by the scientific community, and its validation will therefore be useful for many working groups scattered among various universities and research centers. Different relevant aspects of the numerical scheme can be found in reference [3]. This approach assumes that the shock wave is generated by an instantaneous release of energy concentrated in one point (spherical shock wave), or along a line (cylindrical shock wave). The verification processes of the KNP scheme on simulating the unsteady reflection processes of shock waves were performed by considering results from the pseudo-steady theory (There are no unsteady theoretical results.) [10,19] (The pseudo-steady result in Ben-Dor’s book [10] is found in p. 299.) and using other numerical studies [17]. Finally, a validation procedure was performed by using the experimental data for the triple point trajectory from [25].
This paper is composed of eight sections. Section 2 briefly describes the physical model and the implemented numerical scheme. In Section 3, a theoretical description of the reflections over straight surfaces generated by spherical and cylindrical shock waves is introduced. Section 4 presents the numerical results comparing them with the pseudo-steady results from [10,19]. Section 5 introduces a comparison between the OpenFOAMTM simulations with numerical results from [17], while a comparison with experimental data from [25] is performed in Section 6. Finally, a brief description and a discussion of the results are presented in Section 7, while the main conclusions are described in Section 8.

2. Physical Model and Numerical Scheme

This section briefly introduces the fundamental equations of the gas dynamics and the finite volume method implemented.

2.1. Gas Dynamics Equations

To simulate supersonic and non-viscous flows of gases, the Euler equations can be applied. These are a system of non–linear hyperbolic conservation laws that govern the dynamics of gases for which the effects of viscous stresses, body forces, and heat flow by conduction and radiation are neglected (p. 3 in [33]). This system can be written in a Cartesian coordinate system as
u t + F c x i = 0 , i = 1 , 2 , 3
where the vector of conservative variables is
u = ρ ρ U 1 ρ U 2 ρ U 3 E ,
and the flux vector is given by
F c = ρ V c ρ U 1 V c + n 1 p ρ U 2 V c + n 2 p ρ U 3 V c + n 3 p ρ E + p V c .
where U = U 1 , U 2 , U 3 T is the velocity vector, p the pressure, E the total energy, ρ the density, V c = U · n ^ the contravariant velocity [34], and n is the unit normal vector to the cell face.
The total energy for a perfect gas is
E = T 0 T c v d T R T 0 + 1 2 U 2
where c v is the specific heat for constant volume, and R is the gas constant.
To close the system, the state equation for perfect gases is utilized (thermally ideal gases, p. 7 in [33])
p = ρ R T
Physically, the system (1)–(3) arises from the application of the fundamental laws of conservation. The first equation corresponds to the mass conservation or continuity equation, the Second Law of Newton or momentum conservation is given by the second, third, and fourth equations, and finally, the fifth equation is the First Law of Thermodynamics or energy conservation.

2.2. Finite Volume Formulation

The system of equations described in the previous subsection, Equations (1)–(5), are solved using the finite volume method in a non-structured framework [35]. The physical domain is discretized in cells over which an approximated solution for transport equations will be obtained. All geometrical information related to the employed cell-centered discretization framework is presented in Figure 1, in which all the flow variables and the thermo-physical properties are stored in each cell’s centroid.
Then, utilizing the finite volume method, the system (1) is written as:
t V u dV + s F c d S = 0
where the surface integral is approximated by
s F c d S f ϕ f u f
in which ϕ = V c S f is the volumetric flux through the face, u f denotes the value of any field (i.e., pressure, velocity, density, etc.), the subscript f refers to the variable values at the face centroid, and V is the volume of the cell.
In high-speed flow, the flow properties are also transported by waves. This transport needs to be considered properly to obtain an accurate description of the flow behavior. Thus, the flux interpolation procedures are constructed using biased information to be more “upwind” of the query point. Here, convective terms are evaluated using the second-order central-upwind of Kurganov, Noelle and Petrova (KNP) [27]. The solver employs the adapted version to unstructured meshes of the KNP scheme. It is worth noting that this scheme has shown to be accurate for several supersonic flow patterns [29,30,31,36,37]. Therefore, the convective terms are evaluated as
f ϕ f u f = f α + ϕ f + u f + + α ϕ f u f + ω f u f u f +
where the first two terms correspond to flux evaluated in outward (+) and inward (-) directions of the cell’s face. The weighting factors, α ± , confer to this method an upwinding degree accounting for the local propagation rates of discontinuities [28]. The diffusive volumetric flux (third term) only arises in cases where the convective term to be calculated is part of a substantial derivative, which occurs in most cases. However, it is interesting to note by way of example that in the energy equation when calculating · U p this term will not be necessary. For the KNP scheme, the diffusive volumetric flux ( ω f ) and the weighting factors ( α ± ) are as follows [28]
ω f = ζ f + ζ f ζ f + + ζ f 1 α ± = ζ f ± ζ f + + ζ f 1 ,
where ζ f ± are the local volumetric fluxes, which are determined by the propagation velocities of the discontinuities at the interfaces. These fluxes in the + and − directions, respectively, are
ζ f + = max c f + S f + ϕ f + , c f S f + ϕ f , 0 ζ f = max c f + S f ϕ f + , c f S f ϕ f , 0
where c f ± = γ R T f ± is the sound velocity in incoming and outgoing directions of the interface, respectively.
Note that α ± is evaluated based on one-sided local propagation velocities (biased in the upwind direction), while ω f is based on the maximum of the propagation velocity of any discontinuity that may exist at a face between interpolated values in the + and − directions.
The reconstruction in unstructured meshes, used in the OpenFOAMTM software, in the outward (+) direction is given by:
u f + = 1 g f + u P + g f + u N
where the subscripts P and N refer to the centroids of the owning and neighboring cells, respectively, (Figure 1). The function
g f + = β 1 w f
where
w f = | S f · d f N | / | S f · d |
defines the shape of the reconstruction of β , allowing a wide range of schemes to be covered in a compact form. Where N is the neighboring cell, d is the vector between the centroids of the local and neighboring cells, d f N is the vector between the centroid of the face and that of the neighboring cell and S f is the face surface vector (Figure 1).
If β = β ( r ) , that is, the limiting functions are defined by successive gradients r of the reconstructed variable, which for an unstructured mesh results
r = 2 d · u P d u f 1
where u P is the evaluated gradient at P and d u f is the normal to the face gradient projection of u scaled by the magnitude of d . As noted, it is possible to switch between schemes of different orders by selecting the appropriate limiting function β ( r ) with r 0 .
For the evaluation of the gradients, to maintain consistency, the Kurganov schemes are also used [28]:
V u d V = s u f d s f u f s f = f α + u f + s f + α u f s f
If all the approximations previously introduced for convective terms and gradients are replaced in the finite volume discretization, a semi-discrete form of the transport equations is achieved:
t t + Δ t u t p V p + f u f s f + f ϕ f u f d t = 0
To complete the discretization, an evaluation of the temporal term is needed. If the time derivative is integrated, and all the spatial terms are denoted f ( t , u ( x , t ) ) , the last equation can be written:
u t + Δ t = f ( t , u ( x , t ) ) + u t V p
According to the information in time used to evaluate f ( t , u ( x , t ) ) , a particular time integration scheme is attained. Here, the explicit Euler scheme has been used.
The explicit Euler scheme was implemented and studied in [3] to simulate the propagation of intense shock waves generated by explosions. In that paper, it was established that:
  • To obtain better stability of the solutions, upwind reconstruction functions were recommended. Upwind reconstructions produce accurate results, and they do not generate non-physical oscillations.
  • Time-accurate solutions are obtained when CFL ≤ 0.5.

3. Spherical and Cylindrical Shock Wave Reflections over Straight Surfaces

For steady and pseudo-steady flows, the flow fields generated by the shock wave reflection phenomena depend on two independent variables, ( x , y ) in steady flows, and ( x / t , y / t ) in pseudo-steady flows. However, the reflection of a cylindrical and spherical shock wave on a straight solid wall produces unsteady flow fields, which are dependent on three variables ( x , y , t ) , where x and y are the space coordinates and t denotes the time (Cartesian coordinate system). Therefore, the analytical and numerical evaluation of the reflection phenomenon in unsteady flow conditions is more difficult than in steady or pseudo-steady flows (p. 246 in [10]).
There are three ways to generate unsteady shock wave reflections (p. 246 in [10]):
(1)
Reflecting on a non-straight surface of a shock wave moving with constant velocity.
(2)
Reflecting on a straight surface of a shock wave moving with variable velocity.
(3)
Reflecting on a non-straight surface of a shock wave moving with variable velocity.
Here, case (2) is studied using the finite volume method described in the previous section. The shock waves with non-constant velocity are spherical and cylindrical. They are generated by the sudden release of energy.
In general, the shock wave reflection configuration can either present as a regular reflection RR or an irregular reflection IR. The RR ⇄ IR transition occurs when the flow Mach number behind the reflected shock wave, with respect to the reflection point, verifies (p. 159 in [10]):
M 2 R = 1 ,
where M R 2 is the flow Mach number in the state (2), behind the reflected shock wave (see Figure 2), concerning the reflection point R.
The irregular reflection can either be a von-Neumann reflection vNR or a Mach reflection MR. Courant and Friedrichs [38] theoretically described three different types of Mach reflection. These depend on the propagation direction of the triple point of the Mach reflection concerning the solid reflecting surface:
  • If the triple point moves away from the reflective surface, then the MR is a direct Mach reflection: DiMR;
  • If the triple point moves parallel to the solid surface, then it is a stationary Mach reflection: StMR;
  • If the triple point moves closer to the surface, then it is an inverse Mach reflection: InMR.
Note that the three configurations for the Mach reflection (DiMR, StMR, and InMR) were theoretically described by Courant and Friedrichs in 1948 [38]. However, its experimental verification occurred almost 40 years later through the work of Ben-Dor and Takayama [39]. Furthermore, the StMR and InMR reflections only occur in unsteady flows. However, the DiMR reflections can be generated in pseudo-steady and unsteady flows and the DiMR are divided into three different types:
  • A single-Mach reflection (SMR);
  • A transitional-Mach reflection (TMR);
  • A double-Mach reflection (DMR).
Additional wave configurations may occur in the DMR depending on the direction of the triple point. Further discussion can be found in the book of Ben-Dor (p. 167 in [10]). In this paper, we are interested when a cylindrical or spherical shock wave reflects over a straight surface. Initially, a regular reflection (RR) is produced, which then, depending upon the strength of the shock wave, changes into a DMR, a TMR, and/or an SMR as it propagates outwards. Furthermore, as the wave propagates outward, then the point where the wave touches the reflecting surface encounters an ever-decreasing effective reflecting wedge angle. Thus, the instantaneous shock wave Mach number, M s , decreases with time.
To predict the transition between these reflections, different criteria were established in [10] (To better understand the shock reflections and their transition criterion, see Chapter 3 (p. 136) of [10].):
  • The first transition is from RR to IR. In this study, the IR can be a DMR, TMR or SMR depending on the value of H O B (height of burst). According to [19], for example:
    -
    If HOB = 0.8 , the transition will be from RR to DMR;
    -
    If HOB = 1.25 , the transition will be from RR to TMR;
    -
    If HOB = 2.0 , the transition will be from RR to SMR.
    The transition criterion for RR ⇄ IR was established in Equation (18).
  • The SMR ⇄ TMR criterion transition is (p. 161 in [10])
    M 2 T = 1
    where M 2 T is the flow Mach number in region (2) in a frame of reference attached to the first triple point T (see Figure 2).
  • The TMR ⇄ DMR transition criterion is (p. 176 in [10])
    M 2 T = 1
    where M 2 T is the flow Mach number in region (2) in a frame of reference attached to the second triple point T (see Figure 2).
Hu and Glass [19] showed different types of reflections for air (in pseudo-steady flows) in the ( M s , θ w ) plane. They presented a figure which relates the Mach number of the shock wave, M s , the angle between the shock and the surface, the non-dimensional height of burst, H O B , and the non-dimensional distance between the first contact point and the local instantaneous contact point, x ^ . The non-dimensional height of burst, H O B , is defined as
H ^ = H O B W W r p a r p a 1 / 3
where W is the weight (equivalent to TNT) of the explosive charge used to generate the blast (spherical shock) wave, W r = 1   kg of TNT, p a r = 1 atm , and p a the atmospheric pressure. Furthermore, the x ^ is
x ^ = x W W r p a r p a 1 / 3
Hu and Glass also showed the position of the triple point Y T / H O B . However, these results should be considered only as an approximation, since the pseudo-steady transition lines are not entirely applicable for unsteady flows. Differences of up to 10° have been verified in experiments (p. 300 in [10]).

4. Numerical Simulations

In this section, the test of study, the method to obtain the meshes, the boundary conditions, the mesh sensibility analysis and the comparison with the theoretical results of Ref. [19] are described.

4.1. Test of Study

Consider Figure 3 where a spherical (or cylindrical) shock wave at two different times is drawn. The shock wave propagates outwards to the point where it touches the reflecting surface, and it encounters an ever decreasing effective reflecting wedge angle θ w 2 < θ w 1 .
For simplicity, the denominator of Equation (21) is equal to one if the pressure outside of the explosion is chosen as p a = 1 atm and the energy is W = 1 kg of TNT, so H ^ = H O B in Equation (21). Air is used in all numerical simulations in this paper.
In addition, a cylindrical shock wave with an initial radius of R ( 0 ) = 10 cm is considered. The released energy is equal to 1.1936 kg of TNT ( E 0 i = 4.186 MJ), and the initial density of the undisturbed field is 1 kg/m3. The initial pressure is computed in the same way as [3]:
p 0 i = 3 ( γ 1 ) E 0 i ( β + 1 ) π R ( 0 ) β
where β = 2 for a cylindrical geometry and β = 3 for a spherical shock wave. Consequently, the pressure inside the high energy zone is p i = 53.3 MPa and the temperature is 155,524.02 K (considering that ρ 0 i = ρ 0 e = 1.1936 kg/m3 and the temperature outside the high energy zone is T i = 294.84 K).
To validate the solver, several simulations are performed varying the value of the height of burst. The following parameters are compared with the results of refs. [10,19] (The pseudo-steady result in Ben-Dor’s book [10] is found in p. 299.): the reflecting angle θ w , the distance of the reflecting shock x from the initial point of the explosion, the value of the Mach shock wave, and the position of the triple point ( Y T ).

4.2. Meshes and Boundaries

OpenFOAMTM always implements three-dimensional meshes. However, two dimensional simulations can be developed by using only one element in the flow normal direction and applying the empty boundary condition [35]. This approach was applied in all 2D cases here studied. Because of physical symmetry conditions, only 1/4 of the full domain was simulated as Figure 4 shows. An explosive 2D region (liberation zone) is centered on the origin of coordinates with an initial radius of R ( 0 ) = 0.1 m. The size of the domain varies according to the value of the H O B from 0.3 m to 2.0 m. All the numerical simulations are performed using the blockMesh utility included in the OpenFOAMTM software. For the grid independence analysis, performed in Section 4.3, structured meshes ranging from 87,600 to 2,180,000 cells were used. All the other numerical tests, described in Section 4.5, were performed with a fine mesh of 1.4 million of elements.
Figure 4 shows the geometrical domain configurations and the boundary conditions. Due to the symmetry of the problem, symmetryPlane condition is used on the bottom and left patches. Since a supersonic flow is guaranteed after the blast passes, only a simple extrapolated boundary condition is allowed in the top and right boundaries. Therefore, the zeroGradient condition is imposed on all flow variables for the pressure and temperature. For the velocity, the boundary condition in the right patch is fixedValue equal to zero.
The values of pressure p and temperature T for the high energy zone are determined as a function of the liberated energy, keeping the densities outside and inside equal. The utility setFields is used for setting these initial conditions.

4.3. Grid Independence

Firstly, a grid independence study is performed. The reflection angle θ w for the three different transitions described in Section 3 is compared with the pseudo-steady results of Refs. [10,19] for H O B = 0.8 mx. Figure 5 shows this study for the three transitions when the H O B = 0.8 m. The solid line shows the pseudo-steady result of [10], while the circles indicate the numerical simulations. Figure 5a shows the reflection angle from the RR ⇄ DMR. Figure 5b exhibits the same information but for the DMR ⇄ TMR. Finally, Figure 5c display the the TMR ⇄ SMR. In all cases, as the number of elements is increased, differences between the pseudo-steady arise and the numerical results decreases.
Therefore, a mesh of 1,440,000 ( 800 × 1800 ) cells will be used in the following sections. This mesh presents an adequate balance between grid sensibility and CPU time.

4.4. Time Step Sensitivity Analysis

The reflection angle θ w for the three different transitions described in Section 3 is compared with the pseudo-steady results of refs. [10,19] for H O B = 0.8 m. Instead of fixing a time step, we defined in OpenFOAMTM a maximum CFL number and the software automatically computes the Δ t . A mesh of 1,440,000 ( 800 × 1800 ) elements is used for all the simulations. Figure 6 shows the reflection angle for the three transitions RR ⇄ DMR, DRR ⇄ TMR, and TMR ⇄ SMR. Note that there is convergence for CFL < 0.5.

4.5. 2D Case. Comparisons with Pseudo-Steady Results

Let us start with the analysis by studying the case where H O B = 0.8 m. This H O B is interesting because there are three transitions: RR ⇄ DMR ⇄ TMR ⇄ SMR.

4.5.1. H O B = 0.8 m

Transition from RR to DMR

According to the results of Hu and Glass [19], the three different reflections described in Section 3 will occur. None of them are instantaneous as they occur in very short time periods.
Figure 7 shows the different times when the first transition occurs. Figure 7a,b shows the incident and reflected shock wave. Figure 7c–f displays the phenomenon very close to the wall. As soon as the shock wave arrives at the wall, an RR occurs. Graphically, it is possible to observe that the triple point T from Figure 2 has begun to set up together with the Mach stem, which may indicate that the reflection process has begun. However, the transition criterion from [10] states that the Mach number in the reflection point near the wall must be M 2 R = 1 . Table 1 shows the M 2 R = 1 for different times. Considering this, the RR ⇄ DMR transition takes place in t 550 μs.
Table 1 shows the M 2 R = 1 for different times. Considering this, the RR ⇄ DMR transition takes place in t 550 μs.
Table 2 shows the comparison between the pseudo-steady results and the numerical simulation for t = 550 μs and t = 600 μs. The difference between numerical and theoretical results, except for the Mach of the shock wave, is lower for the case t = 550 μs. Note that Table 1 shows similar results. A difference of 2.25 % exists between the pseudo-steady result and the numerical result for the reflecting angle θ w . This percentage is higher when the position of the triple point Y T / H O B is compared. This analysis should be considered descriptive from a phenomenological point of view since the pseudo-steady results are not completely applicable for unsteady flows. However, for the tests here studied the velocity of the shock wave at the straight surface changes only up to 12 % . Therefore, the pseudo-steady results could approximate the unsteady ones.

Transition from DMR to TMR

As the shock wave continues colliding with the reflecting surface, the DMR transitions to a TMR. From Figure 7f, it is possible to observe that near the triple point the reflected wave has changed its slope, which denotes that the second triple point is in the train to be formed. Once again, the transition criteria will give us the approximate time (Mach number in the second triple point is M 2 T = 1 ) where the transition takes place. Figure 8 shows the reflection for three different times. The right side of Table 3 depicts the values of the Mach number at the second triple point T . Therefore, the transition occurs at t 620 μs.
Table 3 also shows the difference between the pseudo-steady results and the numerical simulation. Again, a good agreement can be noted in the reflecting angle and position x from the center of the shock wave. However, major differences exist in the position of the triple point.

Transition from TMR to SMR

The last transition is the TMR ⇄ SMR transition. As mentioned in Section 3, the transition occurs when the Mach number in the triple point verifies M 2 T = 1 . The Mach numbers for different times are listed in Table 4 and the transitions occur in t 950 μs. Figure 9 shows the wave structure at these times and the position of the triple point.

4.5.2. H O B = 2.0 m

According to [10,19], for H O B = 2.0 m, only one transition takes place: from RR to SMR. This transition occurs when M 2 R = 1 . Figure 10 shows the triple-point formation between t = 2900 μs and t = 3100 μs, and the Mach stem is observed. Figure 10a–c shows the incident and reflected shock wave before the formation of the triple point. In Figure 10d ( t = 3600 μs), the triple point has already been formed.
Values of the M 2 R are summarized in Table 5, so it is considered that in t 3000 μs the transition takes place. The comparison between the pseudo-steady results and the numerical ones is also shown in Table 5. A good agreement in the reflection angle and distance x is observed.

4.5.3. H O B = 1.5 m

In this case, the phenomenon is similar to H O B = 2 m. Only exists one transition: RR ⇄ SMR. Table 6 summarizes the results of this case and the Mach number in the reflection point.

4.5.4. H O B = 1.0 m

This case is similar to H O B = 0.8 m (Section 4.5.1), but with the difference that the RR ⇄ DMR does not exists, and a RR ⇄ TMR transition takes place. In Figure 11, the Mach stem can be observed in t = 850 μs, which indicates that the transition is taking place. The Mach number in the reflection point M 2 R of Table 7 shows that the transition occurs in t 825 μs.
The differences between spherical and pseudo-steady with cylindrical and numerical results for the reflection angle and the distance are around 5 % . However, it is higher in the position of the triple point, as observed in the other cases studied.
Table 8 shows the results for the TMR ⇄ SMR transition. Note that the transition occurs between t = 1300 μs and t = 1325 μs. To compare with the spherical and pseudo-steady results, it is assumed that the transition occurs in t = 1325 μs. Differences of approximately 10 % are observed for θ w , x, and M s .
Figure 12 shows the numerical schlieren images for the pressure and Mach number. It allows us to observe the classic contact discontinuity near the triple point for some irregular reflections [10].

4.5.5. H O B = 0.4 m

The transitions and the wave structures are similar to those depicted for the H O B = 0.4 m case where the three transitions occur. The results for the three transitions are in Table 9.

4.6. 3D Case. Spherical Shock Waves

In this section, the results obtained for spherical configurations are presented. For the 3D computational domains, symmetry boundary conditions were employed, which allows simulating an octave of the physical domain. Dimensions of the computational domain were defined by considering the results from Ref. [19] in the following way: for a given height of burst and the transition angles, the lateral size of the computational domain can be determined by considering the H O B as the hypotenuse, and then the opposite side is determined as L = HOB tan θ , where θ is the reflecting wedge angle.
A spherical region with high-pressure and high-temperature gas was used as an igniter. The explosive charge energy is equivalent to 1 kg of TNT (4.184 MJ), and the pressure and temperature for non-perturbed regions are 101, 325 Pa and 300 K, respectively. The simulation is performed for a HOB of 0.8 m. This HOB with the non-perturbed conditions produces a blast wave with a strength equivalent to M s ≈ 3.7.
According to theory, the initial reflection is an RR which changes to a DMR, then to a TMR, and finally, to an SMR [10]. To obtain all possible reflections, L = 1.75 m was employed. Figure 13 shows the computational domain. For the 3D simulation, an octave of the domain is employed, as previously mentioned. Therefore, in the surfaces given by the vertices ache, abfh and cabd the symmetryPlane boundary condition was utilized. The reflecting conditions were imposed on the other three surfaces.
The initial conditions considered here generate an incident wave that first impacts the surface defined by the vertices cdih. Thus, given the domain dimensions, all the possible reflections mentioned have to be observed.
First, to examine the grid convergence in 3D configurations we consider three meshes: M 1 (2.7 × 10 6 cells), M 2 (5.4 × 10 6 cells) and M 3 (7.3 × 10 6 cells). As a first comparison, we present the reflected wave structure on the plane ache at t = 550 μs in Figure 14. At this time, the RR ⇄ DMR transition has already begun. As expected, a coarser grid (M1) gives the most diffusive results, and on the finer one, the structure of the reflected wave is better captured.
As mentioned earlier, the transition process is not instantaneous. In Figure 15, the transition from a RR to a DMR can be observed, and the triple point is clearly defined. The transition criterion states that the Mach number in the reflection point near the wall must be M 2 R = 1 . Values from the 3D simulation are presented in Table 10, according to this criterion the RR ⇄ DMR transition takes place at t = 550 μ s. As in the 2D cases, comparisons with the pseudo-steady results are presented in Table 11. It should be noted that the numerical results from the 3D simulation are in good agreement with the pseudo-steady predictions.
The DMR ⇄ TMR transition process can be observed in Figure 16. The figure presents a 3D contour plot of the Mach number distribution. From this figure can be appreciated the formation of a second triple point, also it can be observed that at the triple point location the Mach number approaches unity. The second triple point is wholly defined at t = 640 μ s. The right side of Table 12 shows the values of the Mach number at the second triple point T . Therefore, the transition occurs at t 620 μs. Table 12 also shows the difference between the pseudo-steady results and the numerical 3D simulation. Again, it is observed a good agreement with the reference data for the position of the triple point ( Y T / H O B ) and the position of the shock wave (x). However, major differences exist in the position of the reflecting angle.
Finally, results for the transition TMR ⇄ SMR are presented. This transition takes place when the Mach number in the triple point verifies M 2 T = 1 . The mach numbers for different times are listed in Table 13. From this information, it can be accepted that the transition occurs when M = 1.032 at t = 950 μs. Figure 17 shows the wave reflection configuration at four time steps. It is observed that, at this stage there is a single triple point that is completely defined. The results obtained are in well agreement with those provides by the pseudo-steady theory.

5. Comparison with Other Numerical Simulations

Liang et al. [17] solve the compressible and two-dimensional Euler/Navier–Stokes equations in a finite volume approach using a fifth-order weighted essentially non-oscillatory scheme with a fourth-order Runge–Kutta method. The numerical method was tested in [18] showing that the solver could resolve the four types of shock wave reflection: RR, SMR, TMR, and DMR.
The variation of the triple-point trajectory is shown in Figure 18 for two different initially incident shock Mach number M s = 4 and M s = 5 . The black dots are the results of ref. [17] and the red crosses are the OpenFOAMTM results. If we compare both figures for a given HOB, the higher the shock Mach number, the higher the triple point. The differences of the transition point (from RR ⇄ IR) have a difference of less than 8 % between our simulation and [17].
To develop the OpenFOAMTM simulations, we have used the same methodology described in the previous section. The C F L 0.5 , and a mesh of 1,440,000 ( 800 × 1800 ) elements is used.

6. Comparison with Experimental Data

In this section, we present the comparison between the experimental data given in Ref. [25] and the numerical results calculated using OpenFOAMTM. Ridoux et al. in [25] studied the dynamics of shock wave propagation (blast waves) and their interaction with straight surfaces.
Figure 19 shows the numerical results here obtained and the experimental data [25] for an explosion of 1 kg of TNT at a HOB of 1.59 m. The black dots represent the experimental data and the red crosses are the numerical results. After the transition RR ⇄ IR, which is well predicted by the numerical simulation, the difference between the numerical results and the experimental results increases from differences up to 10 % (in x 4 m) to a maximum of 20 % in x 4.5 m. These differences seem to continue growing in time as the shock wave continues with the reflection. A mesh of 1,440,000 ( 800 × 1800 ) elements is used for the simulation.

7. Discussion of Results

The reflections on straight surfaces of cylindrical and spherical shock waves are studied in this paper. The shock waves are generated by the sudden release of energy (blast waves). These types of reflections produce unsteady flow fields. To validate and analyze the numerical results here calculated, they are compared with theoretical studies on pseudo-steady reflections, numerical data obtained by other authors, and experimental results.
To contrast with the pseudo-steady analysis, several numerical tests for cylindrical shock waves were developed. These tests consider different distances between the shock waves origin and the reflecting surface (different height of burst): H O B = 0.4 , 0.8 , 1.0 , 1.5 , 2.0 . Accordingly, different types of transitions between reflections were analyzed. To develop these numerical tests, the mesh and CFL number convergences were previously established (see Figure 5 and Figure 6). From this comparison, we can conclude the following for the 2D cases:
  • There was accuracy for the reflecting angle θ w where the maximum difference is up to 5 % .
  • The greatest differences for the Mach number of the shock wave, M s , occur for the TMR ⇄ SMR ( 13 % , 15 % , and 20 % ). For other transitions, the differences are less significant (between 3 % and 6 % ) .
  • The numerical simulations could not properly predict the position of the triple point (from the reflecting wall). The differences between the pseudo-steady and numerical results measure approximately between 20 % and 30 % .
For spherical shock waves, the results show fair agreement with the pseudo-steady theory and the two-dimensional simulations. However, it is questionable to compare the 3D simulations with those in 2D and with the pseudo-steady results. This is because the 3D simulation can capture the wave structure and its evolution in greater detail. It was observed that for all the considered transitions that the simulation can correctly capture them in a descriptive way. It is worth noting that the 3D contour plot provides an effective way to analyze the complex reflection wave structure in detail.
For 3D simulations, we can conclude the following from the comparisons with the pseudo-steady theory:
  • In all transitions, the poorest approximation was observed for the triple point position (from the reflecting wall).
  • The numerical predictions for the Mach of the incident wave are in good agreement, but the highest difference was observed on the TMR ⇄ SMR ( 6.03 % ).
  • The predictions of the wedge angle can be considered in fair agreement, with the poor prediction for the RR ⇄ DMR ( 6.1 % ).
Another aspect of the 3D simulations is the involved data volume. The finer grid (7.2 million cells) produces an output of around 24 Gb. The data quantity increases as is needed to compute other derived variables such as pressure gradient. This implies that the required computational resources for post-processing activities must be appropriate configurations for proper data manipulation.
The comparison between the theoretical pseudo-steady study of refs. [10,19] and the numerical results should be considered only as descriptive from a physical perspective since the pseudo-steady analysis is not completely applicable to unsteady flows. However, the theoretical pseudo-steady results help to determine the sequence of transitions. Furthermore, for the cases numerically studied in this paper, the velocity of the shock wave at the straight surface is approximately constant (variation of less than 12 % ). Hence, the pseudo-steady analysis could approximate the unsteady phenomenon.
The numerical results published in Ref. [17] were used to compare with those calculated here via OpenFOAMTM. We can conclude that both numerical analyses show similar results and indicate the same physical behavior. Furthermore, both studies show for the same HOB that the height of the triple point increases as the Mach number of the incident shock wave increases. The relative difference in the position of the triple point between both numerical studies does not exceed 12 % .
Finally, the comparison with the experimental data from [25] shows a good prediction of the transition point. Minimal differences exist near the transition point between the simulations and the experiment. However, a maximum difference of 20 % is observed for x 4.5 m. Furthermore, we note that the numerical simulations correctly describe the physics involved in the shock wave reflection phenomenon.

8. Conclusions

In this paper, the verification of the proper behavior of the Kurganov, Noelle and Petrova numerical scheme included in the rhoCentralFoam solver of the OpenFOAMTM software for simulating the unsteady reflection of shock waves over straight surfaces was achieved. Several comparisons have been developed: with pseudo-steady theoretical results (there are no unsteady theoretical results), with other numerical simulations, and with experimental data. From these comparisons, the KNP scheme in the form used in this work has shown the ability to capture the phenomena involved in the unsteady reflections.
Good accuracy between the numerical and the pseudo-steady results was obtained for the 2D and 3D simulations. However, the numerical simulations could not properly predict the position of the triple point. The comparison with other numerical results shows that the difference in the position of the triple point between both numerical studies does not exceed 12 % . Similar result differences were observed in the position of the triple point between the numerical simulation and experimental result. Therefore, this scheme implemented in the OpenFOAMTM software is a robust tool for simulating the reflection on a straight surface of a variable velocity shock wave and for capturing the transitions between regular and irregular wave reflections involved in this phenomenon.

Author Contributions

Conceptualization, L.M., S.E. and L.G.M.; methodology, L.M., S.E. and L.G.M.; software, L.M. and L.G.M.; validation, L.M.; formal analysis, L.M.; investigation, L.M. and L.G.M.; resources, L.M., L.G.M. and S.E.; writing—original draft preparation, L.M.; writing—review and editing, L.M., S.E. and L.G.M.; visualization, L.M.; supervision, S.E. project administration, S.E. and L.G.M.; funding acquisition, S.E. and L.G.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the following projects CONICET-PUE-IDIT, “Vulnerabilidad de la infraestructura y del medio físico asociado al transporte de combustibles y almacenamiento”; FONCyT-PICT-2017 “Estudio de vulnerabilidad estructural de tanques y tuberías de almacenamiento de combustibles frente a cargas generadas por viento y explosiones”; Universidad Nacional de Córdoba “Desarrollo y aplicación de conocimientos teóricos, numéricos, experimentales y códigos computacionales en mecánica de fluidos e intermitencia caótica” and CONICET PIP 2021–2023 GI “Efectos de viento, explosiones y fuego en tanques de almacenamiento de combustibles”. Elaskar is director of these projects, Gutierrez Marcantoni participates as a researcher in the projects, and Monaldi participates within the framework of a doctoral fellowship.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank José Tamagno for their teachings.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
HOBHeight of burst
RRRegular reflection
IRIrregular reflection
SMRSingle-Mach reflection
TMRTransitional-Mach reflection
DMRDouble-Mach reflection
vNRvon-Neumann reflection
MRMach reflection
DiMRDirect Mach reflection
SiMRStationary Mach reflection
IniMRInverse Mach reflection
TTriple point
T’Double triple point
M s Mach number of the shock wave
θ w Reflected angle
Y T Position of the triple point
V p Proprietary cell volume
F c Convective flux
uTransported conservative field
f Face value
s Face surface vector

References

  1. Radchenko, P.A.; Batuev, S.P.; Radchenko, A.V. Numerical analysis of concrete fracture under shock wave loading. Phys. Mesomech. 2021, 24, 40–45. [Google Scholar] [CrossRef]
  2. Figuli, L.; Zvaková, Z.; Kavický, V.; Loveček, T. Dependency of the Blast Wave Pressure on the Amount of Used Booster. Symmetry 2021, 13, 1813. [Google Scholar] [CrossRef]
  3. Marcantoni, L.F.; Elaskar, S.; Tamagno, J.; Saldía, J.; Krause, G. An assessment of the OpenFOAM implementation of the KNP scheme to simulate strong explosions. Shock Waves 2021, 31, 193–202. [Google Scholar] [CrossRef]
  4. Chauhan, A.; Arora, R.; Siddiqui, M.J. Propagation of blast waves in a non-ideal magnetogasdynamics. Symmetry 2019, 11, 458. [Google Scholar] [CrossRef] [Green Version]
  5. Lechat, T.; Emmanuelli, A.; Dragna, D.; Ollivier, S. Propagation of spherical weak blast waves over rough periodic surfaces. Shock Waves 2021, 31, 379–398. [Google Scholar] [CrossRef]
  6. Cullis, I.G. Blast waves and how they interact with structures. BMJ Mil. Health 2001, 147, 16–26. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Draganić, H.; Sigmund, V. Blast loading on structures. Tehnički Vjesnik 2012, 19, 643–652. [Google Scholar]
  8. Nartu, M.K.; Kumar, M.; Ramisetti, S.B. Improved Methodology for Accurate Prediction of Blast Wave Clearing on a Finite Target. J. Eng. Mech. 2022, 148, 04022049. [Google Scholar] [CrossRef]
  9. Anas, S.M.; Alam, M. Comparison of existing empirical equations for blast peak positive overpressure from spherical free air and hemispherical surface bursts. Iran. J. Sci. Technol. 2022, 46, 965–984. [Google Scholar] [CrossRef]
  10. Ben-Dor, G. Shock Wave Reflection Phenomena, 2nd ed.; Springer: Berlin/Germany, Germany, 2007. [Google Scholar]
  11. Brode, H.L. Numerical solutions of spherical blast waves. J. Appl. Phys. 1955, 6, 766–775. [Google Scholar] [CrossRef]
  12. Brode, H.L. Blast wave from a spherical charge. Phys. Fluids 1959, 2, 217–229. [Google Scholar] [CrossRef]
  13. Dewey, J.; McMillin, D.; Classen, D. Photogrammetry of spherical shocks reflected from real and ideal surfaces. J. Fluid Mech. 1977, 81, 701–717. [Google Scholar] [CrossRef]
  14. Dewey, J.; McMillin, D. An analysis of the particle trajectories in spherical blast waves reflected from real and ideal surfaces. Can. J. Phys. 1981, 59, 1380–1390. [Google Scholar] [CrossRef]
  15. Takayama, K.; Sekiguchi, H. Formation and diffraction of spherical shock waves in a shock tube. Rep. Inst. High Speed Mech. Tohoku Univ. 1981, 43, 89–119. [Google Scholar]
  16. Colella, P.; Ferguson, R.E.; Glaz, H.M.; Kuhl, A.L. Mach reflection from an HE-driven blast wave. In Dynamics of Explosions; International Colloquium on Dynamics of Explosions and Reactive Systems: New York, NY, USA, 1986; pp. 388–421. [Google Scholar]
  17. Liang, S.; Hsu, J.; Wang, J. Numerical study of cylindrical blast-wave propagation and reflection. AIAA J. 2001, 39, 1152–1158. [Google Scholar] [CrossRef]
  18. Liang, S.; Hsu, J.; Chen, H. Numerical Study of Spherical Blast-Wave Propagation and Reflection. Shock Waves 2002, 12, 59–68. [Google Scholar] [CrossRef]
  19. Hu, T.; Glass, C.J.I. Blast wave reflection trajectories from a height of burst. AIAA J. 1986, 24, 607–610. [Google Scholar] [CrossRef]
  20. Bazhenova, T.V.; Gvozdeva, L.G.; Nettleton, M.A. Unsteady interactions of shock waves. Prog. Aerosp. Sci. 1984, 21, 249–331. [Google Scholar] [CrossRef]
  21. Gvozdeva, L.G.; Predvoditeleva, O.A.; Fokeev, V.P. Double Mach reflection of strong shock waves. Fluid Dyn. 1968, 3, 6–11. [Google Scholar] [CrossRef]
  22. Bazhenova, T.V.; Fokeev, V.P.; Gvozdeva, L.G. Regions of various forms of Mach reflection and its transition to regular reflection. Acta Astronaut. 1976, 3, 131–140. [Google Scholar] [CrossRef]
  23. Dixon-Hiester, L.; Reisler, R.; Keefer, J.; Ethridge, N. Shock enhancement at transition from regular to Mach reflection. Am. Inst. Phys. 1990, 208, 204–209. [Google Scholar]
  24. Kleine, H.; Timofeev, E.; Takayama, K. Reflection of blast waves from straight surfaces. In Shock Waves; Springer: Berlin/Heidelberg, Germany, 2005; pp. 1019–1024. [Google Scholar]
  25. Ridoux, J.; Lardjane, N.; Monasse, L.; Coulouvrat, F. Extension of geometrical shock dynamics for blast wave propagation. Shock Waves 2020, 30, 563–583. [Google Scholar] [CrossRef]
  26. Kurganov, A.; Tadmor, E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. J. Comput. Phys. 2000, 160, 241–282. [Google Scholar] [CrossRef] [Green Version]
  27. Kurganov, A.; Noelle, S.; Petrova, G. Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton—Jacobi equations. SIAM J. Sci. Comput. 2001, 23, 707–740. [Google Scholar] [CrossRef] [Green Version]
  28. Greenshields, C.; Weller, H.; Gasparini, H.; Reese, J. Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows. Int. J. Numer. Methods Fluids 2010, 63, 1–21. [Google Scholar] [CrossRef] [Green Version]
  29. Marcantoni, L.F.G.; Tamagno, J.P.; Elaskar, S.A. Two-dimensional numerical simulation of detonation celluler structures in H2-O2-Ar mixture with OpenFOAM. Int. J. Hydrogen Energy 2017, 42, 26102–26113. [Google Scholar] [CrossRef]
  30. Marcantoni, L.F.G.; Tamagno, J.P.; Elaskar, S.A. rhocentralRffoam: An Openfoam solver for high speed chemically active flows—Simulation of planar detonations. Comput. Phys. Commun. 2017, 219, 209–222. [Google Scholar]
  31. Marcantoni, L.G.; Tamagno, J.; Elaskar, S. A numerical study on the impact of chemical modeling on simulating methane-air detonations. Fuel 2019, 240, 289–298. [Google Scholar] [CrossRef]
  32. Azadboni, R.K.; Wen, J.X.; Heidari, A.; Wang, C. Numerical modeling of deflagration to detonation transition in inhomogeneous hydrogen/air mixtures. J. Loss Prev. Process. Ind. 2017, 49, 722–730. [Google Scholar] [CrossRef] [Green Version]
  33. Toro, E. Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed.; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
  34. Hung, C. Definition of contravariant velocity components. Theor. Fluid Mech. Meet. 2002, 3, 3202. [Google Scholar]
  35. Jasak, H.; Jemcov, A.; Tukovic, A. OpenFOAM: A C++ library for complex physics simulations. Int. Workshop Coupled Methods Numer. Dyn. 2007, 1000, 1–20. [Google Scholar]
  36. Espinoza, D.; Casseau, V.; Scanlon, T.; Brown, R. An open-source hybrid CDF-DNSC solver for high speed flows. AIP Conf. Proc. 2016, 1786, 50–70. [Google Scholar]
  37. Zang, B.; Vevek, U.; Lim, H.; Wei, X.; New, T. An assessment of OpenFOAM solver on RANS simulations of round supersonic free jets. Comput. Sci. 2018, 28, 18–31. [Google Scholar] [CrossRef]
  38. Courant, R.; Friedrichs, K. Supersonic Flow and Shock Waves; Wiley Interscience: New York, NY, USA, 1948. [Google Scholar]
  39. Ben-Dor, G.; Takayama, K. The dynamics of the transition from Mach to regular reflection over concave cylinders. Isr. J. Technol. 1986, 23, 71–74. [Google Scholar]
Figure 1. Finite volume discretization framework.
Figure 1. Finite volume discretization framework.
Symmetry 14 02048 g001
Figure 2. Case of pseudo-steady transitional Mach-reflection. The position of the reflection point R, the triple point T and the second triple point T’ is shown. Adapted from [10], p. 148.
Figure 2. Case of pseudo-steady transitional Mach-reflection. The position of the reflection point R, the triple point T and the second triple point T’ is shown. Adapted from [10], p. 148.
Symmetry 14 02048 g002
Figure 3. Schematic illustration of the instantaneous Mach number for the incident shock wave and the wedges angles of the spherical shock wave propagating over a straight surface.
Figure 3. Schematic illustration of the instantaneous Mach number for the incident shock wave and the wedges angles of the spherical shock wave propagating over a straight surface.
Symmetry 14 02048 g003
Figure 4. In the left figure, the domain configurations and boundary conditions for all the numerical simulations are sketched. The figure on the right shows the initial pressure for the case of H O B = 0.8 m.
Figure 4. In the left figure, the domain configurations and boundary conditions for all the numerical simulations are sketched. The figure on the right shows the initial pressure for the case of H O B = 0.8 m.
Symmetry 14 02048 g004
Figure 5. Values of the reflection angle θ w for different numbers of elements in the mesh for H O B = 0.8 m. The solid line shows the value of [10], while the circles indicate the numerical simulations.
Figure 5. Values of the reflection angle θ w for different numbers of elements in the mesh for H O B = 0.8 m. The solid line shows the value of [10], while the circles indicate the numerical simulations.
Symmetry 14 02048 g005
Figure 6. Values of the reflection angle θ w for different C F L maximum in the mesh for H O B = 0.8 m. The solid line shows the value of [10], while the circles indicate the numerical simulations.
Figure 6. Values of the reflection angle θ w for different C F L maximum in the mesh for H O B = 0.8 m. The solid line shows the value of [10], while the circles indicate the numerical simulations.
Symmetry 14 02048 g006
Figure 7. Numerical schlieren images of the pressure in the RR ⇄ DMR for H O B = 0.8 m.
Figure 7. Numerical schlieren images of the pressure in the RR ⇄ DMR for H O B = 0.8 m.
Symmetry 14 02048 g007
Figure 8. Wave structure obtained as numerical schlieren images of the pressure in the DMR ⇄ TMR transition for H O B = 0.8 m.
Figure 8. Wave structure obtained as numerical schlieren images of the pressure in the DMR ⇄ TMR transition for H O B = 0.8 m.
Symmetry 14 02048 g008
Figure 9. Wave structure obtained as numerical schlieren images of the pressure in the TMR ⇄ SMR transition for H O B = 0.8 m.
Figure 9. Wave structure obtained as numerical schlieren images of the pressure in the TMR ⇄ SMR transition for H O B = 0.8 m.
Symmetry 14 02048 g009
Figure 10. Wave structure obtained as numerical schlieren images of the pressure in the RR ⇄ SMR transition for H O B = 2.0 m.
Figure 10. Wave structure obtained as numerical schlieren images of the pressure in the RR ⇄ SMR transition for H O B = 2.0 m.
Symmetry 14 02048 g010
Figure 11. Wave structure obtained as numerical schlieren images of the pressure in the RR ⇄ TMR transition for H O B = 1.0 m.
Figure 11. Wave structure obtained as numerical schlieren images of the pressure in the RR ⇄ TMR transition for H O B = 1.0 m.
Symmetry 14 02048 g011
Figure 12. Wave structure obtained as numerical schlieren images in t = 1325 μs in the RR ⇄ TMR transition for H O B = 1.0 m of (a) pressure and (b) Mach number.
Figure 12. Wave structure obtained as numerical schlieren images in t = 1325 μs in the RR ⇄ TMR transition for H O B = 1.0 m of (a) pressure and (b) Mach number.
Symmetry 14 02048 g012
Figure 13. 3D domain for the case of H O B = 0.8 m.
Figure 13. 3D domain for the case of H O B = 0.8 m.
Symmetry 14 02048 g013
Figure 14. 3D domain for the case of H O B = 0.8 m.
Figure 14. 3D domain for the case of H O B = 0.8 m.
Symmetry 14 02048 g014
Figure 15. RR ⇄ DMR transition on 3D configuration H O B = 0.8 m.
Figure 15. RR ⇄ DMR transition on 3D configuration H O B = 0.8 m.
Symmetry 14 02048 g015
Figure 16. DMR ⇄ TMR transition on 3D configuration H O B = 0.8 m.
Figure 16. DMR ⇄ TMR transition on 3D configuration H O B = 0.8 m.
Symmetry 14 02048 g016
Figure 17. TMR ⇄ SMR transition on 3D configuration H O B = 0.8 m.
Figure 17. TMR ⇄ SMR transition on 3D configuration H O B = 0.8 m.
Symmetry 14 02048 g017
Figure 18. Height of triple point for H O B = 3 m. The black dots show the value of [17], while the red crosses are the numerical simulations with OpenFOAMTM.
Figure 18. Height of triple point for H O B = 3 m. The black dots show the value of [17], while the red crosses are the numerical simulations with OpenFOAMTM.
Symmetry 14 02048 g018
Figure 19. Comparison of the triple point trajectory between the experiment presented in [25] and the simulations with OpenFOAMTM. The black dots represent the experimental data. The red crosses are the numerical results.
Figure 19. Comparison of the triple point trajectory between the experiment presented in [25] and the simulations with OpenFOAMTM. The black dots represent the experimental data. The red crosses are the numerical results.
Symmetry 14 02048 g019
Table 1. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ DMR transition.
Table 1. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ DMR transition.
Time [μs] M 2 R = 1
5401.073
5501.071
5601.079
6001.085
Table 2. Comparison between the pseudo-steady results and the numerical simulation for the RR ⇄ DMR ( H O B = 0.8 m) transition.
Table 2. Comparison between the pseudo-steady results and the numerical simulation for the RR ⇄ DMR ( H O B = 0.8 m) transition.
Ben-Dor [10] t = 550 μs t = 600 μs
θ w [deg]50.7849.63546.85
M s 2.772.882.77
x [m]0.650.680.75
Y T / H O B 0.00050.00060.0065
Table 3. Values of the Mach Number in the double triple point M 2 T = 1 for different times to determine the DMR ⇄ TMR transition.
Table 3. Values of the Mach Number in the double triple point M 2 T = 1 for different times to determine the DMR ⇄ TMR transition.
Time [μs] M 2 T = 1 Ben-Dor [10] t = 620 μs
6101.0191 θ w [deg]45.0544.12
6201.0082 M s 2.492.62
6301.0258x [m]0.80.825
6401.0264 Y T / H O B 0.010.01625
Table 4. Values of the Mach Number in the triple point M 2 T = 1 for different times to determine the TMR ⇄ SMR transition.
Table 4. Values of the Mach Number in the triple point M 2 T = 1 for different times to determine the TMR ⇄ SMR transition.
Time [μs] M 2 T = 1 Ben-Dor [10] t = 950 μs
9401.0464 θ w [deg]33.1933.25
9501.0264 M s 1.992.306
9601.0313x [m]1.251.22
9701.0345 Y T / H O B 0.050.0687
Table 5. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ SMR transition ( H O B = 2.0 m).
Table 5. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ SMR transition ( H O B = 2.0 m).
Time [μs] M 2 R = 1 Ben-Dor [10] t = 3000 μs
28001.00817 θ w [deg]46.0147.231
29001.00523 M s 1.3211.6
30001.02218x [m]1.951.85
31001.0673 Y T / H O B 0.0010.0015
Table 6. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ SMR transition ( H O B = 1.5 m).
Table 6. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ SMR transition ( H O B = 1.5 m).
Time [μs] M 2 R = 1 Ben-Dor [10] t = 3000 μs
17001.0186 θ w [deg]48.9648.118
17251.0095 M s 1.61.872
17501.0032x [m]1.41.345
17751.0086 Y T / H O B 0.00050.0003
Table 7. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ TMR transition ( H O B = 1.0 m).
Table 7. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the RR ⇄ TMR transition ( H O B = 1.0 m).
Time [μs] M 2 R = 1 Ben-Dor [10] t = 825 μs
8001.02084 θ w [deg]50.9848.65
8251.0084 M s 2.132.235
8501.01458x [m]0.810.88
8751.02164 Y T / H O B 0.00060.0005
Table 8. Values of the Mach Number in the reflection point M 2 T = 1 for different times to determine the RR ⇄ TMR transition ( H O B = 1.0 m).
Table 8. Values of the Mach Number in the reflection point M 2 T = 1 for different times to determine the RR ⇄ TMR transition ( H O B = 1.0 m).
Time [μs] M 2 T = 1 Ben-Dor [10] t = 1325 μs
12751.01899 θ w [deg]39.079135.34
13001.00382 M s 1.8052.051
13250.9996x [m]1.31.051
13500.9933 Y T / H O B 0.050.06
Table 9. Values of the Mach Number in the reflection point M 2 T = 1 for different times to determine the RR ⇄ TMR transition ( H O B = 0.4 m).
Table 9. Values of the Mach Number in the reflection point M 2 T = 1 for different times to determine the RR ⇄ TMR transition ( H O B = 0.4 m).
RRDMRBen-Dor [10] t = 95 μs
θ w [deg]31.1631.6
M s 5.735.2
x [m]0.350.36
Y T / H O B 0.00050.0006
DMRTMRBen-Dor [10] t = 275 μs
θ w [deg]49.6148.01
M s 4.133.85
x [m]0.620.65
Y T / H O B 0.070.1
TMRSMRBen-Dor [10] t = 500 μs
θ w [deg]19.2521.22
M s 3.024.058
x [m]0.951.03
Y T / H O B 0.20.425
Table 10. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the transition RR ⇄ DMR on the 3D configuration.
Table 10. Values of the Mach Number in the reflection point M 2 R = 1 for different times to determine the transition RR ⇄ DMR on the 3D configuration.
Time [μs] M 2 R = 1
5400.9993
5501.004
6001.073
Table 11. Comparison between the pseudo-steady results and the numerical simulation for the RR ⇄ DMR on the 3D configuration ( H O B = 0.8 m) transition.
Table 11. Comparison between the pseudo-steady results and the numerical simulation for the RR ⇄ DMR on the 3D configuration ( H O B = 0.8 m) transition.
Ben-Dor [10] t = 550 μs t = 600 μs
θ w [deg]50.7847.8645.16
M s 2.772.632.70
x [m]0.650.6930.797
Y T / H O B 0.00050.00850.009
Table 12. Values of the Mach Number in the double triple point M 2 T = 1 for different times to determine the DMR ⇄ TMR transition.
Table 12. Values of the Mach Number in the double triple point M 2 T = 1 for different times to determine the DMR ⇄ TMR transition.
Time [μs] M 2 T = 1 Ben-Dor [10] t = 620 μs
6100.9930 θ w [deg]45.0543.89
6201.039 M s 2.492.1
6301.052x [m]0.80.77
6401.046 Y T / H O B 0.010.0157
Table 13. Values of the Mach Number in the triple point M 2 T = 1 for different times to determine the TMR ⇄ SMR transition on 3D configuration.
Table 13. Values of the Mach Number in the triple point M 2 T = 1 for different times to determine the TMR ⇄ SMR transition on 3D configuration.
Time [μs] M 2 T = 1 Ben-Dor [10] t = 6950 μs
9401.064 θ w [deg]33.1934.75
9501.032 M s 1.992.11
9601.075x [m]1.251.026
9701.076 Y T / H O B 0.050.0817
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Monaldi, L.; Marcantoni, L.G.; Elaskar, S. OpenFOAMTM Simulation of the Shock Wave Reflection in Unsteady Flow. Symmetry 2022, 14, 2048. https://doi.org/10.3390/sym14102048

AMA Style

Monaldi L, Marcantoni LG, Elaskar S. OpenFOAMTM Simulation of the Shock Wave Reflection in Unsteady Flow. Symmetry. 2022; 14(10):2048. https://doi.org/10.3390/sym14102048

Chicago/Turabian Style

Monaldi, Lucas, Luis Gutiérrez Marcantoni, and Sergio Elaskar. 2022. "OpenFOAMTM Simulation of the Shock Wave Reflection in Unsteady Flow" Symmetry 14, no. 10: 2048. https://doi.org/10.3390/sym14102048

APA Style

Monaldi, L., Marcantoni, L. G., & Elaskar, S. (2022). OpenFOAMTM Simulation of the Shock Wave Reflection in Unsteady Flow. Symmetry, 14(10), 2048. https://doi.org/10.3390/sym14102048

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop