1 Introduction

The pulse detonation engine (PDE) has been the focus of propulsion research efforts in the last few decades, due to its potential to drastically increase the efficiency when compared to conventional gas turbines [1,2,3]. In a hybrid-PDE configuration, an annular array of pulse detonation combustors (PDCs) replaces conventional isobaric combustion chambers. One of the main challenges for implementing a hybrid PDE is to maintain reliable and efficient operation of the downstream turbine. Poor coupling of the PDC with the turbine could eliminate any potential gain in cycle efficiency provided from pressure gain combustion. The inherent cyclic operation of a PDC results in highly transient pressure and temperature fields at the PDC exhaust, which are highly undesirable for conventional turbines.

Efficient attenuation of the leading shock wave in the PDC exhaust flow is essential for turbine work extraction. Numerous studies of PDC–turbine applications confirm the occurrence of strong reflected shock waves at the turbine inlet resulting in significant losses [4,5,6,7,8]. Therefore, it is desirable to attenuate the strong leading shock wave transmitted from the detonation wave before it enters the turbine.

Fig. 1
figure 1

Schematic diagram of the shock tube experimental facility with the shock divider assembly attached

Full size image

There are a number of studies on devices influencing the PDE exhaust flow. While most of the studies focus on the impact of different nozzles and ejectors on the performance of a single PDE engine [9,10,11,12,13,14], limited research has focused on the mitigation of the leading shock pressure in the PDE exhaust flow [15], despite its relevance for the integration of the PDC with a downstream turbine. One approach is to divide the leading shock wave in multiple weaker shock waves. In this manner, the energy produced in a single PDE cycle will be redistributed temporally. The temporal redistribution of the pressure across one PDE cycle results inherently in lower pressure fluctuations in a multi-cycle operation mode. Hence, a more suitable turbine-inlet flow can be obtained by dividing the leading shock wave into weaker shock waves.

A shock divider (thereafter simply termed divider) is proposed here as a method to spread the leading shock energy. This divider consists of a bifurcating section and a recombination section. In the bifurcating section, the incident strong shock wave is split and guided through two different ducts. The separated shock waves travel with different velocities along two pathways with different cross sections. The two shock waves are then transmitted into a single duct, after leaving the device. The key requirement for a divider design is a large temporal separation of the shocks at a minimum addition of losses.

Numerous studies on the propagation of shock waves through ducts with area changes were conducted in the last decades [16,17,18,19]. It was shown that the shock wave Mach number and the area ratio significantly impact the flow evolution. Moreover, the shock wave propagation in branched ducts was subject of various studies [20,21,22,23]. It was shown that losses due to reflection and diffraction of the shock wave can result in significant shock attenuation. However, to the best of the authors’ knowledge, there is only very limited research on the recombination of separated shock waves.

For an efficient divider design, it is crucial to understand the determining fundamental flow dynamical mechanisms. Therefore, the aim of the current study is to provide a better understanding of the flow physics inside a generic divider. We employ numerical simulations and experimental methods to study the flow inside a divider with an incident shock Mach number of \(M_\mathrm{s} = 1.61\). For this purpose, the divider is mounted at the exit of an open-end shock tube, where high-speed schlieren images along with total pressure measurements are conducted. We use a second-order dimensional split finite volume MUSCL-scheme to solve the compressible Euler equations. Furthermore, the propagation of the leading shock waves inside the divider is modelled by using a simplified model called geometrical shock dynamics (GSD) [18]. We first validate our numerical results based on schlieren images and pressure recordings. Different divider geometries are then evaluated numerically to gain a better understanding of their impact on the flow evolution.

2 Methodology

2.1 Shock tube facility

Experiments are undertaken using an open-end shock tube facility. The facility, manufactured from stainless steel, consists of a 35.1-mm-diameter driver section and a 12.5-mm-diameter driven section. A diagram of the facility is shown in Fig. 1. The driver and driven sections are separated by a diaphragm made of polyester film with a 0.1-mm thickness. The operating gas for the driver section is a mixture of atmospheric air and helium, while the operating gas for the driven section is atmospheric air. Bursting of the diaphragm occurs through actuation of a linear solenoid and plunger. The driver section is 250 mm long with the solenoid located in a 90-mm section, unsealed from the primary driver section. The driven section is 675 mm long and contains a 25-mm conical converging section beginning 50 mm from the diaphragm. The driver and driven sections are held together under pressure using a high-pressure quick release clamp. Further details on the facility are given in [14].

Gauge pressure is measured in the driver section using a Gems Series 3100 Pressure Transducer to give the diaphragm pressure ratio. The signals from the pressure transducer are captured on a 16-bit National Instruments DAQ. The trigger signal for the Ledex linear solenoid is provided as an output from the DAQ.

2.2 Shock divider

To separate a strong shock wave in multiple weaker shock waves, the assembly shown in Fig. 2 is utilised. This assembly is termed the shock divider. The shock divider connects to the end of the shock tube facility, as shown in Fig. 1. A shock wave is generated by the shock tube, which then travels into the shock divider assembly. The shock first travels through the round-to-square transition section (Fig. 2). There, the cross section of the geometry is transformed from a circular cross section to a square cross section using a seventh-order spline, as proposed by Wilson [24]. After the transition, the cross section of the geometry remains square throughout the divider assembly. The origin of the Cartesian coordinate system is the inlet of the divider centre body (Fig. 2).

Fig. 2
figure 2

Schematic diagram of the shock divider assembly

Full size image

After the shock wave passes through the round-to-square transition section, it enters the central divider assembly. The divider assembly, indicated in Fig. 2, separates the initial incident shock wave into two shock waves. Each shock is allowed to travel along different paths with different duct width and path length. The central divider assembly consists of two flat transparent perspex sides surrounded by a central machined aluminium section. The perspex walls allow for optical access to the internal flow within the divider. The separated shock waves then enter the same exit pathway before being diffracted out of the open end of the divider.

Three different shock dividers are investigated in the current study, as shown in Fig. 3. The investigated shock dividers consist of two branches: the upper branch and the lower branch. The lower branch in all dividers is a simple constant-width straight duct, whereas the upper branch is characterised by a curvature.

Fig. 3
figure 3

Divider assemblies. Dimensions are normalised by the divider square side length D

Full size image

The leading edge of the divider centre body is located at the divider inlet prior to the separation of the divider branches (Fig. 3a). Therefore, the shock wave does not diffract before it enters the divider. Consequently, the shock waves transmitted to the divider branches are expected to have nearly the same Mach number at the divider inlet. To induce different arrival times of the shock waves at the divider outlet, the path of the upper branch is longer than that of the lower branch. While the duct width of the lower and upper branches is equal in divider I (Fig. 3b), this is not the case for dividers II and III (Fig. 3c, d). In a parametric study, the width of the upper branch is increased from 0.5 to \(2\ D\) from dividers I–II (Fig. 3). The upper branches are modified simply by changing the centre body of the divider. Hence, the angle \(\alpha = 19\) shown in Fig. 3a remains unchanged for all dividers, whereas \(\beta \) decreases from dividers I–III. The angle \(\beta \) is 19, 13, and 7 for dividers I, II, and III, respectively. Accordingly, the transmitted shock wave into the upper channel of divider III faces a distinctive increase in cross section, a smaller increase in divider II, and no change in cross section in dividers I. The length of the channel centreline is 5, 4, and 2% longer than the lower branch for dividers I–III, respectively. The Mach number of the incident shock wave at the divider inlet is \(M_\mathrm{s} = 1.61\) for all investigated configurations in the current study.

2.3 Schlieren setup

To complement the numerical results, schlieren measurements of the shock divider flow are taken using a Toepler Z-Type schlieren system [25], as shown in Fig. 4. The 8” schlieren mirrors have a focal length of 1219 mm. Images are acquired using a Photron FASTCAM SA-Z 2100K up to a frame rate of 210 kHz. This provides ultra-high-speed images of the motion of the shock wave within the divider. The flow motion is illuminated using a pulsed LED light source with an exposure time of 1 \(\upmu \)s [26]. The signal from PCB pressure transducer located 50 mm from the divider inlet is used as a trigger for the camera. A precise timing control is provided by a Beaglebone Black processor, developed by Fedrizzi and Soria [27].

Fig. 4
figure 4

Schematic diagram of the shock tube experimental facility with attached shock divider assembly

Full size image

2.4 High-frequency total pressure probe

Time resolved measurement of total pressure behind transient shock waves is very challenging [28,29,30]. In this work, an in-house-made high-frequency total pressure probe is used to measure the total pressure at the exit of the divider. Figure 5a shows two photographs of the probe. The probe design is based on the probe used previously by Paxson and Dougherty [31] to measure the total pressure behind a pulsejet. To allow for high-frequency pressure measurements, the transducer is placed directly at the head of the probe (Fig. 5). A Kulite XCE-062 miniature transducer is mounted in a L-shaped tube, allowing the sensor to be placed within the divider (Fig. 5b). A Kulite KSC-2 signal conditioner is used to amplify the measured signal. The relatively high resonance frequency of the sensor allows for a nearly non-oscillating signal behind the shock wave [29, 32]. The probe is mounted at \(x/D = 12.5\) inside the divider. Simultaneous schlieren measurements of the probe head are used to ensure the absence of a bow shock in front of the sensor.

Fig. 5
figure 5

a Photographs of the total pressure probe and b the measurement position inside the divider

Full size image

2.5 Numerical methods

2.5.1 Geometrical shock dynamics (GSD)

Geometrical shock dynamics (GSD) is a simplified approach to predict shock wave propagation. It was first introduced by Whitham [33]. An illustration of the approach is presented in Fig. 6. The method is based on the decomposition of the shock front into elementary parts propagating independently along ray tubes. A ray tube is treated as a duct of cross-sectional area A with rigid walls. The approach is based on the assumption that the motion of the shock only depends on the variation of the local ray tube. Hence, the motion of the shock wave is determined without calculating the flow field downstream of the shock wave.

Fig. 6
figure 6

Rectangular (x, y) and curvilinear \((\alpha , \beta )\) coordinate systems. Shock front given by \(\alpha = \mathrm{const}\) and rays given by \(\beta =\mathrm{const}\)

Full size image

The motion of the shock wave is determined by a relation between the cross-sectional area A and the local Mach number M, using the A–M relation:

$$\begin{aligned} \frac{A_i(t)}{A_i(0)} = \frac{f(M_i(t))}{f(M_i(0))}, \end{aligned}$$
(1)

where the subscript i denotes the index of each segment and A and M are their area and shock Mach number, respectively. The function f(M) is given as

$$\begin{aligned}&f(M) = \exp (-f_e(M)), \end{aligned}$$
(2)
$$\begin{aligned}&f_e(M) =\int _{M_\mathrm{min}}^{M} \frac{M\lambda (M)}{M^2 -1}dM, \end{aligned}$$
(3)
$$\begin{aligned}&\lambda (M) = \left( 1+\frac{2}{\gamma +1}\frac{1-\mu ^{2}}{\mu }\right) \left( 1+2\mu +\frac{1}{M^2}\right) , \end{aligned}$$
(4)
$$\begin{aligned}&\mu ^2 = \frac{(\gamma -1)M^2+2}{2\gamma M^2-(\gamma -1)}, \end{aligned}$$
(5)

where \(\gamma \) is the specific heat ratio.

Numerous methods are used for numerical implementation of GSD [34,35,36,37,38]. The implementation used in the current study is based on the front tracking method introduced by Henshaw et al. [34]. A short overview of the applied numerical implementation is given in the following. Detailed information regarding the numerical scheme can be found in [34].

The discretised elements of the shock wave propagate along rays normal to the shock front \(\mathbf {n_i}\), with a local speed depending on the local Mach number \(M_i\) (Fig. 6). As shown by Henshaw et al. [34], the transformation from the curvilinear coordinate system to the rectangular system along a ray leads to

$$\begin{aligned}&\frac{\partial x(\beta , t)}{\partial \alpha }=M(\beta , t) \cos \theta , \quad \end{aligned}$$
(6)
$$\begin{aligned}&\frac{\partial y(\beta , t)}{\partial \alpha }=M(\beta , t) \sin \theta , \end{aligned}$$
(7)
$$\begin{aligned}&\frac{\mathrm {d}}{\mathrm {d} t} {\mathbf {x}}_{i}(t)=a_0 M_{i}(t) {\mathbf {n}}_{i}(t), \quad i=1, \ldots , N. \end{aligned}$$
(8)

Equation (8) is the vector form of (6)–(7), while \(\alpha \) is eliminated in favour of time t using \(\alpha = a_0 t\), where \(a_0\) is the speed of sound upstream of the shock. This results in a nonlinear system of ordinary differential equations, which is closed by (2)–(5).

Fig. 7
figure 7

a CFD domain and b grid with 46 cells in y-direction discretising the inlet channel of divider II (\(N_\mathrm{D} = 46 \))

Full size image

An explicit second-order two-step leap-frog scheme is used for the numerical integration, while the time step \(\varDelta t\) is adapted in every step to maintain the CFL condition:

$$\begin{aligned}&\mathbf {x_i}(t+\varDelta t)=\mathbf {x_i}(t-\varDelta t)+2\varDelta t a_0 M_i(t) \mathbf {n_i}(t). \end{aligned}$$
(9)
$$\begin{aligned}&A_{i}(t)=\frac{1}{2}\left\{ \begin{array}{ll}{s_{i+1}(t)-s_{i}(t),} &{} { \text{ if } i=1;} \\ {s_{i+1}(t)-s_{i-1}(t),} &{} { \text{ if } i=2, \ldots , N-1;} \\ {s_{i}(t)-s_{i-1}(t),} &{} { \text{ if } i=N,}\end{array}\right. \end{aligned}$$
(10)

where the arclength s(t) represents the geometry of the shock and is given by

$$\begin{aligned} s_{i}(t)=\left\{ \begin{array}{ll}{0,} &{} { \text{ if } i=1;} \\ {s_{i-1}(t)+\left| x_{i}(t)-x_{i-1}(t)\right| ,} &{} { \text{ if } i=2, \ldots , N.}\end{array}\right. \end{aligned}$$
(11)

The Mach number in (9) is determined by using (1). As given in (10), the cross-sectional area \(A_i(t)\) is determined by a one-sided scheme at the endpoints and a centred scheme about the point \(x_i(t)\).

2.5.2 Numerical simulation scheme (CFD)

Numerical simulation of the flow inside the shock divider is conducted based on the compressible Euler equations. The equations are given in their two-dimensional conservative form as

$$\begin{aligned} \frac{\partial }{\partial t} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho E \end{pmatrix} + \frac{\partial }{\partial x} \begin{pmatrix} \rho u \\ \rho u^2 + p \\ \rho v u \\ u(\rho E + p) \end{pmatrix} + \frac{\partial }{\partial y} \begin{pmatrix} \rho u \\ \rho u v \\ \rho v^2 + p \\ v(\rho E + p) \end{pmatrix} = 0, \end{aligned}$$
(12)

where \(\rho \) is the density, \({\mathbf {u}} = (u, v)\) the particle velocity, \(E = e + \Vert {\mathbf {u}} \Vert ^2 / 2\) the total energy, e the internal energy, and p the pressure. Equations (12) are closed by the caloric perfect gas approximation \(p = \rho (\gamma - 1) e\), where the adiabatic index \(\gamma = C_p/ C_v = 1.4\) is assumed to be constant. To solve these equations, a fully conservative second-order dimensional split finite volume MUSCL-scheme is used [39, 40]. To prevent artificial oscillations in the numerical solution, the slopes of the reconstruction step are limited with the van-Leer limiter. The solution is discretised on a structured grid. A level set is used to embed the boundary of the divider into the structured grid and results into irregular cut cells, which need special treatment. A conservative cut-cell method is used to ensure stability for arbitrarily small cells such that no mass flows through the boundary [41, 42].

In the simulations, the shock tube and the divider are initially filled with air at atmospheric pressure and room temperature. An initial Riemann problem is set up at the exit of the shock tube. Using the Rankine–Hugoniot conditions, the left state of the Riemann problem is set to be the post-shock state of a shock with a Mach number of \(M_\mathrm{s} = 1.61\). Transmissive boundary conditions are set at the upstream and downstream domain boundaries.

The measurement domain used for the CFD simulations is shown exemplary for divider II in Fig 7a. Figure 7b illustrates the Cartesian structured grid using 46 cells in y-direction to discretise the divider channel width (\(N_\mathrm{D} = 46\)). While \(N_\mathrm{D} = 46\) is used for all simulations conducted in this work, the total number of cells increases slightly from dividers I to III. The total number of cells for the entire domain shown in Fig 7a is 198,603, 204,846, and 216,297 for dividers I–III, respectively.

Five simulations using a Cartesian structured grid with different resolutions from \(N_\mathrm{D} = 11\)–184 are conducted to investigate the grid dependency. For the grid dependency study, divider III is used, which represents presumably the most critical configuration, as it results in the strongest diffraction of the shock wave in its upper branch. Figure 8a illustrates the shock wave inside the upper branch of divider III for the time when the shock wave reaches the centre of the divider’s centre body; lower end of the blue line in Fig. 8a is at \(x/D = 6\). A pink box in Fig. 8a represents the range chosen for the data shown in Fig. 8b. The pressure shown in Fig. 8b represents the pressure from \(x/D = 2\) to 6.4 and \(y/D = 0.4\) for different grid resolutions at \(t = 0.156\) ms. Downstream of the shock wave \(2< x/D < 6\), the pressure profiles differ only slightly, except for \(N_\mathrm{D} = 11\). The pressure in the vicinity of the shock wave is presented in Fig. 8c. The results show that regardless of grid resolution the shock wave is captured approximately across four to six cells. Hence, the pressure rise across the shock wave is captured more distinctively for finer grid resolution. However, even the coarsest grid using only 11 cells to discretise the divider’s channel width (\(N_\mathrm{D} = 11\)) resolves the position of the shock wave remarkably well when compared to the finer grids (Fig. 8b, c). While the pressure downstream of the shock wave is about 5% larger for \(N_\mathrm{D} = 11\) compared to the others, the results of \(N_\mathrm{D} = 23\)–184 do not differ significantly regarding both the position and the pressure downstream of the shock.

Fig. 8
figure 8

a Schematic illustration of the incident shock wave inside the upper channel of divider III. b Pressure in the cell next to the centre body in the vicinity and far downstream of the incident shock at \(t = 0.156\) ms. c Zoomed view of the data shown in b

Full size image

3 Results and discussion

We start with the validation of the numerical schemes by comparison with the experimental results. The validated numerical results are then used to analyse the flow evolution inside the shock dividers. Finally, we evaluate the impact of the divider width ratio on the separation of the shock waves.

3.1 Validation of numerical schemes

Figure 9 presents a series of experimental and numerical snapshots spanning the early evolution of the flow inside divider I. Experimental (EXP) schlieren images are compared with CFD schlieren images and GSD results. Figure 9a shows the incident shock wave just before it enters the shock divider. The time given above the images is given relative to the moment the incident shock wave passes the divider inlet at \(x/D = 0\).

As seen in Fig. 9b, at \(t = 0.008\) ms, the incident shock wave is divided in two separated shock waves, which propagate into the upper and lower branches of the divider. The figure shows that the position and shape of the separated incident shock waves in both branches of the divider are very well captured by both numerical schemes. As mentioned in Sect. 2.5.1, the GSD approach considers only the evolution of the incident shock waves. Therefore, the post-shock structures in the flow are not captured in the GSD results. However, the CFD results show an excellent agreement with the experimental schlieren images in terms of separated shock waves as well as flow structures upstream of the incident shock waves. This agreement is well demonstrated in Fig. 9b–d. These flow structures will be discussed in more detail in Sect. 3.2. The qualitative comparison of the results based on Fig. 9 shows an excellent agreement between the numerical and experimental results.

As mentioned before, GSD neglects any influence of the post-shock flow on the propagation of the leading shock. Therefore, the good agreement between GSD and experiment leads to the conclusion that the post-flow does not significantly affect the shock propagation inside the divider. Furthermore, no viscosity is considered in the CFD approach, yet the agreement between CFD and experiment is very good. Consequently, the impact of turbulence on the flow inside the divider must be marginal.

Fig. 9
figure 9

Experimental schlieren (left row), CFD schlieren (mid row), GSD results (right row) showing the flow evolution inside the divider I

Full size image
Fig. 10
figure 10

X–t diagram of the leading shock at the upper and lower branches of the divider based on experiment, CFD and GSD simulations for divider I

Full size image

To compare the results in a quantitative manner, the position of the separated shock wave in the lower and upper branches is captured at the outer wall of the divider: at the upper wall of the upper branch and at the lower wall of the lower branch. Figure 10 presents the corresponding x–t diagram of the transient shock waves based on the CFD, GSD, and experimental data. The comparison demonstrates a very good agreement between the CFD and GSD results regarding the evolution of the leading shock wave in both divider branches. However, both numerical schemes slightly overestimate the shock propagation velocity compared to the experiment. The spatial displacements between the shock waves in the upper and lower branches \(\xi \) differ slightly between the numerical and experimental data, as shown in Fig. 10. The quantity \(\xi \) at the divider exit is in the CFD case 4% and in the GSD case 20% smaller compared to the experiment (Fig. 10).

The numerical approaches slightly overestimate the shock strength as they do not account for all mechanisms, leading to entropy generation. The discrepancy between the experimental and numerical results is larger in the upper divider branch compared to the lower branch. This is to be expected, as multiple mechanisms in the upper branch lead to entropy generation, whereas in the lower branch the shock wave simply propagates through a constant cross-sectional straight duct. The entropy generation mechanisms will be discussed to some extent in the next section. The comparison between the numerical and experimental data shows the capability of both GSD and CFD approaches to resolve the propagation of the incident shock waves in the divider branches. However, unsurprisingly the CFD results replicate the experiments more accurately compared to the GSD.

3.2 Divider flow evolution

The CFD results are used to analyse the flow evolution inside the shock divider. Figure 11 presents four snapshots of contour plot pairs showing the pressure and Mach number distributions inside divider II. In the Mach number plots, the subsonic and supersonic regions of the divider flow are colour-coded with blue and red colours, respectively.

Fig. 11
figure 11

Time series of pressure and Mach number contours showing the evolution of the flow inside the divider II

Full size image

Once the shock wave reaches the divider inlet at \(x/D = 0\), a number of different events take place. The planar incident shock separates in two different shock waves. One propagates along the lower branch and the other along the upper branch. As shown in Fig. 11a, the part in the lower branch is transmitted as a planar shock wave, whereas the separated shock in the upper branch is slightly curved. Here, the upper part of the shock wave diffracts at the convex corner of the divider at \(x/D = 0\). Furthermore, a small part of the incident shock wave reflects at the leading edge of the centre body and propagates upstream (Fig. 11a-i). This reflected shock wave is linked at the triple point with the Mach stem close to the wall, and the diffracted shock (Fig. 11a-iv). This triple point configuration is known as a single Mach reflection [43]. As shown in the Mach number contour plots, the flow inside the entire divider is subsonic at this stage. This is to be expected as the shock strength of \(M_\mathrm{s} = 1.61\) planar shock results in subsonic post-shock flow conditions based on 1D gas dynamic equations.

The flow inside the divider accelerates due to various mechanisms. The diffraction of the shock wave at the convex corner of the upper branch (\(x/D = 0\)) results in an upstream propagating expansion wave (Fig. 11a-i). Hence, the flow velocity increases, resulting in supersonic flow at the convex corner of the divider inlet. This supersonic region and a steady, but growing, Prandtl–Meyer (PM) expansion wave are highlighted in the close-up view, shown in Fig. 11b-iii. A quasi-steady normal shock adjusts the pressure mismatch between the expanded flow through the PM expansion wave and the compressed flow by the transient shock in the upper branch (Fig. 11b-iii). Furthermore, an additional PM expansion wave in the upper branch at \(x/D = 1\) can be observed in Fig. 11b-iv. This is initiated by a barely visible convex corner in the divider centre body at \(x/D = 1\) (Fig. 11b-ii). The duct width in the upper branch of the divider is constant from \(x/D = 0\) to 1, while for \(x/D > 1\), the duct width increases linearly up to the centre of the divider at x\(/D = 6\), resulting in a convex corner at \(x/D = 1\) (Fig. 11b). The convex corner leads to further acceleration of the supersonic flow through a steady PM expansion wave (Fig. 11b-iv).

The flow inside the upper and lower branches differs considerably from one another. Within the diverging section of the upper branch, the leading shock wave is a single shock with a nearly planar shape. However, the shock structure changes significantly as the planar shock reaches the turning point of the divider at \(x/D = 6\). As the shock wave transmits to the concave section, a relatively strong Mach reflection occurs at the divider outer wall, resulting in a non-planar shock wave (Fig. 11c-i). The triple point of the Mach reflection reflects several times from both the inner and outer walls of the duct as the shock wave propagates through the divider (not shown here). Consequently, the pressure distribution in the converging section of the upper branch (\(x/D = 6\) to 11) is highly non-uniform (Fig. 11d-i). Furthermore, the flow Mach number decreases from the convex corner at \(x/D = 1\) up to the leading shock wave at \(x/D \approx \) 5.5, as shown in Fig. 11b-ii. This is mainly caused by the diverging shape of the upper duct, which decelerates the transient shock wave. On the contrary, the shock wave becomes gradually stronger as it propagates along the converging section of the upper branch. Hence, the flow Mach number in the converging section (\(x/D = 6\) to 11) increases along the narrowing duct (Fig. 11d-ii).

In contrast to the upper branch, the flow evolution inside the straight lower branch is inherently simple. The transient shock wave maintains its planar shape up to the trailing edge of the centre body at \(x/D = 12\) (Fig. 11a–d). The shock velocity barely changes in the constant cross-sectional duct. Hence, the shock wave propagates at almost the same speed as the undivided incident shock wave (\(M_\mathrm{s} = 1.61\)). Therefore, the pressure and Mach number distributions are almost uniform upstream of the transient shock in the lower branch of the divider (Fig. 11a–c).

Fig. 12
figure 12

X–t diagram of both divider branches based on CFD simulations for divider I–III showing the shock propagation in the a entire divider and in its b inlet and c outlet sections

Full size image

While the lack of noticeable shock attenuation in the lower branch is expected for the numerical results due to neglection of viscous effects, the nearly constant shock velocity in the experiments may be surprising. Hence, the expected shock attenuation throughout the lower branch is estimated based on the literature on shock attenuation inside shock tubes. Non-ideal effects resulting in attenuation of the incident shock wave of shock tubes are of significant interest in particular for chemical kinetic studies [44, 45]. These non-ideal effects mainly depend on the boundary layer, experimental conditions, and the shock tube geometry, particularly the inner diameter of the tube [46,47,48]. To estimate the shock attenuation for the divider, the recently determined empirical relations given in [47] are used. Nativel [47] measured the shock attenuation across a range of pressures and incident shock Mach numbers for four different shock tube facilities. The driven section diameters of the shock tubes varied from 5 to 16 cm. For each diameter, an empirical relation was given based on a least-squares method to obtain a linear fit. Accordingly, the empirical relations for different shock tube diameters given in [47] are extrapolated to the divider’s lower branch channel width. The resulting shock attenuation for the lower branch is

$$\begin{aligned} \text {attenuation } ( \%/m) = 0.73+ 4.12\ (p_1^{-0.14} \sqrt{M_\mathrm{s}}), \end{aligned}$$

where \(p_1\) is the shock upstream pressure. For the considered configuration in this study, the shock attenuation is 7.9%/m. Taking into account the relatively short length of the divider, the shock wave is estimated to attenuate by only 0.8% throughout the divider’s lower branch.

The upper and lower separated shock waves reach the divider outlet at different times, as shown in Fig. 11d. The shock wave in the upper branch reaches the divider outlet at a later time due to a number of reasons. The path length of the upper branch is higher than that of the lower branch. Furthermore, the upper shock wave propagates at a lower average velocity. This is mainly due to two different mechanisms. Firstly, the divergent section of the divider weakens the shock wave due to diffraction. Secondly, the reflection of the shock wave from the divider walls, in particular at the turning point of the divider, weakens the shock wave. However, the shock wave in the lower branch propagates at a nearly constant velocity along the shortest path from the divider inlet to the outlet through a straight duct. Therefore, the shock wave from the lower branch reaches the divider exit first.

When the shock wave in the lower branch reaches the end of the centre body, it encounters an increase in duct width from D/2 to D. The area expansion results in diffraction of the shock wave (Fig. 11d). As shown in Fig. 11d-iv, a recirculation region occurs at the corner of the centre body (\(x/D = 12\)). The resulting vortex separates from the corner at a later time and propagates further downstream while interacting with the following shock wave from the upper branch (Fig. 13d). These results show the dynamic evolution of the wave patterns inside the divider and the capability of the divider to separate shock waves into multiple, weaker, consecutive shock waves.

3.3 Impact of divider width ratio on shock separation

The impact of the divider width ratio on the separation of the shock waves is investigated in this section. The divider width ratio is defined as the ratio of the upper branch maximum channel width to the lower branch. The width ratios of the investigated dividers are 1, 2, and 4 for dividers I, II, and III, respectively.

An x–t diagram, based on the CFD simulations, is presented in Fig. 12a. The results show a nearly identical shock velocity in the lower branch for all three dividers, as all x–t lines for the lower branches lie on top of each other. As presented in a close-up view in Fig. 12b, the separated shock wave in the lower branch propagates with nearly the same velocity as the initial incident shock wave upstream of the divider. However, the velocity of the shock in the upper branches varies between the dividers. As shown in Fig. 12c, the shock velocity decreases as the width ratio increases. Hence, the shock separation distance \(\xi \) at divider outlet increases with a larger width ratio, as indicated in Fig. 12c. This is to be expected, as a larger width ratio results in further deceleration of the shock velocity due to diffraction.

3.4 Recombination of the separated shocks

A straight exit duct is attached to the outlet of the dividers to investigate the further evolution of the separated shock waves and their interaction with each other (Fig. 13). As shown in Fig. 12c, for all dividers, the velocity of the shock wave from the upper branch increases continuously after the shock wave is transmitted to the exit duct. The flow inside the exit duct is mainly driven by two consecutive shock waves propagating in the same direction. The leading shock increases the speed of sound and sets the flow into motion as it propagates through the air, initially at rest. Therefore, the subsequent shock wave velocity increases gradually. Given a long enough exit duct, two shock waves travelling in the same direction must collide [49].

Fig. 13
figure 13

Numerical schlieren images showing the interaction of the separated shock wave after exiting the divider for dividers I–III at two different times

Full size image
Fig. 14
figure 14

Total pressure at the divider outlet at \(x/D = 12.5\) showing two separated shock waves

Full size image

Figure 13 shows numerical schlieren images of the three dividers and their exit ducts. The images in Fig. 13a–c show the very moment the shock wave from the upper branch reaches the trailing edge of the divider centre body. The shock waves from the lower branch are already diffracting into the exit ducts. The images in Fig. 13a–c clearly show that a larger width ratio results in a larger separation of the shock waves at the divider outlet. The separation distance \(\xi \) between the shock waves at this stage for divider III is 2.6 times larger than divider I. The shock wave from the lower branch propagates in two different directions after it is diffracted around the trailing edge of the divider centre part. A part of the shock wave is transmitted to the exit duct. This part of the shock propagates as a strong Mach reflection in downstream direction, whereas the remaining part of the shock propagates back into the upper branch.

The schlieren images in Fig. 13d–f show the flow evolution for a short time later for all three dividers. As shown in Fig. 13d for divider I, the shock wave from the upper branch has diffracted around the corner. A part of this shock propagates back into the lower branch. Another part of this shock wave is reflected from the lower wall of the exit duct. The remaining part of the shock wave propagates downstream along the exit duct. In Fig. 13d, the part which is moving downstream is behind the leading shock wave from the lower branch. In contrast, in the case of divider II the rearward shock wave in the exit duct is separated from the preceding shock (Fig. 13e). In Fig. 13f, the shock waves in the exit duct of divider III are separated by a larger distance. In addition, the schlieren images show that only a part of the lower and upper branch shock waves is transmitted to the exit duct, while the rest propagates back towards the divider inlet. Furthermore, the separation of the shock waves is larger and lasts longer for larger divider width ratio.

3.5 Temporal redistribution of the incident shock wave energy

The temporal redistribution of the energy of a single shock wave is further examined by analysing the total pressure. Figure 14 presents the total pressure in the exit duct for all dividers based on the CFD simulations. For the sake of comparability, the time is set to zero when the first shock wave reaches the sensor position. As shown in Fig. 13, the diffracted shock from the lower branch reaches the exit duct and therefore the pressure sensor, first. In addition to the pressure obtained from CFD simulations, the experimentally measured pressure for divider II is shown in Fig. 14. The signal is filtered with a 50-Hz–150-kHz band-pass filter to remove noise. The delay of nearly 9 µs between the shock arrival and the plateau pressure corresponds to the rise time of the sensor [32, 50]. The measured plateau pressure after the passage of the first shock wave is 8% less than the CFD predicted pressure. The numerically overestimated pressure complies with our observations in Sect. 3.1 that the numerical shock strengths are slightly overestimated compared to the experiments. Furthermore, there is very good agreement between CFD and experiments for the second shock wave and the proceeding oscillations. The oscillations of the pressure are mainly caused by multiple reflections of the shock wave from the upper and lower walls of the exit duct. The oscillations converged towards a plateau pressure of nearly 3.7 bar for divider II. These results support the ability of the high-frequency total pressure probe to measure the transient total pressure.

In Fig. 14, the total pressure behind a \(M_\mathrm{s} = 1.61\) shock wave is given as a reference based on normal shock expression. It allows for estimation of the entropy generation produced by the divider. The plateau total pressure behind the secondary shock is 5, 12, and 20% smaller than the reference configuration, without a divider, for dividers I to III. Hence, the pressure loss increases with increasing separation of the shock waves. These results show that any attempt to reduce impulsive loading on the turbines through the use of shock divider will come at the cost of total pressure loss.

As shown in Fig. 14, the time \(\tau \) is taken as a measure for the temporal redistribution of the shock energy. The time \(\tau \) is the interval between the moment when the pressure first increases and the moment at which the maximum pressure is reached (Fig. 14). For the configuration without a divider \(\tau = 0\), as the maximum pressure occurs directly behind the incident shock wave. However, \(\tau > 0\) when the divider separates the shock waves and increases with larger divider width ratios. The time \(\tau \) based on the CFD simulations is 15, 19, and 27 \(\upmu \)s for the dividers I, II, and III, respectively. A good agreement for \(\tau \) between CFD and experiment is evident, as \(\tau _\mathrm{EXP}\) = 25 \(\upmu \)s and \(\tau _\mathrm{CFD}\) = 19 \(\upmu \)s for divider II (Fig. 14). The discrepancy between these quantities is mainly attributed to the sensor rise time. The temporal evolution of the total pressure at the divider exit demonstrates the ability of the divider to redistribute the energy of a shock wave. Furthermore, the comparison between the dividers in Fig. 14 shows that the temporal redistribution of the energy can be further increased through the divider width ratio, but at the expense of total pressure.

4 Conclusion

The redistribution of transient incident shocks into multiple shock waves is investigated by utilising a device termed a shock divider. The shock divider consists of two pathways with different path lengths and channel widths. A transient shock wave is generated using an open-end shock tube. The divider is connected to the end of the shock tube. As the shock wave enters the divider, it separates into two primary shock waves propagating through the two divider branches. The path length of the upper branch is higher than that of the lower branch, allowing the separated shock waves to exit the divider one after each other. The separated shock waves transmit to a single duct at the downstream end of the divider.

Numerical and experimental methods are used to evaluate the flow evolution inside the divider for an initial shock wave with a Mach number of \(M_\mathrm{s} = 1.61\). A second-order dimensional split finite volume MUSCL-scheme (CFD) is used to solve the compressible Euler equations for a 2D Riemann problem. Furthermore, low-cost geometrical shock dynamics (GSD) simulations are used to estimate its eligibility to predict the wave dynamics inside the divider. High-speed schlieren images of the divider flow are used to validate the numerical approaches. In addition, a high-frequency pressure probe is used to measure the transient total pressure at the divider exit.

The comparison between the CFD results with experimental schlieren and total pressure measurements shows very good agreement. Moreover, the GSD approach predicts the propagation of the leading shock waves inside the divider with high accuracy. The CFD approach is shown to be more accurate while requiring relatively low computational resources. However, the numerical approaches slightly overpredict the shock velocity in the upper divider branch, as they do not account for the boundary layer growth due to the lack of viscosity. The computational time for the simulation of the divider flow is in the order of 5 minutes for the GSD and 10 minutes for the CFD using a 2.9-GHz i7 work station. We conclude that both schemes are appropriate for efficient design studies due to the very low computational cost.

Supersonic regions, quasi-steady shock waves, and Prandtl–Meyer expansion waves occur inside the divider, even though the initial flow behind a \(M_\mathrm{s} = 1.61\) planar shock is subsonic. The flow transitions from subsonic to supersonic as the shock wave diffracts at the convex corner of the divider inlet. Induced by the shock diffraction, an upstream-propagating expansion wave increases the flow velocity. Furthermore, the planar incident shock wave separates into a planar shock and a slightly curved shock in the lower and upper branches, respectively. While the shock wave in the lower branch propagates at a nearly constant velocity, the leading shock in the upper branch is exposed to various loss mechanisms. A relatively strong Mach reflection occurs at the turning point of the upper branch. Its triple point reflects a number of times from the inner walls of the divider before it transmits to the exit duct downstream of the divider. Furthermore, the gradual diffraction of the incident shock within the diverging section decelerates the shock wave further. Therefore, the shock wave from the upper branch reaches the divider exit after the lower branch.

The impact of the divider channel width ratio on the separation of the divided shock waves is studied. For this purpose, three different divider width ratios of 1, 2, and 4 are investigated. The time interval between the arrival time of the two shock waves at the divider exit increases with increasing width ratio. This is mainly caused by the diffraction of the shock wave, as further diverging of the upper channel results in further deceleration of the shock wave. Therefore, the leading shock wave from the upper channel exits the tube at a later time. Consequently, the separation of the shock waves downstream of the divider lasts longer with increasing divider width ratio.

Total pressure is measured at the exit of the divider using a high-frequency pressure probe. The comparison with the numerical results shows the capability of the probe for measurement of highly transient flow. Furthermore, the temporal redistribution of the initial shock wave energy is investigated by analysing the temporal evolution of the total pressure. It is shown that the temporal redistribution of energy increases with increasing divider width ratio, but at the expense of total pressure.

The divider may be an approach for mitigating the load on the turbine in future PDC–turbine applications. The results of the current study show promise towards the ability of the divider to redistribute the energy of a shock wave. However, future studies are needed to minimise the entropy generation induced by the divider.