self-gravitating orbiting cylinder test for spectral solver
This part is adapted (and translated) from the PhD-thesis from Jung (2016, [jung2016] ). It contains two different types of tests with cylinders. Both originally come from Chan (2006, [chan2006] ):
- A static simulation of the gravitational potential in a
- razor-thin disk.
- gaussian distribution in vertical direction.
- A dynamic evolution of the cylinders.
Density-Potential Pairs
This part compares the solution of the potential by the gravitational spectral solver with an existing analytical result. Since the module only calculates the potential within the simulation area, the potential produced by the boundaries need to be neglectable. However, the potential from outside the potential can calculated with additional gravitational modules.
Assuming a density distributions around a central point
\[ \Sigma_{(r_k,\varphi_k)}(r,\varphi) = \frac{1}{2 \pi \sigma^2} \exp{\left( - \frac{R_k^2}{\sigma} \right)}, \]
with \( \sigma \) a measurement for the width of and \( R_k \) the distance of a coordinate point form the center of the distribution. For a razor thin disk, this yields the potential (see [chan2006]):
\[ \Phi_{(r_k, \varphi_k)} = - \frac{G}{\sigma} \left(I_0(y_k)K_1(y_k) - I_1(y_k)K_0(y_k)\right), \]
where \( y_k = \frac{R_k}{2 \sigma} \). \( I_n \) and \( K_n \) are the modified Bessel functions of first and second kind with integral order (see functions ).
The density distribution in our example is
\[ \Sigma(r,\varphi) = 2 \Sigma_{(1,10^{-3})}(r,\varphi) + 0.5 \Sigma_{(1,\pi+10^{-3})} (r,\varphi) + \Sigma_{(0.9,\frac{3}{4}\pi)} (r,\varphi), \]
with the according potential
\[ \Phi(r,\varphi) = 2 \Phi_{(1,10^{-3})}(r,\varphi) + 0.5 \Phi_{(1,\pi+10^{-3})} (r,\varphi) + \Phi_{(0.9,\frac{3}{4}\pi)} (r,\varphi). \]
Additionally, the implementation is tested for a vertical gaussian distribution, with a scale height similar to \( \sigma =0.1 \). with the same setup as above. This gives us a solution for the potential
\[ \Phi(r_k, \varphi_k) = - \frac{1}{R_k} \mathrm{erf}\left( \frac{R_k}{\sqrt{2} \sigma} \right). \]
Below the relative error for the potential is shown, solved for a razor-thin disk (left) and a vertical gaussian distribution (right). For the first case we see a maximum error of \( 10^{-3} \) at steep gradients and \( 10^{-5} \) for the second case.
Self-Gravitating Rotating Cylinders
A time-dependent test from [chan2006], which generates two cylinders with gaussian density distribution in polar grid placed opposite towards each other. A background field adds a rigid rotation. The thereby arising centrifugal force has to be balanced from the self-gravitation of the mass-distribution. This is a very sophisticated setup for the accuracy of the acceleration by self-gravitation and the general angular-momentum conservation. The setup is
\[ \varrho(r,\varphi) = \frac{10^{-2}}{\pi \left( r_{\mathrm{max}}^2 - r_{\mathrm{min}}^2 \right)} + 0.99 \left( \frac{1}{2 \pi \sigma^2} \exp{\left(- \frac{R_1}{2 \sigma^2} \right)} + \frac{1}{2 \pi \sigma} \exp{\left(- \frac{R_2}{2 \sigma^2} \right)}\right) \\ p(r,\varphi) = \frac{10^{-2}}{\pi \left( r_{\mathrm{max}}^2 - r_{\mathrm{min}}^2 \right)} + \frac{G}{2 \pi \sigma^2} \left(\mathrm{Ei}\left(- \frac{R_1}{\sigma^2} \right) - \mathrm{Ei}\left(- \frac{R_1}{2 \sigma^2} \right) + \mathrm{Ei}\left(- \frac{R_2}{\sigma^2} \right) - \mathrm{Ei}\left(- \frac{R_2}{2 \sigma^2} \right)\right), \]
and \( \mathbf{v} = r \mathbf{e}_{\varphi} \). The function \( \mathrm{Ei}(x) := \int_{-\infty}^x \frac{\exp{t}}{2} \mathrm{d}t \) is the exponential integral, \( R_1 \), \( R_2 \) are the radii from the center of the cylinders and \( \sigma = 0.1 \) a measure for spreading of the cylinders.
| Additional Settings: |
|---|
| angular velocity: \( \Omega = 1.0 \) |
| heat capacity ratio: \( \gamma = \frac{5}{3} \) |
| resolution: \( N_r \times N_{\varphi} = 256 \times 1024 \) |
| simulation time: \( t_{\mathrm{sim}} = 100 \) |
| gravitational constant: \( G = 1.0 \) |
| \( r_{\mathrm{min}}, r_{\mathrm{max}} = 0.2, 1.8 \) |
The images below show the results at \( t = 0, 100 \). The cylinders retain their shape very well. Within \( 100 \) orbits of the cylinders in the rotating frame there is a loss of \( 0.06 \% \) of total angular momentum. Compared with [chan2006], we have a \( 200 \times \) better angular momentum conservation in one time step.
References:
- [chan2006] Chan, Chi-kwan; Psaltis, Dimitrios; Özel, Feryal, 2006 "Spectral Methods for Time-dependent Studies of Accretion Flows. II. Two-dimensional Hydrodynamic Disks with Self-Gravity" http://adsabs.harvard.edu/abs/2006ApJ...645..506C
- [jung2016] Jung, Manuel. “Multiskalensimulationen von Schwarzlochakkretion in Balkengalaxien.” PhD-thesis, Christian-Albrechts Universität Kiel, CAU Kiel, 2016. http://macau.uni-kiel.de/receive/dissertation_diss_00019153.
