Problem Setup

The simulations evolve the equations of ideal MHD,

\begin{displaymath} {\partial\rho\over \partial t} + \nabla\cdot (\rho {\bf v}) = 0, \end{displaymath} (1)


\begin{displaymath} \rho {\partial{\bf v} \over \partial t} + (\rho {\bf v}\cdot... ...a \Phi + \left( {{\bf B}\over 4\pi}\cdot \nabla\right){\bf B}, \end{displaymath} (2)


\begin{displaymath} {\partial\rho\epsilon\over \partial t} + \nabla\cdot (\rho\epsilon {\bf v}) = -(P+{\mathcal Q}) \nabla \cdot {\bf v}, \end{displaymath} (3)


\begin{displaymath} {\partial{\bf B}\over \partial t} = \nabla\times\left( {\bf v} \times {\bf B} \right), \end{displaymath} (4)

where $\rho$ is the mass density, $\epsilon$ is the specific internal energy, ${\bf v}$ is the fluid velocity, P is the pressure, ${\mathcal Q}$ is an explicit artificial viscosity (Stone & Norman 1992a), ${\bf B}$ is the magnetic field vector, and $\Phi$ is the gravitational potential. We employ cylindrical coordinates, $(R,\phi,z)$, and work in the ``cylindrical disk'' limit which assumes a cylindrical gravitational potential and ignores vertical stratification. The gravitational potential $\Phi$ is a cylindrically-symmetric form of the Paczynski & Wiita (1980) pseudo-Newtonian potential $\Phi = - GM/(R-R_g)$. The equation of state is either isothermal, $P = c_s^2 \rho$, or adiabatic, $P=\rho\epsilon(\Gamma -1)$. Radiation transport and losses are omitted. Since there is no explicit resistivity or physical viscosity, the adiabatic gas can heat only by compression, or in nonadiabatic heating through the action of the artificial viscosity ${\mathcal Q}$.

These equations are solved using time-explicit Eulerian finite differencing. The numerical algorithm is that employed by the ZEUS code for hydrodynamics (Stone & Norman 1992a) and MHD (Stone & Norman 1992b; Hawley & Stone 1995). Time and length units are established by setting GM = Rg = 1. To ensure a substantial reservoir of matter and to minimize the influence of the outer boundary, the radial grid runs from R=1.5 to R=61.5 using 256 radial grid zones. These grid zones are either evenly spaced, or graded so as to concentrate resolution in the inner regions. The azimuthal angle $\phi $ runs from 0 to some fraction of $2\pi$ depending upon the simulation. In the vertical direction z runs over a somewhat arbitrary periodic length. Periodic boundary conditions are employed in $\phi $ and z. The radial boundary conditions are simple zero-gradient outflow conditions; no flow into the computational domain is permitted. The radial magnetic field boundary condition is set by requiring the transverse components of the field to be zero outside the computational domain, while the perpendicular component satisfies the divergence-free constraint.

In the pseudo-Newtonian potential the angular momentum of a circular orbit (here referred to as the Keplerian angular momentum) is

\begin{displaymath} \ell_{kep} = (GMR)^{1/2} {R\over R-R_g} = \Omega R^2. \end{displaymath} (5)

With GM = Rg = 1 the innermost marginally stable circular orbit rms is located at R=3; at this point $\Omega = 0.29$, and the orbital period is $P_{\rm orb} = 2\pi/\Omega = 21.8$. At the grid outer boundary the orbital period is almost 3000. The advantage of the large radius of the outer boundary is that many orbits can elapse in the inner regions of the disk before the outer boundary conditions become important to the simulation. This offers the possibility of establishing an accretion flow at the inner edge of the disk independent of the outer boundary condition (although still dependent of the initial conditions chosen for bulk of the disk).

In reporting results from the simulations many of the values will be averaged over both height z and angle $\phi $. For example, the averaged mass density is

\begin{displaymath} \langle \rho\rangle = {\int \rho R d\phi dz \over \int R d\phi dz}. \end{displaymath} (6)

The net local mass flux is

\begin{displaymath} \langle \dot M\rangle = \int -\rho v_R R d\phi dz, \end{displaymath} (7)

here defined to be positive for net accretion (inflow). Similar averages are constructed for pressure, magnetic energies, stresses and angular momentum. These averaged values are computed at specific time intervals. Where desired one can further average these quantities over time.

A drawback to reporting space and time averaged values is that they blur a very real and important property of the flow: it is highly variable. One way to measure the azimuthal fluctuation level at a given time is with the quantity

\begin{displaymath} {\delta \rho \over \rho}\left(R\right) = {1 \over \langle \r... ... \left[\rho - \langle \rho \rangle \right]^2 \right\}^{1/2} , \end{displaymath} (8)

where ${\langle \rho \rangle}$ is defined in (6).

The azimuthal structure in the disk is examined with a fourier transform of a vertically averaged value such as density, e.g.,

\begin{displaymath} \tilde \rho (R,m) = {1\over 2\pi} \int_{0}^{2\pi} \int_{-z_{max}}^{z_{max}} \rho(R,\phi,z) e^{i m\phi} dz d\phi . \end{displaymath} (9)

The power, $\vert\tilde\rho \vert^2$, is further averaged over radius within the interior of the active region in the disk. Similarly, the time-dependent behavior of a disk is examined with a fourier transform over time of azimuthally and vertically averaged values.

Table 1 lists the simulations carried out as part of this study. The models are identified by labels: CK stands for cylindrical Keplerian, HK for hydrodynamic Keplerian, and NK for Newtonian Keplerian, two models that use the standard 1/R Newtonian potential rather than the pseudo-Newtonian potential. The table also lists the computational domain, grid resolution, equation of state, initial field topology, and end time in code units.

Title Title Page   |   Introduction 1. Introduction   |   Results 3. Results