Becker–Morduchow–Libby solution is an exact solution of the compressible Navier–Stokes equations, that describes the structure of one-dimensional shock waves. The solution was discovered in a restrictive form by Richard Becker in 1922, which was generalized by Morris Morduchow and Paul A. Libby in 1949. ^{[1] [2]} The solution was also discovered independently by M. Roy and L. H. Thomas in 1944 ^{[3] [4]}The solution showed that there is a non-monotonic variation of the entropy across the shock wave. Before these works, Lord Rayleigh obtained solutions in 1910 for fluids with viscosity but without heat conductivity and for fluids with heat conductivity but without viscosity. ^[5] Following this, in the same year G. I. Taylor solved the whole problem for weak shock waves by taking both viscosity and heat conductivity into account. ^{[6] [7]}

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

${\displaystyle {\begin{aligned}{\frac {d}{dx}}(\rho u)&=0,\\\rho u{\frac {du}{dx}}+{\frac {dp}{dx}}-{\frac {4}{3}}{\frac {d}{dx}}\left(\mu '{\frac {du}{dx}}\right)&=0,\\\rho u{\frac {d\varepsilon }{dx}}+p{\frac {du}{dx}}-{\frac {d}{dx}}\left(\lambda {\frac {dT}{dx}}\right)-{\frac {4}{3}}\mu '\left({\frac {du}{dx}}\right)^{2}&=0,\end{aligned}}}$

where ${\displaystyle \rho }$ is the density, ${\displaystyle u}$ is the velocity, ${\displaystyle p}$ is the pressure, ${\displaystyle \varepsilon }$ is the internal energy per unit mass, ${\displaystyle T}$ is the temperature, ${\displaystyle \mu '=\mu +3\zeta /4}$ is an effective coefficient of viscosity, ${\displaystyle \mu }$ is the coefficient of viscosity, ${\displaystyle \zeta }$ is the second viscosity and ${\displaystyle \lambda }$ is the thermal conductivity. Instead of ${\displaystyle \varepsilon }$, it is convenient to work with the specific enthalpy ${\displaystyle h=\varepsilon +p/\rho .}$

Let us denote properties pertaining upstream of the shock with the subscript "${\displaystyle 0}$" and downstream with "${\displaystyle 1}$". The shock wave speed itself is denoted by ${\displaystyle D=u_{0}}$. The first integral of the governing equations, after imposing the condition that all gradients vanish upstream, are found to be

${\displaystyle {\begin{aligned}\rho u&=\rho _{0}D,\\p+\rho u^{2}-{\frac {4}{3}}\mu '{\frac {du}{dx}}&=p_{0}+\rho _{0}D^{2},\\h+{\frac {u^{2}}{2}}-{\frac {1}{\rho _{0}D}}\left(\lambda {\frac {dT}{dx}}+{\frac {4}{3}}\mu 'u{\frac {du}{dx}}\right)&=h_{0}+{\frac {D^{2}}{2}},\end{aligned}}}$

By evaluating these on the downstream side where all gradients vanish, one recovers the familiar Rankine–Hugoniot conditions, ${\displaystyle \rho _{1}u_{1}=\rho _{0}D}$, ${\displaystyle p_{1}+\rho _{1}u_{1}^{2}=p_{0}+\rho _{0}D^{2}}$ and ${\displaystyle 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

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 we can write ${\displaystyle h=c_{p}T}$, where ${\displaystyle c_{p}}$ is the specific heat at constant pressure. The third equation then becomes

${\displaystyle h+{\frac {u^{2}}{2}}-{\frac {4\mu '}{3\rho _{0}D}}\left({\frac {3}{4Pr'}}{\frac {dh}{dx}}+{\frac {1}{2}}{\frac {du^{2}}{dx}}\right)=h_{0}+{\frac {D^{2}}{2}}}$

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

${\displaystyle h+{\frac {u^{2}}{2}}=h_{0}+{\frac {D^{2}}{2}}.}$

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

${\displaystyle {\frac {p}{p_{0}}}={\frac {1}{\eta }}\left[1+{\frac {\gamma -1}{2}}M_{0}^{2}(1-\eta ^{2})\right]={\frac {\eta _{1}(\gamma +1)/(\gamma -1)-\eta ^{2}}{[\eta _{1}(\gamma +1)/(\gamma -1)-1]\eta }},\quad \eta _{1}={\frac {\gamma -1}{\gamma +1}}+{\frac {2}{(\gamma +1)M_{0}^{2}}},}$

where ${\displaystyle M_{0}=D/c_{0}}$ is the Mach number of the wave with respect to upstream. Combining this with momentum and continuity integrals, we obtain the equation for ${\displaystyle \eta (x)}$ as follows

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

We can introduce the coordinate

${\displaystyle \xi =\rho _{0}D{\frac {3(\gamma +1)}{8\gamma }}\int _{0}^{x}{\frac {dx'}{\mu (x')}}}$

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

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

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

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

As ${\displaystyle \xi \to -\infty }$ (or, ${\displaystyle x\to -\infty }$), we have ${\displaystyle \eta \to 1}$ and as ${\displaystyle \xi \to +\infty }$ (or, ${\displaystyle x\to +\infty }$), we have ${\displaystyle \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 ${\displaystyle T/T_{0}}$ is esaily to found to given by

${\displaystyle {\frac {T}{T_{0}}}=1+{\frac {\gamma -1}{2}}M_{0}^{2}(1-\eta ^{2})}$

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

${\displaystyle {\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 ${\displaystyle \gamma =1.4}$ and ${\displaystyle M_{0}=2}$. The notable feature is that the entropy does not monotonically increae 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 becuase of heat conduction. While the viscous dissipation associated with the term ${\displaystyle (du/dx)^{2}}$ always increases the entropy, heat conduction increases the entropy in the colder layers where ${\displaystyle ds/dx\sim -d^{2}T/dx^{2}<0}$, whereas it decreases the entropy in the hotter layers where ${\displaystyle ds/dx\sim -d^{2}T/dx^{2}>0}$.

Taylor's solution: Weak shock waves

When ${\displaystyle Pr'\neq 3/4}$, analytical solution is possible only in the weak shock-wave limit, as first shown by G. I. Taylor in 1910. ^[6] In the weak shock-wave limit, all terms such as ${\displaystyle p-p_{0}}$, ${\displaystyle \rho -\rho _{0}}$ etc., will be small. The thickness ${\displaystyle \delta }$ of the shock wave is of the order ${\displaystyle \delta \sim p_{1}-p_{0}}$ so that diffentiation with respect to ${\displaystyle x}$ increases the order smallness by one; e.g. ${\displaystyle 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 ${\displaystyle p(x)}$ is found to be related to the steady travelling-wave solution of the Burgers' equation and is given by ^[9]

${\displaystyle p(x)={\frac {1}{2}}(p_{1}+p_{0})+{\frac {1}{2}}(p_{1}-p_{0})\tanh {\frac {x}{\delta }}}$

where

${\displaystyle \delta ={\frac {8a\rho c^{4}}{(p_{2}-p_{1})\Gamma }},\quad a={\frac {\lambda }{2\rho c^{3}c_{p}}}\left({\frac {4}{3}}Pr'+\gamma -1\right),}$

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

${\displaystyle s(x)=s_{0}+{\frac {\lambda \Gamma }{8a\rho ^{2}c^{5}T}}\left({\frac {\partial T}{\partial p}}\right)_{s}(p_{1}-p_{0})^{2}{\frac {1}{\cosh ^{2}(x/\delta )}}.}$

Note that ${\displaystyle s(x)-s_{0}}$ is a second-order small quantity, although ${\displaystyle s_{1}-s_{0}}$ is a third-order small quantity as can be inferred from the above expression which shows that ${\displaystyle s=s_{0}}$ for both ${\displaystyle x\to \pm \infty }$. This is allowed since ${\displaystyle s(x)}$, unlike ${\displaystyle 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 ${\displaystyle c}$ and the kinematic viscosity is of the order ${\displaystyle lc}$, where ${\displaystyle k}$ is the mena free path of the gas molecules,, we have ${\displaystyle a\sim l/c^{2}}$; an estimation based on heat condcution gives the same result. Combining this with the relation ${\displaystyle p/\rho \sim c^{2}}$, shows that

${\displaystyle \delta \sim l,}$

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

Rayleigh's solution

Two problems that were origially considered by Lord Rayleigh is given here. ^[5]

Fluids with heat conduction and without viscosity

The problem when viscosity is neglected but heat conduction is allowed is of signifcant innterest in astrophysical context due to presence of other heat exhange mechasnisms such as radiative heat transfer, electron heat trasnfer in plasmas, etc. ^[8] 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

${\displaystyle {\begin{aligned}\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{aligned}}}$

All the required ratios can be expreses in terms of ${\displaystyle \eta }$ immediately,

${\displaystyle {\begin{aligned}{\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{aligned}}}$

By eliminating ${\displaystyle \eta }$ from the last two equations, one can obtaina equation ${\displaystyle dT/dx=f(T)}$, which can be integrated. It turns out there is no continuous solution for strog shock waves, precisely when ^[8]

${\displaystyle M_{0}^{2}>{\frac {3\gamma -1}{\gamma (3-\gamma )}};}$

for ${\displaystyle \gamma =1.4}$ this condition becomes ${\displaystyle M_{0}>1.2.}$

Fluids with viscosity and without heat conduction

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.

See also

• Taylor–von Neumann–Sedov blast wave

References