A Posteriori Error Estimation for Incompressible Viscous Fluid With a New Boundary Condition

abstract: This paper describes numerical solutions of incompressible NavierStokes equations with a new boundary condition. To solve this problem, we use the discretization by mixed finite element method. We use a vector extrapolation method for computing numerical solutions of the steady-state Navier-Stokes equations. In addition, two types of a posteriori error indicator are introduced and are shown to give global error estimates that are equivalent to the true error. A numerical experiment on the driven cavity flow is given to demonstrate the effectiveness of the vector extrapolation method. We compare the result with the solution from commercial code like ADINA system as well as with values from other simulations.


Introduction
The idea of mixed finite element is to approximate simultaneously the piezometric head and the velocity.This approximation gives velocity throughout the field and the normal component of the velocity is continuous across the inter-element boundaries.Moreover, with the mixed formulation, the velocity is defined with the help of Raviart Thomas basis functions [1,2,3,45] and, therefore, a simple integration over the element gives the corresponding streamlines.This method was widely used for the prediction of the behavior of fluid in the hydrocarbons tank.
In computational mechanics one usually faces the problem of increasing the accuracy of a solution without adding unnecessary degrees of freedom.It is, therefore, necessary to update the mesh to ensure that it becomes fine enough in the critical regions while remaining reasonably coarse in the rest of the domain.Local a posteriopri error estimator is the adequate tool for identifying these critical regions automatically, using an input information only the given data and the numerical solution itself.A large amount of work has already been done in the construction of error estimators for several problems of computational fluid dynamics.In the conforming case there are several ways to define error estimators by using the residual equation.Ainsworth and Oden [5] and Verfurth [6] give a general overview.In the specific case of the Stokes and Navier-Stokes equations governing the steady flow of a viscous incompressible fluid, the work of Bank and Welfert [7], Verfurth [8] and Oden and Ainsworth [9] laid the basic foundation for the mathematical analysis of practical methods.Other works for the stationary Navier-Stokes problem have been introduced in [10,11,12,13,14,15,16].
Supported by . . .(Grant No. . . . ) In this paper we explore the potential of vector extrapolation in the context of computing steady state solutions of incompressible Navier-Stokes equations.Spatial discretization using mixed approximation of the velocity and pressure variables naturally lead us to consider highdimensional nonlinear algebraic systems-so effective iteration methods are crucial for efficient computation.Existing vector extrapolation methods can be broadly classified into two categories: polynomial methods and ε-algorithms.The first family includes three approaches: minimal polynomial extrapolation (MPE) of Cabay and Jackson [17]; reduced rank extrapolation (RRE) of Eddy [18] and Mesina [19], and the modified minimal polynomial extrapolation (MMPE) derived in [20,21,22].The second family includes the topological ε-algorithm (TEA) and the scalar and vector ε-algorithms (SEA and VEA).We will restrict our attention to the RRE methodology in this work.Some different recursive algorithms for implementing these methods have been presented in [23,31,33].We note that, when applied to linearly generated vector sequences, the MPE, the RRE and the TEA methods are related to Krylov subspace methods.
The plan of the paper is as follows.The model problem is described in section 2, followed by the discretization by mixed finite element method in section 3. A vector extrapolation methods is described in section 4 .Section 5 shows the methods of a posteriori error estimator of the computed solution and numerical experiment is described in section 6.

Governing equations
We consider the steady-state Navier-Stokes equations for a positive constant viscosity ν, where − → u is the fluid velocity, p is the pressure field, ∇ is the gradient, ∇. is the divergence operator, − → n denote the outward pointing normal to the boundary and α is a nonzero bounded continuous function defined on ∂Ω.This system is the basis for computational modeling of the flow of an incompressible Newtonian fluid such as air or water.The presence of the nonlinear convection term − → u .∇− → u means that boundary value problems associated with the Navier-Stokes equations can have more than one solution.
We set and Let the bilinear forms a : These inner products induce norms on V and W denoted by ∥.∥ V and ∥.∥ W respectively. ) Given the continuous functional l : Then the standard weak formulation of the Navier-Stokes flow problem (2.1) is the following: Let the subspace of divergence-free velocities be given by We obtain (2.12) The problem (2.9)-(2.10) is known [35] to possess a unique solution whenever the data is sufficiently small.In particular, if there exists a constant C such that for some fixed ω ∈ [0, 1).Then, there is a unique solution (2.14)

Mixed finite element discretization
Let P be a regular partitioning of the domain Ω into the union of N subdomains K such that • each K is a convex Lipschitzian domain with piecewise smooth boundary ∂K.
The common boundary between subdomains K and J is denoted by: Γ KJ = ∂K ∩ ∂J.
The finite element subspaces X h and M h are constructed in the usual manner so that the inclusion The finite element approximation to (2.9)-(2.10) is then Find for all ( − → v h , q h ) ∈ X h × M h .We define the appropriate bases for the finite element spaces, leading to a non linear system of algebraic equations.Linearization of this system using Newton iteration gives the finite dimensional system: for all − → v h ∈ X h and q h ∈ M h , where R k ( − → v h ) and r k (q h ) are the non linear residuals associated with the discrete formulations To define the corresponding linear algebra problem, we use a set of vector-valued basis functions { − → φ j }, so that and we fix the coefficients u j : j = n u + 1, . . ., n u + n ∂ , so that the second term interpolates the boundary data on ∂Ω D .
We introduce a set of pressure basis functions {Ψ k } and set where n u and n p are the numbers of velocity and pressure basis functions, respectively.We obtain a system of linear equations This system is referred to as the discrete Newton problem.The matrix A is the vector Laplacian matrix and B is the divergence matrix for i and j = 1, . . ., n u and k = 1, . . ., n p .The vector-convection matrix N and the Newton derivative matrix W are given by for i and j = 1, . . ., n u .The Newton derivative matrix is symmetric.
The right-hand side vectors in (3.7) are the non linear residuals associated with the discrete solution − → u h , p h and the function α , for i = 1, ..., n u and k = 1, . . ., n p .For Picard iteration, we give the discrete problem The lowest order mixed approximations like Q 1 − P 0 and Q 1 − Q 1 are unstable.We use a stabilized element pair Q 1 − P 0 , this is the most famous example of an unstable element pair, using bilinear approximation for velocity and a constant approximation for the pressure.

Vector extrapolation methods
We use a vector extrapolation method for computing numerical solutions of the steady-state Navier-Stokes equations.Extrapolation methods are of interest whenever an iteration process converges slowly.For a survey of these methods see for example the papers [31,32] and the book [33].The most popular vector extrapolation methods are the minimal polynomial extrapolation (MPE) [17], the reduced rank extrapolation (RRE) [18,19], and the modified minimal polynomial extrapolation (MMPE) [20,21,22].Convergence analysis of these methods can be found in [23,24].In particular, A. Sidi [25] shows that the MPE and the RRE approaches are mathematically equivalent to Arnoldis method [26] and to the generalized minimal residual method (GMRES) [27], respectively.Vector extrapolation methods are considered to be most effective when applied to nonlinear systems of equations [28,29,30,31].We use the reduced rank extrapolation method [18,19].
To define the RRE method, we will consider the solution of the linear or nonlinear system of equations: Let l be the solution of the system (4.1),l 0 be a given initial approximation, and generate the sequence of vectors l 1 , l 2 ,. .., according to the fixed-point iterative method where x − H(x) = 0 represents a preconditioned form of (4.1).Let (l n ) be a given sequence of N-dimensional column vectors formed by (4.2), and lim n→∞ l n = l , we set 3) The RRE method, when applied to the sequence l n , can be shown to generate an approximation ξ RRE n,k of the limit or the antilimit of (l n ).This approximation can be expressed in the form where and the scalars θ ij defined by the l 2 inner product: is given by the ratio of two determinants ) the matrices whose columns are ∆ i l n , ., ., .∆i l n+k−1 .
Using Schurs formula, we have where △ 2 L + n,k−1 is defined by exists and is unique if and only if det(△ 2 L T n,k−1 △ 2 L n,k−1 ) ̸ = 0. We uses the algorithms proposed in [23].

A posteriori error Analysis
In this section, we propose two types of a posteriori error indicator : the local Poisson problem estimator and the residual error estimator, which are shown to give global error estimates.Theorem 5. 1 We have sup ) and we have sup We gather (5.2) and ( 5. ) Theorem 5.2 Let (2.14) hold.Then there exist positive constants K 1 and K 2 such that Proof.See T.J. Oden, W. Wu, and M. Ainsworth [5].
The local velocity space on each subdomain K ∈ P is and the pressure space is (5.7) Let the bilinear forms a (5.9) Given the continuous functional l K : (5.10) (5.11) (5.12) (5.13) (5.14) The velocity space V(P) is defined by and the broken pressure space W (P ) is defined by (5.16) Examining the previous notations reveals that (5.17) We consider the space of continuous linear functional τ on V (P ) × W (P ) that vanish on the space V × W .Therefore, let H(div, Ω) denote the space equipped with norm (5.19)

Theorem 5.3 A continuous linear functional τ on the space V (P )) × W (P ) vanishes on the space V × W if and only if there exists
where − → n K denotes the unit outward normal vector on the boundary of K.
It will be useful to introduce the stresslike tensor σ( − → v , q) formally defined to be Where δ ij is the Kronecker symbol.
In order to define the value of the normal component of the stress on the interelement boundaries it is convenient to introduce notations for the jump on Γ KJ : An averaged normal stress on Γ KJ is defined by where α (i) KJ : Γ KJ −→ R are smooth polynomial functions.Naturally, the stress should be continuous then it is required that the averaged stress coincide with this value.On Γ KJ , we have The notation [[ .]] is used to define jumps in the elements of V(P) between subdomains.We define and (5.26) For − → v ∈ V (P ), we have for all ( − → w , q) ∈ V (P ) × W (P ).
We define the linear functional R : (5.29) for all ( − → w , q) ∈ V (P ) × W (P ).For ( − → w , q) ∈ V × W , we obtain (5.30) Let the lagrangian functional L : (5.31) So that and, for ( Using (5.34), we obtain: We have the problems on each subdomain (5.37) Suppose that the minimum exists, then the minimising element is characterized by finding The necessary and sufficient conditions for the existence of a minimum are that the data satisfy the following equilibration condition: When the subdomain K lies on the boundary ∂Ω the local problem (5.38) will be subject to a homogeneous Dirichlet condition on a portion of their boundaries and thus will be automatically well posed.However, elements away from the boundary are subject to pure Neumann conditions and the null space of the operator a ( ., . ) will contain the rigid motions where We construct data which satisfy the condition (5.39).We define Using (5.21), we obtain The averaged interelement stress may be rewritten where ⟨ − → n K .σ(− → v h , q h )⟩ 1 2 denotes the interelement averaged stress obtained using the symmetrical weighting corresponding to α = 1 2 .Then (5.47) For example, one might choose the piecewise bilinear pyramid functions associated with interior nodes in the partition.The relation (5.47) must hold at all points x contained in elements which do not interest the boundary of the domain.
The functions λ (k) KJ : Γ KJ −→ R are chosen to be of the form λ KJ,A X A (s), (5.48)where λ k KJ,A are constants to be determined.Owing the constraint (5.44), it is required that JK,A = 0, (5.49) for each A. Lemma 2. Suppose that for each X A the constants {λ and Proof.The result follows immediately by using (5.50), (5.48) and 5.44).
Summarizing and incorporating the results of section 5 we have Theorem 5.5 Let the conditions of Theorem 1 hold.Then there exists a constant C > 0 such that where (5.55) We define the global error estimator η by η = ( (5.56) We define the stress jump across edge or face E adjoining elements T and K where − → n E,K is the outward pointing normal.
We define the equidistributed stress jump operator and the interior residuals and (5.59) The element contribution η r,K of the residual error estimator is given by and the global residual error estimator η r is given by Theorem 5. 6 The estimator η r,K is equivalent to the η K estimator : there exist positive constants c 1 and C 2 such that (5.61) Proof.Same steps of the prof of Theorem 3.9 in [44].
Theorem 5.7 There exist positive constant C ′ such that (5.62)

Numerical examples
In this section some numerical results of calculations with finite element method and ADINA system will be presented.Using our solver, we run the test problem driven cavity flow [37, 38, 39, 40, 41, 42, 43, ].This is a classic test problem used in fluid dynamics, known as driven-cavity flow.It is a model of the flow in a square cavity with the lid moving from left to right.Let the computational model:        The solution shown in Figure 2 corresponds to a Reynolds number of 2000.The particles in the body of the fluid move in a circular trajectory.Steady flow in a two dimensional cavity is not stable for Reynolds number much greater than 10 4 .Indeed, we have made calculations for Reynolds number 10 4 .In addition, our code does not converge because the turbulence phenomena is not taken into account in our model.At a critical Reynolds number (approximately 13,000) the flow pattern develops into a time-periodic state with "waves" running around the cavity walls.The profiles of the u-velocity component along the vertical center line and the v-velocity component along the horizontal center line are shown in Figures 3 and 4 for Re=100 and Re=1000, respectively.In these figures, we have also included numerical predictions from [42] and ADINA system.There is an excellent agreement between the computed results, those published in [42] and the results computed with ADINA system.Figure 5 show the evolution of the nonlinear residual norm, using a logarithmic scale, for the simple Picard method, the Picard-Newton method and reduced rank extrapolation method (RRE), for two types of finite element discretization: a high order Q 2 − P −1 solution and a low order Q 1 − P 0 solution.We observes that the RRE method can be seen to accelerate the the fixed-point (Picard) iteration.The computational results of Figure 6 and Table 1 suggest that all two estimators seem to be able to correctly indicate the structure of the error.

Conclusion
We were interested in this work in the numeric solution for two dimensional partial differential equations modelling (or arising from) model steady incompressible fluid flow with a new boundary condition.It includes algorithms for discretization by finite element methods and a posteriori error estimation of the computed solutions.We use a reduced rank extrapolation method for computing numerical solutions of the steady-state Navier-Stokes equations.RRE is an effective techniques that have been used in accelerating the convergence of vector sequences.Two types of a posteriori error indicator are introduced and are shown to give global error estimates that are equivalent to the true discretization error.The computational results suggest that all two estimators seem to be able to correctly indicate the structure of the error.

3 )
to get (5.1).Let ( − → e , E) ∈ V × W be the error in the finite element approximation, − → e = − → u − − → u h and E = p − p h and define ( − → ϕ , ψ) ∈ V × W to be the Ritz projection of the modified residuals a leaky cavity.The streamlines are computed from the velocity solution by solving the Poisson equation numerically subject to a zero Dirichlet boundary condition.

Fig. 3 .
Fig. 3. Velocity component u at vertical center line (left plot), and the velocity component v at horizontal center line (right plot) with a 129 × 129 grid and Re=100.

Fig. 4 .
Fig. 4. Velocity component u at vertical center line (left plot), and the velocity component v at horizontal center line (right plot) with a 129 × 129 grid and Re=1000.

Table 1 .
η r is the residual error estimator and η is the local Poisson problem error estimator for leaky driven cavity, with Reynolds number Re =600.