Piston driven shock waves in non-homogeneous planar media

In this work, we analyze in detail the problem of piston driven shock waves in planar media. Similarity solutions to the compressible hydrodynamics equations are developed, for a strong shock wave, generated by a time dependent pressure piston, propagating in a non-homogeneous planar medium consisting of an ideal gas. Power law temporal and spatial dependency is assumed for the piston pressure and initial medium density, respectively. The similarity solutions are written in both Lagrangian and Eulerian coordinates. It is shown that the solutions take various qualitatively different forms according to the value of the pressure and density exponents. We show that there exist different families of solutions for which the shock propagates at constant speed, accelerates or slows down. Similarly, we show that there exist different types of solutions for which the density near the piston is either finite, vanishes or diverges. Finally, we perform a comprehensive comparison between the planar shock solutions and Lagrangian hydrodynamic simulations, by setting proper initial and boundary conditions. A very good agreement is reached, which demonstrates the usefulness of the analytic solutions as a code verification test problem.

In Ref. [24], self-similar solutions for ablation driven shock waves [25][26][27][28][29][30][31] in homogeneous media, are presented.The ablation pressure serves as a pressure piston boundary condition for the shock wave, with a given temporal power-law dependency of the form p 0 t τ .The separate treatment of the ablation region and the shock region enables the use of a binary-equation of state which is essential for an accurate modeling of some real solid materials experiments [32].For example, this solution enables an analytic theoretical modeling of subsonic radiative heat flow, used for evaluating the drive temperature in hohlraum experiments via a measurement of the shock wave velocity in aluminum wedges [11].
The ablation driven shock solutions [24], assume an ideal gas equation of state (EOS) and a medium with an initial homogeneous density.In this work we extend this problem and develop similarity solutions to the compressible hydrodynamics equations for a non-homogeneous medium with an initial density profile which is a spatial power law of the form ρ 0 x −ω , as described in Fig. 1.Power law density profiles are widely used in many fields.For example, in astrophysical modeling, as the interior and atmospheres of stars [33][34][35][36][37][38] and galaxies [39][40][41].
The generalized shock solutions are presented and analyzed in detail.Different qualitative behavior of the propagation and structure of the piston driven shock wave is shown to exist for different ranges of the temporal ex- * menahem.krief@mail.huji.ac.il ponent τ and spatial exponent ω.It is shown that the shock speed could be constant, accelerate and decelerate and that the density near the piston can be finite, vanish or diverge.The solutions are studied in both Lagrangian mass coordinates and in real space.We show that the transformation between these two representations is given by a closed formula, which demonstrates the behavior of the shocked fluid in the x − t plane.We study the temporal dependence of the energy supplied to the fluid by the pressure load, and how it is divided into kinetic and internal fluid energy.
Finally, we perform a comprehensive comparison between the semi-analytical planar shock solutions to numerical hydrodynamic simulations.As shown, a very good agreement is reached, which highlights the use of the new generalized solutions for the purpose of verification and validation of numerical hydrodynamics simulation codes.In addition, the solutions presented in this work can be used directly to generalize the ablation driven shock wave solutions [11,24], to the case of an initially non-homogeneous media.

II. STATEMENT OF THE PROBLEM
In this section we will list the governing equations of motion and boundary conditions for the hydrodynamic problem consisting of a piston which exerts a time dependent external pressure on a planar medium which is initially at rest, with an initial spatial power law density profile.The time dependent pressure is assumed to be a temporal power law.This problem is described schematically in Fig. 1.
Since the system boundary changes over time, it is customary to work in Lagrangian mass coordinates: where ρ (x, t) is the planar mass density.The hydrodynamic equations in planar symmetry and Lagrangian mass coordinates are given by (conservation of mass, momentum and energy, respectively): where v (m, t) = 1/ρ (m, t) is the specific volume, u (m, t) is the fluid velocity, p (m, t) is the pressure and e (m, t) is the specific internal energy (energy per unit mass).In this work we will assume a polytropic ideal gas equation of state with an adiabatic index γ = r + 1 so that: The boundary condition consists of a pressure piston with a power law time dependence in the form: The initial conditions are of a medium at rest: with a density profile which is a spatial power law of the form: The resulting initial specific volume is: where v 0 = 1/ρ 0 .In order to write the initial specific volume profile in terms of Lagrangian mass coordinates, we note that: where we have used the fact that the total mass is finite, which leads to the constraint: Eq. 10 gives the initial Eulerian to Lagrangian relation: and the initial specific volume in Lagrangian coordinates:

III. SELF-SIMILAR REPRESENTATION
The hydrodynamic problem defined in the previous section is defined in terms of N = 7 dimensional quantities: two independent variables m, t, three dependent hydrodynamic profiles p (m, t) , u (m, t) , v (m, t) and two dimensional parameters v 0 , p 0 .The dimensions of these quantities are listed in table I, in terms of K = 3 units: mass, length and time.Therefore, from the central theorem of dimensional analysis (a.k.a. the Pi theorem) [42][43][44], the problem can be fully characterized in terms of N − K = 4 dimensionless variables, which are given in terms of power laws of the dimensional quantities.Therefore, it is customary to use an independent dimensionless similarity variable ξ and three dependent dimensionless similarity profiles V (ξ) , U (ξ) , P (ξ) for the specific volume, velocity and pressure fields, respectively.Writing these similarity variables in terms of power laws of the dimensional quantities, yields the following self-similar representation: The dimensional quantities in the problem, and their dimensions in terms of mass (M ) length (L) and time (T ) units.
We note that the piston boundary condition 6 is written in terms of the similarity pressure profile as: By substitution of this self-similar representation in the hydrodynamics equations 2-4 and using the relations 2−ω ξ t , all dimensional quantities factor out, and a set of nonlinear ordinary differential equations for the similarity profiles is obtained: It is seen that for the special case of a spatially constant initial density (ω = 0), the equations are reduced to the form given in equations 43a-43c in Ref. [24].

IV. JUMP CONDITIONS
Since the medium is initially stationary with zero pressure, a strong shock is formed due to the operation of the pressure piston.The hydrodynamic profiles in the unshocked region are given simply by the initial conditions (zero pressure and velocity and a power law density profile, see equations 7,13).In the shocked region, the hydrodynamic profiles are obtained from the solution of the similarity equations 19-21.The transition between the shocked and unshocked regions, is obtained by Rankine-Hugoniot jump conditions.In this section we will derive expressions for these jump conditions in terms of the values of the similarity profiles at the shock surface.
First, we note that the general form of the Rankine-Hugoniot jump conditions for a set of conserved quantities u (y, t) which obey a set of (nonlinear) hyperbolic conservation laws in the form ∂u ∂t + ∂ ∂y F (u) = 0, are written as ẏ∆u = ∆F , where ẏ is the local discontinuity velocity and ∆u, ∆F are, respectively, the differences in the conserved quantities u and fluxes F (u) across the discontinuity [45,46].We want to derive the jump conditions for the hydrodynamic equations 2-4, which are written in Lagrangian form.We note that the internal energy equation 4, is not written in conservative form.Instead, we use the equation for the conservation of total energy: which is in conservative form.The Rankine-Hugoniot jump conditions for the conservation of mass (equation 2), momentum (equation 3) and total energy (equation 22) are given by: where 'ref' represents quantities in the (reference) unshocked region and 's' represents quantities in the shocked region.ṁ is the shock mass flux (the shock velocity in mass space).Since in the unshocked region u ref = p ref = 0, the jump conditions describe a strong shock [43,47]: In order to write these equations in terms of the similarity profiles, we note that the shocked mass and mass flux are given by: where ξ s is the similarity variable at the shock front, which is not known at this point.The temporal similarity exponent for the mass, (2+τ , is plotted in Figure 2 as a function of τ, ω.We see that the requirement for the shock to propagate, gives the following constraint on the piston pressure temporal exponent: The specific volume in the unshocked side is given by (see equation 13): while in the shocked side: where we have defined . By substituting equations 30-35 in the jump conditions 26-28, all dimensional terms factor out, and the following dimensionless jump conditions are obtained: Solving these equations for V s , U s , P s gives the following expressions for the similarity profiles at the shock in terms of ξ s : As before, we see that for ω = 0, these equations are reduced to equations 46a-46d in Ref. [24].

V. SOLUTION OF THE SIMILARITY EQUATIONS
In order to solve the ordinary differential equations (ODE) 19-21 for the similarity profiles, we write them in the canonical form y ′ = f (y) by employing Kramers' law: where: with the following coefficients: The initial conditions for equation 42 are the jump conditions 39-41 at ξ = ξ s .The piston boundary condition 18 gives a constraint from which the value of ξ s is determined.This can be performed numerically by finding the root of the function F (ξ s ) = P ξs (0) − 1 where P ξs (0) is the similarity pressure profile at ξ = 0 that was obtained by a numerical integration of the ODE 42 from a trial value ξ = ξ s to ξ = 0.An example of the root finding procedure for the calculation of ξ s is presented in Fig. 3.In Figures 4-5 we present calculations of ξ s as a function of τ, ω in a wide range.In Figure 6 the similarity profiles V (ξ), U (ξ), P (ξ) are plotted for various values of τ, ω.It is seen that for all cases shown, except the special case τ = ω = 0 (a constant piston pressure and a uniform initial density), the specific volume either diverges or vanishes near ξ = 0 (the location of the piston).This behavior can be characterized by performing an asymptotic analysis near ξ → 0, of equation 42 for the specific volume.Using the fact that U (ξ), P (ξ) are finite near ξ → 0 and, as will be shown subsequently in section VI, we also havelim ξ→0 ξV (ξ) = 0, the equation for the specific volume near ξ → 0 can be written as: where: This equation has a power law solution of the form V (ξ) = Aξ q .This asymptotic power is plotted in Figure 7 for a wide range of τ, ω.The line q (τ, ω) = 0, which is given by the relation: represents a family of solutions which have a finite density near the piston, and divides the τ − ω plane into regions for which the density near the piston diverges (q > 0), or vanishes (q < 0).We note that for the trivial case τ = ω = 0, we have q = 0, and the density is finite near the piston.In this case, the coefficients in equations 48-49 are zero, which results in the vanishing of the determinants 44-46, so that the ODE 42 has a trivial analytical solution for which the similarity profiles are constant (as also seen in Figure 6).Since the pressure at the piston (ξ = 0) is known from the boundary condition 18, the value of V, U and ξ s can be obtained analytically by inverting the strong shock jump conditions 39-41, resulting in the trivial solution: For other values of τ, ω for which q = 0 (according to equation 52), even though the density is finite near the origin, the determinants 44-46 are nonzero, and the solution is not constant, so a numerical procedure for the determination of ξ s , as described above is necessary.In summary, other than the trivial solution for τ = ω = 0,    A color plot for the asymptotic power q (see equation 51) of the specific volume which behaves as V (ξ) ≈ Aξ q near the piston, as a function of τ, ω and γ = 1.25.The line q (τ, ω) = 0 (given in equation 52) represents a family of solutions with a finite density near the piston, and divides the τ − ω plane into regions for which the density near the piston diverges, (q > 0) or vanishes (q < 0).
we are not aware of any closed form analytical solutions of equation 42, which therefore must be solved numerically, yielding a semi-analytical self similar solution.
Finally, we note that the ODE in the form 42 has a potential singularity if ∆ (V, P ) = 0.The emergence of singularities is a typical behavior which appears in the construction of similarity solutions of the second kind [43,44], and specifically in the work on converging shock waves [17,18,[48][49][50].Regarding the planar piston driven shock wave solution, if the determinant in equation 43 is zero, we must have the relation This relation results in a vanishing pressure near the piston (ξ = 0), and therefore cannot hold there, since the pressure is finite according to the piston boundary condition 18.Moreover, if a singular behavior existed in some finite range 0 < ξ 1 < ξ < ξ 2 , this relation can be inserted to the energy ODE 21, resulting in a non-singular power law solution of the form V (ξ) = Bξ l , where B is a constant and In practice, we did not encounter any singular behavior in the numerical integration of eq.42, for a wide range of τ, ω (see results in figures 4,9,14), which shows that the solutions did not pass through any singular points.We conclude that the piston driven shock wave equations 42 do not seem to posses any singularities, which is a common property of self-similar solutions of the first kind, such as the Noh stagnation shock or Sedov-Taylor-von Neumann blast wave problems.

VI. THE SOLUTION IN REAL SPACE
In order to write the solution in real space (that is, in Eulerian coordinates), we will now derive exact an expression for the temporal evolution of the positions of Lagrangian mass elements.For an element m > m s (t) which the shock did not reach at time t, the position is given by the initial position (via equation 12).For a shocked mass element m < m s (t), we have: where x p (t) is the piston position.By substituting the similarity volume profile (equation 15), one obtains the relation: where the upper integration limit ξ = ξ (m) is given from equation 14.We note that the integral in equation 58 can be calculated analytically, by integrating equation 19 one finds: (59) where we have used the fact that in order for the integral ξ 0 V (ξ ′ ) dξ ′ , which represents position, to be finite, we must have lim ξ→0 [V (ξ) ξ] = 0.
The piston position is obtained by a temporal integration of the piston velocity u p (t) ≡ u (m = 0, t), and using equation 16: In summary, the position of a Lagrangian element m as a function of time is given by the following formula: We note that for the shock position x s (t) ≡ x (m s (t) , t), both expressions in equation 61 coincide, since ξ = ξ s and due to the jump conditions 39-40.This gives a simple expression for the shock position in term of ξ s : The position temporal exponent 2+τ 2−ω , which governs the temporal behavior of any mass element (and specifically the shock and piston position), is plotted in Figure 8 as a function of τ, ω.The line 2+τ  2−ω = 0 represents a family of solutions with a constant shock velocity, and divides the τ − ω plane into regions for which the shock accelerates or slows down over time.
The quantity U (0), which sets the piston velocity, is plotted as a function of τ, ω in Figures 9-10.It was obtained from a numerical integration of the similarity ODE, as described in the previous section.
In Figure 11 several x − t diagrams are plotted for various values of τ, ω.These diagrams describe the exact temporal evolution of the locations of Lagrangian mass elements (via equation 61), through which the shock is propagating.One can see that different values of τ, ω result in constant speed, accelerating and decelerating motion (according to Figure 8).In addition, it is seen that for some combinations of τ, ω, the cells near the piston are compressed/expanded substantially, due to the fact that the density near the piston boundary is either constant, vanishing or diverging (according to Figure 7).
In Figures 12-13 we present the evolution of hydrodynamics profiles in both Lagrangian in Eulerian coordinates.The profiles were by obtained from the self-similar solution (via equations [15][16][17] for two cases: ω = τ = −0.5 and ω = τ = 0.5.As can be seen from the Figures, in the first case the density vanishes near the piston (according to Figure 7) and the shock propagation slows down over time (also evident in Figure 11), while for the latter case the density diverges near the piston and the shock accelerates.

VII. THE ENERGY IN THE SYSTEM
The total internal and kinetic energy in the system are given by, respectively: By substituting the self-similar representation 15-17 we find: from which it is evident that the ratio between kinetic and internal energy is time independent, and given by the numerical constant: 2 dξ ξs 0 This ratio is plotted as a function of τ, ω in Figure 14, where the integrals were calculated numerically.The total energy in the system is given by: We note that the integral expression can be calculated analytically, in a manner similar the derivation of equation 59.To that end, we use the similarity ODE for the conservation of total energy (equation 22): By a direct integration of this equation, we obtain an exact formula for the total energy from the piston to a mass element m: Using this identity, the total energy integral for the entire system can be obtained by using the jump condition 38: (71) It is not surprising that U (0) sets the total energy in the problem, as it is the only quantity that detriments the piston velocity (see equation 60).

VIII. COMPARISON TO HYDRODYNAMIC SIMULATIONS
In this section we compare the semi-analytical planar shock solutions to numerical hydrodynamic simulations.We employ a standard staggered-mesh artificialviscosity one-dimensional fully Lagrangian code (which  The locations of the piston (blue) and the shock (red) are also presented.The diagrams describe the exact locations of Lagrangian mass elements (via equation 61), through which the shock is propagating.It is evident that different values of τ, ω result in constant speed, accelerating and decelerating shocks (according to Figure 8).In addition, the behavior of the density (either constant, vanishing or diverging) near the piston boundary is evident (according to Figure 7).
is described in detail in appendix D of Ref. [17]).The pressure piston boundary condition is applied by setting a ghost cell adjacent to the leftmost zone with a temporally applied pressure p ghost (t) = p 0 t τ .The initial pressure and velocity are zero throughout the system, while the initial density is set according to the spatial power-law 8.
We compare the analytical solutions and the simulated results for the density, velocity, pressure and specific internal energy profiles, which are represented in mass coordinates as well as in real space.1000 computational zones are used in these simulations.In addition, in Fig. 21 we compare the analytical x − t diagrams (see section VI), to simulations with a lower number of 100 computational zones.
All simulation are performed with the parameters γ = 1.25, ρ 0 = 1g/cm and p 0 = 1dyn/cm 2 /sec τ .The final simulation time is defined by the time for which the analytical shock position is at 0.9cm (see eq. 62).We define 3 test cases.
1. τ = 0, ω = 0 -this is the trivial case of a constant applied pressure, which generates a shock which propagates into a homogeneous cold medium.All hydrodynamic profiles are constant in space (via the trivial shock jump conditions 39-41), as shown Figures 17-18.The shock velocity is constant in time (see equation 62).The total energy increases linearly in time, and the internal and kinetic energy are equal (see equations 65-67 and Figure 14).
2. τ = −0.5, ω = −0.5 -in this case the applied pressure decreases over time, and the generated shock propagates into a cold medium with a decreasing density profile.According to Figure 7, the density vanishes near the piston, as shown in Figures 17-18.According to equation 62 (see also Figure 8), the shock propagation slows down over time as x s (t) ∝ t 0.6 (also evident in Figure 21).In order to properly resolve the diminishing density near the piston, in this case the simulations were performed with non uniform (geometrically) spaced zones where the initial zone widths were set according to ∆x n = q∆x n−1 with q = 1.002 for the fine simulations (with 1000 zones), and q = 1.05 for the coarser simulations (with 100 zones).This geometrically spaced grid is evident in the x − t plots in Figure 21.
3. τ = 0.5, ω = 0.5 -in this case the applied pressure increases over time, and the generated shock propagates into a cold medium with an increasing density profile.According to Figure 7, the density diverges near the piston, as shown in Figures 19-20.According to equation 62 (see also Figure 8), the shock propagation accelerates over time as x s (t) ∝ t 5/2 (also evident in Figure 21).
As can be seen in Figures 15-21, in all cases we observe a very good agreement between the analytical solutions and the numerical simulations.In Figure 22, we show the spatial convergence of the numerical simulations to the analytic results for the final shock position (which is 0.9cm by construction).The resulting order of convergence, which is obtained by the slope of the relative error in log-log scale, is about 1.0 (that is, a linear convergence) for the cases ω = 0, τ = 0 and ω = −0.5, τ = −0.5, and about 1.3 (slightly better than linear) for the case ω = 0.5, τ = 0.5.

IX. SUMMARY
In this work we have generalized the solutions to the compressible hydrodynamics problem of piston driven shock waves in homogeneous planar ideal gas media, to non-homogeneous media with an initial density profile in the from of a spatial power law.Similarity solutions were developed systematically in both Lagrangian and Eulerian coordinates.It was shown that the generalized solutions take various qualitatively different forms according to the values the temporal piston pressure exponent τ and the spatial initial density exponent ω.Specifically, it was shown that there exist different families of solutions for which: (i) the shock propagates at constant speed, accelerates or slows down (ii) the density near the piston is either finite, vanishes or diverges and (iii) the total kinetic energy can be larger, equal or smaller than the total internal energy.
An elaborate comparison between the semi-analytical planar shock solutions to numerical hydrodynamic simulations was performed, and a very good agreement was reached.This highlights the use of the new solutions for the purpose of verification and validation of numerical hydrodynamics simulation codes.
Finally, we note that the solutions presented in this work, can be used directly to generalize the ablation driven shock wave solutions [11,24], to media which have a non-homogeneous initial density profile.Figure 12.Hydrodynamic profiles for the density, pressure and velocity, at various times (as written in the legends), for τ = −0.5, ω = −0.5, γ = 1.25, ρ0 = 1g/cm and p0 = 1dyn/cm 2 /sec τ .The profiles are plotted with respect to Lagrangian mass coordinates (left column) as well as in real space (right column), in which the dashed black and blue lines represent the piston and shock positions, respectively.In this case, according to Figure 7, the density vanishes near the piston, as can be seen in the density plots which are also given in log scale.12, but for τ = 0.5, ω = 0.5.In this case, according to Figure 7, the density diverges near the piston, as can be seen in the density plots.A comparison between the analytical solution (dashed blue lines) and a numerical simulation (red lines) for the hydrodynamic profiles: density (top left), specific internal energy (top right), velocity (bottom left) and pressure (bottom right), as a function of the Lagrangian mass coordinate m for the parameters τ = 0, ω = 0, γ = 1.25, ρ0 = 1g/cm and p0 = 1dyn/cm 2 /sec τ .The profiles are shown at various times where the shock positions are at 0.1, 0.2, 0.4, 0.5, 0.7, 0.9 cm (as can be seen in Fig. 16).The simulation was performed with 1000 zones.

Figure 1 .
Figure 1.A schematic description of the hydrodynamics problem: a time dependent pressure piston drives a strong shock wave, into a medium which is initially at rest with a power law density profile.

Figure 2 .
Figure 2. A color plot for the mass similarity exponent (the temporal power of the shocked mass, see equation 29), as a function of τ, ω.

25 Figure 3 .
Figure 3.The function F (ξs) = P ξs (0) − 1 which has a root at the correct value of the similarity variable at the shock, ξs.The calculation was performed for ω = 0.5, τ = −0.45 and γ = 1.25, resulting in ξs = 2.0189 (as shown by the black dashed vertical line).

Figure 4 .
Figure 4.A color plot for the shock front similarity variable ξs as a function of the piston pressure temporal power τ and the initial density spatial power ω and an adiabatic index γ = 1.25.

)Figure 5 .
Figure 5. Shock front similarity variable ξs as a function of τ for selected values of ω and as a function of τ, ω and γ = 1.25.

Figure 7 .
Figure 7.A color plot for the asymptotic power q (see equation 51) of the specific volume which behaves as V (ξ) ≈ Aξ q near the piston, as a function of τ, ω and γ = 1.25.The line q (τ, ω) = 0 (given in equation 52) represents a family of solutions with a finite density near the piston, and divides the τ − ω plane into regions for which the density near the piston diverges, (q > 0) or vanishes (q < 0).

Figure 8 .
Figure 8.A color plot for the position similarity exponent2+τ 2−ω (which is the temporal power of the position of any mass element, including the shock and piston positions, see text), as a function of τ, ω.The line 2+τ 2−ω = 1 represents solutions with a constant shock speed, and divides the τ − ω plane into regions of solutions with accelerating (above the line) or decelerating (below the line) shocks.

Figure 11 .
Figure 11.x − t diagrams for γ = 1.25, ρ0 = 1g/cm, p0 = 1dyn/cm 2 /sec τ and various values of τ, ω, as listed in the title of each subplot.The locations of the piston (blue) and the shock (red) are also presented.The diagrams describe the exact locations of Lagrangian mass elements (via equation 61), through which the shock is propagating.It is evident that different values of τ, ω result in constant speed, accelerating and decelerating shocks (according to Figure8).In addition, the behavior of the density (either constant, vanishing or diverging) near the piston boundary is evident (according to Figure7).

Figure 14 .
Figure 14.A color plot of the ratio between the total kinetic and internal energies in the system (equation 67), as a function of τ, ω for γ = 1.25.

Figure 15 .
Figure 15.A comparison between the analytical solution (dashed blue lines) and a numerical simulation (red lines) for the hydrodynamic profiles: density (top left), specific internal energy (top right), velocity (bottom left) and pressure (bottom right), as a function of the Lagrangian mass coordinate m for the parameters τ = 0, ω = 0, γ = 1.25, ρ0 = 1g/cm and p0 = 1dyn/cm 2 /sec τ .The profiles are shown at various times where the shock positions are at 0.1, 0.2, 0.4, 0.5, 0.7, 0.9 cm (as can be seen in Fig.16).The simulation was performed with 1000 zones.

Figure 16 .Figure 17 .Figure 18 .
Figure16.Same comparison as in Fig.15, with the hydrodynamic profiles as a function of position (real space representation).The piston positions are indicated by dashed-dotted vertical lines.

Figure 22 .
Figure 22.Spatial convergence plots for the shock position at the end of the three simulations defined in section VIII.The values of ω, τ are listed in the titles.Shown are the relative error between the final shock positions as obtained in the simulations and the analytic result, as a function of the number of computational zones N cells .The (fitted) convergence rate, which is given by the slope of the convergence line in log-log scale, is written in the legends of each figure. 5.