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 gridsSamarskii1989 , 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 usedHundsdorferVerwer2003 ; 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 processesMarchuk1990 ; 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
where is the density and is the velocity. The momentum equation is written in the conservative form
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:
Initial conditions for the density and velocity are also specified:
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
Integration over the domain , in view of (3), leads to
The second term in (7) is expressed from the renormalized equation of continuity. Define the pressure potential from the equation
In particular, for an ideal fluid, we have
From the continuity equation (1), we have
Integration of the renormalized equation of continuity (10) results in the expression
Adding this equality to (7), we get
We arrive to the conservation law for the total mechanical energy:
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:
Assuming that the boundary condition (3) is satisfied for the velocity , we obtain
The continuity equation (1) can be written in the form of an operator–differential equation:
where the notation is used. Similarly, equation (2) is written in the form
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
Define the subspace of finite elements and the discrete operator as
Similarly to (14), we have
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
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.
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.
At each time level , .
which has been obtained by integrating equation (23).
For the first term, we have
In view of this, from (26), we obtain
From (22) and the definition of the operator , we have
This makes possible to rewrite the inequality (27) as
The following equality takes place:
Under these natural assumptions, we get
In view of this, integration of equation (29) results in
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.
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).
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:
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
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.
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.
The publication was financially supported by the Ministry of Education and Science of the Russian Federation (the Agreement # 02.a03.21.0008).
- (1) G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, 2000.
- (2) L. D. Landau, E. Lifshitz, Fluid Mechanics, Butterworth-Heinemann, 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, McGraw-Hill, 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 Time-Dependent Advection-Diffusion-Reaction 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 Convection-Diffusion Problems, Chapman & Hall London, 1996.
- (14) A. A. Samarskii, P. N. Vabishchevich, Numerical Methods for Solving Convection-Diffusion 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, Butterworth-Heinemann, 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: Steady-State and Time-Dependent 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, North-Holland, 1990, pp. 197–462.
- (25) P. N. Vabishchevich, Additive Operator-Difference 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.