User:Lothar.brendel/Gravity driven flow into a perfect spherical sink

Isothermal Bondi flow

$$

The stationary continuity equation for spherical symmetry and $\overrightarrow{u}=-u{\overrightarrow{e}}_r$ reads $$\begin{array}{rl}0& =\mathrm{\nabla }\cdot (\rho \overrightarrow{u})=-{\partial }_r(r^2\rho u)/r^2\\ \iff \frac{\dot{M}}{4\pi }& =r^2\rho u=\text{const.},\end{array}$$ $\dot{M}$ being the "anihilation rate" of the spherical sink. It complements the stationary Euler equation (which is invariant wrt $u\to -u$) $$\begin{array}{rl}u{\partial }_{r}u+\frac{GM}{r^2}& =-\frac{{\partial }_{r}P}{\rho }=-\frac{\partial P}{\partial \rho }\frac{{\partial }_r\rho }{\rho }\\ & =-c^2{r}^{2}u{\partial }_r(r^{2}u{)}^{-1}=c^2\frac{2ru+r^2{\partial }_{r}u}{r^{2}u}\\ & =c^2(\frac{2}{r}+\frac{{\partial }_{r}u}{u})\\ \iff \frac{u}{c}\frac{{\partial }_{r}u}{c}-\frac{{\partial }_{r}u}{u}& =\frac{2}{r}-\frac{2R_{\text{Bo}}}{r^2},\end{array}$$ which suggests $R_{\text{Bo}}=GM/(2c^2)$ as natural length unit and $R_{\text{Bo}}/c$ as natural time unit, i.e. we continue with $$\begin{array}{rl}(u-1/u){\partial }_{r}u& =\frac{2}{r}-\frac{2}{r^2}\\ \iff {\partial }_{r}u& =\frac{2u(r-1)}{r^2(u^2-1)}\underset{r\to \mathrm{\infty }}{\overset{u\to 0}{\to }}-\frac{2u}{r},(1)\end{array}$$ where in the last step we have put in the "natural" boundary condition of $u$ vanishing for $r\to \mathrm{\infty }$.

The subsonic solution

The asymptotic solution of (1) is $$u_{\mathrm{\infty }}(r)=\frac{u_{\wedge }}{r^2},$$ from which the real solution deviates strongly for $r\gg ̸1$. As long as $u<1$, ${\partial }_{r}u$ will be negative for $r>1$, i.e. $u$ will increase for decreasing $r$, but slower than $u_{\mathrm{\infty }}(r)$. We assume $u<1$ for $r\ge 1$, which implies ${\partial }_{r}u>0$ for $r<1$, i.e. $u$ decreases again for decreasing $r$, making $u(1)=u_{\text{max}}<1$ a local maximum. For $u_{\text{max}}\to 1$, the slope switching around $r=1$ gets more localised, rendering the maximum more and more cusp-like.

For $u\to 0$, (1) gets linear again: $$\begin{array}{rl}{\partial }_{r}u& =\frac{2u(r-1)}{r^2}\end{array}$$ It yields a solution vanishing at zero together with all its derivatives: $$\begin{array}{rl}{u}_0(r)& =\frac{u_{\vee }}{r^2}\exp (-2/r)\end{array}$$ The amplitudes $u_{\vee }$ and $u_{\wedge }$ are non-linear functions of $u_{\text{max}}$, varying linearly with $u_{\text{max}}\ll 1$ and saturating both at $u_{\vee }\approx u_{\wedge }\approx 4.48$.

Larger amplitudes would lead to $u$ reaching unity already away from $r=1$, which corresponds to a diverging slope ${\partial }_{r}u$. Hence, faster stationary, subsonic solutions with spherical symmetry do not exist.

Further reading:

• Bondi 1952
• McCrea 1956 (solutions with shocks)
• Keto 2020