1 Introduction
Applied models of continuum mechanics Batchelor ; LandauLifshic1986 are based on conservation laws for the mass, momentum and energy. The transport of scalar and vector quantities due to advection determines a mathematical form of conservation laws Godunov ; LeVeque . In addition, some parameters of a flow have the positivity property (monotonicity). Such important properties of the differential problem of continuum mechanics must be inherited in a discrete problem Anderson ; Wesseling .
Flows of ideal fluids are governed by the Euler equations, whereas the Navier–Stokes equations are applied to describing viscous flows. Mathematical problems of validation of such models are considered, for example, in the books lions1 ; lions2 . When discussing the existence of solutions in various Sobolev spaces, the principal problems of the positivity (non–negativity) of the fluid density are also should be highlighted. Such a consideration is also carried out (see, for instance, feireisl2016mathematical ) at the discrete level for various approximations in time and space.
In computational fluid dynamics, the most important problems are associated with two contradictory requirements. Namely, it is necessary to construct monotone approximations for advective terms and to fulfil conservation laws. The construction of monotone approximations is discussed in many papers (see, e.g., kulikovskii2000mathematical ; HundsdorferVerwer2003 ; Kuzmin ). In MortonKellogg1996 ; SamarskiiVabischevich1999a , standard linear approximations are considered for the basic problems of continuum mechanics (convection–diffusion problems).
For discretization in space, conservative approximations are constructed on the basis of using the conservative (divergent) formulation of continuum mechanics equations. This approach is most naturally implemented using integro–interpolation method (balance method) for regular and irregular grids
Samarskii1989 , and in the control method volume LeVeque ; Versteeg . Nowadays, the main numerical technique to solve applied problems is the finite element method Guermond ; Larson . It is widely used in computational fluid dynamics Donea ; Zienkiewicz , too.Discretizations in time for computational fluid dynamics are often constructed using explicit schemes that have strong restrictions on a time step in sense of stability. Moreover, explicit schemes have similar restrictions on the monotonicity of an approximate solution. So, it is more natural to focus on implicit schemes. To solve boundary value problems for partial differential equations, two–level schemes are widely used
HundsdorferVerwer2003 ; Ascher2008 ; LeVeque2007 (–method, schemes with weights). For linear problems, a study of discretizations in time can be based on the general theory of stability (well–posedness) for operator–difference schemes Samarskii1989 ; SamarskiiMatusVabischevich2002 . In particular, it is possible to apply unimprovable (coinciding necessary and sufficient) stability conditions, which are formulated as operator inequalities in finite–dimensional Hilbert spaces.In the present work, an initial–boundary value problem is considered for the Euler equations describing barotropic fluid flows (Section 2), which are conservation laws for the mass, momentum, and total mechanical energy. Discretization in space is performed (Section 3) using standard Lagrange elements for the density and cartesian velocity components. To evaluate an approximate solution at a new time level, the fully implicit scheme is employed. For the approximate solution, the mass conservation law holds and an estimate for the dissipation of the total mechanical energy is fulfilled. The fully implicit scheme is not convenient for numerical implementation. The solution at the new time level is determined from a system of coupled nonlinear equations for the density and velocity. A decoupling scheme is proposed in Section 4, which refers to the class of linearized schemes of splitting by physical processes
Marchuk1990 ; Vabishchevich2014 . Linearization is carried out over the field of advective transport in such a way that at each time level we solve individual problems for the density and velocity. Possibilities of the proposed schemes are illustrated by the results of numerical solving a model two–dimensional problem with a perturbation of the fluid density being initially at rest (Section 5). To solve numerically the nonlinear discrete problem at the new time level, the Newton method is used. In the above calculations, a small number of iterations (two or three) is sufficient for the process convergence. The influence of the grid size in space and time is investigated. It was observed that decreasing of the time step results in the monotonization of the numerical solution. The main result of the paper is the proof of the robustness of the linearized decoupling scheme. The scheme involves separate solving standard advection problems for the density and velocity and demonstrates high iteration convergence. Such an approach can be used for other problems of continuum mechanics, e.g., for numerical solving initial–boundary value problems for the Navier–Stokes equations.2 Mathematical models
An initial–boundary value problem is considered for describing barotropic fluid flows. The system of equations includes the scalar advection equation for the density and the vector advection equation for the velocity with a given pressure dependence on the density. The conservation laws for the mass, momentum, and total mechanical energy are discussed.
2.1 Barotropic fluid
The continuity equation in a bounded domain has the form
(1) 
where is the density and is the velocity. The momentum equation is written in the conservative form
(2) 
where is the pressure. The considered fluid is assumed to be barotropic, i.e., we have a known dependence of the pressure on the density , .
Assume that the domain boundaries are rigid and so, the impermeability condition is imposed:
(3) 
Initial conditions for the density and velocity are also specified:
(4) 
The initial–boundary value problem (1)–(4) describes transient flows of an ideal barotropic fluid.
The direct integration of the continuity equation (1) over the domain taking into account the boundary condition (3) results in the mass conservation law
(5) 
In the Hilbert space , we define the scalar product and norm in the standard way:
In a similar way, the space of vector functions is defined. If the density is non–negative, the conservation law for the mass can be written as
This relation can be treated as an a priori estimate for in .
The equation (2) directly expresses the conservation law for the momentum. Integrating this equation over , we obtain
Thus, we have
(6) 
Multiplying by and taking into account equation (1), rewrite equation (2) as
Integration over the domain , in view of (3), leads to
(7) 
The second term in (7) is expressed from the renormalized equation of continuity. Define the pressure potential from the equation
(8) 
In particular, for an ideal fluid, we have
(9) 
2.2 Operator–differential formulation
For the convenience of consideration, we introduce operators of advective (convective) transport for the system of Euler equations. The advection operator in the divergent form is written as follows:
(13) 
Assuming that the boundary condition (3) is satisfied for the velocity , we obtain
(14) 
The continuity equation (1) can be written in the form of an operator–differential equation:
(15) 
where the notation is used. Similarly, equation (2) is written in the form
(16) 
For the system of equations (15), (16) with a prescribed dependence , we consider the Cauchy problem, where the initial conditions (see (4)) have the form
(17) 
3 Implicit two–level scheme
To solve numerically the initial–boundary value problem for the Euler equations, we use the fully implicit (backward Euler) scheme for time–stepping with finite element discretizations in space. The problems of fulfilment of the conservation laws at the discrete level are discussed.
3.1 Discretization in space
To solve numerically the problem (15)–(17) (or (15)–(19)), we employ finite element discretizations in space (see, e.g., Thomee2006 ; brenner2008mathematical ). For (13), we define the bilinear form
Define the subspace of finite elements and the discrete operator as
Similarly to (14), we have
(18) 
For two– and three–dimensional vector quantities, the coordinate–wise representation is employed: . For a simple specification of the boundary conditions (3), assume that separate parts of the boundary of the computational domain are parallel to the coordinate axes. A finite element approximation is used for the individual components of the vector .
After constructing discretizations in space, we arrive at the Cauchy problem for the system of semi–discrete operator equations in the corresponding finite–dimensional space, namely, we have the Cauchy problem for the system of ordinary differential equations. For instance, for (
15)–(17), we put into the correspondence the problem(19) 
(20) 
(21) 
Here , with denoting –projection onto .
3.2 Discretization in time
Let be a step of a uniform, for simplicity, grid in time such that , . To construct and study time–stepping schemes, the main attention is given to the fulfillment of the corresponding conservation laws (a priori estimates). Such an important problem as the positivity (non–negativity) of the density at each time level requires a more in–depth study and so, it is not considered in the present work.
To solve numerically the problem (19)–(21), the fully implicit scheme is applied. In this case, the approximate solution at the new time level is determined from
(22) 
(23) 
using the prescribed (see (21)) value . The basic properties of the approximate solution are related to the fulfillment of the conservation laws for the mass and total energy. To simplify our investigation, assume that the density is positive.
Assumption 1
At each time level , .
In view of (18), integration of equation (22) over the domain leads to
(24) 
The equality (24) is a discrete analog the mass conservation law (5). For the momentum conservation law (6), we put into the correspondence the equality
(25) 
which has been obtained by integrating equation (23).
An estimate for the total mechanical energy can be established, e.g., following the work feireisl2016mathematical . Multiplying equation (23) by and integrating it over , we arrive at
(26) 
For the first term, we have
In view of this, from (26), we obtain
(27) 
From (22) and the definition of the operator , we have
This makes possible to rewrite the inequality (27) as
(28) 
To estimate the second term on the left–hand side of the inequality (28), we apply the discrete analogue of the renormalized equation of continuity (10). Multiply the continuity equation (22) by :
(29) 
The following equality takes place:
where
Assumption 2
Assume that
(30) 
Under these natural assumptions, we get
For the second term in (29), taking into account (8), we have
In view of this, integration of equation (29) results in
(31) 
Combining (28) and (31), we obtain the inequality
(32) 
Comparing (32) with (12), we can conclude that at the discrete level, instead of fulfillment of the conservation law for the total energy, we observe a decrease of the energy. It should be noted that this property has been established under the additional assumption (30). The result of our consideration can be expressed in the following statement.
Proposition 1
The fully implicit scheme (22), (23) produces the approximate solution of the problem (19)–(21) that satisfies the mass conservation law in the form ( ref 24) and the momentum conservation law (25). Moreover, if the assumptions 1 and 2 hold, the estimate (32) for the total mechanical energy is also fulfilled.
4 Decoupling schemes
A linearized scheme is used, where the solution at a new time level is evaluated by advective transport taken from the previous time level. Using such a linearization, an iterative process is constructed for the numerical implementation of the fully implicit scheme. The approximate solution at the new time level is determined by sequential solving, first, the linear problem of advection for the density and, secondly, the linear problem for the velocity.
4.1 Linearized scheme
We focus on the use of such time–stepping techniques that demonstrate the following properties:

the transition to a new time level is implemented by solving linear problems;

splitting with respect to physical processes is employed, namely, the problems for the density and velocity are solved separately (with individual problems for the velocity components).
An example of the simplest decoupling scheme for the Euler equations system (19), (20) is the linearized scheme, where the advective transport involves the velocity from the previous time level.
Instead of (22), (23), we employ the scheme
(33) 
(34) 
First, from the linear transport equation (33), we evaluate the density at the new time level. Next, from the linear decoupled system (34) for the velocity components, we calculate the velocity .
Remark 0
The system of equations (34) with a given density is, in the general case, coupled for individual cartesian velocity components. In the case, where all parts of the boundary of the computational domain are parallel to the axes of the cartesian coordinate system, the system of equations (34) is decoupled and so, we can evaluate independently the individual components of the velocity.
4.2 Iterative decoupling scheme
On the basis of the linearized scheme (33), (34), it is possible to construct an iterative algorithm for the numerical implementation of the fully implicit scheme (22), (23). The approximate solution for , at the –th iteration is denoted by , , with the initial approximation from the previous time level:
(35) 
Assume that the new approximation at the new time level is calculated, when the previous iterations have been done. Similarly to (33), (34), we use the system of equations:
(36) 
(37) 
where
(38) 
Thus, at each iteration, we firstly solve the linear problem for the density and only then calculate the linear problem for the velocity.
5 Numerical results
The possibilities of the fully implicit scheme and decoupling schemes are illustrated by numerical results for a model two–dimensional problem with a density perturbation of an initially resting fluid.
5.1 Test problem
Here, we present the results of numerical solving a model problem obtained using different time–stepping techniques. The problem (1)–(4) is considered in the square
Assume that the dependence of the density on the pressure has the form ((9) with , . We simulate the motion of the initially resting fluid ( in (4)) with the initial density (see (4)) specified in the form
where and .
5.2 Fully implicit scheme
The problem is solved using the standard uniform triangulation on segments in each direction. The piecewise–linear finite elements are employed for discretization in space. To implement the fully implicit scheme (22), (23) for the nonlinear discrete problem at the new time level, the Newton method with a direct solver is applied for the corresponding system of linear algebraic equations.
Time–evolution of the compressible fluid flow is shown in Fig. 1
, which presents the density at various time moments. In this calculation, we used the spatial grid with
, the time step was . Time–histories of the density in the center of the computational domain (), maximum () and minimal () values of the density over the entire domain are given in Fig. 2.Newton’s iterative method for solving the discrete problem at each new time moment converges very quickly (two or three iterations are enough). Table 1 demonstrates convergence of the iterative process for the first step in time. Here, we present the relative error for the first three iterations for the model problem obtained with and various time steps.
iteration  

1  2.219e03  5.602e04  1.404e04 
2  5.104e08  9.100e10  1.458e11 
3  5.566e15  9.994e15  1.788e14 
The accuracy of the approximate solution of the test problem will be illustrated by the data on the density in the section . The solution calculated on the grid with using various grids in time is shown in Fig. 3. Similar data for and are shown in Fig. 4, 5, respectively. It is easy to see a good accuracy in the reconstruction of the leading edge of the wave, when the spatial grid is refined. Also, we observe the effect of smoothing, namely, elimination of non–monotonicity with increasing of time step.
5.3 Decoupling scheme
In using the decoupling scheme (35)–(38), the greatest interest is related to the convergence rate of the iterative process. The time step in the calculations was equal to . We present numerical results for the model problem under consideration obtained on different grids in space. The dependence of the solution on the number of iterations for the grid with is shown in Fig. 6. Figure 7 and 8 presents similar data for grids with and , respectively. It is easy to see see that on the finest grid (see Fig. 8) the linearized scheme (34), (34)) ( in (35)–(38)) yields a substantially non–monotonic solution, which is monotonized on subsequent iterations.
The main conclusion of our study is the demonstration of high computational efficiency of the iterative decoupling scheme (35)–(38). Namely, for the problems under consideration it is sufficient to do only two iterations using (35)–(38).
For the above time–stepping methods, the key point is a violation of the conservation law for the total energy. For the fully implicit scheme (22), (23), instead of conservation of the energy (see the estimate (32)), decreasing of the total energy is observed.
The dynamics of the total mechanical energy using a linearized scheme (33), (34) () and iterative decoupling schemes (35)–(38)) for on various grids is shown in Fig. 9–11. Here, according to (12), is calculated at each time moment
For , the solution obtained using the decoupling scheme (35)–(38) practically coincides with the solution derived from the fully implicit scheme (22), (23). The above data indicate that the conservation law for the total energy is satisfied with a good accuracy. Decreasing of the time step results in increasing of the accuracy of the conservation law fulfillment.
Acknowledgements
The publication was financially supported by the Ministry of Education and Science of the Russian Federation (the Agreement # 02.a03.21.0008).
References
 (1) G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, 2000.
 (2) L. D. Landau, E. Lifshitz, Fluid Mechanics, ButterworthHeinemann, 1987.
 (3) S. K. Godunov, E. I. Romenskii, Elements of Continuum Mechanics and Conservation Laws, Springer, 2003.
 (4) R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.
 (5) J. Anderson, Computational Fluid Dynamics: The Basics with Applications, McGrawHill, 1995.
 (6) P. Wesseling, Principles of Computational Fluid Dynamics, Springer, 2010.
 (7) P.L. Lions, Mathematical Topics in Fluid Mechanics: Incompressible models, Oxford University Press, 1996.
 (8) P.L. Lions, Mathematical Topics in Fluid Mechanics: Compressible Models, Oxford University Press, 1998.
 (9) E. Feireisl, T. G. Karper, M. Pokornỳ, Mathematical Theory of Compressible Viscous Fluids: Analysis and Numerics, Springer, 2016.
 (10) A. G. Kulikovskii, N. V. Pogorelov, A. Y. Semenov, Mathematical Aspects of Numerical Solution of Hyperbolic Systems, Taylor & Francis, 2000.
 (11) W. H. Hundsdorfer, J. G. Verwer, Numerical Solution of TimeDependent AdvectionDiffusionReaction Equations, Springer Verlag, 2003.
 (12) D. Kuzmin, A Guide to Numerical Methods for Transport Equations, University Erlangen–Nuremberg, 2010.
 (13) K. W. Morton, R. B. Kellogg, Numerical Solution of ConvectionDiffusion Problems, Chapman & Hall London, 1996.
 (14) A. A. Samarskii, P. N. Vabishchevich, Numerical Methods for Solving ConvectionDiffusion Problems, URSS, Moscow, 1999, in Russian.
 (15) A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001.
 (16) H. Versteeg, W. Malalasekra, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Prentice Hall, 2007.
 (17) A. Ern, J.L. Guermond, Theory and Practice of Finite Elements, Springer, 2004.
 (18) M. G. Larson, F. Bengzon, The Finite Element Method: Theory, Implementation and Applications, Springer, 2013.
 (19) J. Donea, A. Huerta, Finite Element Methods for Flow Problems, Wiley, 2003.
 (20) O. C. Zienkiewicz, R. L. Taylor, N. P., The Finite Element Method for Fluid Dynamics, ButterworthHeinemann, 2013.
 (21) U. M. Ascher, Numerical Methods for Evolutionary Differential Equations, Society for Industrial and Applied Mathematics, 2008.
 (22) R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: SteadyState and TimeDependent Problems, Society for Industrial and Applied Mathematics, 2007.
 (23) A. A. Samarskii, P. P. Matus, P. N. Vabishchevich, Difference Schemes with Operator Factors, Kluwer Academic, Dordrecht, 2002.
 (24) G. I. Marchuk, Splitting and alternating direction methods, in: P. G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, Vol. I, NorthHolland, 1990, pp. 197–462.
 (25) P. N. Vabishchevich, Additive OperatorDifference Schemes: Splitting Schemes, de Gruyter, Berlin, 2014.
 (26) V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Verlag, Berlin, 2006.
 (27) S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York, 2008.