We study the stability and convergence properties of a semi-implicit time stepping scheme for the incompressible Navier–Stokes equations with variable density and viscosity. The density is assumed to be approximated in a way that conserves the minimum-maximum principle. The scheme uses a fractional time-stepping method and the momentum, which is equal to the product of the density and velocity, as a primary unknown. The semi-implicit algorithm for the coupled momentum-pressure is shown to be conditionally stable and the velocity is shown to converge in
This paper focuses on the stability and error estimates for a numerical approximation of incompressible viscous Newtonian fluids with variable density and viscosity in an open bounded domain Ω ⊂ R d , with d = 2 or 3, up to a given time T . These flows are governed by the time-dependent incompressible Navier-Stokes equations ∂ t ρ + u • ∇ρ = 0, in Ω × (0, T ],
(1.1a)
where ρ is the density, u the velocity, p the pressure, ε(u) = 1 2 (∇u + (∇u) ⊤ ) the strain rate tensor, and η(ρ) is the dynamical viscosity that we assume to be a function of the density. We note that this assumption, i.e. η a function of the density, is consistent for Newtonian multi-fluids where temperature effects are disregarded. For non-Newtonian fluid, the dynamical viscosity can be computed using viscosity models like the power-law model [10,25]. However, such models are beyond the scope of this paper. The function f represents a given source term such as gravity force. The above system is supplemented with initial conditions u(x, 0) = u 0 and ρ(x, 0) = ρ 0 in Ω,
and homogeneous Dirichlet condition on the velocity
such that no boundary condition needs to be enforced on the density as there is no inflow on the boundary. The approximation of the solution of the system (1.1a) -(1.1c), whose existence and uniqueness have been established for two-dimensional domains with constant viscosity in [13,28,31,36], presents many challenges such as enforcing the incompressiblity constraint, the tracking of the interface between fluids with different density, or more generally, the approximation of density functions that present sharp gradients. The difficulties associated to the incompressibility constraint (1.1c) can be handled using projection type methods [2,5,6,14,17,39,46] and also artificial compression techniques [33] which were both originally developed for constant density flow [8,9,29,41,42]. The tracking of the interface between fluids of different density has also been the focus of extensive studies, which led to the development of various approximation methods such as front tracking methods [43], phase field method [1,7,32] and level set techniques [37,38].
While the error analysis of methods for incompressible flows with constant density has been extensively studied, for instance see [19] and references therein, the case of incompressible flows with variable density presents more difficulties and fewer results are available. In [18], the space-time error estimates for the velocity equations are established under the assumption that the density is approached via a method that preserves the minimum-maximum principle. Recent work in [6] shows the time-convergence properties for the full set of equations (1.1a)-(1.1c) using a convective Gauge Uzawa scheme, and [4] establishes space-time error for a fully discretized algorithm that couples the velocity-pressure. We note that all the above works provide error estimates for schemes that involve a constant positive dynamical viscosity η. Moreover, these schemes involve the discretization of the quantity ρ∂ t u which leads to algorithms with time-dependent stiffness matrices that need to be reassembled at each time iteration.
To improve the computational efficiency, in [5], one of the authors introduced a time independent stiffness matrix algorithm to approximate incompressible flows with variable density and viscosity. One of the main ideas consists of using the momentum m, equal to the product of the density and the velocity, as a primary unknown, so the time derivative in equation (1.1b) can be rewritten as ∂ t m. As the diffusion operator involves the velocity and a variable viscosity, following [2,14,15], it is rewritten as -∇•(ν max ε(m)) + ∇•(ν max ε(m) -ηε(u)), allowing us to treat the first term implicitly and the correction explicitly, so the resulting stiffness matrices become time-independent when advective terms are made explicit. While this method has been studied numerically on various settings, see [23,24], no theoretical result on its stability and convergence properties has been established yet. This paper aims to fill this gap by establishing stability and temporal error estimates of a semi-implicit discretization version of the velocity-pressure coupling algorithm where we assume that the density is approximated in a way that satisfies the minimum-maximum principle.
The paper is organized as follows: in the next section, we introduce some notations and hypotheses on the dynamical viscosity and the regularity of the solutions of (1.1a)- (1.1c). In Section 3, we state our hypothesis on the approximated density, and we introduce a semi-implicit time stepping scheme for the velocity-pressure system described by the equations (1.1b)- (1.1c). The conditional stability and convergence of the algorithm are established in Sections 4-5. Numerical results are reported in Section 6, which also includes numerical investigations of an explicit version of the algorithm presented in Section 3. Concluding remarks are given in Section 7.
In this section, we define some notations that will be used in the rest of the paper. We also recall some classic results, such as Korn and Poincaré inequalities, and we introduce some hypotheses that are used later in the paper to establish the stability and convergence of the scheme proposed in Section 3.
We consider the incompressible Navier-Stokes system (1.1a)-(1.1c) with variable density and viscosity on the time interval [0, T ] and in an open, connected bounded domain, Ω ⊂ R d , with d = 2 or 3.
We first denote the time step τ := T /N , where N is the number of time iterations. For any timedependent function f , we set f (t n ) := f (nτ ), 0 ≤ n ≤ N , and we denote by f n an approximation of f at time t n . The sequence f 0 , f 1 , ..., f N is denoted by f τ .
The space W m,p (Ω) stands for the standard Sobolev spaces, with 0 ≤ m ≤ ∞, 1 ≤ p ≤ ∞. We use the simplified notation H s (Ω) for the Hilbert spaces W s,2 , and set L 2 (Ω) := H 0 (Ω). The closure of the space of infinitely differentiable functions compactly supported in Ω,
In the rest of the paper, we use bold fonts for vector-valued functions and spaces, and regular fonts for scalar-valued functions and spaces. We also remind the readers of the Korn inequality,
and Poincaré inequality,
that hold for positive constants c 1 and c 2 that are independent of the space and time discretization. In addition, the proofs in this paper use repeatedly, without reminder, the polarization identity
, and the Young inequality
that holds for any constant > 0.
In the remainder of the paper, we make the following assumptions. First, we assume that the solution of (1.1a)-(1.3) possesses the following regularity:
where the space L p (X), with X a Banach space, is a Bochner space with norm ‖u‖ L p (X) = ( T 0 ‖u(t)‖ p X dt) 1/p . Second, we assume that the dynamical viscosity η is a Lipschitz function of the density, meaning that there exists a positive constant L, such that for any reals ρ, ρ, we have
This hypothesis will allow us to bound the error in dynamical viscosity by the error in density in the convergence analysis established in Section 5. We note that the above relation is consistent with immiscible two phase field models where the dynamical viscosity varies from η 1 to η 2 while the density varies from ρ 1 to ρ 2 . In these models, the dynamical viscosity can be defined as
which is a Lipschitz function of the density with L = η2-η1 ρ2-ρ1 .
We introduce in this section the semi-implicit time-stepping algorithm used to approximate the velocitypressure couple, and we describe the hypotheses made on the given approximation of the density.
We assume that the density equation (1.1a) is approximated by a sequence ρ τ that satisfies an equation of the form:
where R n is a consistent term at each time step n, meaning that it is assumed to be small. We also assume that R n is defined is a way such that the approximated density ρ τ satisfies for all integer 0 ≤ n ≤ N , the minimum-maximum principle
and the following relations: 3) present many challenges to enforce as the mass conservation is an hyperbolic equation. Describing in details methods that can satisfy the above hypothesis is out of the scope of our paper, but we note that such properties can be achieved using monotone scheme, flux transport corrected techniques [3,7,21,27,47], discontinuous Galerkin methods with limiters technologies [11,26,45], or even entropy viscosity techniques [16]. Before describing the numerical scheme we use to approximate the couple velocity-pressure, we note that for immiscible incompressible flows, e.g. two phase flows, the condition (3.3) cannot be enforced as the gradient of the density scales with inverse of the mesh size. As shown in Sections 4-5, our stability and error analysis still holds for such problems by relaxing the hypothesis (3.3) into a condition of the form:
with γ 1 a positive constant. Introducing a spatial discretization with a mesh size h, we note that the above hypothesis becomes very restrictive as it implies τ ∝ h 2 when the density's gradient scales in h -1 . However, we will show with numerical illustrations in Section 6 that such condition does not need to be enforced, even for problems with discontinuous density, meaning that the proposed algorithm remains stable and convergent under classic CFL condition.
We introduce a semi-implicit scheme for the incompressible Navier-Stokes equations where the incompressiblity condition is enforced using a pressure-correction projection method. We note that traditional pressure splitting methods lead to the introduction of an auxiliary scalar pressure φ that is solution of the following Poisson problem with variable coefficient:
Such formulation yields two disadvantages: (1) as the density depends of space and time, the matrix associated to the above problem needs to be reassembled at each time iteration, increasing computation costs;
(2) the matrix can become ill-conditioned when the ratio of density, ρ max /ρ min , is large. Following [17,18], we overcome these difficulties by replacing ρ n+1 by a constant χ in the above equation which was proven to not impact the stability and accuracy of the resulting projection-like method. The algorithm described below is inspired from [5] where the momentum, m, is used as a primary unknown. The main difference here is that we consider a semi-implicit version of the algorithm so we can establish the stability and convergence of the algorithm in Sections 4-5. After initialization, the algorithm reads as follows: for 0 ≤ n ≤ N ,
Step 1. Update the density ρ n+1 using an algorithm that satisfies the conditions (3.1)-(3.4).
Step 2. The momentum m n+1 is the solution of
where m *,n := ρ n+1 u n and ν is a positive constant defined in the following section. The velocity u n+1 , that does not satisfy the incompressiblity condition (1.1c), is then computed as follows
(3.7)
Step 3. The auxiliary pressure φ n+1 is computed by solving
where χ is a positive constant.
Step 4. The pressure p n+1 and the incompressible velocity field ũn+1 are updated as follows:
(3.10)
The analyses performed in Sections 4-5 show that the positive constant ν and χ, used in (3.6)-(3.8), can be defined as follows:
We note that for Newtonian fluid, bounds on the density and dynamical viscosity can be derived using the properties of the fluid considered under classic condition (e.g. given temperature) such that the above constants can be computed with this knowledge before approximating the solutions of (1.1a)-(1.1c). Also, to show the stability and convergence of the algorithm in the following sections, we assume the time step τ and positive constants γ 0 , γ 1 and γ 2 satisfy the following conditions:
As mentioned in Section 3.1, we remind that the condition (3.14) is heuristic, meaning that in practice, the resulting algorithm is stable and convergent under classic CFL condition even for problems with density function that either presents large gradient or is discontinuous, see Section 6.2.
This section is dedicated to establishing the conditional stability of the velocity-pressure semi-implicit scheme (3.6)-(3.10) described in Section 3.2. The main result is summarized in the following theorem.
where C is a positive constant that depends on the geometry of the domain Ω, the fluid's properties, e.g. η min , and the positive constant γ 1 defined in Section 3.3.
Proof. We start by multiplying equation (3.6) by 2τ u n+1 and integrate over Ω to obtain
We note that the following identities holds
Multiplying the mass conservation (3.1) by τ |u n+1 | 2 yields the following
Moreover, using Cauchy-Schwarz inequality combined with the definition of m *,n = ρ n+1 u n and the fact that the the density and viscosity are bounded away from zero, we get the following inequality
Substituting the above identities into equation (4.2) yields the following inequality:
3)
The following aims to bound the terms S i . Using the hypothesis (3.11) on ν and the minimum-maximum principle (3.2), we get
.
Using the assumptions (3.11) and (3.15), the term S 2 can be bounded by
.
Following [17], see equations (3.11)-(3.13) therein, we get the following bound for S 3 :
Applying Young inequality, with a constant 1 > 0 whose value is set later in the proof, we get
, where we apply Korn and Poincaré inequalities (2.1)-(2.2) to obtain the second inequality. Using the hypothesis (3.14) and applying Young inequality to the term S 5 reads
Substituting the above bounds for the terms S i back in (4.3) we obtain the resulting inequality
where we drop the positive term
in the above left handside. Using the hypothesis (3.12) and setting 1 = ηmin c1c2 (1 -8ν 2 γ 2 1 ), a positive term thanks to the hypothesis 3.14, we can drop the terms of the second line in the above inequality. We conclude the proof by applying a standard telescopic argument which yields the stability estimate (4.1) with C = -1 1 .
Remark 1. We note that the hypothesis (3.14) can be replaced by the hypothesis (3.3), i.e. assuming that the gradient of the approximated density can be bounded in L ∞ -norm independently of the time step. In such case, the above result still holds if we assume that the time step satisfies
As discussed in the Section 3, the assumption (3.3) cannot be considered for general setups of immiscible flows where the gradient of the approximated density, ρ τ can become very sharp around the fluids interface, i.e. scaling with the inverse of the mesh size which led us to use the assumption (3.14). We note that in such case, this hypothesis is very restrictive as it can be read as the time step is proportional to the square of the mesh size. We will show in the Section 6 that it is not necessary to enforce such restrictions even for problem with discontinuous density, meaning that the algorithm remains stable under classic CFL condition.
In this section, we carry out the time error analysis of the scheme introduced in Section 3 for the couple velocity-pressure. The main result is described in Theorem 3. First, we introduce the following quantities that will be used to represent the errors in density, pressure, velocity, and auxiliary velocity, respectively:
Using hypothesis (2.4), we note that the error for the dynamical viscosity η, defined by e n η := η(t n ) -η n , is bounded by the error in density, e ρ , as we can write e n η ≤ Le n ρ .
(5.5)
Before performing the error analysis of the scheme (3.6)-(3.10) in Section 5.3, we introduce an equivalent scheme that uses the velocity as primary unknown, and we describe our hypothesis on the density error in the following sections.
First, we show that the steps (3.6)-(3.7) of the algorithm can be replaced by the following update: Find u n+1 solution of
(5.6)
The main difference with the algorithm introduced in Section 3.2 is that the momentum is not used as primary unknown. While rewriting the algorithm in this form allows us to perform the following analysis, we note that in practice, using the momentum as primary unknown is more efficient as it results in a mass matrix that does not depend on the density and can be assembled at initialization. Proof. We note that equation (3.1) can be rewritten as
Thanks to the above identity, the definition of m n+1 := ρ n+1 u n+1 , and the product rule, equation (3.6) is equivalent to (5.6). Thus, both schemes are equivalent.
In addition to the hypothesis (3.1)-(3.5), we assume that the approximated density ρ τ and the residual R n , defined in (3.1), satisfy the following inequalities
(5.9)
These assumptions are consistent with the use of a first order algorithm for the density equation that preserves the density's lower and upper bounds.
Our purpose in the rest of this section is to establish a priori estimates in L 2 (Ω) for the solutions of the scheme (5.6)-(3.8)-(3.9).
The proof is split into three lemmas, where we introduce error equations for the velocity, for the pressure, and bounds for some of the resulting terms, before we establish the final error estimate in Theorem 3. In the following, we assume that the approximated density results from an algorithm of the form (3.1) and that the velocity, pressure and auxiliary pressure approximation, (u τ , p τ , φ τ ), are computed using the algorithm (5.6)-(3.8)-(3.9).
First, the equation that controls the error on the velocity is given in the following lemma.
Lemma 1. Assume that the solution to (1.1a)-(1.3) satisfies the regularity hypothesis (2.3) and the minimummaximum principle (3.2). Then the following holds for all integers 0
where
Proof. Using the mass conservation (1.1a) and the incompressibility condition (1.1c), we rewrite equation (1.1b) as follows
Taking the difference of the above equation at time t n+1 and equation (5.6), multiplying the difference by 2τ e n+1 u and integrating over Ω, we get
(5.
Combining the above identity with equation (5.11) yields the error equation (5.10) for e u .
The following lemma establishes the error equation for the pressure.
Lemma 2. Under the same assumptions as Lemma 1, we get
(5.12)
Proof. The proof follows closely some arguments from the analysis of pressure correction projection methods such as the one in [18] that we report them here for completeness.
Combining the equations (3.8)-(3.9), we write Δ(e n+1 p -e n p ) = χ τ ∇•e n+1 u + Δ(p(t n+1 ) -p(t n )). (5.13) Taking the inner product with 2τ 2 χ (e n+1 p -2e n p + e n-1 p
), we obtain
(5.14)
We then multiply (5.13) with 2τ 2 χ e n+1 p . It reads
(5.15)
Taking the difference of (5.14) and (5.15), we obtain the following equality
) .
(5.16)
Next, we multiply the difference of (5.13) at times t n+1 and t n with τ (e n+1 p -2e n p + e n-1 p
which implies
(5.17)
Combining equations (5.16) and (5.17), we get the desired inequality (5.12).
We now outline the bounds for the residual terms R n+1 η , R n+1 N L1 , R n+1 N L2 defined in (5.10) in the following Lemma.
Lemma 3. Under the assumptions (2.3)-(2.4), (3.11)-(3.12), and (5.8)-(5.9), there exists a positive constants C, independent on the time step, such that the following holds for any arbitrarily small positive constants 1 , 2 , 3 , 4 , 5 that are also not dependent on the time step.
(5.20)
Proof. First, we use integration by parts to rewrite 2τ
The first term, A 1 , is bounded as follows
where we use the Young inequality and the assumption (3.3) to obtain the last inequality. Applying polarization identity on the term A 2 gives
) .
Combining the hypothesis (3.11) and (3.15), we get
.
Using the regularity hypothesis (2.3), the term A 2 can be bounded as follows
In a similar use of polarization identity and regularity (2.3), we get:
Using Young inequality, the regularity hypothesis (2.3), and the Lipschitz condition (2.4), we get
Summing up these A i terms, we obtain the inequality (5.18).
The first nonlinear residual term is reformulated as
then applying Cauchy-Schwarz inequality and regularity (2.3), we get
Using the definition of ũn , see (3.10), we get the following identity
and thus, by regularity (2.3), we obtain
which yields the inequality (5.19).
Eventually, the second nonlinear residual term is bounded as follows
thanks to minimum-maximum principle (3.2) and regularity (2.3).
We are now ready to state the main result of this section. 12), and (5.8)-(5.9) hold, and that the initial conditions are chosen such that e 0 u = e 0 p = 0. Assume also that the time step is small enough, meaning that it satisfies
where C is a positive constant that is independent of the time step τ .
Proof. In the following, we denote by C and C 1 some generic constants that only depend on the geometry of the domain, the final time T , the minimum and maximum of the density and dynamical viscosity, and the regularity of the strong solutions (ρ, u, p). Thus, these constants are independent of the time discretization. Adding the relations (5.10) and (5.12), and using the hypothesis (3.2), yields the following
(
The following shows how the terms of the right hand side of the above inequality are bounded. First, we note that the term R n+1 u can be rewritten as
Using regularity assumptions (2.3) and minimum-maximum principle (3.2), we bound this term as follows
.24) Next, we note that the identity below holds R n+1 p , e n+1 u = -2τ ⟨e n u , ∇(p(t n+1 ) -2p(t n ) + p(t n-1 ))⟩ -2τ e n+1 u , ∇(2e n p -e n-1 p ) + 2τ e n+1 u e n u , ∇(p(t n+1 ) -2p(t n ) + p(t n-1 )) . Combining all the pressure terms of the right hand side, applying Poincaré and Korn inequality, we obtain 2τ R n+1 p , e n+1 u + 2τ e n+1 u , ∇(2e n p -e n-1 p ) + 2τ e n+1 u e n u , ∇(p(t n+1 ) -2p(t n ) + p(t n-1 )) + 2τ 2 χ ∇(p(t n+1 ) -p(t n )), ∇(2e n p -e n-1 p
(5.26)
We introduce the variable A n defined by
and we set the positive constant i as follows
We then drop the second term of the first and third line of the equation (5.26) as they can be shown to be positive using (5.21). The above inequality reads
where C 1 is a positive constant independent of the time step. Using the assumptions (5.8)-(5.9), and summing over n, from 0 to N -1, yields
where we use the relation N τ = T with T being the final time. Reminding that the time step τ is chosen small enough, so that 2C 1 τ < 1, the proof is concluded using discrete Grönwall lemma.
Proof. The proof is similar to that of Theorem 3. The only difference arises when bounding the term A 1 in the proof of Lemma 3. The bound now reads as follows
where we use the Young inequality and the assumption (3.14) to obtain the last inequality. As the bound now involves a term in O(τ 2 ), versus a O(τ 3 ) for the proof of Theorem 3, the resulting order of convergence in L 2 norm is reduced by one half.
Remark 2. We state the above result to show that the above analysis still holds for more complex problems where the density is allowed to present very sharp gradients, or discontinuity. While one could expect a better rate of convergence for the velocity, i.e one, we note that in such settings the density, and so the momentum, becomes singular so the best order of convergence rate one can expect for these quantities is half in L 1 norm. We refer to Section 6.2.1 where numerical tests are performed with a discontinuous density to show that the algorithm is stable under classic CFL condition and converges with order one for the velocity-pressure couple in L 2 norm and for the density in L 1 norm.
In this section, we study the space-time convergence properties of the algorithm described in Section 3. First, we introduce a spatial discretization of the algorithm (3.1)-(3.9) using finite element methods. Second, we test the accuracy of our semi-implicit algorithm using smooth and nonsmooth manufactured solutions. We conclude by introducing an explicit version of our scheme and by comparing its convergence and stability properties to the above semi-implicit scheme. All tests are performed on two dimensional domains, denoted by Ω, using the code FreeFEM++. More information on this software can be found in [22].
The spatial discretization is done using a continuous Galerkin finite element method. We introduce a conforming, shape regular mesh sequence E h of the domain Ω that consists of simplex elements. The mesh size is defined by h = max K∈E h h K , where h K represents the diameter of a cell K. The unknowns (ρ, m, p) are approximated using the following spaces respectively:
)
where P k represents the vector space of polynomials with total degree of at most k. We note that the pair (X h , M h ) used to approximate the couple velocity-pressure is the P 2 -P 1 Taylor-Hood approximation space.
The fully discrete problem, using the time discretization algorithm (3.1)-(3.9), reads as follow. Find ρ n+1 ∈ X h such that the following holds for all ψ ∈ X h
holds for all v ∈ X h where m *,n = ρ n+1 u n . The velocity u n+1 is then updated using (3.7). Find φ n+1 ∈ M h such that:
The pressure p n+1 and the velocity ũn+1 can then be updated using (3.9)-(3.10). Notice that the above formulation of the mass equation has been stabilized with an artificial viscosity ν h as the density can present large gradients, which represents the term R n+1 introduced in Section 3.1. There exist many stabilization techniques in literature such as entropy-based viscosity [5,16], residual viscosity [34,40], and first order viscosity (i.e. ν h is made proportional to O(h) with h the mesh size). In the following, we use an h-viscosity to stabilize the equation. While such artificial viscosity has a strong diffusive effect, it is one of the most robust. Moreover, it will not impact the convergence properties of the algorithm under classic CFL condition, i.e. τ = O(h), as we show in the following that the algorithm converges with order one in time for smooth density functions. We also note that the above formulation does not guarantee minimum/maximum principle of the density due to the discretization in space. Maximumpreserving methods can be construct using flux transport corrected techniques [3,7,21,27,47]. However, the simulations reported in the following conserve the density's bound up to a few percents so we did not implement such techniques.
In the following of the paper, the constants ν and χ are defined as follows:
which satisfies the hypotheses (3.11)-(3.12) used in Sections 4-5 to establish the stability and convergence properties of the semi-implicit scheme.
In this section, we use manufactured solutions to check the accuracy of the above semi-implicit method. We perform convergence tests on four different setups with various choices of dynamical viscosity functions and ratios of magnitude of fluid parameters. The first test, referred to as Test 1, considers a nonsmooth density function and a constant dynamical viscosity. The goal of the test is to show that the condition (3. 14) is heuristic so that, in practice, the proposed scheme remains stable and convergent under classic CFL condition τ ∼ h. The subsequent tests aim to study the convergence property of our scheme for various dynamical viscosity functions where the density is smooth. The second test, Test 2, uses a density and velocity that satisfy the mass conservation equation (1.1a). The dynamical viscosity is made linearly dependent on the density, and the kinematic viscosity, ν = η/ρ, is constant. The third test, Test 3, considers a more complex set of solutions which leads us to add a non-zero source term in the mass conservation equation (1.1a). It allows us to study the convergence of our method on problems where the kinematic viscosity varies in space and time and where the diffusion term in the Navier-Stokes equations is non-zero. Finally, in Test 4, the dynamical viscosity is chosen to be a nonlinear Lipschitz function of the density by setting η = ρ -1 .
We start our numerical investigations section by considering a problem where the density function is discontinuous. The domain is set to the unit disk and the time domain is set to [0, T ] = [0, 1]. We consider the following exact solutions:
1 else, where θ = arctan( y x ). The dynamical viscosity η is set to one so the source term f in the momentum equation can be computed as follows f = -0.25xρ(x, y, t) + cos(x) sin(y) sin(t) -0.25yρ(x, y, t) + sin(x) cos(y) sin(t) .
We note that although the velocity is time-independent, the momentum, equal to ρu, is not. Moreover, both the density and momentum are discontinuous functions. Thus, we compute the density errors using L 1 norm, as the best rate of convergence in space we can expect is one in L 1 norm. We also note that the dynamical viscosity is set to one so the source term does not involve a Dirac delta function. As shown in [5], stabilizing the mass conservation equation with a h-viscosity can lead to a poor approximation when the density presents discontinuity or a sharp gradient around an interface Σ. To obtain a density profile that is not overly diffused, we combine this stabilization method with a compression technique introduced in [12,35,44]. In a nutshell, an additional term, equal to
, is added in the left hand side of the mass conservation equation.
Our results, summarized in Table 1, show that the algorithm converges with order one when the time step is set to τ = h/10, i.e. the algorithm is stable and convergent under a CFL condition τ ∼ h even when the gradient of the density scales with the inverse of the mesh size. The main conclusion of this investigation is that the condition (3.14), that reads τ ≤ Ch 2 when the density gradient is sharp, is mostly heuristic and does not need to be enforced in practice.
τ Velocity Pressure Density L 2 error Order L 2 error Order L 1 error Order 0.05 4.21E-3 -8.74E-2 -6.71E-2 -0.025 1.04E-3 2.02 2.87E-2 1.60 3.50E-2 0.94 0.0125 2.04E-4 2.35 8.54E-3 1.75 1.82E-2 0.95 0.00625 5.97E-4 1.77 3.09E-3 1.47 8.91E-3 1.03 0.003125 1.65E-5 1.86 1.55E-3 1.00 4.46E-3 1.00
Tab. 1: Results of convergence test for nonsmooth density with τ = h/10. Semi-implicit scheme with ν h = 0.125h.
In this test, the domain is set to Ω = {(x, y) ∈ R 2 : x 2 + y 2 < 1} and Dirichlet boundary conditions are applied for the density and the momentum. We consider the following exact solutions u(x, y, t) = 1 2 -y cos(0.5t)
x cos(0.5t) , p(x, y, t) = sin(x) sin(y) sin(t),
where θ = arctan( y x ). Moreover, ρ 1 and ρ 2 are constants satisfying 0 < ρ 1 < ρ 2 that allow us to control the ratio of density. The source term f in the momentum equation can then be computed as follows f = 0.25ρ(x, y, t) y sin(0.5t) -x cos 2 (0.5t) + cos(x) sin(y) sin(t) -0.25ρ(x, y, t) x sin(0.5t) + y cos 2 (0.5t) + sin(x) cos(y) sin(t) .
We note that although the velocity u can be represented exactly in X h , the momentum m = ρu, one of the primary unknowns of the algorithm, cannot. We first perform a series of tests using a small ratio of density and viscosity's magnitudes by setting (ρ 1 , ρ 2 ) = (η 1 , η 2 ) = (1, 3). We use five different grids of respective mesh sizes h = 0.1, 0.5, 0.25, 0.125, and 0.0625. To run simulations with a constant CFL of order 1, we define the time step as τ = 0.5h and set the final time to T = 1. The L 2 relative errors at the final time T of the velocity, pressure, and density are displayed in Table 2. We recover a rate of convergence close to or larger than one for all the unknowns. We note that initially the density convergence rates are smaller than one, which could be improved by using a less diffusive stabilization method than the h-viscosity. These results are compatible with the theoretical order of convergence O(τ ) for smooth density functions as the error in space is bounded by O(h), due to the use of h-viscosity in the mass equation, and that the mesh size is proportional to the time step τ .
τ Velocity Pressure Density L 2 error Rate L 2 error Rate L 2 error Rate 0.05 3.01E-3 -8.52E-2 -8.67E-3 -0.025 6.62E-4 2.18 2.27E-2 1.91 4.71E-3 0.88 0.0125 1.81E-4 1.87 8.98E-3 1.34 2.32E-3 1.02 0.00625 5.69E-5 1.67 4.01E-3 1.16 1.04E-3 1.17 0.003125 1.99E-5 1.52 2.11E-3 0.93 5.35E-4 0.95 Tab. 2: Results of Test 2 for smooth density with (ρ 1 , ρ 2 ) = (1, 3) and τ = h/2. Semi-implicit scheme with ν h = 0.125h and η = ρ.
For the second series of tests, we consider solutions with a larger ratio of density by setting (ρ 1 , ρ 2 ) = (η 1 , η 2 ) = (1, 100). We perform five runs using the above mesh sizes and time steps. The L 2 relative errors at the final time T are displayed in Table 3. Both series of tests confirm that the algorithm converges with order one in time. We note that the errors in pressure are initially larger when a large ratio of density is considered due to the modification of the projection step (3.8) where ρ is replaced by a constant χ set to 0.99ρ min . We also note that, although the dynamical viscosity is variable, the above tests are performed
τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 1.92E-3 -1.13E0 -1.57E-2 -0.025 6.29E-4 1.61 9.09E-1 0.31 8.54E-3 0.88 0.0125 2.13E-4 1.56 4.28E-1 1.09 4.23E-3 1.01 0.00625 7.12E-5 1.58 1.97E-1 1.12 1.97E-3 1.11 0.003125 3.25E-5 1.13 9.49E-2 1.05 9.84E-4 1.00 Tab. 3: Results of Test 2 for smooth density with (ρ 1 , ρ 2 ) = (1, 100) and τ = h/2. Semi-implicit scheme with ν h = 0.125h and η = ρ.
with η = ρ, meaning the kinematic viscosity ν is constant. Thus, the following setup focuses on problems with variable kinematic viscosity.
This test aims to study the impact of variable dynamical and kinematic viscosity of the convergence properties of the proposed semi-implicit scheme. A setup inspired from [30] is considered. The manufactured solutions are given by u(x, y, t) = 3 4 + 1 4 sin(t) -sin 2 (x) sin(y) cos(y) sin(x) cos(x) sin 2 (y) , p(x, y, t) = sin(x) sin(y) sin(t),
We note that unlike the problem considered in the previous section, the couple (ρ, u) is not solution of the mass conservation (1.1a). Thus a source term, denoted by f ρ , is added to the right hand-side of the mass conservation. The source terms f ρ and f are computed accordingly such that (ρ, u, p) is a solution of the system (1.1a)-(1.1b). We note that this setup is more challenging than in the previous section as the right hand-side of the Navier-Stokes equations now depends on the dynamical viscosity η as the strain-rate tensor, ε(u), is non zero. The computational domain is kept to the unit disk and the time interval is set to [0, T ] = [0, 1]. To study the impact of variable kinematic viscosity, we set (ρ 1 , ρ 2 ) = (1, 100) and consider three set of parameters (η 1 , η 2 ) for the dynamical viscosity: (1, 50), (1,10), and (1, 2). As the dynamical viscosity goes from η 1 to η 2 when the density varies from ρ 1 to ρ 2 , the effective ratios of kinematic viscosity considered here are 2, 10 and 50. Our simulations are performed with the same mesh size and time step as those in the previous section, meaning that τ = h/2 and h = 0.1, 0.5, 0.25, 0.125, and 0.0625. Our results are summarized in Tables 4, 5, and 6. They show that the algorithm converges with order 1 under a CFL of order 1 for problems with variable dynamical and kinematic viscosities.
τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 2.18E-2 -3.54E0 -1.04E-2 -0.025 5.60E-3 1.96 1.04E0 1.77 6.14E-3 0.75 0.0125 1.72E-3 1.70 5.23E-1 0.99 3.21E-3 0.94 0.00625 6.35E-4 1.44 2.07E-1 1.34 1.53E-3 1.07 0.003125 2.68E-4 1.24 8.38E-2 1.31 7.55E-4 1.01 Tab. 4: Results of Test 3 for smooth density with (ρ 1 , ρ 2 ) = (1, 100), (η 1 , η 2 ) = (1, 50), and τ = h/2. Semi-implicit scheme with ν h = 0.125h. 6.2.4 Test 4: Nonlinear dynamical viscosity This last test uses the same setup than Test 3 described in section 6.2.3 modulo that the dynamical viscosity is now chosen to be a nonlinear function of the density. Specifically, we set η(x, y, t) = 1 ρ(x, y, t) . τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 3.02E-2 -2.38E0 -1.05E-2 -0.025 1.09E-2 1.47 6.22E-1 1.94 6.10E-3 0.78 0.0125 4.59E-3 1.25 1.91E-1 1.70 3.20E-3 0.93 0.00625 2.24E-3 1.03 7.21E-2 1.41 1.45E-3 1.14 0.003125 1.09E-3 1.05 2.86E-2 1.33 7.51E-4 0.95 Tab. 5: Results of Test 3 for smooth density with (ρ 1 , ρ 2 ) = (1, 100), (η 1 , η 2 ) = (1, 10), and τ = h/2. Semi-implicit scheme with ν h = 0.125h. τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 3.90E-2 -2.30E0 -1.05E-2 -0.025 1.61E-2 1.27 4.34E-1 2.41 6.09E-3 0.78 0.0125 7.68E-3 1.07 1.33E-1 1.71 3.20E-3 0.93 0.00625 4.01E-3 0.94 5.39E-2 1.30 1.45E-3 1.14 0.003125 2.01E-3 1.00 2.54E-2 1.09 7.49E-4 0.95 Tab. 6: Results of Test 3 for smooth density with (ρ 1 , ρ 2 ) = (1, 100), (η 1 , η 2 ) = (1, 2), and τ = h/2. Semi-implicit scheme with ν h = 0.125h.
We note that the definition of the density yields ρ(x, y, t) ∈ [ρ 1 , ρ 2 ] for all (x, y) ∈ Ω and t ∈ [0, T ]. Thus, the above dynamical viscosity is a Lipschitz function of the density with Lipschitz constant L = ρ -2 1 . The computational domain is set to Ω = {(x, y) ∈ R 2 : x 2 + y 2 < 1} and the time interval is set to [0, T ] = [0, 1]. We perform three series of tests with (ρ 1 , ρ 2 ) = (1, 3), (ρ 1 , ρ 2 ) = (1, 50) and (ρ 1 , ρ 2 ) = (1, 100). As a result, the kinematic viscosity ν = η/ρ lives in the intervals [0.11, 1], [4 × 10 -4 , 1], and [10 -4 , 1] respectively. The L 2 relative errors at final time are displayed in Tables 7-8-9. They confirm that the algorithm converges with order one even for problems where the dynamical viscosity is a nonlinear function of the density and where the kinematic viscosity presents large ratio of magnitude.
τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 3.52E-2 -1.86E-1 -5.78E-3 -0.025 6.37E-3 2.47 4.11E-2 2.18 3.37E-3 0.78 0.0125 2.40E-3 1.40 1.32E-2 1.64 1.76E-3 0.93 0.00625 1.12E-3 1.10 3.97E-3 1.73 8.41E-4 1.07 0.003125 5.27E-4 1.09 1.22E-3 1.70 4.20E-4 1.00 Tab. 7: Results of Test 4 for smooth density with (ρ 1 , ρ 2 ) = (1, 3), η = 1/ρ and τ = h/2. Semi-implicit scheme with ν h = 0.125h.
In this section, we study the accuracy of an explicit version of the semi-implicit scheme (6.4)-(6.6). A similar scheme has originally been introduced by one of the authors in [5] and has been tested with a pseudo-spectral code named SFEMaNS [20]. The main difference here is that we do not use a level set technique to approximate the mass conservation equation and that the space discretization is done using finite elements, see Section 6.1. In addition to be suitable for spectral codes, another motivation for investigating the accuracy of explicit schemes is that the stiffness matrices of the mass and momentum
τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 4.04E-2 -6.48E-1 -1.00E-2 -0.025 1.75E-2 1.21 1.90E-1 1.77 5.97E-3 0.75 0.0125 9.71E-3 0.85 5.36E-2 1.82 3.13E-3 0.93 0.00625 5.14E-3 0.92 2.51E-2 1.10 1.49E-3 1.07 0.003125 2.62E-3 0.97 1.23E-2 1.02 7.45E-4 1.00 Tab. 8: Results of Test 4 for smooth density with (ρ 1 , ρ 2 ) = (1, 50), η = 1/ρ and τ = h/2. Semi-implicit scheme with ν h = 0.125h. τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 3.97E-2 -8.97E-1 -1.02E-2 -0.025 1.84E-2 1.11 4.33E-1 1.05 6.07E-3 0.74 0.0125 1.03E-2 0.84 1.31E-1 1.73 3.18E-3 0.94 0.00625 5.33E-3 0.95 5.04E-2 1.37 1.51E-3 1.07 0.003125 2.73E-3 0.97 2.49E-2 1.01 7.57E-4 1.00 Tab. 9: Results of Test 4 for smooth density with (ρ 1 , ρ 2 ) = (1, 100), η = 1/ρ and τ = h/2. Semi-implicit scheme with ν h = 0.125h.
equations become time-independent. As a consequence, these matrices only need to be assembled and preconditioned at initialization. On the other hand, the semi-implicit scheme studied earlier requires a reassembling of these matrices at each time iteration which can hinder the computational performance of the algorithm. The following compares the convergence and stability properties of the explicit scheme with the semi-implicit scheme on various setups.
The time discretization of the explicit version of our algorithm reads as follows. After initialization of the unknowns, solve the following sequential scheme for n ≥ 0.
)
We note that in addition to making explicit the nonlinear terms ũ • ∇ρ and (u • ∇)m, we also replace the projected velocity ũ by u in the mass conservation equation. This replacement allows us to avoid computing ũ using equation (3.10) without impacting the accuracy of the resulting algorithm as shown in the following.
We test the explicit scheme with the setup of Test 2 described in Section 6.2.2. Similarly to the semi-implicit scheme, we perform simulations with the sets of parameters, (ρ 1 , ρ 2 ) = (η 1 , η 2 ) = (1, 3) and (1, 100). Our numerical results are summarized in Tables 10-11. They recover a similar behavior compared to the semi-implicit algorithm in the sense that both algorithms converge with first order under a CFL condition of order 1 and that the respective L 2 relative errors are of similar order of magnitude.
τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 1.64E-3 -8.53E-2 -2.27E-2 -0.025 4.33E-4 1.92 3.04E-2 1.49 1.40E-2 0.70 0.0125 1.62E-4 1.42 9.04E-3 1.75 7.48E-3 0.90 0.00625 6.69E-5 1.27 2.59E-3 1.80 3.59E-3 1.06 0.003125 3.43E-5 0.97 1.04E-3 1.32 1.89E-3 0.93 Tab. 10: Results of Test 2 with (ρ 1 , ρ 2 ) = (η 1 , η 2 ) = (1, 3) and τ = h/2. Explicit scheme with ν h = 0.125h. τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 2.72E-3 -1.18E0 -4.05E-2 -0.025 8.70E-4 1.65 4.84E-1 1.29 2.53E-2 0.68 0.0125 3.42E-4 1.35 1.64E-1 1.56 1.36E-2 0.90 0.00625 1.60E-4 1.10 6.94E-2 1.24 6.52E-3 1.06 0.003125 8.01E-5 1.00 3.42E-2 1.02 3.34E-3 0.96 Tab. 11: Results of Test 2 with (ρ 1 , ρ 2 ) = (η 1 , η 2 ) = (1, 100) and τ = h/2. Explicit scheme with ν h = 0.125h.
The convergence of the explicit scheme on problems with variable kinematic viscosity and nonlinear dynamical viscosity is also studied using the setups described in Sections 6.2.3 (Test 3) and 6.2.4 (Test 4). First, we use the explicit scheme to approximate the solution of Test 3 with (ρ 1 , ρ 2 ) = (1, 100) and (η 1 , η 2 ) = (1, 2) so that the kinematic viscosity presents a ratio of magnitude equal to 50. Then, we consider the setup of Test 4, where η = ρ -1 , with (ρ 1 , ρ 2 ) = (1, 100). Our results are summarized in Tables 12 and 13. To facilitate the comparisons with the semi-implicit algorithm, see Tables 6 and 9, all the simulations are performed with the same temporal and spatial discretization, i.e. we set τ = h/2. We observe that, similarly to our simulations with Test 2, the explicit algorithm displays a rate of convergence akin to that of the semi-implicit scheme for simulations with a CFL of order 1. The following section investigates the possible advantage of using the semi-implicit for simulations performed with large time step, meaning CFL larger than one, over long-time integration.
τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 3.90E-2 -2.31E0 -1.12E-2 -0.025 1.62E-2 1.27 4.25E-1 2.44 6.54E-3 0.78 0.0125 7.69E-3 1.07 1.32E-1 1.69 3.47E-3 0.91 0.00625 4.01E-3 0.94 5.35E-2 1.30 1.59E-3 1.13 0.003125 2.01E-3 1.00 2.52E-2 1.09 8.21E-4 0.95
Tab. 12: Results of Test 3 for smooth density with (ρ 1 , ρ 2 ) = (1, 100), (η 1 , η 2 ) = (1, 2) and τ = h/2. Explicit scheme with ν h = 0.125h.
To conclude our numerical illustrations, we propose to study the behavior of the semi-implicit and explicit algorithms for long-time integration simulations with various time steps. Our goals are to determine if the semi-implicit scheme offers a better stability property, in the sense that if it allows the use of larger time
τ Velocity Pressure Density L 2 error Order L 2 error Order L 2 error Order 0.05 3.96E-2 -9.00E-1 -1.10E-2 -0.025 1.85E-2 1.10 4.25E-1 1.08 6.58E-3 0.74 0.0125 1.03E-2 0.84 1.29E-1 1.72 3.43E-3 0.94 0.00625 5.33E-3 0.95 4.95E-2 1.38 1.65E-3 1.06 0.003125 2.73E-3 0.97 2.46E-2 1.01 8.28E-4 1.00 Tab. 13: Results of Test 4 for smooth density with (ρ 1 , ρ 2 ) = (1, 100), η = 1/ρ and τ = h/2. Explicit scheme with ν h = 0.125h.
steps compared to the explicit scheme, and to check if the semi-implicit algorithm offers a better accuracy for long-time integration simulations. We consider two setups: Test 2 with (ρ 1 , ρ 2 ) = (η 1 , η 2 ) = (1, 100), and Test 4 with (ρ 1 , ρ 2 ) = (1, 100). While both setups involve large ratio of density, we note that they still present different challenges, e.g. large variation of dynamical viscosity for the Test 2, and large variation of kinematic viscosity with a nonlinear dynamical viscosity η(ρ) (Test 4). We set the mesh size to h = 1/40 and the final time to T = 100. L 2 relative error at final time obtained with the semi-implicit and explicit algorithms for a large range of time steps are displayed in Tables 1415. Our results for the Tests 2, where the dynamical viscosity is a linear function of the density, show that the semi-implicit scheme remains stable for larger time steps compared to the explicit scheme, with an observed factor of 10. We also note that the errors obtained with the explicit scheme are always larger than the one obtained with the semi-implicit algorithm (with a factor of around 2). For the Test 4, where the dynamical viscosity is a nonlinear function of the density, both the semi-implicit and the explicit schemes behave in a similar way. We note that simulations are stabilized with a first order artificial viscosity. A less diffusive stabilization method, such as an entropy or residual based viscosity [16,34], could be used to improve the accuracy of both schemes. Overall, we conclude from our numerical investigations that the semi-implicit scheme always outperforms or at least matches the stability and accuracy of the the explicit scheme. While both schemes behave similarly for simulations with CFL of order one, the semi-implicit scheme seems to be more robust than the explicit scheme as it also allows the use of larger time-step for some setups.
We introduced a semi-implicit scheme for the incompressible Navier-Stokes equations with variable density and viscosity inspired from an explicit scheme introduced in [5]. Although the explicit scheme has been successfully applied to complex problems in engineering such as aluminum reduction cell [23] and liquid metal batteries [24], no theoretical analysis of the convergence and stability of the coupled density-velocity-pressure
τ Semi-Implicit Explicit Velocity Pressure Density Velocity Pressure Density 0.1 NaN NaN NaN NaN NaN NaN 0.05 1.27E-1 3.62E0 7.70E-2 1.26E-1 3.71E0 7.43E-2 0.025 4.30E-2 8.83E-1 5.94E-2 4.34E-2 8.73E-1 5.97E-2 0.01 1.69E-2 2.45E-1 3.64E-2 1.70E-2 3.43E-1 3.64E-2 Tab. 15: Stability semi-implicit vs explicit scheme for Test 4 with (ρ 1 , ρ 2 ) = (1, 100). L 2 error at final time T = 100 with mesh size set to h = 1/40 and ν h = 0.125h.
problem has been established. This paper fills this gap by showing the stability of its semi-implicit version and by establishing error estimates for the time discretized algorithm. The order of convergence is shown to be either one when the gradient of density is bounded or half when the density gradient is controlled by the time step, meaning the density's gradient scales with the inverse of the mesh size for simulations under CFL condition. Numerical investigations with manufactured solutions show that the resulting algorithm converges with order one in time under classic CFL condition. Numerical comparisons between the semi-implicit and explicit schemes hint that the semi-implicit scheme can yield a better stability and accuracy properties for some setups, while the explicit scheme has the advantages to be more suitable for spectral methods and to involve time-independent stiffness matrices. However, both algorithms behave similarly when working with a CFL of order one. For future work, we plan to extend this model to problems coupled with the temperature equation, where the heat capacity and heat diffusion depend on space and time, and to extend this numerical analysis for a fully discretized algorithm using finite element methods.