Bondi Accretion

Am working on a 3D simulation of Bondi Flow onto a central sink particle and studying the accretion properties using the Krumholz accretion algorithm. From ​Bondi (1952), we have the following quantities (in order: Bondi radius, accretion rate, sonic radius, nondimensional density and velocity):

R_{BH} = \frac{GM_*}{C_\infty^2}

\dot{M}_{BH} = 4 \pi \lambda_{crit}(G M_*)^2c_\infty^{-3}\rho_{\infty}

R_s=\frac{5-3\gamma}{4}\frac{GM_*}{C_{\infty}^2}

z\_=\frac{\rho}{\rho_{\infty}}

y\_=\frac{v_r}{C_\infty}

In the above, \lambda is a numerical coefficient which controls the strength of the accretion flow. At its critical value, \lambda=\lambda_{crit}, \dot{M} is strongest. In the code, \lambda_{crit} is solved for using Eqn. 18 in Bondi's paper. With respect to these 'Bondi' variables, the sims used the following (in computational units), unless otherwise noted:

T_{\infty}	~1
\rho_{\infty}	~1
M_*	.14318E+13

with the computational domain set up using:

Lx=Ly=Lz	20
Mx=My=Mz	64
MaxLevel (AMR)	2
Effective resolution (dx)	.15625
Sink dx = 4dx	.625

Note, the sink dx is equivalent to the accretion volume radius (r_{acc}).

Research Log:

3/18/18, \gamma=7/5 case

Am running the \gamma=7/5 case because the sonic point in this case is located far beyond the accretion radius (R_s=6>>R_{acc}=.625). This will enable us to see the generation of accretion shocks in the sim and check the behavior of the accretion algorithm under conditions of supersonic infall. We expect that the accretion algorithm should be fairly 'well-behaved' in the supersonic regime, as any pressure errors generated in the accretion kernel will be unable to propagate upstream and effect the hydrodynamical solution beyond the accretion volume.

For this run, I made the box bigger than the fiducial case above, but used the same effective resolution. The bigger box now encompasses the bondi radius (R_{BH}=30.5). The following table lists this sims params:

\gamma	1.4
\lambda_{crit}	.625
C_{\infty}	1.183
R_{BH}	30.5
R_{s}	6.1
\dot{M}_{BH}	8646.9
t_{sim}	60
L_x=L_y=L_z	80
mx=my=mz	128
dx_{min}	.15625

The final time of the simulation is ~2 sound crossing times, using the sound speed of the gas at the edge of the domain (C_{amb}=1.283) and the distance between this edge and the sink particle (r=40).

This case no longer forces an inner spherical boundary in the sim where the Bondi solution is copied. Instead the solution is sampled for each cell within the computational volume (including ghost zones) at t=0. For each subsequent cycle, the Bondi solution is pasted everywhere r>20.

Results

The non-dimensional profiles for this run line up nicely with Bondi's Fig. 3 (curves III and IV):

(Note — I did this lineout from the center of the first zone in the octant surrounding the sink, to the zone center at the edge of the boundary. Aside from doing a spherical average over shells, I think this will provide the most accurate results from Visit, as I was getting somewhat wonky results sampling along cell boundaries — probably due to Visit's interpolatation method)

Calculating the mass flux through spherical shells surrounding the sink seems to agree well with the predicted value of \dot{M}_{BH}=8646.9. I tried at 3 different radii: r=5<R_s, r=8>R_s, r=25>OBS, where OBS is the spherical boundary in the grid beyond which the Bondi solution is pasted into cells:

Despite these agreements with theory, and the fact that the sonic point lies well beyond the accretion kernel, there is a spurious outflow being generated within the accretion kernel. Contrary to intuition, this outflow propagates beyond the accretion kernel into the grid, disrupting the upstream hydrodynamic solution.. This happens very quickly in the simulation (\sim1/30 ~t_{sc, amb}). Here is a snapshot of density illustrating this effect:

A time sequence of the first few frames shows that this occurs due to an increase in pressure inside of the accretion kernel, which blows outward into the grid:

The upwind hydro solution is intact until contact with the outward traveling shock generated in the accretion volume:

A zoom-in of the pressure inside the accretion kernel with overlaid velocity vectors shows the generation of an outward traveling shock:

It is also notable how symmetric the pressure is within the kernel (in this plane as well as the other coordinate planes). The relative orientation of the velocity vectors is also interesting. Despite the appearance of a subsonic, post-shock flow field by frame 1, within the kernel (central 8x8 zones) the velocity vectors point inward to the sink particle predominantly along the coordinate axes.. so much so that the flow field somewhat resembles a convective instability.. A similar plot of density:

The various profiles within the kernel seem to be identical regardless of viewing plane, until the 6th frame:

The mass flux through a shell with r=r_{acc}, as well as onto the particle is much lower (and variable) than the predicted \dot{M} (note the last time query corresponds to the 0th frame, and the first the 8th — the points between are equally spaced between these):

A movie of the pressure with overlaid velocity vectors and contours of the sonic surfaces​ shows the generation of multiple shock layers, both at the surface of the kernel, as well as inside the kernel.

The same movie, but zoomed out​, shows that the sonic radius of the Bondi solution is correctly located at R_s=6, and doesn't move over time.

Finally, we have lineout movies of pressure​, density​, and vrad​. As these movies show, the Bondi solution is unaltered beyond r>20, as to be expected from the boundary conditions of the sim. Surprising, however, is that despite the density falling below its initial value in the kernel (due to accretion), the pressure spikes there seemingly artificially (rather than due to a build of material in the kernel).

3/18/18, \gamma=1.66, steady-state case

The case of \gamma=5/3 flow is an extreme limit for the solution space. This value of \gamma has a few interesting features:

In other words, when \gamma=5/3 exactly, we have no accretion flow.

In addition to the Bondi/flow params cited in the top section of this page, the sim was set up using the following (in CU):

\gamma	1.66
\lambda_{crit}	.2526
C_{\infty}	1.291
R_{BH}	25.63
R_{s}	.0128
\dot{M}_{BH}	2692.17

Note, the sonic radius is not resolved by ~ factor of 20 in this simulation, as seen when checking dx quoted above.

In this sim, the Bondi solution was initialized everywhere in the grid at t=0 (except within a small inner radius, described below). In subsequent timesteps, the Bondi solution is re-pasted into the ghost zones and a small inner spherical region at the origin. A sink particle is initialized at the origin with mass given above (M_{sink}=M_*). The particle can accrete gas following the Krumholz prescription. The only source of gravity in the sims is the point gravity object associated with the particle (self-gravity is turned off).

The simulation setup reproduces the correct nondimensional profiles. Here's a comparison of the t=0 nondimensional density and velocity profiles astrobear calculates for the bondi module (left), compared to the Bondi solutions for \gamma=5/3 (curves II and III for the case of \lambda=\lambda_{crit}).

The profiles are cut-off within an inner radius of r=.875 to avoid extremely high speeds there (am going to get rid of this in the next round of sims). Sampling the mass flux across a spherical shell less than this radius doesn't match up with the theoretical prediction of \dot{M}=2692 since the solution is getting stepped on there. Sampling the mass flux across a shell larger than this, however, produces agreement (2689 compared to 2692).

Since the solution is getting stepped on within some small inner radius, can't meaningfully check the behavior of the accretion algorithm and any spurious waves it might be generating there. Instead, will be removing this inner boundary in the next run. These next sims will then be able to test the effect of the sonic radius being inside the accretion radius as opposed to outside. To test this will do a resolution study on \gamma=7/5 flow, where the sonic point is \sim .2R_{BH}.

As the following images show, this set of ICs produces a steady-state solution.. This is due to the solution being stepped on within the accretion volume, which would likely be the first place any deviation from the similarity solutions would occur due to spurious pressure waves that the accretion algorithm might produce in such strongly subsonic flow as this. The attached module files that produce these files are thus called "*_steadystate*"

Mesh
Radial velocity
Radial mach
Velocity field
Isosurface (r=rsink)
Mass flux (r=rsink)

Library