Hydrodynamic modeling of self-gravitating astrophysical objects on tetrahedral meshes

The paper proposes a new numerical method for solving the equations of gravitational hydrodynamics on a tetrahedral mesh. The proposed numerical method is focused on modeling the evolution of astrophysical objects of spherical shape, which is appropriate for gravitational collapse and star formation, and also for supernova explosion. The construction of tetrahedral grids is carried out in three stages. At the first stage, a geodesic grid methodology is used to construct a triangular grid on the surface of the sphere, which encompasses the computational domain. At the second stage, the resulting triangular mesh is serialized from the surface of the sphere to its center, and at the third stage, the obtained prisms are divided into tetrahedra. This approach allows us to simulate spherical objects without singularities that occur when using spherical or cylindrical coordinates. The paper describes numerical methods for solving the equations of hydrodynamics and the Poisson equation. Numerical examples are given that verify the developed numerical methods.


Introduction
Mathematical modeling is one of the most important theoretical tools for studying astrophysical objects. The main trends in hydrodynamic modeling are the Lagrangian method of smoothed particles [1,2] and Euler methods based on adaptive grids [3]. Over the past ten years, numerical methods based on moving grids [4,5], as well as on a combination of Lagrangian and Euler approaches, grid-based [6,7] or meshless [8], were actively developed. The main trend of all these approaches is the adaptation of the numerical discretization to the physical features of the problem.
We present an original approach to address such features based on tetrahedral meshes. Given our interest in star formation problems, the construction of numerical meshes is based on the use of geodesic grids [9]. The use of such grids allows us to construct efficient parallel numerical solution methods [10,11,12], including methods of high-order accuracy [13].
The second section describes the procedure for constructing tetrahedral meshes for modeling spherical objects. The third section briefly describes numerical methods for solving hydrodynamic equations on tetrahedral meshes and a method for solving the Poisson equation. The fourth section is devoted to numerical examples. In the fifth section, a conclusion is formulated. The scheme for constructing the tetrahedral mesh is shown in Fig. 1 the main nodes of the prism, the blue color indicates the nodes located at the faces of the prism, the pink color indicates the node that is placed in the center of the prism.

Numerical method
3.1. The method for solving the hydrodynamic equations We will consider the equations of hydrodynamics written in the form of conservation laws for mass, momentum, and total energy ∂ ∂t where ρ is the gas volume density, ⃗ u = (u x , u y , u z ) the gas velocity, p the gas pressure, E = p/(γ−1)+ρ⃗ u 2 /2 the total mechanical energy of gas, γ the adiabatic index. The conservation laws in the vector form are: where u is the vector of conservative variables, functions f (u), g(u), h(u) the fluxes in the corresponding directions. Additioanlly, we will consider the vector of physicsl variables q(u) = (ρ, u x , u y , u z , p). We assume that a simple procedure exists that can recover q(u) from the vector of conservative variables. Let us consider an arbitrary tetrahedral cell with the ordinal number 1. For describing a set of neighboring cells we will use index j. We then define the normal n ij of the cell i in the direction of cells j for the known vector ⃗ M i M ij , M i the center of the sphere inscribed in the tetrahedron i, M ij the center of the circle inscribed in the face between cells i and j. The area of faces and the volume of cells will also be necessary in the calculations. The volume of tetrahedron V i and triangles S ij can easily be written. Integrating the equations and considering the finite volume V i with the boundary Γ i = ∪ j S ij , we obtain: using the Gauss theorem, the integral equations can be written as The finite-volume scheme has the form: for equations: To solve the Riemann problem, we will use the scheme of Rusanov [14,15]: To determine the timestep τ , we define the diameters of spheres d i that are inscribed in every grid cell. Then the timestep can be calculated as: where CF L < 1 ia the the Courant number. To increase the order of accuracy, we employ the apprach analogous to that described in [16].

Method for solving the Poisson equation
The gravitational potential in equations of hydrodynamics is found by solving for the Poisson equation: The computed potential is taken into account in the right-hand-side of equations of hydrodynamics: To solve the Poisson equations, we consider the potential determined at the nodes of the computational grid, and employ the weak-set finite element method for discretization of the Laplace operator. Let us consider an arbitrary tetrahedon ABCD, the vertices of which We define four basis functions, each of which is unique in only one node, and in the others they are set equal to zero: This system can be analytically solved using the Cramer method. In this case, the gradient of each basis function in the tetrahedron has a simple analytical form: the sum of which allows us to build a potential gradient in a tetrahedral cell: where Φ i are the values of the potential in the corresponding nodes. When constructing a highorder accuracy scheme for solving hydrodynamic equations, we projected the density function into nodes, therefore, like the potential, density is also considered at the nodes of the tetrahedron.
In conclusion, we provide the formulas for the stiffness matrix S and the mass matrix M for the finite element formulation for the tetrahedron under consideration: At the outer boundary of the region, a fundamental solution of the Laplace equations in the form of boundary conditions of the first kind is used. To solve the resulting sparse system of equations, the conjugate gradient method is used.

The use of moving meshes
Let us define an analytic smooth function that defines the movement of the computational domain ⃗ w = (w x , w y , w z ). Smoothness is needed to eliminate the possible effect of tetrahydron degeneracy. The following function can be considered d ⃗ w dt ≈ ▽Φ. The equations of hydrodynamics can then be written as where P is the diagonal pressure tensor. After introducing the vector of effective gas velocity with respect to the moving mesh ⃗ v = ⃗ u − ⃗ w, the equations can be rewritten in the form of two systems: The first system can be solved according to the scheme provided earlier with the corresponding modification of the fluxes. The solution of the second system is not actually carried out, the grid nodes move with the conservation of mass, angular momentum, and the total energy in each computational cell. In the current implementation, we retained the possibility of including the movement of the computational grid, but in computational experiments this possibility was not used. In the future, we will study in detail the algorithms for the grid movement for each problem to be solved.

The entropy equation
As an additional equation for writing an overdetermined system of hydrodynamic equations and for regularizing the numerical solution, we can use the equation for entropy, also written in the conservative form. This approach was effectively used in [12,17,18,19] to achieve a non-negative value of pressure in regions with a high Mach number and to guarantee a non-decreasing entropy in the numerical method. The advantages and disadvantages of using the entropy equation are also discussed in detail in [13]. In the current implementation, we have retained the possibility of including the entropy equation in the numerical model. In computational experiments, this equation was not used.

Numerical examples
We consider four test problems as numerical examples: the problem of discontinuity decay to test the ability of the method to treat all hydrodynamic waves, the Sedov problem for testing the treatment of strong shock waves, the problem of solving Poisson's equation with differentiable functions to check the convergence degree of the method, and the collapse of a gas cloud to test the method in a complete mathematical model.

Decay of discontinuity
Let us consider the dynamics of gas on a segment [0; 1] till a time instance of t = 0.2. At x 0 = 0.5 the gas distribution has a discontinuity. The gas parameters on the left are: p L = 2 is the pressure, ρ L = 2 is the density, and v L = 0 is the velocity. On the right from the discontinuity the gas parameters are: p R = 1, ρ R = 1, and v R = 0. The adiabatic index is γ = 1.4. We set the non-reflecting boundary conditions on both sides. The results are presented in Fig. 2. 200 grid cells were used for this test.
The purpose of this test is to demonstrate the correct computation of contact discontinuity, rarefaction waves, and shock waves. Most methods for solving gas dynamics equations give either oscillations or diffusion on the position of the shock. The author's method of high order in accuracy reveals a small (by two cells) smearing of the solution at the shock wave region.

The Sedov test for a point explosion
Sedov's problem of a point explosion in astrophysics is designed to simulate a central supernova explosion. To model the point explosion problem, we consider a region [−0.5; 0.5] 3 , the adiabatic exponent γ = 5/3, the initial density in the region ρ 0 = 1, and the initial pressure p 0 = 10 −5 . At a time instance of t = 0, the internal energy E 0 = 0.6 is released. The explosion region is limited by the radius r central = 0.01. The simulated profile of density and angular momentum at t = 0.05 is shown in Fig. 3.
The Sedov point explosion problem is a standard test that checks the ability of a numerical method to reproduce strong shock waves with large Mach numbers. The sound speed of the background medium is negligible; therefore, the Mach number reaches a value M ≈ 1500. Our numerical method reproduces well the position of the shock wave, as well as the radial profile of the density and momentum.

Solution of the Poisson equation
The solution of the Poisson equation was studied by refining the computational grid and using the known analytical solution of the gravitational potential and the density function: 15r , r > 1.  Table 1. Euclidean norm of the deviation of the solution when refining the mesh. Table 1 provides the values of relative residuals (errors) during the mesh refinement. As can be seen from the Table, the second order of convergence is achieved. The scheme has the second-order of convergence.

Collapse of a gaseous cloud
Let us consider the gravitational collapse of an initially static gaseous cloud with pressure p = ρ/10, where the density profile has the form: 2r 3 − 3r 2 + 1, r ≤ 1, 0, r > 1.
The simulated gas density profile at time t = 0.5 is shown in Fig. 4. As can be seen from the Figure, the maximum density value increases by a factor of six, while most of the mass is concentrated in a sphere that is twice smaller than the initial distribution.

Conclusion
A new numerical method for solving the equations of gravitational hydrodynamics on tetrahedral grids is developed. A procedure for constructing tetrahedral grids for modeling spherical astrophysical objects based on geodetic grids is proposed. The results of computational experiments for verification of the developed numerical method are presented.