 # Decoupling schemes for predicting compressible fluid flows

Numerical simulation of compressible fluid flows is performed using the Euler equations. They include the scalar advection equation for the density, the vector advection equation for the velocity and a given pressure dependence on the density. An approximate solution of an initial--boundary value problem is calculated using the finite element approximation in space. The fully implicit two-level scheme is used for discretization in time. Numerical implementation is based on Newton's method. The main attention is paid to fulfilling conservation laws for the mass and total mechanical energy for the discrete formulation. Two-level schemes of splitting by physical processes are employed for numerical solving problems of barotropic fluid flows. For a transition from one time level to the next one, an iterative process is used, where at each iteration the linearized scheme is implemented via solving individual problems for the density and velocity. Possibilities of the proposed schemes are illustrated by numerical results for a two--dimensional model problem with density perturbations.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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

 ∂ϱ∂t+div(ϱu)=0,x∈Ω,0

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:

 (u⋅n)=0,x∈∂Ω. (3)

Initial conditions for the density and velocity are also specified:

 ϱ(x,0)=ϱ0(x),u(x,0)=u0(x),x∈Ω. (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

 m(t)=m(0),m(t)=∫Ωϱ(x,t)dx. (5)

In the Hilbert space , we define the scalar product and norm in the standard way:

 (w,u)=∫Ωw(x)u(x)dx,∥w∥=(w,w)1/2.

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

 ∥ϱ1/2∥=∥ϱ1/20∥.

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

 ∫Ω∂∂t(ϱu)dx+∫∂Ωp(ϱ)ndx=0.

Thus, we have

 I(t)=I(0)−∫∂Ωp(ϱ)ndx,I(t)=∫Ωϱudx. (6)

Multiplying by and taking into account equation (1), rewrite equation (2) as

Integration over the domain , in view of (3), leads to

 12ddt∫Ωϱ|u|2dx−∫Ωp(ϱ)divudx=0. (7)

The second term in (7) is expressed from the renormalized equation of continuity. Define the pressure potential from the equation

 ϱdΠdϱ−Π(ϱ)=p(ϱ). (8)

In particular, for an ideal fluid, we have

 p(ϱ)=aϱγ,Π(ϱ)=aϱγγ−1,a=const>0,γ>1. (9)

From the continuity equation (1), we have

 ∂Π∂t+div(Πu)+p(ϱ)divu=0,0

Integration of the renormalized equation of continuity (10) results in the expression

 ddt∫ΩΠdx+∫Ωp(ϱ)divudx=0.

Adding this equality to (7), we get

 ddt∫Ω(12ϱ|u|2+Π(ϱ))dx=0. (11)

We arrive to the conservation law for the total mechanical energy:

 E(t)=E(0),E(t)=∫Ω(12ϱ|u|2+Π(ϱ))dx. (12)

The equations (5), (6) and (12) are the basic conservation laws for of the problem (1)–(4).

### 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:

 A(u)φ=div(uφ). (13)

Assuming that the boundary condition (3) is satisfied for the velocity , we obtain

 (A(u)φ,1)=0. (14)

The continuity equation (1) can be written in the form of an operator–differential equation:

 dϱdt+A(u)ϱ=0,0

where the notation is used. Similarly, equation (2) is written in the form

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

 ϱ(0)=ϱ0,u(0)=u0. (17)

For the considered problem (15)–(17), the key point is the property (14) for the advection operator written in the divergent 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

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

 a(φ,ψ)=∫Ωdiv(uφ)ψdx.

Define the subspace of finite elements and the discrete operator as

 (A(u)φ,ψ)=a(φ,ψ),∀ w,u∈Vh.

Similarly to (14), we have

 (A(u)φ,1)=0. (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

 dϱdt+A(u)ϱ=0, (19)
 ϱ(0)=ϱ0,u(0)=u0. (21)

Here , with denoting –projection onto .

The solution of the problem (18)–(21) satisfies the same system of conservation laws as the solution of the problem (14)–(17) (see (5), (6), (12)).

### 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

 ϱn+1−ϱnτ+A(un+1)ϱn+1=0, (22)

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

 (ϱn+1,1)=(ϱn,1),n=0,1,...,N−1. (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

 (ϱn+1un+1,1)=(ϱnun,1)−τ∫∂Ωp(ϱn+1)ndx,n=0,1,...,N−1, (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

For the first term, we have

 ϱn+1un+1−ϱnunτun+1=12ϱn+1|un+1|2−ϱn|un|2τ+12ϱn+1|un+1|2+ϱn|un|2−2ϱnunun+1τ=12ϱn+1|un+1|2−ϱn|un|2τ+12ϱn+1−ϱnτ|un+1|2+12ϱn|un+1−un|2τ≤12ϱn+1|un+1|2−ϱn|un|2τ+12ϱn+1−ϱnτ|un+1|2.

In view of this, from (26), we obtain

From (22) and the definition of the operator , we have

 12(ϱn+1−ϱnτ,|un+1|2)+(A(un+1)(ϱn+1un+1),un+1)=(A(un+1)(ϱn+1un+1),un+1)−12(A(un+1)ϱn+1),|un+1|2)=0.

This makes possible to rewrite the inequality (27) as

 12(ϱn+1|un+1|2−ϱn|un|2τ,1)−(p(ϱn+1),divun+1)≤0. (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 :

 ϱn+1−ϱnτdΠdϱ(ϱn+1)+A(un+1)ϱn+1dΠdϱ(ϱn+1)=0. (29)

The following equality takes place:

 Π(ϱn+1)−Π(ϱn)=dΠdϱ(ϱn+1)(ϱn+1−ϱn)−12d2Πdϱ2(˜ϱn+1)(ϱn+1)(ϱn+1−ϱn)2,

where

 ˜ϱn+1∈[min(ϱn,ϱn+1), max(ϱn,ϱn+1)].
###### Assumption 2

Assume that

 d2Πdϱ2(ϱ)≥0, (30)

Under these natural assumptions, we get

 dΠdϱ(ϱn+1)(ϱn+1−ϱn)≥Π(ϱn+1)−Π(ϱn).

For the second term in (29), taking into account (8), we have

In view of this, integration of equation (29) results in

 (Π(ϱn+1)−Π(ϱn)τ,1)+(p(ϱn+1),divun+1)≤0. (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

 ϱn+1−ϱnτ+A(un)ϱn+1=0, (33)

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.

For the linearized scheme, the discrete analogs of the mass conservation law (see (24)) and the momentum conservation law (see (25)) hold.

### 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:

 ϱ0n+1=ϱn,u0n+1=un. (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:

 ϱk+1n+1−ϱnτ+A(ukn+1)ϱk+1n+1=0, (36)

where

 ϱn+1=ϱKn+1,un+1=uKn+1,n=0,1,...,N−1. (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

 Ω={x | x=(x1,x2),−5

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

 ϱ0(x)=1+αexp(−β|x|2),

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. 45, 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. Figure 3: The solution of the problem at various time moments calculated on the grid M=50: dotted line — τ=0.01, dashed — τ=0.005, solid — τ=0.0025. Figure 4: The solution of the problem at various time moments calculated on the grid M=100: dotted line — τ=0.01, dashed — τ=0.005, solid — τ=0.0025. Figure 5: The solution of the problem at various time moments calculated on the grid M=200: dotted line — τ=0.01, dashed — τ=0.005, solid — τ=0.0025.

### 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). Figure 6: The solution of the problem at various time moments calculated on the grid M=50: dashed line — K=1, dotted — K=2, solid — K=5. Figure 7: The solution of the problem at various time moments calculated on the grid M=100: dashed line — K=1, dotted — K=2, solid — K=5. Figure 8: The solution of the problem at various time moments calculated on the gridM=200: dashed line — K=1, dotted — K=2, solid — K=5.

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. 911. Here, according to (12), is calculated at each time moment

 E(tn)=∫Ω(12ϱn|un|2+Π(ϱn))dx,n=0,1,...,N.

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. Figure 9: Time–history of the total mechanical energy for various time steps obtained on the grid with M=50: dashed line — K=1, solid — K=5. Figure 10: Time–history of the total mechanical energy for various time steps obtained on the grid with M=100: dashed line — K=1, solid — K=5. Figure 11: Time–history of the total mechanical energy for various time steps obtained on the grid with M=200: dashed line — K=1, solid — K=5.

## 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, 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.