Asymptotic Steady State Solution to a Bow Shock with an Infinite Mach Number


The problem of a cold gas flowing past an unmovable object is considered. It is shown that at large distances from the obstacle the shock front forms a parabolic body of revolution. The interior of the shock front is obtained by solution of the hydrodynamic equations in parabolic coordinates. The results are verified with a hydrodynamic simulation. Finally, relations to astrophysical bow shocks and other analytical works on oblique shocks is discussed.


Bow shocks occur when a supersonic flow encounters an obstacle. A prominent example from astrophysics is planetary bow shocks(Treumann 2008). Another common astrophysical scenario where shock waves emerge is when supernova remnants engulf a dense molecular cloud (McKee 1975). A bow shock has even been observed around a star that moves supersonically relative to the ISM (Noriega-Crespo 1997).

Bow shocks have been studied extensively, both theoretically (Wilkin 1996, Farris 1994) and numerically (Mohamed 2012, Miceli 2006). However, all studies were carried out under the assumption of finite Mach number. In this work we rather assume that the incoming matter is cold, so for every finite velocity its Mach number would be infinite. The asymptotic shape far from the obstacle is qualitatively different between the two. In the case of a finite Mach number \(M\), far away from the obstacle the length scale of the obstacle becomes irrelevant, and the shock front coincides with the Mach cone, i.e. a cone with an opening angle \(\alpha=\sin^{-1}\frac{1}{M}\) (LANDAU 1987). However, as we shall later see, in the case of an infinite Mach number, the shock front is parabolic solid of revolution. For a spherical obstacle of radius \(R_{o}\), the shock front loci in cylindrical coordinates is \(z=\xi\frac{r^{2}}{R_{o}}\), where \(\xi\) is a dimensionless constant. This means that no matter how far from the obstacle, the obstacle’s length scale is never negligible. For a very large Mach number, we should expect that the transition from former to the later should occur where the two curves intersect

\begin{equation} r_{t}\approx R_{o}\sqrt{M^{2}-1}\\ \end{equation}
\begin{equation} z_{t}\approx R_{o}\left(M^{2}-1\right)\\ \end{equation}

The solution described here is asymptotic, which means it is only valid at distances much larger than the obstacle. This means that the solution can only manifest itself high Mach number flows, and even then, it would only occupy the region \(zM^{2}>z\gg R_{o}\).

This paper is organized as follows. In section 2 we describe the complete mathematical formulation. In section 3 we present validation of our analytic results with a numerical simulation. Finally, in section 4, we discuss the results.

Mathematica Formulation

Spherical parabolic coordinates are given by

\begin{equation} x=\sigma\tau\cos\phi\\ \end{equation}
\begin{equation} y=\sigma\tau\sin\phi\\ \end{equation}
\begin{equation} z=\frac{1}{2}\left(\tau^{2}-\sigma^{2}\right)\\ \end{equation}

In these coordinates, the steady state, azimuthally symmetric hydrodynamics equations take the following form. The conservation of mass is

\begin{equation} \frac{\partial}{\partial\sigma}\left(\sqrt{\sigma^{2}+\tau^{2}}\sigma\tau v_{\sigma}\right)+\frac{\partial}{\partial\tau}\left(\sqrt{\tau^{2}+\sigma^{2}}\sigma\tau v_{\tau}\right)=0\\ \end{equation}

The conservation of entropy is

\begin{equation} v_{\sigma}\frac{\partial s}{\partial\sigma}+v_{\tau}\frac{\partial s}{\partial\tau}=0\\ \end{equation}

where \(s=\ln p-\gamma\ln\rho\) is the specific entropy. The conservation of momentum in each direction is

\begin{equation} v_{\sigma}\frac{\partial v_{\sigma}}{\partial\sigma}+v_{\tau}\frac{\partial v_{\sigma}}{\partial\tau}+\frac{v_{\tau}\left(v_{\sigma}\tau-v_{\tau}\sigma\right)}{\tau^{2}+\sigma^{2}}+\frac{1}{\rho}\frac{\partial p}{\partial\sigma}=0\\ \end{equation}
\begin{equation} v_{\sigma}\frac{\partial v_{\tau}}{\partial\sigma}+v_{\tau}\frac{\partial v_{\tau}}{\partial\tau}-\frac{v_{\sigma}\left(v_{\sigma}\tau-v_{\tau}\sigma\right)}{\tau^{2}+\sigma^{2}}+\frac{1}{\rho}\frac{\partial p}{\partial\tau}=0\\ \end{equation}

One of the momentum equations can be replaced by Bernoulli’s equation

\begin{equation} \frac{1}{2}v_{in}^{2}=\frac{1}{2}v_{\sigma}^{2}+\frac{1}{2}v_{\tau}^{2}+\frac{1}{\gamma-1}c^{2}\\ \end{equation}

where \(c\) is the speed of sound. We found it most convenient to express the equations in terms of the density \(\rho\), the sound speed \(c\), and the two Mach numbers \(m_{\sigma,\tau}=\frac{v_{\sigma,\tau}}{c}\)

We assume that the shock front coincides with a curve on which \(\sigma=\sigma_{s}\) is constant. The Rankine Hugoniot boundary conditions at the shock fronts are

\begin{equation} \rho_{s}=\frac{\gamma+1}{\gamma-1}\rho_{a}\\ \end{equation}
\begin{equation} c_{s}=\frac{\sqrt{2\gamma\left(\gamma-1\right)}}{\gamma+1}v_{in}\frac{\sigma_{s}}{\sqrt{\tau^{2}+\sigma_{s}^{2}}}\\ \end{equation}
\begin{equation} m_{\sigma s}=-\sqrt{\frac{\gamma-1}{2\gamma}}\\ \end{equation}
\begin{equation} m_{\tau s}=\frac{\gamma+1}{\sqrt{2\gamma\left(\gamma-1\right)}}\frac{\tau}{\sigma_{s}}\\ \end{equation}

So we have a well defined boundary value problem. However, since it involves two coordinates, our approach so far has no benefit over the hydrodynamic equations in cylindrical coordinates. At very high values of \(\tau\gg\sigma_{s}\), the equations can be simplified. We assume that the variables vary with \(\tau\) as the shocked values do, i.e. \(\rho\) is independent of \(\tau\), \(c\propto\frac{1}{\tau}\), \(m_{\sigma}\) is independent of \(\tau\) and \(m_{\tau}\propto\tau\). Using this approximation, \(\tau\) can be eliminated from the hydrodynamic equations, and the problem reduces to a set of ordinary differential equation in \(\sigma\). These equation can be numerically integrated, and we can obtain curves for \(\rho\), \(\tau\cdot c\), \(m_{\sigma}\) and \(\frac{m_{\tau}}{\tau}\). We performed the complete analysis in a mathematica notebook.


To verify our results, we ran a simulation in the RICH code(Yalinewich 2015) (version 2 branch). We have also released a complete listing of the source code for the simulation. We were working in 2D cylindrical coordinates, assuming azimuthal symmetry. The boundaries of the computational domain were: \(1.8>z>-0.2\) and \(0.5>r>0\). The wind velocity was 1, and it was directed in the positive \(z\) direction. The obstacle was a circle of radius 0.01, located at the origin \(\left(0,0\right)\). We ran the simulation to time \(t=10\), and compared the numerical results to our analytic predictions. First, we fit the shock front to a parabola (figure \ref{parabolic_shock_front}). As can be seen from the figure, the deviation from a parabolic is minuscule. Next, we compared the numerical profiles inside the bow shock to the analytical prediction (figure \ref{profiles}). In order to do that, we interpolated the variables on the curve \(\tau=1.5\). The agreement between the two is reasonable. The deviations are due to two main reasons. The first is that since we are working in cylindrical coordinates, the relative error at each radius is equal to the ratio between the cell’s width and the radius. The second is that the slope of the density is very steep next to the shock front, so a very fine resolution is required to resolve it.