On Shock Propagation through Double-Bend Ducts by Entropy-Generation-Based Artificial Viscosity Method

Shock-wave propagation through obstacles or internal ducts involves complex shock dynamics, shock-wave shear layer interactions and shock-wave boundary layer interactions arising from the associated diffraction phenomenon. This work addresses the applicability and effectiveness of the high-order numerical scheme for such complex viscous compressible flows. An explicit Discontinuous Spectral Element Method (DSEM) equipped with entropy-generation-based artificial viscosity method was used to solve compressible Navier–Stokes system of equations for this purpose. The shock-dynamics and viscous interactions associated with a planar moving shock-wave through a double-bend duct were resolved by two-dimensional numerical simulations. The shock-wave diffraction patterns, the large-scale structures of the shock-wave-turbulence interactions, agree very well with previous experimental findings. For shock-wave Mach number Ms=1.3466 and reference Reynolds number Ref=106, the predicted pressure signal at the exit section of the duct is in accordance with the literature. The attenuation in terms of overpressure for Ms=1.53 is found to be ≈0.51. Furthermore, the effect of reference Reynolds number is studied to address the importance of viscous interactions. The shock-shear layer and shock-boundary layer dynamics strongly depend on the Ref while the principal shock-wave patterns are generally independent of Ref.


Introduction
Shock/blast wave propagation involves complex wave interactions with the media and surface boundaries owing to several phenomena such as shock reflection, shock focusing, shock diffraction and shock-turbulence interaction. Understanding of these phenomena is crucial for a wide range of engineering applications in bio-medicine, disaster management, detonation, mining, aviation/transport industry and others. Knowledge of such complex dynamics is integral part of the design and optimization of devices for shock-wave lithotripsy, shock/blast-wave attenuation, suppression of tunnel sonic boom, etc. The existence of a wide range of flow scales together with unsteady flow discontinuities and the coupled shock-turbulence interactions make it a challenging and active research field. Numerical prediction is a cost effective option compared to its experimental counterpart accounting for several restrictions. However, resolving these unsteady dynamics by numerical techniques demands development of high-fidelity numerical tools. The inspiration behind this work is particularly related to the applicability and effectivity of the high-order numerical tools resolving flow dynamics associated with shock-wave propagation and attenuation.
The literature shows various approaches to attenuate shock-waves, e.g., foams, textiles, porous materials, granular filters, metallic grids, perforated plates/walls, rigid barriers, branched/bend duct, duct with rough walls, etc., as mentioned in [1]. In this regard, a comprehensive review of various methods to attenuate shock/blast waves was reported by Igra et al. [2] addressing both experimental and numerical approaches. Essentially, shock-wave attenuation by geometrical means such as rigid barriers or sudden changes in the flow geometries is governed by compression, rarefaction regions arising from shock diffraction and intense shock-turbulence interactions [1][2][3][4][5][6][7][8][9][10][11]. These involve multiple wave interactions, reflecting from surrounding boundaries and complex shock-shock, shock-vortex and shock-boundary layer interactions. The dynamics may involve transition between regular reflection (RR) and Mach reflection (MR) with triple point. However, depending upon the geometrical complexities and flow conditions, there exists numerous type of irregular reflections [12] and the evolution may consists of complex shock systems with one or more triple point [13,14].
Owing to the enormous advent of computational power and numerical methodologies in recent time, detailed flow dynamics can be retrieved via high-resolution numerical approach compared to experimental measurements. The numerical prediction of complex viscous shocked flow requires stable, high-order numerical schemes to capture shocks as well wide range of flow scales. Often, a competing opposing goal is to be optimized by numerical tools to sustain the accuracy and stability of the scheme. For example, high-order numerical schemes need the addition of dissipation to reduce inherent oscillatory behavior near discontinuities and at the same time one has to reduce numerical dissipation for resolving fine turbulent scales. Successful utilization of explicit high-order Weighted Essentially Non-oscillatory (WENO) has been reported to resolve unsteady shock dynamics related to shock diffraction, shock-reflection and shock focusing [1,5,[13][14][15][16][17] problems. On the other hand, high-order Discontinuous Galerkin Methods (DGM) or high-order Discontinuous Spectral Element Methods (DSEM) are comparatively more flexible dealing with complex flow geometries exploiting unstructured mesh-based formulation. Discontinuous spectral element methods (DSEM) can be considered as a nodal DG method [18][19][20][21][22], belonging to the class of weighted residual methods having Dirac delta function as a test function. Spectral element based studies cover a wide range of applications areas related to smooth, non-smooth and complex fluid flow problems [23][24][25][26][27][28][29][30][31][32][33]. Nevertheless, dealing with Gibb's oscillations for non-smooth fluid flow problems requires implantation of suitable slope limiters or artificial viscosity based approaches. In this regard, Chaudhuri et al. [34] reported the stable shock capturing capabilities of explicit high-order DSEM in conjunction with entropy-generation-based artificial viscosity (AV) method. Although inviscid simulations quite accurately predict the shock wave diffraction wave patterns over different geometrical shapes, it is evident that viscous effects are important for resolving long-time behavior of shock-vortex evolution, shock-shear layer and shock-boundary layer interactions [35][36][37][38][39][40][41]. These issues are also discussed in the recent works [42,43] using high-order DSEM with AV method.
In their works, Igra et al. studied shock propagation and diffraction in ducts with cavity [44,45] and in branched ducts [46]. The effect of attenuation of these configurations compared to uniform cross-section duct is arising from the multiple shock-wave reflections during the propagation of the moving shock. Numerical predictions solving inviscid Euler equations of these configurations agree well with the experimental findings related to general wave patterns and pressure profiles. Subsequently, Biamino et al. investigated the effect of the length of the branched segment revealing the fact that it is questionable to achieve any protection at the end of a branched duct [47]. Nevertheless, favorable shock attenuation by employing abrupt changes in tunnel geometries was reported by Igra et al. [48], highlighting the shock attenuation associated with smooth-walled, rough-walled double-bend ducts. In addition, a detailed experimental findings with respect to varying volume of the double-bend ducts are presented in their work.
The motivation of the present work is to resolve such complex flow dynamics with high-fidelity viscous simulations, which have not been reported in the literature to the best of our knowledge. We aim to present an analysis of shock propagation through the double-bend duct addressing flow behavior with varying reference flow Reynolds number Re f . The M s s (shock-wave Mach numbers) are chosen similar to those presented in the experimental work of Igra et al. [48]. A high-order artificial viscosity based DSEM [34] is used for this purpose. The flow analysis sheds light on numerical perspective about the robustness and accuracy of the scheme, as well as the shock diffraction and shock attenuation aspects in double-bend ducts. The paper is organized as follows. A brief description of the governing equations followed by the numerical procedure is given in Section 2. The problem setup is illustrated in the Section 3. The flow analysis and the comparison of the results with the experimental are discussed in Section 4. Finally, conclusions are drawn in Section 5.

Governing Equations & Numerical Approach
In this section, we briefly present the governing equations and numerical approach. Detail of the numerical procedure is reported in [34] and is not repeated here for brevity.

Governing Equations
Gaseous system involving moving shock-wave is governed by the compressible Navier-Stokes (NS) system of equations. The non-dimensional form of the governing equations with artificial transfer coefficients is given by, and where U is the conservative solution vector, and F a and F v are the inviscid and viscous flux, respectively. ρ is the density, v is the velocity vector, and E t is the total internal energy. p is the static pressure and T is the temperature. γ is the ratio of specific heats. δ is the Kronecker delta tensor.
Here, M f , Re f and Pr f are the reference Mach number, Reynolds number, and Prandtl number, respectively. Assuming Stokes' hypothesis with zero bulk viscosity, the viscous shear stress tensor can be (∇v) T + ∇v is the symmetric part of the velocity gradient tensor. The superscript "T" designates a transpose. Similarly, κ eff = κ h Re f + κ is the effective thermal conductivity. Note that µ h and κ h are yet to be determined artificial transfer coefficients. For NS system of equations, µ = κ = 1.
In this study, γ = 1.4, M f = 1 and Pr f = 0.72 were prescribed. The ideal gas equation of state,

Numerical Approach
We used the staggered Chebyshev collocation method to approximate the compressible Navier-Stokes system of equations together with the explicit marching algorithm in time.
In three-dimensional nodal collocation formulation of DSEM, the physical domain is subdivided with hexahedral physical elements. An iso-parametric transformation is then used to map each physical element to a unit cubic computational element in the computational domain. This reduces to quadrilateral physical elements and square computational elements in the two-dimensional formulation. The solution vector is collocated at the Chebyshev-Gauss quadrature points and the fluxes are collocated at the Chebyshev-Lobatto quadrature points. The detail of the discretization procedure can be found in [19][20][21][22]. Several methodologies can be adopted to suppress the Gibb's oscillation near the flow discontinuities. For this, a cost effective approach is the usage of AV compared to slope limiters. A short review of different AV approaches in conjunction with high-order methods is reported in [34]. The basic idea of an AV method is to explicitly add even-order dissipation term to stabilize the numerical scheme. Nevertheless, this method requires to define and assign arbitrary model constants as flow-dependent tuning parameters to achieve optimal solution. Additionally, time step restriction arising from artificial viscous terms should also be taken care of. We used entropy-generation-based AV coefficients (µ h and κ h ) and set a suitable upper bound of these coefficients to judicially address these issues. Here, we present, again, a very brief account of the method of calculation of µ h and κ h to facilitate the understanding the essential features of the numerical procedure.
We consider the non-negative entropy generation terms of the entropy conservation equation in-order to scale the AV coefficients. The artificial momentum and thermal conductivity are scaled with the viscous and conductive, entropy generating terms Φ and Γ, respectively. In non-dimensional form, this yields the following expression for the artificial viscosity coefficients: ln p ρ γ and ρ s is the spatial average of ρ s. Here, ||ρ s − ρ s|| ∞ is the globally computed supremum based on the global average entropy. The artificial viscosity coefficients defined above ensure the positivity (thus dissipative behavior) and µ h and κ h scale with the grid spacing and vanish as ∆h → 0.
We further used a shock sensor θ [49] to control the artificial viscosity coefficients, so that the modified coefficients can be expressed as µ h θH(−∇ · v), and κ h θ. This reduces the artificial dissipation in rotation-dominated regions of the flow-field. The purpose of the Heaviside function H(−∇ · v) is to ensure that the dissipation is small in regions of isentropic expansion fans and contact discontinuities. The coefficients are kept below an upper bound so that the inviscid time step ∆t inv = CFL inv ∆h/(|v| + √ T), is smaller than the viscous time step. From this constraint, the upper bound of the artificial coefficients becomes µ h,m = C m ρ∆h(|v| + √ T), where C m ∝ CFL vis /CFL inv represents another model parameter. Note that κ h is also kept bounded by the µ h,m . Interested readers are suggested to find the detailed description of the overall scheme in [34].

Problem Setup
The geometry of the double-bend duct configuration is taken similar to that presented in the experimental and inviscid study of Igra et al. [48]. Figure 1 shows the setup for the non-dimensional physical domain in the two-dimensional (2D) x − y plane. A shock wave with a suitable Mach number M s is allowed to pass through the duct geometry. The initial shock location is at x s = 0.75 and the conditions are prescribed using the Rankine-Hugoniot relations for stagnant state (1) and shocked gas state (2). The left and right boundaries are set by the initial states. The top and bottom boundaries are assigned with no-slip wall conditions. The reference state φ f , is taken as the stagnant state conditions φ 1 , to satisfy M f = 1 and the reference length is taken as the height of the entrance of the double-bend duct in the experiment [48] L f = H. The length of volume inside the double-bend duct is L and the chosen geometry satisfies L/H = 16 (see Figure 1). Table 1 summarizes the flow properties of stagnant state (1) and shocked gas state (2) for respective values of M s .
In total, 68,000 P3 (fourth order) elements with 1.088 × 10 6 degrees of freedom were considered in the domain. In the present study, the value of the CFL number was taken as 0.9 and the model constants for artificial viscosity were set as: (C µ , C κ , C m ) = (0.5, 0.25, 0.15), unless stated otherwise. Additionally, to enhance the stability of the method, we used adaptive spectral filter with filter order of 16 near the regions where max(µ h , κ h ) attains its local upper bound and a filter order of 64 otherwise [34]. The simulations were performed on Supermicro X9DRT compute nodes (dual Intel E5-2670). A typical test case utilized 18 nodes consuming about 3000 CPU hours.

Results and Discussion
We simulated several cases by setting M s = 1.3466 and M s = 1.53 with Re f = 10 3 , 10 4 , 10 5 and 10 6 . Simulations were performed until non-dimensional time t = 15. We first discuss the basic flow features for the case (M s = 1.53, Re f = 10 6 ) and then subsequently present the comparison for several other cases towards the shock attenuation aspect in double-bend ducts.

General Flow Evolution
The moving shock-wave after the passage through the inlet section gets diffracted over the top left corner of the domain. It can be realized that the viscous interactions are important to account the shock-wave diffraction over such 90°convex corners. This issue has been discussed with respect to DSEM-AV based numerical scheme in [42]. The flow evolution here, in the double-bend duct, remains similar to that of 90°convex corner until the upward moving reflected shock-wave from the bottom-wall interacts with the diffracted flow field around the corner (see density contour at t = 7.2 in Figure 1). In addition, the shear layer interacts with the boundary layer at the top wall through expansion and secondary shock-wave. The density contours in Figure 2 illustrate the subsequent complex flow evolution with multiple transverse wave interactions. The shear layer gets intensely perturbed by the transverse reflected shock-waves, yielding a complex shock-vortex interaction. The unstable shear layer is associated with vortex shedding and vortices become deformed and convected forward in the flow-field. This flow dynamics is also illustrated by the numerical interferograms computed from the density field. The interferograms were computed using the expression: I = β 1 + cos 2π ρ − ρ 1 ∆ρ . Here, ∆ρ = (ρ max − ρ min )/N, and N is the number of interferential fringes. We chose β = 1, as recommended in [50], and N = 10 to estimate I. The overall shock-wave patterns, shock-boundary layer interaction, shock-shear layer interactions and the secondary viscous vortex interactions at the left vertical wall interacting at the convex corner are clearly visible in these interferograms. Note that the simulation resolved the RR → MR transition of the shock-wave reflection at the bottom wall (see contours at t = 8.

Comparison with Experimental Results
We present the time snaps of the flow evolution for the case (M s = 1.53, Re f = 10 6 ), to compare with the experimental work of Igra et al. [48]. The numerical Schlieren pictures were computed from the density field using the expression: S = β exp −λ |∇ρ| |∇ρ| max . We chose β = 1 and λ = 15, as recommended in [50], to estimate S. The Schlieren pictures shown in Figure 4 can be compared with the experimental shadowgraphs [48]. The shock-wave patterns from the present simulation are in excellent agreement with the experiment. These resemble to the experimental Shadowgraphs at 220 µs, 260 µs, 300 µs, 380 µs, 400 µs and 420 µs (see Figure 15k,m,o,s,t,u of Igra et al. [48]). The predicted large scale flow structures in rotational dominated regions, and the shear layer dynamics are in very good agreement with the experiment. The global structures of the shock-wave boundary layer interactions at the top wall is also in accordance with the experimental structures. Note that the initial location of the shock-wave for the experimental findings are not known. The flow-field time snaps of the simulation are not exactly at the same instant. The numerical tool equipped with entropy generation based AV together with adaptive spectral filter clearly provides stable and accurate prediction of the flow dynamics. The strain-enstrophy angle can be defined as Ψ = tan −1 S · S A · A .
where, S = 1 2 (∇v) T + ∇v and A = 1 2 ∇v − (∇v) T . The contours of Ψ illustrated in Figure 5 at the interaction zone of window section (0.8 ≤ x ≤ 9, −2.5 ≤ y ≤ 1). The patches of Ψ contours below 45°signify the rotation dominated regions in the domain. For a 2D flow field, the two eigenvalues of the velocity gradient tensor (∇v) T can be expressed as λ 1 = λ r1 + iλ i1 and λ 2 = λ r2 + iλ i2 , where i 2 = −1. The instantaneous contours of these eigenvalues are shown in Figure 6 at different time instants. The non-zero high values of λ i1 are associated with the vortex core regions of shear layer and in the boundary layer regions. From the contours of λ r1 , it can be seen that the high values are adhering to the outer edges of vortices in the shear layer. In addition, in the boundary layer regions, its values are lower where λ i1 assumes higher values. On the other hand, λ r2 shows larger values near high dilatation shock regions. These characteristics are similar to that reported in [42]. Nevertheless, it can be realized that three-dimensional (3D) simulations are necessary to further resolve the fine turbulent scales of the flow evolution. An insightful attempt resolving such 3D structures was presented by Chaudhuri et al. [17] via large eddy simulation.

Shock-Wave Attenuation and Effect of Re f
In the review by Igra et al. [2], it is reported that the height of the physical domain is H = 60 mm. The experimental condition of the stagnant state is P 1 = 0.987 bar, T 1 = 23.4°C for M s = 1.3466 and P 1 = 0.982 bar, T 1 = 23.7°C for M s = 1.53 [48]. These give rise to a reference Reynolds number in the range of ≈ 10 6 , based on the length parameter H and the speed of sound at stagnant state. We performed simulations for both M s = 1.3466 and M s = 1.53 with varying Re f from 10 3 to 10 6 . For Re f = 10 3 , the model constants for artificial viscosity were set as: (C µ , C κ , C m ) = (0.2, 0.2, 0.1) and the model constants for remaining cases were assigned to the values mentioned in Section 3. Figure 7 shows the numerical Schlieren pictures for various test cases at t = 14. At this time instant, for M s = 1.3466, the incident shock-wave reaches near the right bottom corner of the duct while, for M s = 1.53, it enters the exit section of the domain, after experiencing the second diffraction at this corner of the domain. The basic shock-wave patterns remain unaffected with the variation of Re f . However, the incident shock wave is lagging behind relatively for Re f = 10 3 compared to the higher Re f s for both M s s. This effect is larger for higher M s . For each M s , the shock-wave interaction with the main vortex and large scale structures of the shear layer are very much identical for Re f = 10 6 and 10 5 . Evidently, dominant viscous effects on the shear layer structures and vortex shedding can be clearly seen for lower Re f s. The shear layer becomes stable for Re f = 10 3 . The entropy generation based AV methodology used in this study is designed to add optimal dissipations to get stable solution. Figure 8 depicts the representative contours of µ h . Note that the intense shock regions of wave patterns are highlighted by the non-zero values of µ h . The contours of the ratio µ h /µ h,m also show that only few small patches of region are associated with values close to unity. The entropy generation can be expresses as is related to thermal contribution (see details of the entropy transport equation in [34]). Figure 9 illustrates the time evolution of area weighted average (defined as ψ = ψdA dA ) of these generation terms. This is consistent with the evolution of estimated µ h and κ h . Evidently, the values of µ h remain higher than κ h . In addition, note that C µ = 0.5, is set higher than C κ = 0.25 for these cases. It is interesting to compare the entropy generation (G S,eff ), estimated with effective transport coefficients µ eff , κ eff (see Section 2.1), and the entropy generation (G S,h ), estimated by artificial coefficients µ h Re f , κ h Re f . Figure 10 shows the contours of G S,h /G S,eff . It is clear from these contours that the higher values of the ratio are associated with the shock dominated regions and these corroborate with the effectiveness of the shock sensor and the dilatation based Heaviside function for estimating AV coefficients mentioned in Section 2.2.  We further analyzed the pressure signals at the bottom wall of the double-bend duct for different cases to highlight the attenuation aspect of the flow configuration. The pressure profiles at the bottom wall for both M s = 1.3466 and 1.53 with Re f = 10 6 are shown in Figure 12. Note that L/H = 16 for the double-bend duct considered in this study. For M s = 1.3466, the value of P/P 1 predicted from the present simulation near the exit section of the domain is ≈ 1.5. This value closely matches with the findings of Igra et al. [48] for this configuration. On the other hand, we observed P/P 1 ≈ 1.8 for M s = 1.53. The over pressure can be defined as Π = P − P 1 P 2 − P 1 . The outcomes yield the over pressure Π ≈ 0.53 for M s = 1.3466 and Π ≈ 0.51 for M s = 1.53. Note that, for M s = 1.53, the pressure signal also shows the signature of the diffraction from the right bottom corner of the domain. The reflected shock-wave from that corner interacts with bottom boundary and this is clearly seen in Figure 12 at t ≈ 15 with a second peak of P/P 1 ≈ 2. In Figure 13, the effect of low Re f is clearly visible. For Re f = 10 3 , we observed a shock-wave retardation apart from attenuation. The attenuation features, however, remain unaffected with low Re f when compared with the higher Re f . The transverse wave reflection at upstream locations are also evident in these pressure signals. Note that, there exists an additional hump in pressure signal for the case of M s = 1.3466. On the other hand, this is similar for M s = 1.53 at the early stage, while at the later stage, one can notice other humps in the pressure signals (see profiles at t = 14, 15). Evidently, these are in accordance with the contours presented in Figure 7.

Conclusions
In this work, we studied the planar shock-wave propagation through a double-bend duct having L/H = 16 via numerical simulations. The physics of the flow involves complex shock wave dynamics with shock-wave diffraction and multiple wave interactions together with the shock-shear layer and shock-boundary layer interactions. A high-order DSEM equipped with entropy-generation-based AV method was used to solve Navier-Stokes system of equations for this purpose. The flow evolution for the case (M s = 1.53, Re f = 10 6 ) was found to be in excellent agreement with the previous experimental findings of the literature. In addition, for the case (M s = 1.3466, Re f = 10 6 ), the predicted pressure signal agrees very well with the literature. Additionally, several simulations were performed for two M s s with varying reference Re f . The principal shock-wave patterns were found to be generally independent of Re f . On the other hand, the shock-shear layer and shock-boundary layer dynamics strongly depend on the Re f . At Re f = 10 3 , we observed marginal retardation of the incident shock-wave, owing to the important role of viscous effects. The results show the applicability and effectiveness of the AV based methodology in resolving complex flow physics associated with the shock propagation and attenuation through double-bend ducts. Three dimensional Detached Eddy Simulation (DES), Delayed Detached Eddy Simulation (DDES) or Large Eddy Simulation (LES) could be performed in the future to resolve later stage fine turbulent flow scales, observed in the previous experiments.
Funding: This research received no external funding.