Solutions of the converging and diverging shock problem in a medium with varying density

The Guderley problem is a well-known hydrodynamic problem of a strong shock propagating radially in an ideal gas medium. The shock originates at infinity and collapses to a central point or an axis, where it is reflected and diverges toward infinity.^1–9 The problem was first formulated with an initially uniform density profile and was analyzed comprehensively and solved for this case by Lazarus and Richtmyer.^2,3 The theory of converging and diverging compressible shocks is an active research area, where, for example, the Guderley problem was generalized to include non-ideal equations of state,^10–15 magneto-hydrodynamics,^16,17 and gravitation.¹⁸ Recently, a general solution for continuous converging waves was developed.¹⁹ 

Imploding shocks are often invoked as a possible ignition mechanism in reactive material. These include chemical compounds (see Ref. 20 for a review on laboratory experiments), thermonuclear in inertial confinement fusion experiments,^21,22 and detonation in white dwarfs in astrophysics.^23,24 A fundamental feature in such ignition scenarios is that both the converging and diverging shocks may lead to ignition, since the particulars depend on a competition between the hydrodynamical time and the heating time due to the exothermal burning processes. Interestingly, the critical strength of the shock required to reach ignition is usually determined by the diverging shock. While this shock is typically not strong (as opposed to the earlier converging shock), the post-shock material in this stage has higher densities and temperatures, having being shocked twice; see the analytical and numerical analysis by Ref. 23. Hence, an accurate calculation of the entire converging–diverging shock sequence is crucial to understand (and designing) an ignition candidate system.

In Ref. 25, we presented generalized solutions for the converging shock for a power law initial density profile $ρ ( r ) ∼ r μ$⁠, with a general spatial exponent μ (see also Refs. 6, 7, 26, and 27). In this work, we extend this generalization of the initial density profile, to the reflected diverging shock solution. By combining the solutions for the converging and diverging flows, we obtain a solution of the complete Guderley problem for an initial power law density profile.

The Guderley problem is a self-similar problem of the second kind, in which the similarity is found by requiring the solution to pass through a singularity, such as a sonic point. As a result, its solution is given in terms of similarity profiles as the hydrodynamic Euler equations are transformed into a system of ordinary differential equations (ODEs). An extension of the problem to a general power law initial density is possible, since this generalization does not introduce additional dimensional parameters to the problem.^28–31 The equations describing the similarity profiles contain two unknown parameters that need be determined, either approximately^32,33 using simple analytical methods, or exactly using more sophisticated numerical methods. An exact numerical method for the calculation of the first parameter, the similarity exponent λ, was developed for general μ, in our previous work.²⁵ The second parameter, which is the diverging shock constant B, was calculated in the literature for μ = 0 by Lazarus et al.^2–4 and Ramsey et al.⁴ 

The Guderley problem can be used as a nontrivial test problem for hydrodynamic simulation codes.^{4,8,10,25,34–37} Using the generalized semi-analytic solutions, we define two types of test problems for hydrodynamic simulation codes, by properly defining the initial and boundary conditions. The tests represent the dynamics of the diverging shock as well as the dynamics of both the converging and diverging shock in the same simulation. We perform a detailed comparison between the simulations and the analytic solutions, for a variety of γ (the ideal gas adiabatic index) and μ.

In this work, we focus on the calculation of the diverging shock. We generalize the method of Lazarus³ and present a robust method for the calculation of the diverging shock constant B, for a general value of the initial density exponent μ. The structure of the article is as follows: In Sec. II, we review the generalized Guderley problem and its self-similar solutions for the converging and diverging shocks are presented. In Sec. III, a robust method for the calculation of the diverging shock constant and the resulting similarity profiles, for general values of μ, is developed. The solutions are calculated for a wide range of values of γ and μ and are compared to the results of previously published works when such exist. In Sec. IV, we define two types of test problems for hydrodynamic simulation codes, by properly defining the initial and boundary condition of a piston with velocity given from the analytic solution.³⁶ We perform detailed hydrodynamic simulations and compare the results to the similarity solutions. We conclude in Sec. V.

In this work, we consider the hydrodynamic problem of a shock wave that diverges to infinity in an ideal gas medium. This diverging shock is originated by a previously converging shock wave according to the Guderley problem, which has reached the origin and reflected from a central point (spherical symmetry) or an axis (cylindrical symmetry), as demonstrated in Figs. 1 and 2.

In this section, we describe a robust method for calculating the parameter B using iterative algorithms. We generalize the method of Lazarus^2,3 which dealt with the specific case of an initially homogeneous density (μ = 0), to a general power law density [Eq. (6)]. In Sec. III A, we describe the proper integration path for Eq. (16), which in turn gives a constraint that enables the calculation of the constant B. This method resembles the calculation of the similarity exponent λ, which is determined by a constraint on the integration path for the incoming shock, as it crosses a singular point.²⁵ 

The key to the full solution of the diverging shock is to determine the correct value of the parameter B, which we do iteratively. As described in Fig. 3 that depicts the integration path in the C–V plane, for a given trial value of B, the ODEs 16 are integrated from x = −1 to x = B, and then using the jump conditions, we obtain V ₂ and C₂ from Eqs. (25) and (27) that we denote as $( V a , C a )$⁠.

In the next step, the ODEs are integrated from the other direction (⁠ $x → ∞$⁠) toward x = B. This is done by switching the integration variable to w [Eq. (33)] and integrating Eqs. (34) and (35) from some point $w 0 ≪ 1$ with the initial values $( V , C ) = ( V 0 , 1 w 0 )$⁠, according to the asymptotic relations (30) and (39). In practice, we take $w 0 = 10 − 7$⁠. The integration is carried out until the value $V ( w ) = V a$ is reached. The value of $C ( V a )$ obtained at the end of the integration on w is denoted by C_b. If the correct value of B is used, the integration from both sides will meet at the same point in the C–V plane (Fig. 3), that is, we must have C_a = C_b.

In order to set the initial bracketing range $[ B min , B max ]$ such that $f ( B min ) · f ( B max ) < 0$⁠, we linearly scan the segment $[ 0 , 2 γ + 1 γ − 1 ]$ until we find $B min and B max$ which result in different signs of $f ( B )$⁠. After finding $B max$ and $B min$⁠, we use a standard bisection method to bracket the root B. Some caution is necessary, since $C a − C b$ may diverge for certain values of B, as shown in Fig. 4. If such divergence occurs, we reduce the scanning increment $Δ B$ and repeat the process from the last successful iteration, $B *$⁠, until once again we find an appropriate range of $[ B min , B max ]$⁠.

In Fig. 5, we plot the curve $C ( V )$ for several cases, as obtained from the numerical integration described in Fig. 3. To the best of our knowledge, numerical results for the values of B appear in the literature only for μ = 0. In this case, the values published by Lazarus³ are used as a benchmark to evaluate the algorithm—as was done, for example, in Ref. 4. Thus, we compare our results for the calculation of $B ( γ , μ = 0 , n )$ in Fig. 6, for a wide range of γ and for spherical (n = 3) and cylindrical symmetries (n = 2). Except for a few cases, the agreement is better than $10 − 5$⁠, mostly to all digits given in Ref. 3. In Fig. 7, we present the value $γ − 1 γ + 1 B$ as a function of γ for various values of μ. In Table I, we list the numerical values of B for $γ = 1.4 , 5 3$ for several values of μ. Finally, in Fig. 8, the self-similar profiles $V ( x ) , C ( x ) , and R ( x )$ are shown for several examples.

In this section, we compare the semi-analytical diverging shock solutions to numerical hydrodynamic simulations. As in Ref. 25, we employ a standard staggered-mesh artificial-viscosity one-dimensional Lagrangian code (for details, see Appendix D in Ref. 25).

We define two types of setups: (i) an initialization of a diverging shock flow that propagates outward and (ii) an initialization of an incoming shock flow that first converges to the origin and reflects outward as a diverging shock that propagates outward into the converging shock flow. The latter case tests the ability of the simulation to accurately describe both the focusing of the shock at the center and the generation of the reflected shock, which then propagates into the incoming shock profile. The former case is more simple, as it does not include the shock focusing and reflection, but is still non-trivial as it describes the propagation of the outgoing shock.

We define three test cases as follows:

1. $γ = 1.4 , μ = 0$—this is the canonical Guderley problem with a uniform initial density profile.

2. $γ = 2 , μ = 0.5$—a test case with a monotonically outward-increasing density profile. While the converging shock accelerates inward in such a profile, we expect the diverging shock to decelerate, as it advances into a denser material.

3. $γ = 1.2 , μ = − 0.8$—a test case with a monotonically decreasing density profile. In this case, the opposite occurs: The converging shock decelerates on its way inward, and the diverging shock accelerates, as it advances into an increasingly rarefied material.

Notably, a sharp discontinuity appears in the hydrodynamic profiles in all three cases (Figs. 9–11 and 15–17). As discussed in Ref. 25, these discontinuities exist over a single computational cell and were originated due to the presence of discontinuities in the initial hydrodynamic profiles.³⁸ As seen in the figures, they do not affect the overall stability and accuracy of the simulations. A detailed analysis of such initialization errors can be found in Ref. 4.

As in Ref. 25, we see that while there exists some variability in the accuracy of the simulations regarding the different physical quantities, we observe that a grid of 1000 computational cells will generally suffice to ensure an accuracy of a few percent or even less than one percent in all physical quantities at the end of the simulation. The spatial convergence rate (reduction in error w.r.t. the analytical solution as a function of increase in the grid size) is similar for all physical quantities in all simulations. These rates are all of order unity (as seen explicitly in Figs. 14 and 20), which is to be expected for an artificial viscosity numerical scheme, which is of first order accurate in the presence of shocks.

As in the converging shock setup defined in Refs. 25, 34, and 36, we initialize the simulation with density, velocity, and pressure profiles given by the analytical diverging solution at time t = B, where the diverging shock is at r = 1. We let the flow evolve until the time $t = 2 B$⁠.

For each case, we present the following plots: hydrodynamic profiles at the final time $t = 2 B$ (Figs. 9–11), the shock position (Fig. 12), the applied piston boundary condition (Fig. 13), the convergence of L₁ norms for $ρ , p , u , e$ (Fig. 14), and the convergence of $λ , B$ that are calculated from the simulation using a fit of the shock position to the form $( t B ) 1 λ$⁠.

In this setup, we initialize the simulation with density, velocity, and pressure profiles given by the analytical converging solution at time t = −1, where the converging shock is at r = 1. We let the flow evolve as the shock reaches the center at time t = 0 and reflects outward, until it reaches the position r = 1 at time t = B. For each case, we present the following plots: density and pressure profiles at the final time t = B (Figs. 15–17), the converging and diverging shock position as function of time (Fig. 18), the applied piston boundary condition (Fig. 19), the convergence of L₁ norms for $ρ , p , u , e$ (Fig. 20), and the convergence of the diverging shock parameters $λ out , B$⁠, which are calculated from the simulation using a fit of the diverging shock position to the form $( t B ) 1 λ$ for t > 0, and the converging shock exponent $λ in$⁠, which is calculated from the simulation using a fit of the converging shock position to the form $( − t ) 1 λ$ for t < 0.

In this work, we studied in detail the complete Guderley problem of a converging and then diverging shock in an ideal gas medium with an initial power-law density profile, $ρ ( r ) ∼ r μ$⁠. We extend our previous work in Ref. 25, which was dedicated to a detailed study of the converging shock. We presented the theoretical framework required to construct exact self-similar solutions for the converging and diverging flows. A generalized algorithm for a numerical calculation of the semi-analytical similarity solutions was developed. The diverging shock solution was studied in a wide range of parameters, namely, the adiabatic constant γ, the initial density exponent μ, and for spherical and cylindrical symmetries. The results are in excellent agreement with previous published work, which, to the best of our knowledge, are only given for the case of an initially constant density profile (μ = 0).

The semi-analytic solutions were used to define nontrivial compressible flow problems, which could serve for code verification. The solutions were compared to hydrodynamic simulations employing a one-dimensional Lagrangian hydrodynamic code, by applying appropriate initial and boundary conditions. A very good agreement was reached between numerical hydrodynamic simulations and the semi-analytical solutions. This success highlights the use of the solutions for the purpose of verification and validation of numerical hydrodynamics simulation codes.

The authors have no conflicts to disclose.

Itamar Giron: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Shmuel Balberg: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Menahem Krief: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In this appendix, we use the similarity representation [Eqs. (12)–(14)], in order to write the dimensional Rankine–Hugoniot jump relations [Eqs. (7)–(9)] for the diverging shock, in terms of the similarity profiles. We note that the location of the diverging shock varies according to the power law $r s ( t ) = ( t B ) 1 λ$⁠. Hence, the shock velocity in the jump relations (7)–(9) is $D = d r s d t = r s λ t$⁠.