Analysis of an augmented mixed-FEM for the Navier-Stokes problem

Analysis of an augmented mixed-FEM for the Navier-Stokes problem∗ Jessika Cama˜no† , Ricardo Oyarz´ua‡and Giordano Tierra§ Abstract In this paper we propose and analyze a new augmented mixed finite element method for the Navier-Stokes problem. Our approach is based on the introduction of a “nonlinearpseudostress” tensor linking the pseudostress tensor with the convective term, which leads to a mixed formulation with the nonlinear-pseudostress tensor and the velocity as the main unknowns of the system. Further variables of interest, such as the fluid pressure, the fluid vorticity and the fluid velocity gradient, can be easily approximated as a simple postprocess of the finite element solutions with the same rate of convergence. The resulting mixed formulation is augmented by introducing Galerkin least-squares type terms arising from the constitutive and equilibrium equations of the Navier-Stokes equations and from the Dirichlet boundary condition, which are multiplied by stabilization parameters that are chosen in such a way that the resulting continuous formulation becomes well-posed. Then, the classical Banach’s fixed point Theorem and Lax-Milgram’s Lemma are applied to prove well-posedness of the continuous problem. Similarly, we establish well-posedness and the corresponding Cea’s estimate of the associated Galerkin scheme considering any conforming finite element subspace for each unknown. In particular, the associated Galerkin scheme can be defined by employing Raviart-Thomas elements of degree kfor the nonlinear-pseudostress tensor, and continuous piecewise polynomial elements of degree k+ 1 for the velocity, which leads to an optimal convergent scheme. In addition, we provide two iterative methods to solve the corresponding nonlinear system of equations and analyze their convergence. Finally, several numerical results illustrating the good performance of the method are provided. Key words: Navier-Stokes, mixed finite element method, augmented formulation, Raviart- Thomas elements Mathematics Subject Classifications (1991): 65N15, 65N30, 76D05, 76M10 ∗This research was partially supported by CONICYT-Chile through project Inserci´on de Capital Humano Avanzado en la Academia 79130048; project Fondecyt 11121347, project Anillo ACT1118 (ANANUM); by DIUBB project 120808 GI/EF; and by Ministry of Education, Youth and Sports of the Czech Republic through the ERCCZ project LL1202. †Departamento de Matem´atica y F´ısica Aplicadas, Universidad Cat´olica de la Sant´ısima Concepci´on, Casilla 297, Concepci´on, Chile, and CI2MA, Universidad de Concepci´on, Casilla 160-C, Concepci´on, Chile, email: [email protected] ‡Departamento de Matem´atica, Universidad del B´ıo-B´ıo, Casilla 5-C, Concepci´on, Chile, and CI2MA, Universidad de Concepci´on, Casilla 160-C, Concepci´on, Chile, email: [email protected] §Mathematical Institute, Faculty of Mathematics and Physics, Charles University, Prague 8, 186 75, Czech Republic, email: [email protected] 1 1 Introduction This paper is concerned with the numerical approximation of the incompressible Navier-Stokes problem: −ν∆u+ (u·∇)u+∇p=fin Ω, div u=0in Ω, u=0on Γ, ZΩ p= 0, (1.1) where Ω ⊂Rn,n∈ {2,3}, is a polyhedral domain, ν > 0 is the fluid viscosity, uis the velocity, pis the pressure, and fis a source term. In particular, we are interested in pseudostress-type mixed variational formulations to solve problem (1.1). Pseudostress-based formulations, applied to the Stokes problem, have been extensively studied during the last decades (see e.g. [5], [14], [20], [22], and the references therein). This study has been extended to solve important problems in engineering, such as the Stokes-Darcy coupled problem and transport problems (see, for instance, [1], [21], [24], [25]). There are two main advantages to utilize this kind of formulations. On the one hand, pseudostress-based formulations have a natural applicability to non-Newtonian flows. Indeed, since in this case the constitutive equation is nonlinear, the stress can not be eliminated, and hence it becomes an unavoidable unknown in the corresponding solvability analysis. Actually, the main advantage of this kind of formulations is that their allows for a unified analysis for linear and nonlinear flows (see [11], [19], [23], [28]). On the other hand, when using pseudostressbased formulations, further variables of interest can be computed as a simple postprocess of the pseudostress, and then, they can be approximated without loosing any accuracy. In fact, after computing the pseudostress tensor σ:= ν∇u−pI one can easily obtain the fluid stress tensor ˜ σ:= ν(∇u+∇ut)−pI, the fluid pressure p, the fluid vorticity ω:= 1 2(∇u−(∇u)t), and the fluid velocity gradient ∇u, respectively, as follows: ˜ σ=σ−1 ntr (σ)I + σt, p =−1 ntr (σ),ω=1 2ν(σ−σt) and ∇u=1 ν(σ−1 ntr (σ)I). Now, concerning mixed formulations for the Navier-Stokes equations, we first mention the works of Farhloul et al. (see [12] and [13]), where the authors extend the analysis of dual-mixed formulations for the Stokes equations to the Navier-Stokes problem. They propose quasi-optimal convergent numerical methods for the fluid flow problem considering the strain tensor (in [12]) and the velocity gradient tensor (in [13]) as the main unknowns of the corresponding systems. In [8] (see also [6] and [7]), Cai et al. have extended the analysis of pseudostress-based mixed methods for the Stokes problem to the Navier-Stokes equations. They introduce and analyze a conforming H(div ) method for a pseudostress-based mixed formulation which turns to be of accuracy O(hk+1−n/6) (n= 2,3) in the L3norm. Finally, in [29], Howell and Walkington, present a dual-mixed finite element method for the Navier-Stokes problem, where the main unknowns are the velocity-gradient (in L2), the velocity (in L2) and a modified version of the pseudostress (or stress) tensor (in H(div )) that links the pseudostress tensor and the convective term. The authors prove well-posedness of the continuous formulation utilizing the classical Babuˇska-Brezzi theory, the fact that the velocity can be seen as a function in H1(Lemma 2.4), and a fixed-point argument. The first tool involves inf-sup conditions which are not easy to prove for the classical 2 H(div )-conforming discrete spaces. In order to deal with this issue, the authors introduce new finite element subspaces for the corresponding unknowns, which turn to be expensive. In this work we introduce a new augmented mixed finite element method for problem (1.1). Similarly to [29], we introduce a “nonlinear-pseudostress” tensor (in H(div )), linking the pseudostress tensor and the convective term, which together to the velocity (in H1), are the main unknowns of the system. The pressure is eliminated by using the incompressibility property, and can be recovered as a simple postprocess of the nonlinear-pseudostress tensor. Since the velocity is kept as a function in H1, we utilize similar techniques developed in [16] and [17] for the linear elasticity problem (see also [14] for similar results applied to the Stokes problem), and augment the mixed formulation with Galerkin least-squares type terms arising from the constitutive and equilibrium equations of the Navier-Stokes equations and from the Dirichlet boundary condition, which are multiplied by stabilization parameters that are chosen in such a way that the resulting continuous formulation becomes well-posed. The introduction of these terms allows us to circumvent the necessity of proving inf-sup conditions, and as a result, to relax the hypotheses on the corresponding discrete subspaces. In this way, the classical Banach’s fixed point Theorem and Lax-Milgram’s Lemma can be applied to prove well-posedness of the continuous problem. Analogously, and considering any pair of conforming subspaces for the nonlinear-pseudostress tensor and the fluid velocity, we establish well-posedness and the corresponding Cea’s estimate of the associated Galerkin scheme. In particular, making use of Raviart-Thomas elements of degree kfor the nonlinear-pseudostress tensor, and continuous piecewise polynomial elements of degree k+1 for the velocity, we obtain optimal rate of convergence, which is one of the principal advantages of our method. The rest of this paper is organized as follows. In Section 2 we present the main aspects of the continuous problem. We reformulate the fluid flow problem as an equivalent first-order set of equations. Then, we multiply the corresponding equations by suitable test functions and integrate by parts to derive our augmented mixed variational formulation. In Section 3 we choose appropriate stabilization parameters, and apply the classical Banach’s fixed point Theorem and Lax-Milgram’s Lemma to prove unique solvability and stability of the continuous formulation. In Section 4 we introduce and analyze the associated Galerkin scheme. We proceed analogously to the continuous case, and establish well-posedness of the corresponding discrete problem, considering any subspace for the nonlinear-pseudostress tensor and the velocity. In addition, we establish the corresponding Cea’s estimate and provide a suitable choice of finite element spaces ensuring optimal rate of convergence. In Section 5 we present two iterative strategies to compute the solutions of the aforementioned augmented scheme, and analyze their convergence. Finally, several numerical results illustrating the good performance of the method, are reported in Section 6. 2 Continuous problem 2.1 Preliminaries Let us first introduce some general functional spaces. If Ois a domain and r∈R, we define Hr(O) := [Hr(O)]n,Hr(O) := [Hr(O)]n×n. In the particular case r= 0 we usually write L2(O), and L2(O), instead of H0(O), and H0(O), respectively. The corresponding norms are denoted by k·kr,O(for Hr(O), Hr(O), and Hr(O)). 3 Also, the Hilbert space H(div ; O) := w∈L2(O) : div w∈L2(O), is standard in the realm of mixed problems (see [4], [18] [15] for instance), where div is the usual divergence operator acting on vector fields. The space of matrix valued functions whose rows belong to H(div ; O) will be denoted H(div ;O). The Hilbert norms of H(div ; O) and H(div ;O) are denoted by k·kdiv ;Oand k·kdiv ;O, respectively. Note that if S∈H(div ;O), then div S ∈L2(O). Note also that H(div ;O) can be characterized as the space of matrix valued functions Ssuch that ctS∈H(div ; O) for any constant column vector c. In addition, we have the following decomposition: H(div ;O) = H0(div ;O)⊕P0(O) I ,(2.1) where H0(div ;O) := S∈H(div ;O) : ZO tr S= 0 ,(2.2) I is the n×nidentity matrix, P0(O) is the space of constant polynomials on O, and tr stands for the usual trace of tensor. More precisely, each S∈H(div ;O) can be decomposed uniquely as: S=S0+cI,with S0∈H0(div ;O) and c:= 1 n|O| ZO tr S∈R.(2.3) Finally, for the sake of simplicity, in what follows we denote (u, v)Ω:= ZΩ uv, (u,v)Ω:= ZΩ u·v,(u,v)Γ:= ZΓ u·vand (T,S)Ω:= ZΩ T:S. 2.2 The first-order set of equations As we already mentioned, in order to describe our mixed variational formulation, we first rewrite problem (1.1) as an equivalent first-order set of equations. We begin by introducing the “nonlinear-pseudostress” tensor T:= ν∇u−pI−u⊗uin Ω.(2.4) From the incompressibility condition div u= tr (∇u) = 0 in Ω, it is not difficult to see that div (u⊗u)=(u·∇)uin Ω and tr (T) = −np −tr (u⊗u) in Ω.(2.5) In particular, the second equation in (2.5) allows us to write the pressure pin terms of the tensor Tand the velocity uas p=−1 n(tr (T) + tr (u⊗u)) in Ω,(2.6) which in turns, together with (2.4), leads us to the equation Td=ν∇u−(u⊗u)din Ω,(2.7) where Sd:= S−1 n(tr S) I, 4 is the deviatoric part of tensor S. On the other hand, from (2.4), and the first equations of (1.1) and (2.5), we easily get the equilibrium equation −div T =fin Ω.(2.8) Finally, from (2.6) we observe that the condition RΩp= 0, ensuring the uniqueness of solution of problem (1.1), is equivalent to (tr (T) + tr (u⊗u),1)Ω= 0.(2.9) According to the above, we rewrite equations (1.1) equivalently as follows Td=ν∇u−(u⊗u)din Ω,−div T =fin Ω, u=0on Γ,(tr (T),1)Ω=−(tr (u⊗u),1)Ω. (2.10) where the unknowns of the system are the tensor Tand the velocity u. The pressure pcan be easily computed as a postprocess of the solution by using (2.6). Here, ∇u=∂ui ∂xj1≤i,j≤n ,and div T =div (ti1, . . . , tin)1≤i≤n. Remark 2.1 Notice that, utilizing the solution of (2.10), besides recovering the pressure p, we can also compute other variables of interest, such as the stress tensor ˜ σ=ν(∇u+ (∇u)t)−pI, the velocity gradient ∇uand the vorticity ω:= 1 2(∇u−(∇u)t), respectively, as follows ˜ σ=Td+ (u⊗u)d+Tt+u⊗u, ∇u=1 ν(Td+ (u⊗u)d), ω=1 2ν(T−Tt). (2.11) These formulae will be employed later to approximate this variables as a simple postprocess of the finite element solutions, conserving the same rate of convergence. 2.3 The augmented mixed variational problem We now turn to the derivation of the weak formulation of (2.10). To do this we first test the first equation of (2.10) with arbitrary S∈H(div ; Ω), integrate by parts, utilize the Dirichlet boundary condition u=0on Γ, and the identity Td:S=Td:Sd, to obtain (Td,Sd)Ω+ν(div S,u)Ω+ (u⊗u,Sd)Ω= 0.(2.12) In addition, the equilibrium equation div T =−fis imposed weakly as follows ν(div T,v)Ω=−ν(f,v)Ω,(2.13) where vis a suitable test function. Observe that the term (u⊗u,Sd)Ωin (2.12) requires the velocity uto be in a space smaller than L2(Ω). Then, as we will see later, in order to obtain an optimal convergent discrete scheme, we propose to consider the velocity uin H1(Ω). 5 Consequently, in order to ensure well-posedness of the resulting variational formulation, we propose to enrich our formulation with the following residual terms arising from the constitutive, equilibrium equations and the Dirichlet boundary condition κ1(div T,div S)Ω=−κ1(f,div S)Ω,(2.14) κ2ν∇u−Td−(u⊗u)d,∇vΩ= 0,(2.15) κ3(u,v)Γ= 0,(2.16) for all S∈H(div ; Ω) and v∈H1(Ω), where κ1,κ2and κ3are positive parameters to be specified later. Then, from (2.12)–(2.16) we obtain the variational problem: Find (T,u)∈H(div ; Ω) × H1(Ω), such that (tr (T) + tr (u⊗u),1)Ω= 0, and (Td,Sd)Ω+κ1(div T,div S)Ω+ν(div S,u)Ω+ (u⊗u,Sd)Ω=−κ1(f,div S)Ω, −ν(div T,v)Ω+κ2ν∇u−Td−(u⊗u)d,∇vΩ+κ3(u,v)Γ=ν(f,v)Ω, (2.17) for all (S,v)∈H(div ; Ω)× ∈ H1(Ω). Let us now define the tensor T0:= T+1 n|Ω|(tr (u⊗u),1)ΩI.(2.18) It is clear that T0∈H0(div ; Ω) if and only if (tr (T) + tr (u⊗u),1)Ω= 0. In this way, thanks to (2.18) and (2.3), problem (2.17) can be reformulated equivalently as: Find (T0,u)∈H0(div ; Ω) ×H1(Ω), such that (Td 0,Sd)Ω+κ1(div T0,div S)Ω+ν(div S,u)Ω+ (u⊗u,Sd)Ω=−κ1(f,div S)Ω, −ν(div T0,v)Ω+κ2ν∇u−Td 0−(u⊗u)d,∇vΩ+κ3(u,v)Γ=ν(f,v)Ω, (2.19) for all (S,v)∈H0(div ; Ω)× ∈ H1(Ω). The following lemma establishes that problems (2.17) and (2.19) are equivalent. Its proof is straightforward. Lemma 2.2 If (T,u)is a solution of (2.17), then (T0,u) = (T+1 n|Ω|(tr (u⊗u),1)ΩI,u)is a solution of (2.19). Conversely, if (T0,u)is a solution of (2.19), then (T,u)=(T0−1 n|Ω|(tr (u⊗ u),1)ΩI,u)is a solution of (2.17). As a consequence of the above, in what follows we focus on analyzing problem (2.19). 3 Analysis of the continuous problem In this section we analyze the well-posedness of problem (2.19). To that end, we first rewrite our problem in a variational setting and state the main properties of the forms involved. 6 3.1 Variational setting First, let us define the global unknown and space: φ:= (T0,u)∈X:= H0(div ; Ω) ×H1(Ω), where Xis endowed with the norm kψk2 X:= kSk2 div ,Ω+kvk2 1,Ω∀ψ= (S,v)∈X. Then, defining the forms A(φ,ψ) := (Td,Sd)Ω+κ1(div T,div S)Ω+κ2ν(∇u,∇v)Ω+κ3(u,v)Γ −ν(div T,v)Ω+ν(div S,u)Ω−κ2(Td,∇v)Ω, C(ζ;φ,ψ) := (u⊗z,Sd)Ω−κ2((u⊗z)d,∇v)Ω, (3.1) for all φ= (T,u),ψ= (S,v),ζ= (R,z) in X, and the functional F(ψ) := −κ1(f,div S)Ω+ν(f,v)Ω∀ψ= (S,v)∈X,(3.2) the variational problem (2.19) reads: Find φ∈X, such that A(φ,ψ) + C(φ;φ,ψ) = F(ψ)∀ψ= (S,v)∈X.(3.3) We now turn to discuss the stability properties of the forms involved. To do this we first establish the following well known inequalities. For (3.4) we refer to Lemma 3.1 in [2] or Chapter IV in [4], whereas for inequality (3.5) we refer to [9] and Chapter I in [15]. CdkSk2 0,Ω≤ kSdk2 0,Ω+kdiv Sk2 0,Ω∀S∈H0(div ; Ω),(3.4) Cpkvk2 1,Ω≤ |v|2 1,Ω+kvk2 0,Γ≤˜ Cpkvk2 1,Ω∀v∈H1(Ω),(3.5) with Cd,Cpand ˜ Cp, only depending on Ω. Notice that the inequality on left hand side of (3.5) is nothing but the well known generalized Poincar´e inequality. Now, for the bilinear form A, we observe that, thanks to (3.5) and the H¨older’s inequality, for any positive constants κ1,κ2, and κ3, there holds |A(φ,ψ)| ≤ CAkφkXkψkX,∀φ,ψ∈X,(3.6) with CA>0 depending only on ν,κ1,κ2,κ3and ˜ Cp. Now, for the form C, recalling that H1is continuously embedded in L4, from H¨older’s inequality, we obtain |C(ζ;φ,ψ)| ≤ (1 + κ2 2)1/2kzkL4(Ω)kukL4(Ω)kψkX ≤(1 + κ2 2)1/2C2 sobkzk1,Ωkuk1,ΩkψkX ≤CCkζkXkφkXkψkX, (3.7) for all ζ= (R,z),φ= (T,u), and ψ∈X, with CC:= C2 sob(1 + κ2 2)1/2. Here, Csob >0 is the constant of the corresponding Sobolev inequality. We now turn to prove the ellipticity of Aon X. This result is established next. 7 Lemma 3.1 Assume that κ1, κ3>0and 0< κ2<2ν. Then there exists αA>0, such that A(ψ,ψ)≥αAkψk2 X∀ψ∈X.(3.8) Proof. Given ψ:= (S,v)∈X, from the definition of A, inequalities (3.4), (3.5), and H¨older’s inequality, we have A(ψ,ψ) = kSdk2 0,Ω+κ1kdiv Sk2 0,Ω+κ2ν|v|2 1,Ω+κ3kvk2 0,Γ−κ2(Sd,∇v)Ω ≥ kSdk2 0,Ω+κ1kdiv Sk2 0,Ω+κ2ν|v|2 1,Ω+κ3kvk2 0,Γ−1 2kSdk2 0,Ω−κ2 2 2|v|2 1,Ω ≥1 2kSdk2 0,Ω+κ1kdiv Sk2 0,Ω+κ2 2(2ν−κ2)|v|2 1,Ω+κ3kvk2 0,Γ ≥1 2min{1, κ1}CdkSk2 0,Ω+κ1 2kdiv Sk2 0,Ω+1 2min {κ2(2ν−κ2),2κ3}Cpkvk2 1,Ω. Then, the result follows with αA=1 2min min{1, κ1}Cd, κ1,min {κ2(2ν−κ2),2κ3}Cp. Remark 3.2 Observe that if κ2∈(0,2ν),κ1=κ2(2ν−κ2)and κ3=1 2κ2(2ν−κ2), then for any viscosity constant ν∈(0,1], the ellipticity constant αAbecomes αA=1 2κ2(2ν−κ2) min{Cd,1, Cp}, which reaches its maximum value when κ2=ν. We end this section by establishing the boundedness of the functional F: |F(ψ)| ≤ kfk0,Ω(κ1kdiv Sk0,Ω+νkvk0,Ω) ≤ kfk0,Ω(κ2 1+ν2)1/2kψkX, (3.9) for all ψ= (S,v)∈X. 3.2 The main result Now, we are in position of providing the main result of this section, namely, existence, uniqueness and stability of problem (3.3). To that end we first introduce the bounded set K:= (ψ∈X:kψkX≤2(κ2 1+ν2)1/2 αAkfk0,Ω),(3.10) and the mapping J:K→K,ζ→ J(ζ) = φ,(3.11) with ζand φsatisfying A(φ,ψ) + C(ζ;φ,ψ) = F(ψ)∀ψ∈X.(3.12) Let us observe that φ∈Xis a solution to (3.3), if and only if, J(φ) = φ. In this way, to prove well-posedness of problem (3.3), in what follows we equivalently prove that Jhas a unique fixed point by means of the classical Banach’s fixed point theorem. Before doing this, we now provide the following lemma establishing the well-definedness of operator J. 8 Lemma 3.3 Assume that 4CC(κ2 1+ν2)1/2 α2 Akfk0,Ω≤1,(3.13) with CCand αAbeing the constants in (3.7) and (3.8), respectively. Assume further that κ1, κ3>0and 0< κ2<2ν. Then, given ζ∈K, there exists a unique φ∈Ksuch that J(ζ) = φ. Proof. Let ζ∈K. From (3.7) and (3.8), there holds A(ψ,ψ) + C(ζ;ψ,ψ)≥A(ψ,ψ)− |C(ζ;ψ,ψ)| ∀ψ∈X, ≥αA−CCkζkXkψk2 X∀ψ∈X.(3.14) In turn, from the definition of set K, and assumption (3.13), we easily get CCkζkX≤2CC(κ2 1+ν2)1/2 αAkfk0,Ω≤αA 2.(3.15) Hence, combining (3.14) and (3.15), we obtain A(ψ,ψ) + C(ζ;ψ,ψ)≥αA 2kψk2 X,ψ∈X,(3.16) that is, the bilinear form A(·,·) + C(ζ;·,·) is elliptic on X. Therefore, applying Lax-Milgram Lemma we obtain that there exists a unique φ∈Xsatisfying (3.12). In addition, from (3.9), (3.12) and (3.16), we get αA 2kφk2 X≤A(φ,φ) + C(ζ;φ,φ) = F(φ)≤ kfk0,Ω(κ2 1+ν2)1/2kφkX, which together to (3.13) implies that φbelongs to Kand concludes the proof.  We now provide the well-posedness of problem (3.3) (or equivalently (2.19)), which corresponds to the main result of this section. Theorem 3.4 Let f∈L2(Ω) such that 4CC(κ2 1+ν2)1/2 α2 Akfk0,Ω<1.(3.17) Assume that κ1, κ3>0and 0< κ2<2ν. Then, there exists a unique φ∈Xsolution to (3.3). In addition, the solution φsatisfies the continuous dependence result kφkX≤2(κ2 1+ν2)1/2 αAkfk0,Ω.(3.18) Proof. First, let us observe that assumption (3.17) ensures that Lemma 3.3 holds true (see assumption (3.13)). Now, as mentioned before, we make use of the classical Banach’s fixed point Theorem to prove that the mapping Jhas a unique fixed point in K. In fact, let ζ1,ζ2,φ1,φ2∈K, such that φ1=J(ζ1) and φ2=J(ζ2). From the definition of Jit follows that A(φ1−φ2,ψ) + C(ζ1;φ1,ψ)−C(ζ2;φ2,ψ) = 0,∀ψ∈X. 9 Proposition 5.1 Let f∈L2(Ω) such that (3.17) holds. Assume that κ1, κ3>0and 0< κ2< 2ν. Then, given an initial guess φ0, satisfying kφ0kX≤2(κ2 1+ν2)1/2 αAkfk0,Ω, the sequence φngiven by (5.1) converges to the unique solution φof problem (3.3) when ntends to infinity, and kφkX≤2(κ2 1+ν2)1/2 αAkfk0,Ω. 5.2 Newton’s method Now we present the well known Newton’s method to solve the nonlinear problem (3.3). Before describing this strategy applied to our problem, we first recall this method in an abstract framework and introduce a classical theorem ensuring its convergence, namely, the Kantorovich Theorem. To do this, we let Xand Ybe two Banach spaces and x0a given element in X. In addition, we let BRand Brbe the bounded sets defined by BR:= {x∈X:kx−x0kX< R} and Br:= {x∈X:kx−x0kX≤r}, where Rand rare positive constants satisfying 0 < r < R. Let us assume that BRcontains a zero of an operator P:BR⊂X→Y, i.e., a point x∗∈BRsuch that P(x∗) = 0, and that P has a continuous derivative in BR. Then the Newton’s method is described as follows: Newton’s method: Given an initial approximation x0∈BRof x∗, compute: xm=xm−1−[P0(xm−1)]−1(P(xm−1)) for m≥1, assuming that [P0(xm−1)]−1exists. The following theorem establishes the convergence of this method (for details, see Theorem 5 in [3] or Theorem 6 (1.XVIII) in [30, p. 708]). Theorem 5.2 (Kantorovich Theorem) Let Pbe defined on BR⊂Xwith continuous second derivative in Br. Moreover assume that (1) there exists the continuous linear operator [P0(x0)]−1, (2) there exists a positive constant K1>0:k[P0(x0)]−1(P(x0))k ≤ K1, (3) there exists a positive constant K2>0:k[P0(x0)]−1P00(x)k ≤ K2, for all x∈Br. If K3=K1K2≤1 2,(5.2) and the radious rof Brsatisfies r≥r0=1−√1−2K3 K3 K1, 16 then, there exists a zero x∗of Pto which the Newton’s iteration converges. In this case, kx∗−x0kX≤r0. Furthermore, if for K3<1/2 r < r1=1 + √1−2K3 K3 K1, or for K3= 1/2 r≤r1, the solution x∗is unique in the set Br. The rate of convergence of the Newton’s method is given by kx∗−xmkX≤1 2m(2K3)2mK1 K3 , m ≥0.(5.3) In what follows, we apply the theory described above to derive and analyze the Newton’s method associated to problem (3.3). To do so, we let ζbe a fixed element in X, and introduce the linear operator N(ζ,·) : X→X0, defined by N(ζ,φ),ψ:= A(φ,ψ) + C(ζ;φ,ψ)∀φ,ψ∈X,(5.4) and the functional F ∈ X0 F,ψ:= F(ψ)∀ψ∈X,(5.5) Then, defining the nonlinear operator P:X→X0, as P(ψ) := N(ψ,ψ)− F ∀ψ∈X,(5.6) we obtain that φ∈Xis the unique solution of (3.3), if and only if, P(φ),ψ=A(φ,ψ) + C(φ;φ,ψ)−F(ψ)=0 ∀ψ∈X, that is, P(φ) = 0∈X0. Now, after simple computations, it is not difficult to see that the Gateaux derivative of P, P0:X→ L(X,X0), is given by (P0(φ))(ζ),ψ=A(ζ,ψ) + C(ζ;φ,ψ) + C(φ;ζ,ψ)∀φ,ζ,ψ∈X,(5.7) where L(X,X0) is the set of linear and bounded operators from Xto X0. According to the above, we obtain that the Newton’s method applied to problem (3.3) reads: Given φ0∈X, for m≥1, find φm∈X, such that (P0(φm−1))(φm−φm−1),ψ=−(P(φm−1)),ψ∀ψ∈X.(5.8) In what follows we apply Theorem 5.2 to prove the convergence of the iterative method (5.8), considering φ0=0as initial guess. This result is stated in the next theorem. 17 Theorem 5.3 Assume that hypotheses of Theorem 3.4 hold and let φ∈Xthe unique solution of problem (3.3). Then, considering φ0=0as initial guess, the sequence {φm}m∈Ngiven by the Newton’s method (5.8) converges to φ. In addition, There exist positive constants K1and K2, independent of the solution and m, such that kφ−φmkX≤1 2m(2K2)2mK1 K2 , m ≥0.(5.9) (Explicit constants K1and K2are provided in (5.12) and (5.14), respectively). Proof. It suffices to verify that assumptions (1),(2) and (3) of Theorem 5.2 hold. First, it is easy to see that P(φ0),ψ=P(0),ψ=−F(ψ)≤ kfk0,Ω(κ2 1+ν2)1/2kψkX∀ψ∈X, from which kP(φ0)kX0≤ kfk0,Ω(κ2 1+ν2)1/2.(5.10) In turn, from Lemma 3.1, we have (P0(φ0))(ψ),ψ=A(ψ,ψ)≥αAkψk2 X∀ψ∈X, which implies that [P0(φ0)]−1exists, and k[P0(φ0)]−1kL(X0,X)≤1 αA .(5.11) In this way, from (5.10) and (5.11) we obtain that assumption (1) and (2) hold with K1=1 αA (κ2 1+ν2)1/2kfk0,Ω.(5.12) Now, for (3) we first notice, after simple computations, that the second derivative of operator Pis given by ((P00(φ)(ϕ))(ψ),ζ=C(ϕ;ψ,ζ) + C(ψ;ϕ,ζ)∀φ,ϕ,ψ,ζ∈X. It follows that ((P00(φ)(ϕ))(ψ),ζ≤2CCkϕkXkψkXkζkX∀φ,ϕ,ψ,ζ∈X, which yields kP00(φ)kL(X,L(X,X0)) ≤2CC.(5.13) Hence, from (5.11) and (5.13) we get that assumption (3) holds with K2=2CC αA .(5.14) Therefore, from (3.17), (5.12) and (5.14) we finally get K3:= K1K2=2CC(κ2 1+ν2)1/2 α2 Akfk0,Ω<1 2, which concludes the proof.  18 Remark 5.4 We observe that, according to the definition of operator Pin (5.6), the Newton’s method described in (5.8) can be rewritten in terms of the forms Aand Cas follows: Find φm∈Xsuch that A(φm,ψ) + C(φm;φm−1,ψ) + C(φm−1;φm,ψ) = C(φm−1;φm−1,ψ) + F(ψ),(5.15) for all ψ∈X. 6 Numerical results In this section we present three examples illustrating the performance of our augmented mixed finite element scheme (4.4) on a set of uniform triangulations of the corresponding domains and considering the finite element spaces introduced in Section 4.5. Our implementation is based on a FreeFem++ code (see [26]), in conjunction with the direct linear solver UMFPACK (see [10]). Regarding the implementation of the iterative methods, the iterations are terminated once the relative error of the entire coefficient vectors between two consecutive iterates is sufficiently small, i.e., kcoeffm+1 −coeffmkl2 kcoeffm+1kl2≤tol, where k · kl2is the standard l2-norm in RN, with Ndenoting the total number of degrees of freedom defining the finite element subspaces Hhand Hh, and tol is a fixed tolerance to be specified on each example. As suggested in Theorem 5.3, for each example shown below we simply take φ0:= (T0 h,0,u0 h) = (0,0) as initial guess. We now introduce some additional notations. The individual errors are denoted by: e(u) := ku−uhk1,Ω,e(T0) := kT0−Th,0kdiv ,Ω,e(p) := kp−phk0,Ω, e(ω) := kω−ωhk0,Ω,e(∇u) := k∇u−Ghk0,Ω,e(˜ σ) := k˜ σ−˜ σhk0,Ω. where ph,ωh,Gh, and ˜ σhare the variables computed by the prostprocess (4.18). Also, we let r(T0), r(u), r(p), r(ω), r(∇u), and r(˜ σ) be the experimental rates of convergence given by r(T0) := log(e(T0)/e0(T0)) log(h/h0), r(u) := log(e(u)/e0(u)) log(h/h0), r(p) := log(e(p)/e0(p)) log(h/h0), r(ω) := log(e(ω)/e0(ω)) log(h/h0), r(∇u) := log(e(∇u)/e0(∇u)) log(h/h0), r(˜ σ) := log(e(˜ σ)/e0(˜ σ)) log(h/h0), where hand h0denote two consecutive meshsizes with errors eand e0. Our first example is carried out in order to assess the influence of the parameters κ1,κ2and κ3on the performance of our augmented mixed formulation (4.4). For this example we consider the domain Ω := (−1,1)2, the viscosity constant ν= 1, and take fso that the exact solution is given by u(x1, x2) =   πsin(πx1)2cos(πx2) −2πcos(πx1) sin(πx1) sin(πx2) , p(x1, x2) = 5x1sin(x2). 19 -4 -2 0 2 4 -4.7 4.7 -4 -2 0 2 4 -4.7 4.7 Figure 6.1: Example 1: ph(left) and p(right) with N= 157339 In order to emphasize the robustness of (4.4) with respect to the parameters κ1,κ2and κ3, as suggested in Remark 3.2, we choose (κ1, κ2, κ3)=(ν2, ν, ν2 2) and (κ1, κ2, κ3)=(3ν2 4,3ν 2,3ν2 8), which, certainly, satisfy the assumptions of Lemma 3.1. In Tables 6.1 and 6.2 we summarize the convergence history for a sequence of uniform triangulations, considering the finite element spaces introduced in Section 4.5 with k= 0 and k= 1, and solving the nonlinear problem with the Newton’s method defined in (5.15). We observe there that the rate of convergence O(hk+1) predicted by Theorem 4.8 (when s=k+ 1) is attained in all the cases. In addition, the specific results displayed in Tables 6.1 and 6.2, showing almost no difference between (κ1, κ2, κ3) = (ν2, ν, ν2 2) and (κ1, κ2, κ3)=(3ν2 4,3ν 2,3ν2 8) illustrate the robustness of scheme (4.4), with respect to the choice of these parameters. Next, in Figure 6.1 we display the approximate pressure (to the left) obtained from the RT0−P1approximation with N= 157339 degrees of freedom, and we compare it with the exact one (to the right). In addition, in Figure 6.2 we display the magnitude and vector field of the approximate velocity (to the left) and we compare it with the exact velocity (to the right). In both cases we observe that the finite element subspaces employed provide very accurate approximations to the unknowns. Our second example focusses on the performance of the iterative methods as a function of the viscosity ν, considering the analytical solution (u, p) obtained by Kovasznay in [31]. For the domain Ω := (−1/2,3/2) ×(0,2) and for a given ν, this solution is given by u(x1, x2) =   1−eλx1cos(2πx2) λ 2πeλx1sin(2πx2) , p(x1, x2) = −1 2e2λx1+ ¯p, 20 Errors and rates of convergence for the mixed RT0−P1approximation N h e(T0)r(T0)e(u)r(u)e(p)r(p) 691 0.4129 90.2251 – 14.0226 – 4.8682 – 2595 0.1901 47.6284 0.8235 7.0122 0.8933 2.7057 0.7571 9955 0.0968 23.9312 1.0192 3.6224 0.9781 1.4606 0.9130 39195 0.0527 12.0739 1.1269 1.8405 1.1154 0.7152 1.1762 157339 0.0307 5.9814 1.2960 0.9051 1.3095 0.3541 1.2970 N h e(ω)r(ω)e(∇u)r(∇u)e(˜ σ)r(˜ σ) Iterations 691 0.4129 6.2193 – 14.069 – 23.0156 – 5 2595 0.1901 2.6606 1.0944 6.6236 0.9710 11.5540 0.8883 4 9955 0.0968 1.2866 0.9580 3.2792 1.0411 5.8339 1.0120 4 39195 0.0527 0.5917 1.2796 1.5895 1.1929 2.8858 1.1595 4 157339 0.0307 0.2887 1.3238 0.7799 1.3138 1.4201 1.3083 4 Errors and rates of convergence for the mixed RT1−P2approximation Nhe(T0)r(T0)e(u)r(u)e(p)r(p) 2315 0.4129 23.6263 – 3.3517 – 1.9654 – 8883 0.1901 5.7915 1.8122 0.8857 1.7154 0.5413 1.6619 34451 0.0968 1.5250 1.9760 0.2188 2.0702 0.1455 1.9456 136407 0.0527 0.3865 2.2611 0.0566 2.2270 0.0397 2.1386 549143 0.0307 0.0939 2.6103 0.0139 2.5957 0.0098 2.5783 N h e(ω)r(ω)e(∇u)r(∇u)e(˜ σ)r(˜ σ) Iterations 231 0.4129 3.3323 – 7.0599 – 10.8747 – 4 8883 0.1901 0.9016 1.6849 1.9217 1.6772 2.9757 1.6704 4 34451 0.0968 0.2165 2.1123 0.4934 2.0135 0.8006 1.9442 4 136407 0.0527 0.0551 2.2520 0.1270 2.2330 0.2087 2.2144 4 549143 0.0307 0.0133 2.6231 0.0300 2.6060 0.0512 2.5946 4 Table 6.1: Example 1: Degrees of freedom, meshsizes, errors, rates of convergence and number of iterations for the mixed RT0−P1and RT1−P2approximations of the Navier-Stokes problem with (κ1, κ2, κ3)=(ν2, ν, ν2 2). 21 Errors and rates of convergence for the mixed RT0−P1approximation N h e(T0)r(T0)e(u)r(u)e(p)r(p) 691 0.4129 90.2374 – 13.9627 – 4.7906 – 2595 0.1901 47.6335 0.8235 7.0124 0.8877 2.7074 0.7356 9955 0.0968 23.9321 1.0193 3.6227 0.9781 1.4593 0.9152 39195 0.0527 12.0741 1.1270 1.8405 1.1155 0.7155 1.1741 157339 0.0307 5.9814 1.2960 0.9051 1.3095 0.3541 1.2978 N h e(ω)r(ω)e(∇u)r(∇u)e(˜ σ)r(˜ σ) Iterations 691 0.4129 6.2267 – 13.9904 – 22.7738 – 4 2595 0.1901 2.6626 1.0950 6.5992 0.9685 11.4952 0.8812 4 9955 0.0968 1.2869 1.0768 3.2748 1.0377 5.8230 1.0072 4 39195 0.0527 0.5917 1.2797 1.5890 1.1912 2.8848 1.1570 4 157339 0.0307 0.2887 1.3239 0.7798 1.3134 1.4199 1.3079 4 Errors and rates of convergence for the mixed RT1−P2approximation Nhe(T0)r(T0)e(u)r(u)e(p)r(p) 2315 0.4129 23.6269 – 3.3457 – 1.9558 – 8883 0.1901 5.7915 1.8122 0.8853 1.7137 0.5402 1.6584 34451 0.0968 1.5250 1.9760 0.2188 2.0697 0.1454 1.9435 136407 0.0527 0.3865 2.2611 0.0566 2.2269 0.0397 2.1377 549143 0.0307 0.0939 2.6103 0.0139 2.5957 0.0098 2.5780 N h e(ω)r(ω)e(∇u)r(∇u)e(˜ σ)r(˜ σ) Iterations 231 0.4129 3.3322 – 7.0533 – 10.8542 – 4 8883 0.1901 0.9016 1.6849 1.9210 1.6765 2.9733 1.6690 4 34451 0.0968 0.2165 2.1124 0.4933 2.0132 0.8003 1.9435 4 136407 0.0527 0.0551 2.2540 0.1272 2.2327 0.2087 2.2140 4 549143 0.0307 0.0133 2.6231 0.0310 2.6060 0.0512 2.5945 4 Table 6.2: Example 1: Degrees of freedom, meshsizes, errors, rates of convergence and number of iterations for the mixed RT0−P1and RT1−P2approximations of the Navier-Stokes problem with (κ1, κ2, κ3)=(3ν2 4,3ν 2,3ν2 8). 22 1 2 3 0 3.15 1 2 3 0 3.15 Figure 6.2: Example 1: velocitiy magnitudes |uh|(left) and |u|(right) and velocity vector fields with N= 157339 where λ:= −8π2 ν−1+√ν−2+ 16π2, and the constant ¯pis such that RΩp= 0. Since, in this case the solution is inhomogeneous on the boundary, we need to modify the functional Fas follows F(ψ) := −κ1(f,div S)Ω+hSn,giΓ+ν(f,v)Ω+κ3(g,v)Γ∀ψ= (S,v)∈X, where h·,·iΓdenotes the duality pairing between H−1/2(Γ) and H1/2(Γ), and g:= u|Γ∈H1/2(Γ). In Table 6.3 we show the behaviour of the iterative methods (5.1) and (5.15) as a function of the viscosity number, considering different meshsizes and a tolerance tol = 1E−10. The parameters κ1, κ2, κ3are simply chosen as (κ1, κ2, κ3) = (1, ν, 1) and we consider the finite element spaces introduced in Section 4.5 with k= 0. Here we observe that the smaller the parameter νthe higher the number of iterations. In addition we see that, as expected, the Newton’s method is faster than the fixed-point iteration. Numerical experiments for smaller values of νare not reported since both iterative methods need too many iterations to converge (more than 300). Next, the numerical results in Table 6.4 show the convergence history of the Newton’s method (5.15) considering the viscosity ν= 0.1. We see there that the rate of convergence O(h) provided by Theorem 4.8 and Corollary 4.9 for s= 1 is attained by the unknowns and all the post-processed variables. Finally, in our third example we assess the capability of a 3D implementation of our augmented scheme. Here, we consider the domain Ω := (0,1)3, the viscosity constant ν= 1, the 23 Number of iterations for the fixed-point method ν h = 0.3801 h= 0.1901 h= 0.0982 h= 0.0530 h= 0.0266 1 15 15 13 11 10 0.1 38 30 22 20 19 0.059 59 59 24 21 20 Number of iterations for the Newton’s method ν h = 0.3801 h= 0.1901 h= 0.0982 h= 0.0530 h= 0.0266 155544 0.1 7 6 6 5 5 0.059 12 7 6 6 6 Table 6.3: Example 2: Convergence behavior of the iterative methods with respect to the parameter ν. N h e(T0)r(T0)e(u)r(u)e(p)r(p) 683 0.3801 5.1122 – 12.9236 – 2.6472 – 2539 0.1901 3.0577 0.7415 6.7165 0.9442 1.2811 1.0471 9923 0.0982 1.3439 1.2448 2.8092 1.3198 0.5635 1.2435 39107 0.0530 0.5933 1.3248 1.1747 1.4129 0.2387 1.3919 157051 0.0266 0.2908 1.0342 0.5554 1.0862 0.1159 1.0476 N h e(ω)r(ω)e(∇u)r(∇u)e(˜ σ)r(˜ σ) 683 0.3801 21.9471 – 43.6479 – 8.9770 – 2539 0.1901 13.7239 0.6773 20.2037 1.1113 3.9804 1.1733 9923 0.0982 7.2897 0.9580 10.1070 1.0488 1.9127 1.1097 39107 0.0530 3.4478 1.2133 4.7831 1.2123 0.8896 1.2405 157051 0.0266 1.7457 0.9868 2.3881 1.0071 0.4405 1.0191 Table 6.4: Example 2: Degrees of freedom, meshsizes, errors, and rates of convergence for the mixed RT0−P1approximation of the Navier-Stokes problem with ν= 0.1, utilizing the Newton’s method. 24 N h e(T0)r(T0)e(u)r(u)e(p)r(p) 442 0.7071 0.2621 – 0.0605 – 0.1578 – 2968 0.3536 0.1400 0.9043 0.0286 1.0824 0.0794 0.9910 21772 0.1768 0.0666 1.0730 0.0114 1.3210 0.0346 1.1991 166804 0.0884 0.0313 1.0901 0.0045 1.3403 0.0151 1.2001 1305892 0.0442 0.0151 1.0479 0.0019 1.2158 0.0070 1.1092 N h e(ω)r(ω)e(∇u)r(∇u)e(˜ σ)r(˜ σ) Iterations 442 0.7071 0.0708 – 0.1259 – 0.3436 – 3 2968 0.3536 0.0413 0.7766 0.0795 0.6628 0.1933 0.8298 3 21772 0.1768 0.0225 0.8745 0.0449 0.8254 0.0980 0.9798 3 166804 0.0884 0.0117 0.9431 0.0236 0.9254 0.0486 1.0121 3 1305892 0.0442 0.0060 0.9748 0.0121 0.9697 0.0242 1.0059 3 Table 6.5: Example 3: Degrees of freedom, meshsizes, errors, rates of convergence and number of iterations for the mixed RT0−P1approximation of the 3D Navier-Stokes problem solved with the Newton’s method. parameters κ1=ν2,κ2=ν κ3=ν2/2, and take fso that the exact solution is given by u(x1, x2, x3) :=   4x2 1x2x3(x3−1)(x2−1)(x2−x3)(x1−1)2 −4x1x2 2x3(x2−1)2(x3−1)(x1−1)(x1−x3) 4x1x2x2 3(x3−1)2(x2−1)(x1−1)(x1−x2) , p(x1, x2, x3) := x1−1 2. For this example we solve the nonlinear problem with the Newton’s method defined in (5.15) with a tolerance tol = 1E−10 and consider the finite element spaces introduced in Section 4.5 with k= 0. In Table 6.5 we summarize the convergence history for a sequence of uniform triangulations. We observe there that the rate of convergence O(h) predicted by Theorem 4.8 and Corollary 4.9 (when s= 1) is attained in all the cases. Next, in Figure 6.3 we display the magnitude and vector field of the approximate velocity (top left and bottom left), and we compare them with their exact counterpart (top right and bottom right). Here, for the magnitude of the velocity we display the section of the cube below the plane x1+x2+x3= 1.5. In addition, in Figure 6.4 we display the second component of the approximate vorticity ω1,h (top left) and the iso-surface of the approximate pressure, and we compare them with their exact counterpart (top right and bottom right). Again, for the vorticity we display the section of the cube below the plane x1+x2+x3= 1.5. All the graphics above were computed with N= 1305892 degrees of freedom. We observe there that the finite element subspaces employed provide very accurate approximations to the unknowns. Acknowledgments. The authors are grateful to Prof. Gabriel N. Gatica and Eligio Colmenares (CI2MA, Universidad de Concepci´on) for their helpful comments that helped us to improve the paper. 25