The Euler equations model the flow of compressible fluids and are a system of non-linear, coupled partial differential equations which mathematically embody the principles of conservation of mass, momentum and energy. In the presence of an external gravitational field, the fluid elements experience an additional force due to gravity which contributes to a source term in the momentum and energy conservation equations. The presence of these source terms leads to non-trivial steady state solutions of the equations. An important class of such solutions are the so calledhydrostatic solutions which are steady solutions with zero velocity. In many problems like weather prediction in the Earth’s atmosphere or in a star, the system may be in a hydrostatic state atleast locally and for very long periods of time. Accurately simulating such stationary solutions with a numerical scheme can be a challenging task since commonly used algorithms may not have the property to maintain hydrostatic solutions. A numerical scheme that can exactly maintain a hydrostatic initial condition for all times and on any grid is said to be well-balanced. Schemes which are not exactly well-balanced have to rely on their consistency to approximately maintain hydrostatic solutions. It becomes necessary to use very high order schemes and/or large grids in order to keep the error to small levels. This becomes especially critical if we are interested in computing small perturbations around the hydrostatic solution, as the truncation errors can dominate the small perturbations leading to very large errors in the solutions.
In the hydrostatic state, the pressure gradient is balanced by the gravitational force. This balance must be achieved in the numerical scheme also for it to be well-balanced. The hydrostatic state is governed by a first order differential equation whose solution requires assuming additional properties on the solution, e.g., isothermal, polytropic, etc. Since there is a wide variety of hydrostatic solutions depending on the equation of state and other conditions of a particular physical problem, it is not possible to construct a scheme that is well-balanced for a wide variety of problems. The well-balanced property is usually achieved by exploiting the structure of the hydrostatic solution which is used in the reconstruction process of finite volume schemes, and also in discretizing the gravitational source terms. Most of the existing schemes have been constructed to be well-balanced for ideal gas equation of state together with either isothermal or polytropic solutions. The scheme in Kappeli2014 is based on isentropic assumption and uses constancy of enthalpy to perform the reconstruction, while Kappeli2017 constructs a scheme for isothermal case by using the constancy of Gibbs free energy. The scheme by Chandrashekar and Klingenberg Chandrashekar2015 was well-balanced for both isothermal and polytropic solutions and did not require á priori knowledge of the temperature (for isothermal case) or polytropic index. Schemes which are well-balanced for arbitrary equations of state have been proposed in Ghosh2015 , Ghosh2016 , Berberich2016 but they require á priori knowledge of the hydrostatic solution which is used to rewrite the source term to achieve the well-balanced property. A few well-balanced discontinuous Galerkin methods have also been proposed in recent years. A well-balanced DG scheme for isothermal case was proposed in Li2016 and extensions were made in Li2018 for polytropic case. Chandrashekar and Zenk Chandrashekar2017 proposed well-balanced nodal discontinuous Galerkin methods for isothermal and polytropic solutions but these schemes require á priori knowledge of the hydrostatic solution which is used in the source term approximation. A finite volume scheme that preserves a discrete hydrostatic solution has been constructed in Kappeli2016 which does not require á priori knowledge of the type of hydrostatic solution that may occur in a particular problem. Note that hydrostatic solutions are not polynomial functions and hence they cannot be exactly represented in a finite element method. What is preserved by the well-balanced finite element method is a projection of the hydrostatic solution onto the finite element space, e.g., Chandrashekar2017
preserves the interpolation of the hydrostatic solution onto the Lagrange polynomials. Many schemes require the use of particular numerical flux functions to achieve well-balanced property, e.g.,Xing2013 , Ghosh2016 , Li2017 make use of a modified local Lax-Friedrich flux where the dissipation operator is altered and Desveaux2014 , Desveaux2016 are built using relaxation solvers, while other schemes like in Chandrashekar2015 , Chandrashekar2017 , Kappeli2014 , Kappeli2017 are well-balanced for any numerical flux function. The latter type of schemes have greater applicability since practitioners can make use of the best available numerical flux functions for their particular problem. Among other approaches, we may mention the wave propagation scheme in LeVeque1999 and central schemes in Chertock2018 .
In finite volume and finite difference schemes, the basic unknowns are values associated with each volume or point. A well-balanced finite volume scheme may be designed to preserve the interpolation of the exact hydrostatic solution on the mesh. In this work, we do not try to exactly preserve the exact hydrostatic solution, which may not even be known in an explicit form for general equations of state and for realistic problems occuring in astrophysics. Instead, we construct a scheme which admits a hydrostatic solution which is a second order accurate approximation to the exact hydrostatic solution. This allows us to make the scheme well-balanced for general equations of state and for any type of hydrostatic solution, and we do not need á priori knowledge of which type of hydrostatic solution may be realized in a particular problem. Note that our scheme is well-balanced for the discrete hydrostatic solution which can be found numerically. This is in contrast to the scheme in Chandrashekar2015 which well-balances the exact hydrostatic solution but only for ideal gas and certain types of hydrostatic solutions, namely isothermal and polytropic. In the Appendix, we give some arguments to show that such a discretely well-balanced scheme can give consistent and accurate approximation to small perturbations around a hydrostatic solution. Any consistent numerical flux function that exactly resolves stationary contact discontinuities can be used which makes this method of interest to practitioners who want to use their own favourite numerical flux function. We construct this scheme on Cartesian grids but the extension to general curvilinear grids can also be made. The well-balanced property is achieved by using two main techniques: (1) a hydrostatic reconstruction scheme and (2) a discretization of gravitational source terms in terms of a local hydrostatic solution. These ingredients are similar to what is used in Chandrashekar2015 but with some modifications that allow us to deal with general hydrostatic solutions, whereas the previous work was restricted to ideal gas equation of state. Both of these techniques are based on exploiting a particular structure of the hydrostatic solution. A restriction of this scheme and also the one in Kappeli2016 , is that the well-balanced property requires the gravitational force to be aligned with the grid lines. In atmospheric or stellar problems where the gravity is predominantly in the radial direction, this requirement can be satisfied with the use of polar or spherical grids. However we show in numerical solutions that in practice, the hydrostatic state is quite well maintained even otherwise.
The rest of the paper is organized as follows. In section (2), we introduce the 1-D Euler equations and some standard equations of state that are of common interest. In section (3), we give a full description of our approach to well-balanced scheme and prove this as a theorem. The extension to two dimensions is made in section (4) where we also discuss the treatment of some boundary conditions. Sections (5) and (6) present several numerical results in one and two dimensions to show the performance of the current scheme. We end the paper with some conclusions in section (7).
2 1-D Euler equations with gravity
Consider the system of compressible Euler equations in one dimension which models conservation of mass, momentum and energy and are given by
Here is the density, is the velocity, is the pressure, is the energy per unit volume excluding the gravitational energy and is the gravitational potential. The energy is given by
where is the internal energy per unit mass. We can write the above set of coupled equations in a compact notation as
where is the set of conserved variables and
is the corresponding flux vector. In the case of a self gravitating system, the gravitational potentialis governed by a Poisson-type equation whose details are not relevant for the present discussion. We will consider the case of static gravitational potential which is assumed to be given as a function of the spatial coordinates. However the scheme we propose can be used for time dependent potential also. To simplify some of the later notation, we also denote the set of primitive variables by .
2.1 Hydrostatic states
Consider the hydrostatic stationary solution, i.e., the solution for which the velocity is zero
We will denote the quantities in a hydrostatic solution with an overbar. In this case, the mass and energy conservation equations are automatically satisfied. The momentum equation becomes an ordinary differential equation given by
We have only one equation but two thermodynamic quantities, , , to be determined. The pressure, density and temperature are related by an equation of state. We will write the equation of state in the form
where is in general a function of the thermodynamic variables, e.g., or , with being the temperature. This form of equation of state is commonly used in astrophysical applications where the function may be given in a tabular form. Integrating the hydrostatic equation (1) we obtain
In the above equation, is the pressure at some reference position . This is a general relation that holds for any hydrostatic solution and any equation of state, and is the crucial building block for our scheme. The above integral can be explicitly evaluated in some simple situations. We next give several examples to illustrate that the equations of state can be put in the form (2) and show the hydrostatic relation in some simple situations.
2.1.1 Ideal gas
For an ideal gas model we have where is the gas constant. Assuming some equilibrium temperature profile , the pressure can be computed using (3) as
If the hydrostatic state is isothermal, i.e., , then
If the hydrostatic solution is polytropic Chandrasekhar1967 , then we have the following relations satisfied
where is some constant. Using these polytropic relations in the hydrostatic equation (1) and performing an integration, we obtain
In this case the temperature profile is determined once the potential is fixed. If then we have an isentropic solution.
2.1.2 van der Waals equation of state
Unlike the ideal gas which is assumed to be made of point particles, the van der Waals equation of state accounts for the volume occupied by the molecules composing the gas and the inter-molecular forces, and is useful in modeling problems with phase transitions. The equation of state is given byVANDERWAALS1873
where is the universal gas constant, and are characteristic values for an individual gas and is its molar mass. We can write this in the form (2) where
If we assume a condition of polytropic van der Waals gas, i.e., the specific heat at constant volume does not change with temperature, the total energy per unit volume for the system can be written as guardone2002
The acoustic speed for a polytropic van der Waals condition is derived to be
which is used in the HLLC numerical flux function. Since we may not have explicitly known hydrostatic solutions, it has to be calculated using a numerical quadrature of the hydrostatic equation.
2.1.3 Ideal gas with radiation pressure
The pressure exerted by radiation becomes significant inside stellar atmospheres due to the presense of high temperatures and hence it must be included in the equation of state. The equation of state is given by Chandrasekhar1967
where is the Stefan-Boltzmann constant. This equation can be written in the form (2) where
and the corresponding internal energy is given by
In this case, we may not have explictly known hydrostatic solutions available to us. The hydrostatic solution has to be computed through some numerical quadrature of the hydrostatic equation.
3 1-D finite volume scheme
Let us divide the domain into finite volumes each of size . The ’th cell is given by the interval and let be the cell center. Consider the semi-discrete finite volume scheme for the ’th cell
where is a consistent approximation of the gravitational source term and is the numerical flux function which we write as a function of primitive variables. The scheme will be completely specified once the source term and numerical flux schemes are explained. For later use, let us define a piecewise linear interpolation of the gravitational potential
and let be a piecewise constant function defined as
Motivated by the structure of the hydrostatic solution given by equation (3), define
where the lower limit of the integral is arbitrary and not relevant for the definition of the scheme since only differences of the function are actually used in the scheme. Note that is a continuous function. This scheme differs from the scheme in Chandrashekar2015 which was specific to ideal gas so that and we used a logarithmic average for the definition of in a piecewise manner over the intervals . The current scheme is much simpler since it does not involve logarithmic averages and the value of is defined to be piecewise constant over each cell.
3.1 Numerical flux
The numerical flux is computed using two states at the interface which are obtained by a reconstruction process. We will assume that the function is consistent in the sense that . Moreover, we assume the following property is also satisfied.
Contact Property The numerical flux is said to satisfy contact property if for any two states and we have
The states , in the above definition correspond to a stationary contact discontinuity. The above property is equivalent to the ability of a numerical flux to exactly support a stationary contact discontinuity. A few important examples of numerical fluxes which satisfy this property are the Roe flux Roe1981 and the HLLC flux Toro1994 .
3.2 Reconstruction scheme
To obtain the states at the cell boundary which are required to calculate the numerical flux , we will reconstruct the following set of variables
which is motivated by the structure of the hydrostatic solution as shown in equation (3). To perform the reconstruction at , define in equation (8) relative to , and compute the quantities , and from the definition of as
The values of the potential at the faces are obtained by the piecewise linear interpolation so that , etc. The variables at the cell centers are defined as
Using these values, we can use any reconstruction scheme to obtain and . E.g., the minmod scheme is given by
where and is the minmod limiter function given by
The primitive variables are obtained from the variables leading to
where . Using these values, we can now compute the numerical flux. Note that by definition and hence the exponential factors can be dropped when we convert to the primitive variables.
We have explained the scheme assuming that the equation of state has the form . In some applications, the function may be given in the form and the above approach can be used in this case also. If is given as some function of , then we can define .
3.3 Approximation of source term
The source term could be calculated in a straight-forward way by using a finite difference approximation of the potential gradient. However for achieving well balanced property, it is necessary to approximate it in a different but equivalent manner. Inside the
’th cell, using the pressure at the cell center, we will estimate the pressure at the faces assuming that the solution is in a hydrostatic state, i.e., using equation (3). Thus the local hydrostatic pressures at the faces inside the ’th cell are given by
Then the source term in the momentum equation is approximated as
In the next theorem we show that this is a consistent approximation. The source term discretization makes use of solution values only within the cell in contrast to the scheme used in Chandrashekar2015 . This scheme is hence expected to be better for discontinuous hydrostatic solutions, e.g., where the density is discontinuous.
The source term discretization is given by
Performing a Taylor expansion of the potential around we get
Hence the source term discretization is second order accurate. ∎
3.4 Well-balanced property
We now state the basic result on the well-balanced property. We essentially show that the finite volume schemes admit a stationary solution which is a discrete analogue of (3).
The finite volume scheme (7) together with a numerical flux which satisfies contact property and reconstruction of variables is well-balanced in the sense that the initial condition given by
is preserved by the numerical scheme.
Let us start the computations with an initial condition that satisfies (11). Since we reconstruct the variables and by our assumption the second and third components are constant, at any interface we have
and hence, since is a continuous function, we get
and similarly at
The above equations together with (9) imply that and . In general, the density may not be continuous, i.e., . But since the numerical flux satisfies contact property, we have and . The flux in mass and energy equations are zero and the gravitational source term in the energy equation is also zero. Hence the mass and energy equations are already well balanced, i.e., and . It remains to check the momentum equation. On the left of the momentum equation, we have
while on the right, the source term takes the form
and hence . This proves that the initial condition satisfying (11) is preserved under any time integration scheme. ∎
Note that the scheme preserves a state that satisfies the conditions given in (11). The exact hydrostatic solution may not satisfy these conditions since the quantity is only known approximately in general. See the next theorem for a special case where we obtain exact well-balanced property. In general, the scheme well-balances an approximate hydrostatic solution which is second order accurate approximation to the true hydrostatic solution. We provide some theoretical justification in the Appendix for why such a scheme can still give good approximations.
It is possible to reconstruct density and still retain the result of the previous theorem. In the isothermal case, the quantity is constant and we can expect the reconstruction of density to be more accurate if we scale the density as in the variables.
Definition 1 (Exact hydrostatic solution).
Let , be a hydrostatic solution which is available in closed form expression. Then the interpolation of this solution on the grid
will be refered to as the exact hydrostatic solution.
The above theorem says that the finite volume scheme admits hydrostatic solutions. This does not mean that it can preserve every exact hydrostatic solution. The next theorem shows that the finite volume scheme is well-balanced for the exact hydrostatic solution in a special situation.
Any hydrostatic solution for ideal gas EOS which is isothermal is exactly preserved by the finite volume scheme (7).
3.5 Computation of discrete hydrostatic solution
The new scheme is well-balanced at a state that satisfies conditions (11) which is essentially a numerical quadrature of the hydrostatic equation (1). Let us discuss how to compute such a solution. We will assume that the equilibrium temperature and the corresponding potential are known at the grid points. We will use equation (11) to determine the pressure and density at all the grid points. The hydrostatic equation is a first order ordinary differential equation. Hence it requires an initial condition; assume that we are given the pressure at the first grid point. If the pressure at is known then the pressure at is given by
By the definition of , we get
The pressure is the root of the function . The root finding problem is solved using a Newton method for which a good initial guess for the root is given by . Once the pressure has been computed the density is obtained from the equation of state where . For equations of state with = such as the van der Waals equation, we define a root-finding problem for which the root is the density of the function
The pressure can be computed using = .
Definition 2 (Discrete hydrostatic solution).
The set of grid point values , obtained from the procedure described in the previous section will be refered to as discrete hydrostatic solution.
It is not necessary to compute the discrete hydrostatic solution in practical computations. We will use this only to demonstrate numerically that the scheme is well-balanced if we start the computations using the discrete hydrostatic solution as the initial condition.
For the ideal gas model and an isothermal hydrostatic solution, theorem (3) shows that and .
4 Two dimensional scheme
The two dimensional Euler equations with gravity are a system of four balance laws which can be written as
where is the set of conserved variables and , are the corresponding flux vector in Cartesian coordinates. Here is the density, is the velocity, is the pressure, is the energy per unit volume excluding the gravitational energy and is the gravitational potential. The energy is given by
where is the internal energy per unit mass. The source term is where and . We also define the set of primitive variables by . The hydrostatic solution satisfies the equation
Similar to the one dimensional case, define
which is independent of the path if we are at hydrostatic solution. Moreover the hydrostatic solution satisfies
In the finite volume scheme, we will choose a path along or axis depending on the flux that is being evaluated. For example, for the computation of the component of the flux, the reconstruction is made using a 1-D scheme along the direction and the path is parallel to the -axis. The discretization of the balance law on Cartesian mesh is obtained in a straight-forward way by using the one-dimensional scheme along the two coordinate directions, and hence we give only a brief description of the 2-D scheme.
4.1 Finite volume scheme
Consider a partition of the computational domain by a Cartesian mesh consisting of and number of points along the coordinate directions, see figure (1a). The coordinates of the points are given by
where and are the cell sizes. Note that we have grid points located on the boundary of the computational domain which is commonly refered to as a cell-vertex scheme. Around each grid point , we construct a cell where and . For the points on a solid wall boundary, the cells are only half the size. For e.g., on the bottom boundary for which , the cells are defined by (see figure (1b)) whereas on the left boundary for which , the cells are defined as . If we are using periodic or transmissive boundary condition in some direction, then the cells on the boundary of the corresponding direction are of full size.
The semi-discrete finite volume scheme is given by
for all interior cells. The source term is computed as where
where the hydrostatic pressures , etc., are computed by applying the one-dimensional scheme along each of the coordinate directions. If the bottom boundary is a solid wall, then the cells are only half the size as shown in figure (1b) and the semi-discrete finite volume scheme for cells on this boundary is of the form
The source term is where