Becker–Morduchow–Libby solution

The Becker–Morduchow–Libby solution is a mathematical description of how the properties of a gas—such as pressure, temperature, and density—change within a one-dimensional shock wave. It is an exact solution of the Navier–Stokes equations for a compressible, viscous, heat-conducting gas. First derived in a simplified form by German physicist Richard Becker in 1922, it was later extended by Morris Morduchow and Paul A. Libby in 1949. Independent versions were also found in 1944 by M. Roy and L. H. Thomas. The solution revealed an unexpected feature: the entropy of the gas can vary in a non-monotonic way across the shock, rather than increasing steadily. Earlier work by Lord Rayleigh and G. I. Taylor had addressed special cases, such as flows with viscosity but no heat conduction, or vice versa.

Mathematical description

In a frame fixed with a planar shock wave, the shock wave is steady. In this frame, the steady Navier–Stokes equations for a viscous and heat conducting gas can be written as $$\begin{array}{cc}\frac{d}{dx}(\rho u)& =0,\\ \rho u\frac{du}{dx}+\frac{dp}{dx}-\frac{4}{3}\frac{d}{dx}\left({\mu }^{\prime }\frac{du}{dx}\right)& =0,\\ \rho u\frac{d\epsilon }{dx}+p\frac{du}{dx}-\frac{d}{dx}\left(\lambda \frac{dT}{dx}\right)-\frac{4}{3}{\mu }^{\prime }{\left(\frac{du}{dx}\right)}^2& =0,\end{array}$$ where $\rho$ is the density, $u$ is the velocity, $p$ is the pressure, $\epsilon$ is the internal energy per unit mass, $T$ is the temperature, ${\mu }^{\prime }=\mu +3\zeta /4$ is an effective coefficient of viscosity, $\mu$ is the coefficient of viscosity, $\zeta$ is the second viscosity and $\lambda$ is the thermal conductivity. To this set of equations, one has to prescribe an equation of state $f(p,\rho ,T)=0$ and an expression for the energy in terms of any two thermodynamics variables, say $\epsilon =\epsilon (p,\rho )$. Instead of $\epsilon$, it is convenient to work with the specific enthalpy $h=\epsilon +p/\rho .$

Let us denote properties pertaining upstream of the shock with the subscript "$0$" and downstream with "$1$". The shock wave speed itself is denoted by $D=u_0$. The first integral of the governing equations, after imposing the condition that all gradients vanish upstream, are found to be $$\begin{array}{cc}\rho u& ={\rho }_{0}D,\\ p+\rho u^2-\frac{4}{3}{\mu }^{\prime }\frac{du}{dx}& =p_0+{\rho }_0{D}^2,\\ h+\frac{u^2}{2}-\frac{1}{{\rho }_{0}D}(\lambda \frac{dT}{dx}+\frac{4}{3}{\mu }^{\prime }u\frac{du}{dx})& =h_0+\frac{D^2}{2}.\end{array}$$

By evaluating these on the downstream side where all gradients vanish, one recovers the familiar Rankine–Hugoniot conditions, ${\rho }_1{u}_1={\rho }_{0}D$, $p_1+{\rho }_1{u}_1^2=p_0+{\rho }_0{D}^2$ and $h_1+u_1^2/2=h_0+D^2/2.$ Further integration of the above equations require numerical computations, except in one special case where integration can be carried out analytically.

Analytical solution

Becker–Morduchow–Libby solution

Two assumptions has to be made to facilitate explicit integration of the third equation. First, assume that the gas is ideal (polytropic since we shall assume constant values for the specific heats) in which case the equation of state is $p/\rho T=c_p(\gamma -1)/\gamma$ and further $h=c_{p}T$, where $c_p$ is the specific heat at constant pressure and $\gamma$ is the specific heat ratio. The third equation then becomes

$$h+\frac{u^2}{2}-\frac{4{\mu }^{\prime }}{3{\rho }_{0}D}(\frac{3}{4P{r}^{\prime }}\frac{dh}{dx}+\frac{1}{2}\frac{d{u}^2}{dx})=h_0+\frac{D^2}{2}$$

where $P{r}^{\prime }={\mu }^{\prime }{c}_p/\lambda$ is the Prandtl number based on ${\mu }^{\prime }$; when $\zeta =0$, say as in monoatomic gases, this Prandtl number is just the ordinary Prandtl number $Pr=\mu c_p/\lambda$. The second assumption made is $P{r}^{\prime }=3/4$ so that the terms inside the parenthesis becomes a total derivative, i.e., $d(h+u^2/2)/dx$. This is a reasonably good approximation since in normal gases, Pradntl number is approximately equal to $0.72$. With this approximation and integrating once more by imposing the condition that $h+u^2/2$ is bounded downstream, we find

$$h+\frac{u^2}{2}=h_0+\frac{D^2}{2}.$$

This above relation indicates that the quantity $h+u^2/2$ is conserved everywhere, not just on the upstream and downstream side. Since for the polytropic gas $h=c_{p}T=\gamma p\upsilon /(\gamma -1)=c^2/(\gamma -1)$, where $\upsilon =1/\rho$ is the specific volume and $c$ is the sound speed, the above equation provides the relation between the ratio $p(x)/p_0$ and the corresponding velocity (or density or specific volume) ratio

$$\eta (x)=\frac{u}{D}=\frac{{\rho }_0}{\rho }=\frac{\upsilon }{{\upsilon }_0},$$

i.e.,Template:R

$$\begin{array}{cc}\frac{p}{p_0}& =\frac{1}{\eta }[1+\frac{\gamma -1}{2}{M}_0^2(1-{\eta }^2)]=\frac{{\eta }_1(\gamma +1)/(\gamma -1)-{\eta }^2}{[{\eta }_1(\gamma +1)/(\gamma -1)-1]\eta },\\ {\eta }_1& =\frac{\gamma -1}{\gamma +1}+\frac{2}{(\gamma +1)M_0^2},\end{array}$$

where $M_0=D/c_0$ is the Mach number of the wave with respect to upstream and ${\eta }_1=u_1/D={\rho }_0/{\rho }_1={\upsilon }_1/{\upsilon }_0$. Combining this with momentum and continuity integrals, we obtain the equation for $\eta (x)$ as follows

$$\frac{8\gamma }{3(\gamma +1)}\frac{{\mu }^{\prime }}{{\rho }_{0}D}\eta \frac{d\eta }{dx}=-(1-\eta )(\eta -{\eta }_1).$$

We can introduce the reciprocal-viscosity-weighted coordinate

$$\xi =\frac{3(\gamma +1)}{8\gamma }{\rho }_{0}D{\int }_0^x\frac{dt}{{\mu }^{\prime }(t)}$$

where ${\mu }^{\prime }(x)={\mu }^{\prime }[T(x)]$, so that

$$\eta \frac{d\eta }{d\xi }=-(1-\eta )(\eta -{\eta }_1).$$

The equation clearly exhibits the translation invariant in the $x$-direction which can be fixed, say, by fixing the origin to be the location where the intermediate value $({\eta }_1+1)/2$ is reached. Using this last condition, the solution to this equation is found to beTemplate:R

$$\frac{1-\eta }{(\eta -{\eta }_1{)}^{{\eta }_1}}={\left(\frac{1-{\eta }_1}{2}\right)}^{1-{\eta }_1}{e}^{(1-{\eta }_1)\xi }.$$

As $\xi \to -\mathrm{\infty }$ (or, $x\to -\mathrm{\infty }$), we have $\eta \to 1$ and as $\xi \to +\mathrm{\infty }$ (or, $x\to +\mathrm{\infty }$), we have $\eta \to {\eta }_1.$ This ends the search for the analytical solution. From here, other thermodynamics variables of interest can be evaluated. For instance, the temperature ratio $T/T_0$ is to found to given by

$$\frac{T}{T_0}=1+\frac{\gamma -1}{2}{M}_0^2(1-{\eta }^2)$$

and the specific entropy $s=c_p\ln (p^{1/\gamma }/\rho )=c_p\{\ln (p/\rho )-[(\gamma -1)/\gamma ]\ln p\}$, by

$$\frac{s-s_0}{c_p}=\ln \frac{T}{T_0}-\frac{\gamma -1}{\gamma }\ln \frac{p}{p_0}.$$

The analytical solution is plotted in the figure for $\gamma =1.4$ and $M_0=2$. The notable feature is that the entropy does not monotonically increase across the shock wave, but it increases to a larger value and then decreases to a constant behind the shock wave. Such scenario is possible because of the heat conduction, as it will become apparent by looking at the entropy equation which is obtained from the original energy equation by substituting the thermodynamic relation $Tds=d\epsilon +pd\upsilon =dh-{\rho }^{-1}dp$, i.e.,

$$\rho uT\frac{ds}{dx}=\frac{4}{3}{\mu }^{\prime }{\left(\frac{du}{dx}\right)}^2+\frac{d}{dx}\left(\lambda \frac{dT}{dx}\right).$$

While the viscous dissipation associated with the term $(du/dx{)}^2$ always increases the entropy, heat conduction increases the entropy in the colder layers where $d(\lambda dT/dx)/dx>0$, whereas it decreases the entropy in the hotter layers where $d(\lambda dT/dx)/dx<0$.

Taylor's solution: Weak shock waves

When $P{r}^{\prime }\ne 3/4$, analytical solution is possible only in the weak shock-wave limit, as first shown by G. I. Taylor in 1910.Template:R In the weak shock-wave limit, all terms such as $p-p_0$, $\rho -{\rho }_0$ etc., will be small. The thickness $\delta$ of the shock wave is of the order $\delta \sim 1/(p_1-p_0)$ so that differentiation with respect to $x$ increases the order smallness by one; e.g. $dp/dx$ is a second-order small quantity. Without going into the details and treating the gas to a generic gas (not just polytropic), the solution for $p(x)$ is found to be related to the steady travelling-wave solution of the Burgers' equation and is given by

$$p(x)=\frac{1}{2}(p_1+p_0)+\frac{1}{2}(p_1-p_0)\tanh \frac{x}{\delta }$$

where

$$\delta =\frac{8a{\rho }_0{c}_0^4}{(p_1-p_0)\mathrm{\Gamma }},a=\frac{{\lambda }_0}{2{\rho }_0{c}_0^3{c}_p}(\frac{4}{3}P{r}^{\prime }+\gamma -1),$$

in which $\mathrm{\Gamma }$ is the Landau derivative (for polytropic gas $\mathrm{\Gamma }=(\gamma +1)/2$) and $a$ is a constant which when multiplied by some characteristic frequency squared provides the acoustic absorption coefficient. The specific entropy is found to be proportional to $(1/T)dp/dx$ and is given by

$$\frac{s(x)-s_0}{c_p}=\frac{{\lambda }_0\mathrm{\Gamma }}{8a{c}_p{c}_0{T}_0}{\left(\frac{\partial T}{\partial p}\right)}_{s_0}\frac{(p_1-p_0{)}^2}{{\rho }_0^2{c}_0^4}\frac{1}{{\cosh }^2(x/\delta )}.$$

Note that $s(x)-s_0$ is a second-order small quantity, although $s_1-s_0$ is a third-order small quantity as can be inferred from the above expression which shows that $s=s_0$ for both $x\to \pm \mathrm{\infty }$. This is allowed since $s(x)$, unlike $p(x)$, passes through a maximum within the shock wave.

Validity of continuum hypothesis: since the thermal velocity of the molecules is of the order $c$ and the kinematic viscosity is of the order $lc$, where $l$ is the mean free path of the gas molecules,, we have $a\sim l/c^2$; an estimation based on heat conduction gives the same result. Combining this with the relation $p/\rho \sim c^2$, shows that

$$\delta \sim l,$$

i.e., the shock-wave thickness is of the order the mean free path of the molecules. However, in the continuum hypothesis, the mean free path is taken to be zero. It follows that the continuum equations alone cannot be strictly used to describe the internal structure of strong shock waves; in weak shock waves, $p_2-p_1$ can be made as small as possible to make $\delta$ large.

Rayleigh's solution

Two problems that were originally considered by Lord Rayleigh are given here.Template:R

Fluids with heat conduction and without viscosity (Pr′ → 0)

The problem when viscosity is neglected but heat conduction is allowed is of significant interest in astrophysical context due to presence of other heat exchange mechanisms such as radiative heat transfer, electron heat transfer in plasmas, etc.Template:R Neglect of viscosity means viscous forces in the momentum equation and the viscous dissipation in the energy equation disappear. Hence the first integral of the governing equations are simply given by

$$\begin{array}{cc}\rho u& ={\rho }_{0}D,\\ p+\rho u^2& =p_0+{\rho }_0{D}^2,\\ h+\frac{u^2}{2}-\frac{\lambda }{{\rho }_{0}D}\frac{dT}{dx}& =h_0+\frac{D^2}{2}.\end{array}$$

All the required ratios can be expressed in terms of $\eta$ immediately,

$$\begin{array}{cc}\frac{p}{p_0}& =1+\gamma M_0^2(1-\eta ),\\ \frac{T}{T_0}& =1+(1-\eta )(\gamma M_0^2\eta -1),\\ \frac{\lambda }{{\rho }_0{D}^3}\frac{dT}{dx}& =\frac{1}{2}\frac{\gamma -1}{\gamma +1}(1-\eta )(\eta -{\eta }_1).\end{array}$$

By eliminating $\eta$ from the last two equations, one can obtain equation $dT/dx=f(T)$, which can be integrated. It turns out there is no continuous solution for strong shock waves, precisely whenTemplate:R

$$M_0^2>\frac{3\gamma -1}{\gamma (3-\gamma )};$$

for $\gamma =1.4$ this condition becomes $M_0>\mathrm{1.2.}$

Fluids with viscosity and without heat conduction (Pr′ → ∞)

Here continuous solutions can be found for all shock wave strengths. Further, here the entropy increases monotonically across the shock wave due to the absence of heat conduction. Here the first integrals are given by

$$\begin{array}{cc}\rho u& ={\rho }_{0}D,\\ p+\rho u^2-\frac{4}{3}{\mu }^{\prime }\frac{du}{dx}& =p_0+{\rho }_0{D}^2,\\ h+\frac{u^2}{2}-\frac{4{\mu }^{\prime }}{3{\rho }_{0}D}u\frac{du}{dx}& =h_0+\frac{D^2}{2}.\end{array}$$

One can eliminate the viscous terms in the last two equations and obtain a relation between $p/p_0$ and $\eta$. Substituting this back in any one of the equations, we obtain an equation for $\eta (x)$, which can be integrated.

See also

• Taylor–von Neumann–Sedov blast wave