Stability of strong ideal-gas shock layers
Extending recent results in the isentropic case, we use a combination of asymptotic ODE estimates and numerical Evans-function computations to examine the spectral stability of shock-wave solutions of the compressible Navier–Stokes equations with ideal gas equation of state. Our main results are that, in appropriately rescaled coordinates, the Evans function associated with the linearized operator about the wave (i) converges in the large-amplitude limit to the Evans function for a limiting shock profile of the same equations, for which internal energy vanishes at one endstate; and (ii) has no unstable (positive real part) zeros outside a uniform ball . Thus, the rescaled eigenvalue ODE for the set of all shock waves, augmented with the (nonphysical) limiting case, form a compact family of boundary-value problems that can be conveniently investigated numerically. An extensive numerical Evans-function study yields one-dimensional spectral stability, independent of amplitude, for gas constant in and ratio of heat conduction to viscosity coefficient within (, for air). Other values may be treated similarly but were not considered. The method of analysis extends also to the multi-dimensional case, a direction that we shall pursue in a future work.
- 1 Introduction
- 2 Preliminaries
- 3 Evans function formulation
- 4 High-frequency bounds
5 Numerical protocol
- 5.1 Approximation of the profile
- 5.2 Bounds on unstable eigenvalues
- 5.3 Approximation of the Evans function
- 5.4 Winding number computation
- 6 Numerical results
- A Gas constants for air
- B Gas constants for other gases
- C Liquids and dense gases
- D Computation of the Mach number
- E Lifted Matrix bounds
- F Temperature-dependent viscosity
A long-standing question in gas dynamics is the stability of viscous shock layers, or traveling-wave solutions
of the compressible Navier–Stokes equations, where is a vector recording specific volume, velocity, and total energy of the fluid at location and time . A closely related question is the relation between Navier–Stokes solutions and solutions of the formally limiting Euler equations in the limit as viscosity and heat conduction coefficients go to zero: more precisely, validity of formal matched asymptotics predicting that the Navier–Stokes solution consists approximately of an Euler solution with smooth viscous shock layers replacing discontinuous Euler shocks.
Recent progress in the form of “Lyapunov-type” theorems established in [38, 54, 22, 23, 24] has reduced both problems to determination of spectral stability of shock layers, i.e., the study of the eigenvalue ODE associated with the linearized operator about the wave: a standard analytically and numerically well-posed (boundary value) problem in ODE that can be attacked by the large body of techniques developed for asymptotic, exact, and numerical study of ODE. Indeed, the cited results hold for a substantially more general class of equations, and in one- or multi-dimensions.
In [29, 44, 15, 16], it has been established in a similarly general context (general equations, one- and multi-dimensions), using asymptotic ODE techniques, that spectral stability holds always in the small-amplitude limit, where amplitude is measured by , i.e., for shocks sufficiently close to a constant solution, thus satisfactorily resolving the long-time stability and small-viscosity limit problems for small-variation solutions.
However, until very recently, the spectral stability of large-amplitude shock waves has remained from a theoretical viewpoint essentially open, the sole exceptions being (i) a result of stability of Navier–Stokes shocks for isentropic gas dynamics with -law gas in the special case , obtained by Matsumura-Nishihara  quite early on using an ingenious energy estimate specific to that case; and (ii) and a result of Zumbrun —again obtained by energy estimates special to the model—which establishes the stability of stationary phase-transitional shocks of an isentropic viscous-capillary van der Waals model introduced by Slemrod .
Progress instead has focused, quite successfully, on the development of efficient and general numerical methods for the study of stability of individual waves, or compact families of waves, of essentially arbitrary systems; see, for example, [6, 7, 8, 5, 30, 4]. These techniques, based on Evans-function computations, effectively resolve the question of spectral stability for waves of large but finite amplitude, but leave open the question of stability in the large-amplitude limit. For discussion of the Evans function and its numerical computation, see [1, 17, 54, 8, 4] or Section 3.4 below.
Quite recently, however, Humpherys, Lafitte, and Zumbrun 
have introduced a new strategy combining asymptotic ODE techniques with
numerical Evans-function computations,
by which they were able to carry out a global analysis
of shock stability in the isentropic -law case,
yielding stability independent of amplitude for
The purpose of the present paper is to extend the approach of  to the full (nonisentropic) Navier–Stokes equations of compressible gas dynamics with ideal gas equation of state, resolving in this fundamental case the long-standing questions of viscous shock stability and behavior in the small-viscosity limit. Specifically, we show, as in the isentropic case, that the Evans function indeed converges in the large-amplitude limit, to a value corresponding to the Evans function of a limiting system. Compactifying the parameter range by adjoining this limiting system, we then carry out systematic numerical Evans-function computations as in [4, 27] to determine stability for gas constant and (rescaled) ratio of heat conduction to viscosity coefficient , well-including the physical values given in Appendices A and B. The result, as in the isentropic case, is unconditional stability, independent of amplitude for an ideal gas equation of state.
1.1. Discussion and open problems
The asymptotic analysis of  is quite delicate; it depends sensitively both on the use of Lagrangian coordinates and on the precise way of writing the eigenvalue ODE as a first-order system. It is thus not immediately clear that the analysis can be extended to the more complicated nonisentropic case. Moreover, since Lagrangian coordinates—specifically, the associated change of spatial variable
where is density—are not available in multi-dimensions, it is likewise, at first glance, unclear how how to extend the analysis beyond one spatial dimension.
Remarkably, we find that the structure of the full, physical equations is much more favorable to the analysis than that of the isentropic model. In particular, whereas in the isentropic case the eigenvalue equations in the large-amplitude limit are a nonstandard singular perturbation of the limiting equations that must be analyzed “by hand”, in the full (nonisentropic) gas case, they are a regular perturbation for which convergence may be concluded by standard theorems on continuous dependence of the Evans function with respect to parameters; see, for example, the basic convergence lemma of .
Indeed, for bounded away from the nonphysical case (see Section 2 for a description of the equations and the physical background), we have the striking difference that, for a fixed left endstate , the density remains uniformly bounded above and below for all viscous shock profiles connecting to a right state , with energy going to infinity in the large-amplitude limit. By contrast, in the isentropic case, the density is artificially tied to energy and thus density goes to infinity in the large-amplitude limit for any ; see, e.g., [46, 47, 48]. This untangling of the large-amplitude behaviors of the density and the energy sets the stage for our analysis. Below, to see this untangling, instead of fixing a left endstate and asking which right endstates may be connected to by a viscous shock profile, we fix the shock speed and all coordinates of the left state except the energy. We find again that the density stays bounded above and below for all possible right states connected by a shock profile to some such .
Since the equations remain regular so long as density is bounded from zero and infinity, one important consequence of this fact is that we need only check a few basic properties such as uniform decay of profiles and continuous extension of stable/unstable subspaces to conclude that the strong-shock limit is in the nonisentropic case a regular perturbation of the limiting system as claimed; see Section 3 for details.
A second important consequence is that Lagrangian and Eulerian coordinates are essentially equivalent in the nonisentropic case so long as remains uniformly bounded from , whereas, in the isentropic case, the equations become singular for Eulerian () coordinates in the large-amplitude limit, by (1.1) together with the fact that density goes to infinity. Here, we have chosen to work with Lagrangian coordinates for comparison with previous one-dimensional analyses in the isentropic case [4, 27, 12]. However, we could just as well have worked in Eulerian coordinates, including full multi-dimensional effects, to obtain by the same arguments that the large-amplitude limit is a regular perturbation of the (unintegrated) limiting eigenvalue equation, and therefore the Evans functions converge in the limit also in this multi-dimensional setting.
Likewise, uniform bounds on unstable eigenvalues may be obtained in multi-dimensions by adapting the asymptotic analysis of  similarly as we have adapted in Section 4 the asymptotic analysis of . Thus, for uniformly bounded from , the analysis of this paper extends with suitable modification to the multi-dimensional case, making possible the resolution of multi-dimensional viscous stability by a systematic numerical Evans-function study as in the one-dimensional case. We shall carry out the multi-dimensional analysis in a following work .
Presumably, the same procedure of compactifying the parameter space after rescaling to bounded domain would work for any gas law with appropriate asymptotic behavior as . Thus, we could in principle investigate also van der Waals gas/fluids, for example, which could yield interesting different behavior: in particular, (as known already from stability index considerations [57, 54]) instability in some regimes. Other interesting areas for investigation include the study of boundary layer stability (see  for an analysis of the isentropic case), and stability of weak and strong detonation solutions (analogous to shock waves) of the compressible Navier–Stokes equations for a reacting gas. A further interesting direction is to investigate the effects of temperature dependence of viscosity and heat conduction on behavior for large amplitudes; see Appendices B.2 and F.
In this work, we have restricted to the parameter range and , where is the gas constant, is a rescaled coefficient of heat conduction ( the Fourier conduction and specific heat), and is the coefficient of viscosity; see equations. (2.1)–(2.3), Section 2. Similar computations may be carried out for arbitrary bounded away from the nonphysical limit . To approach the singular limit would presumably require a nonstandard singular perturbation analysis like that of  in the isentropic case, as the structure is similar; see Remark 2.2. The limits and are more standard singular perturbations with fast/slow structure that should be treatable by the methods of ; this would be a very interesting direction for further study. We note that our results for large do indicate possible further simplification in behavior, as the singular perturbation structure would suggest; see Remark 4.8 and Figure 4. For dry air at normal temperatures, and , well within range; see Appendix A.
Finally, we mention the issue of rigorous verification. Our results, though based on rigorous analysis, do not constitute numerical proof, and are not intended to. In particular, we do not use interval arithmetic. Nonetheless, the numerical evidence for stability appears overwhelming, particularly in view of the fact that the family of Evans contours estimated in the stability computations is analytic in both parameters, yielding extremely strong interpolation estimates by the rigidity of analytic functions.
In any case, our analysis contains all of the elements necessary for numerical proof, the effective realization of which, however, is a separate problem of independent interest. Given the fundamental nature of the problem, we view this as an important area for further investigation.
1.2. Plan of the paper
In Section 2, we set up the problem, describing the equations, rescaling appropriately, and verifying existence and uniform decay of profiles independent of shock strength. In Section 3, we construct the Evans function and establish the key fact that it is continuous in all parameters up to the strong-shock limit. In Section 4, we carry out the main technical work of the paper, establishing an upper bound on the modulus of unstable eigenvalues of the linearized operator about the wave in terms of numerically approximable quantities associated with the traveling-wave profile. In Section 5, we describe our numerical method, first estimating a maximal radius within which unstable eigenvalues are confined, then computing the winding number of the Evans function around the semicircle with that radius to estimate the number of unstable eigenvalues, for (a discretization of) all parameters within the compact parameter domain, including the strong-shock limit. Finally, in Section 6, we perform the numerical computations indicating stability.
In Appendices A and B, we discuss further the dimensionless constants and , and determine their values for air and other common gases. In Appendix C, we discuss equations of state for fluids and dense gases. In Appendix D, we compute a formula for the Mach number, a useful dimensionless quantity measuring shock strength independent of scaling. In Appendix E, we give a general bound on the operator norm of lifted matrices acting on exterior products, useful for analysis of the exterior product method of [6, 5]. In Appendix F, we discuss the changes needed to accommodate temperature-dependence in the coefficients of viscosity and heat conduction, as predicted by the kinetic theory of gases.
In Lagrangian coordinates, the Navier–Stokes equations for compressible gas dynamics take the form
where is the specific volume, is the velocity, is the pressure, and the energy is made up of the internal energy and the kinetic energy:
The constants and represent viscosity and heat conductivity. Finally, is the temperature, and we assume that the internal energy and the pressure are known functions of the specific volume and the temperature:
An important special case occurs when we consider an ideal, polytropic gas. In this case the energy and pressure functions take the specific form
where and are constants that characterize the gas. Alternatively, the pressure may be written as
where , the adiabatic index. Equivalently, in terms of the entropy and specific volume, the pressure reads
In the thermodynamical rarefied gas approximation, is the average over constituent particles of , where is the number of internal degrees of freedom of an individual particle, or, for molecules with “tree” (as opposed to ring, or other more complicated) structure,
where is the number of constituent atoms : for monatomic, for diatomic gas. For fluids or dense gases, is typically determined phenomenologically . In general, is usually taken within in models of gas-dynamical flow, whether phenomenological or derived by statistical mechanics [46, 47, 48].
2.1. Viscous shock profiles
2.2. Rescaled equations
Under the rescaling
where the pressure and internal energy in the (new) rescaled variables are given by
in the ideal gas case, the pressure and internal energy laws remain unchanged
with the same , . Likewise, remains unchanged in (2.6).
2.3. Rescaled profile equations
together with the boundary conditions
Evidently, we can integrate each of the differential equations from to , and using the boundary conditions (in particular and ), we find, after some elementary manipulations, the profile equations:
We note that in the case of an ideal gas, with , these ODEs simplify somewhat, to
where and is as in (2.6).
2.4. Rankine-Hugoniot conditions
where denotes jump between .
, from which we obtain the physicality condition
corresponding to positivity of the denominator, with as . Finally, substituting into and rearranging, we obtain
completing the description of the endstates.
We see from this analysis that the strong-shock limit corresponds, for fixed , to the limit , with all other parameters functions of . In this limit,
At this point, taking without loss of generality , we have reduced to a three-parameter family of problems on compact parameter range, parametrized by , , and .
As , we find from (2.35) that blows up as , i.e., our rescaled coordinates remain bounded only if , constant, as well. (This is reflected in the limiting profile equation for , which admits only profiles from to ; see (2.25), which, for , reduces to .) Thus, our techniques apply for only in the (simultaneous) large-amplitude limit.
2.5. Existence and decay of profiles
Specializing to the ideal gas case, we next study existence and behavior of profiles. Existence and exponential decay of profiles has been established by Gilbarg  for all finite-amplitude shocks . Thus, the question is whether these properties extend to the strong-shock limit, the main issue being to establish uniform exponential decay as , independent of shock strength.
Since the profile equations (2.24)–(2.25) are smooth (polynomial) in , the issue of uniform decay reduces essentially to uniform hyperbolicity of endstates , i.e., nonexistence of purely imaginary linearized growth/decay rates at . Linearizing (2.24)–(2.25) about an equilibrium state, we obtain
Since is , its eigenvalues are
and so hyperbolicity is equivalent to and or .
Computing, we have
so that is equivalent (for , hence ) to
At , this reduces to , or, using (2.34) to
except in the characteristic case , while
except in the characteristic case . Thus, for bounded from zero, hyperbolicity fails at only in the characteristic case .
Next, let us recall the existence proof of , which proceeds by the observation that isoclines and obtained by setting the right-hand sides of (2.24) and (2.25) to zero bound a convex lens-shaped region whose vertices are the unique equilibria , that is invariant under the forward flow of (2.24)–(2.25), and into which enters the unstable manifold of ; recall that at , hence there is a one-dimensional unstable manifold. It follows that the unstable manifold must approach attractor as , determining the unique connecting orbit describing the profile.
By the above-demonstrated hyperbolicity, this argument extends also to the case (), yielding at once existence and uniform boundedness of profiles across the whole parameter range (recall that are uniformly bounded, by the Rankine–Hugoniot analysis of the previous section), in particular for the limiting profile equations at of
Collecting facts, we have the following key result.
Existence, boundedness, and exponential decay of individual profiles follow from the discussion above. Uniform bounds follow by smooth dependence on parameters together with compactness of the parameter range. ∎
3. Evans function formulation
From now on, we specialize to the ideal gas case, setting without loss of generality .
3.1. Linearized integrated eigenvalue equations
Defining integrated variables
The integrated viscous profile,
is a stationary solution of (3.1). Then we write
and we linearize (3.1) about the integrated profile. By an abuse of notation, we denote the perturbation by , , and . Note also that the integrated profile always appears under an -derivative; this explains the appearance of “hats” in place of “tilde-hats” in the expression below. Finally, we use the relationship to simplify some of the expressions. We obtain the linearized integrated equations
Defining , subtracting times the second equation from the third, and rearranging, we obtain, finally, the linearized integrated eigenvalue equations:
3.2. Expression as a first-order system
where, using and (2.24) with ,
Remarkably, similarly as for the profile equations, the entries of are polynomial in . Thus, both profile and linearized eigenvalue equations are perfectly well-behaved for any compact subset of .
3.3. Consistent splitting
Denote by the limiting coefficient matrices at . (These limits exist by exponential convergence of profiles , Lemma 2.3.) Denote by and the stable and unstable subspaces of .
By analytic dependence of on and standard matrix perturbation theory, and are analytic on any domain for which consistent splitting holds.
Consistent splitting and analytic extension to follow by the general results of , except at the strong-shock limit , where
where and .
The matrix is lower block-triangular, with diagonal blocks
corresponding respectively to the scalar convection-diffusion equation
and the isentropic case treated in , The first has eigenvalues , so satisfies consistent splitting on , with analytic continuation (since eigenvalues remain separated) to . The second, as observed in , has eigenvalues , hence likewise satisfies consistent splitting on and (since the single unstable eigenvalue remains separated from the two stable eigenvalues) continues analytically to . Indeed, the unstable manifold has dimension two for all , hence is analytic on that domain. This verifies the proposition at by direct calculation.
At , the computation is more difficult. Here, we refer instead to the abstract results of , which assert that hyperbolic–parabolic systems of the type treated here, including the limiting case, at least for , exhibit consistent splitting on , with analytic extension to , so long as the shock is noncharacteristic, i.e., the flux Jacobian associated with the first-order part of the equations have nonvanishing determinants at . These may be computed in any coordinates, in particular . Neglecting terms originating from diffusion, i.e., including only first-order terms from the left-hand side, we obtain from (3.3) the flux Jacobian
which has determinant , giving for () that and, calculating at that
Thus, we may conclude by the general results of  that consistent splitting holds at both on for , with analytic extension to . ∎
We note that the results of  do not apply at , , where leaves the physical domain. Specifically, at this value the genuine coupling condition of Kawashima , hence the dissipativity condition of  fails, and so we cannot conclude consistent splitting; indeed, the eigenvalue (corresponding to the decoupled hyperbolic mode) is pure imaginary for any pure imaginary .
3.4. Construction of the Evans function
There exist bases , of and , extending analytically in and continuously in to for , , and , determined by Kato’s ODE
where denotes the spectral projection onto , , respectively, and denotes .
The Evans function is analytic in and continuous in on and , , and . Moreover, on , its zeros correspond in location and multiplicity with eigenvalues of the integrated linearized operator , or, equivalently with solutions of (3.3) decaying at .
The specification in (3.12) of initializing bases at infinity is optimal in that it minimizes “action” in a certain sense; see  for further discussion. In particular, for any constant-coefficient system, the Evans function induced by Kato bases (3.12) is identically constant in . For, in this case, bases and are given at by the values prescribed in (3.12), and both evolve according to the same ODE, hence satisfies and by Abel’s Theorem and the fact that , where .
More generally, in the “traveling pulse” case , by the argument of Remark 3.10, whence the Evans function constructed here may be seen to agree with the “natural” Evans function (independent of the choice of )
defined in  for that case. The latter may in turn be seen to agree with a (-modified) characteristic Fredholm determinant of the associated linearized operator about the wave , formally equivalent to