Astron. Astrophys. 319, 547-560 (1997)

2. Initial model and numerical methods

2.1. The initial model

For our test model calculations we have chosen an idealized set of initial conditions corresponding to the standard isothermal test case (Boss & Bodenheimer 1979). The initial model consists of a , uniform density ( g cm^-3), uniformly rotating ( s^-1) sphere of radius cm, composed of pure molecular hydrogen at a temperature of 10 K. An azimuthally varying perturbation of the form

with mode and amplitude , is added to the background uniform density distribution. During the evolution, the cloud is assumed to remain isothermal at a constant temperature equal to the initial value of 10 K. With these initial parameters, the ratios of thermal and rotational to the cloud's self-gravitational energy are, respectively, and .

Previous calculations of the standard test case using first-order accurate FD (Boss & Bodenheimer 1979; Bodenheimer & Boss 1981; de Felice & Sigalotti 1992), SPH (Gingold & Monaghan 1981; ML; Bate et al. 1995) and second-order accurate FD codes (MB; Burkert & Bodenheimer 1993) have shown a fairly good agreement that the outcome of the first evolution is the formation of a binary system. Nevertheless, basic differences between the FD and the SPH methods (ML) still remain to be clarified. In order to address this point we have rerun the standard test model using two new distinct codes. In the sequel we shall refer to them as code FD and code TREESPH.

2.2. Code FD

Code FD is a fully 3D radiation hydrodynamic scheme written in Eulerian spherical coordinates . It is a descendant of a previously implemented FD code (Sigalotti 1993). For the purposes of this paper, we have used a simplified version of the code, which solves the equations of continuity and motion in conservation-law form for a compressible, inviscid, non-magnetic fluid, coupled with the Poisson equation for the gravitational potential. This set of equations is closed by specifying a pressure relation of the form , where is the mass-density. The above equations are solved using FD methods on a radially moving spherical grid along with a multi-step solution procedure for the hydrodynamic equations in which the source and the convective terms are evaluated in separate sub-steps. The time integration of the hydrodynamic equations is made explicitly and is optionally first- or second-order accurate. Temporal second-order accuracy is achieved by means of a predictor-corrector treatment of the source and convective terms. The spatial differences are based on a generalization of the volume-centred discretization method of Mönchmeyer & Müller (1989), due to the incorporation of grid motion. In its standard form, the volume-centred method results in a conservative, second-order preserving scheme in curvilinear coordinates, as was demonstrated through convergence testing. A similar spherical coordinate-based code was previously implemented by Boss & Myhill (1992) who showed that with the use of volume-centred differences, spatial second-order accuracy can be achieved in actual 3D collapse calculations - at least over the range of spatial resolution where 3D runs can be made in practice - without the need of using directional splitting techniques for treatment of the convective terms, as advocated by Finn & Hawley (1989). In collapse calculations with Boss & Myhill's (1992) code, adequate central resolution is maintained by means of a grid redefinition procedure, and so most of the good properties of the volume-centred method are preserved in a direct manner. Here, spatial resolution is maintained by allowing the radial grid to move with matter, and hence in order to preserve the second-order properties of the volume-centred method, a new FD formulation is used to correctly design the convective terms. It consists of appropriate geometrical and temporal corrections which make the hydrodynamical transport algorithm numerically invariant under grid motion, thereby yielding numerical solutions that are free of a class of systematic errors introduced when standard FD replacements are used on moving grids. This new approach is seen to result in a high degree of accuracy on the pressureless collapse test case. Consistent advection (Norman, Wilson & Barton 1980) is also implemented in the definition of the convective fluxes, and transport of mass is performed using a generalization to spherical coordinates of the van Leer (1977) monotonic interpolation method. This results in good local conservation of angular momentum for axisymmetric flows. An artificial viscosity, based on the tensor formulation of Tscharnuter & Winkler (1979), is added to improve stability against shock formation. A more detailed description of the code and results from a variety of tests is under preparation, and will be presented in a forthcoming paper (Sigalotti 1996).

The standard test model calculation is made using a spherical grid consisting of initially uniformly distributed radial points (including the origin ), fairly equidistant points for (including the rotation axis , ), and uniformly spaced points along the azimuthal coordinate for . This resolution corresponds to a total number of effective cells filling the entire computational volume. Reflection symmetry through the equatorial plane is assumed so that only the top hemisphere is represented by the calculation ( cells for ). At the surface of the spherical grid, we enforce a constant-volume boundary condition by setting a zero radial velocity. Thus adequate inner resolution is achieved by allowing all interior radial points to contract during the main collapse phase while keeping the external radius fixed in space and time. For this model, the code was able to handle the gradients which formed in the late evolution without the need of the artificial viscosity and so it was never used throughout the calculation. With this choice of the numerical techniques, the standard test model evolved for time steps corresponding to about 180 hr of CPU time on a CRAY J916/6 supercomputer.

In order to test the code at much lower initial grid resolution, the standard collapse model was rerun with radial points, points for (i.e., points for because of the assumption of equatorial symmetry), and points for . This resolution amounts to a total number of 55121 effective cells which is a factor of lower than the total number of cells allowed in the high-resolution calculation. Hereafter we shall refer to the low- and high-resolution runs as model FDl and model FDh, respectively.

2.3. Code TREESPH

The version of code TREESPH used for the calculations of this paper is a descendant of the scheme originally formulated by Hernquist & Katz (1989). The code combines the method of SPH developed by Lucy (1977) and Gingold & Monaghan (1977) with the hierarchical tree algorithm of Hernquist (1987) for the calculation of the gravitational potential. Basically, the SPH method is a grid-free Lagrangian scheme in which the mass-density as well as all other physical properties of the fluid are determined from the properties of a quasi-random sample of a finite number of fluid elements (particles) through kernel estimation. In TREESPH, the mean value of a field quantity, say the mass-density , within a volume element of extent h, is fully determined by means of a symmetrized smoothing procedure using the spherically symmetric spline kernel proposed by Monaghan & Lattanzio (1985). The parameter h, known as the smoothing length, gives the local spatial resolution. A temporally and spatially varying h is used so that regions of the flow with relatively low and high density of particles can be solved with comparable accuracy. Since the mean density is evaluated through kernel estimates, the continuity equation is satisfied automatically to second order in h and so it need not be solved. The positions and velocities of particles are determined by solving the Euler's equation

where p, , and denote, respectively, the gas pressure, the gravitational potential, and the artificial viscous acceleration associated with the particles. Equations (2) and (3) are complemented by a further equation for the specific internal energy. However, this equation is not used in the calculation of the standard test case because of the assumption of a barotropic equation of state . The form of the artificial viscosity and the symmetrized smoothed representation of the pressure gradient are as given in Hernquist & Katz (1989). The gravitational acceleration term is calculated using the TREECODE of Hernquist (1987), which is an optimized Fortran implementation of the tree method of Barnes & Hut (1986). The details of this code are described in full in Hernquist (1987) and references therein. The time integration of Eqs. (2) and (3) is made explicitly using a second-order accurate leapfrog scheme. In doing so, each particle is allowed to have its own time step which is limited to maintain stability according to a modified form of the Courant condition. In TREESPH, the individual particle time steps are further constrained to be a fraction of the largest system time step , which enters as an input parameter.

In this standard formulation, however, the use of time-varying smoothing lengths can lead to a violation of the fluid conservation laws in certain extreme situations. This was first made evident by Hernquist (1993) who showed that during a head-on collision between two identical polytropes, errors in the conservation of the total energy or entropy may occur at a level of %, depending on whether the internal energy or the entropy equation is integrated. A new SPH formulation, which results in a conservative scheme even with the use of varying smoothing lengths, has been introduced by Nelson & Papaloizou (1993, 1994). Basically, this new formulation differs from that of Hernquist & Katz (1989) in that correction terms that account for the variability of the smoothing lengths, the terms, are added to the usual smoothed representation of the particle equations of motion and internal energy, with h being now prescribed to be a function of the interparticle distances. The inclusion of the terms had the effect of reducing the non-conservation errors to % in a head-on collision test, irrespective of whether the internal energy or the entropy equation is integrated. An alternative SPH formulation, which also results in improved conservation properties for adiabatic flows, was independently proposed by Steinmetz & Müller (1993). In the case of barotropic (isothermal) flows, entropy conservation is no longer required as entropy is not conserved in isothermal systems, and only the particle equations of motion need be corrected. However, isothermal collapse test calculations by Nelson & Papaloizou (1993) showed that including the terms in the equations of motion, although it resulted in improved total energy conservation, had little effect on the qualitative outcome of the results. Therefore, for the purposes of the present study we omit the use of the terms and rely on the standard SPH formulation of Hernquist & Katz (1989). Code TREESPH has been used for a variety of other tests, and results from them together with the description of a modified implemented version of the SPH part of the code will be presented in an oncoming paper (Klapp 1996).

We present two separate model calculations of the standard test case with code TREESPH which differ only in the number of initial particles (N =7123 and N =13997). Similar sets of particles were used by ML in their SPH calculations for the same system. The initial protostellar cloud is set up by first defining a Cartesian box of sides , where R is the initial radius of the cloud. The box is then filled up with regular cubic cells, each of volume , and the origin of the Cartesian grid is made to coincide with the vertex common to the eight central cells (geometrical centre of the square box). An approximate sphere of radius R is then obtained by placing the particles at the cell vertices at distances from the origin. Thus the region surrounding the sphere is a vacuum. The azimuthal density perturbation (1) is added by assigning to each particle q a mass given by

where is the total mass of the sphere, N is the total number of particles inside the spherical volume, and the parameters a and m are as given in § 2.1. The initial cloud is assumed to be at rest except for rotation and so the particle at is assigned an initial velocity given by

where is the initial uniform angular velocity defined in § 2.1. The smoothing length of particle q, say , is adjusted such that it is allowed to interact with a constant number of nearby particles lying within a radius . The largest system time step is chosen to be s for both calculations. With particles the code requires about 7 s/cycle and with it requires about 14 s/cycle on a CRAY YMP-4 supercomputer.

© European Southern Observatory (ESO) 1997

Online publication: July 3, 1998
helpdesk@link.springer.de