Analysis of Estimators for Stokes Problem Using a Mixed Approximation

: In this work, we introduce the steady Stokes equations with a new boundary condition, generalizes the Dirichlet and the Neumann conditions. Then we derive an adequate variational formulation of Stokes equations. It includes algorithms for discretization by mixed ﬁnite element methods. We use a block diagonal preconditioners for Stokes problem. We obtain a faster convergence when applying the preconditioned MINRES method. 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. In order to evaluate the performance of the method, the numerical results are compared with some previously published works or with others coming from commercial code like Adina system.


Introduction
This paper describes a numerical solution of Stokes equations with a new boundary condition.This condition generalizes the known conditions, especially the Dirichlet and the Neumann conditions.We use the discretization by mixed finite element (MFE) method.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 2 A. El Akkad, A. Elkhalfi is defined with the help of Raviart Thomas basis functions [1,2,3] 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.
A wealth of literature on solving saddle point systems exists, much of it related to particular applications.Perhaps the most comprehensive work is the survey by Benzi, Golub and Liesen [10], which considers conditions under which the matrix formulation is solvable and block diagonal preconditioners but which has a focus on linear algebra.The conditions for a unique solution can be found in [23] or, in the substantial area of PDEs, in Babuska [32] and Brezzi [36].Stokes problems also arise in a natural way when the (unsteady) Navier-Stokes equations are simplified using classical operator splitting techniques [39].
A posteriori error analysis in problems related to fluid dynamics is a subject that has received a lot of attention during the last decades.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 [4,5] laid the basic foundation for the mathematical analysis of practical methods.Other works for the stationary Navier-Stokes problem have been introduced in [6,9,20,21,22,24,25].
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 block diagonal preconditioners for Stokes problems 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.

Model problem
We will consider the model of viscous incompressible flow in an idealized, bounded, connected domain in R 2 , The boundary value problem which is posed on two dimensional domains Ω, is defined as: We also assume that Ω has a polygonal boundary Γ , so − → n that is the usual outward pointing normal.The vector field − → u is the fluid velocity, p is the pressure field, ∇ is the gradient, ∇. is the divergence operator, the functional − → f in the space [L 2 (Ω)] 2 , − → g in the space [L 2 (Γ)] 2 , the pressure p in the space L 2 (Ω) and β is a nonzero bounded continuous function defined on ∂Ω.We set and Let the bilinear forms a : These inner products induce norms on V and W denoted by .V and .W respectively. ) ∀q ∈ W. (2.9) Given the continuous functional l : The weak formulation of the Stokes flow problem (2.1) for all ( − → v , q) ∈ V × W .

Finite element approximation
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.
For any K ∈ P , ω K is of rectangles sharing at least one edge with element K.We let ε h = ∪ K∈P ε(K) denotes the set of all edges split into interior and boundary edges, The finite element subspaces X h and M h are constructed in the usual manner so that the inclusion The mixed finite element approximation to (2.11)-(2.12) is then : 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

.5)
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 right-hand side vector f in (3.5) is for i = 1, ..., n u and k = 1, . . ., n p .We use the iterative methods Minimum Residual Method (MINRES) for solving the symmetric system.

Block diagonal preconditioners
We use a fast and robust linear solvers for stabilized mixed approximations of the Stokes equations (2.1)-(2.2)-(2.3).The resulting discrete Stokes system is the saddle-point system [9] Where the vectors U, P are discretized representations of − → u , p, with f , g taking into account the source term − → f as well as nonhomogeneous boundary conditions.The matrix C is the zero matrix when a stable finite element discretization Q 2 − Q 1 Taylor-Hood element is used, and is the stabilization matrix otherwise.
We assume that A is symmetric positive definite, which is the case when a Dirichlet condition is imposed on at least part of the boundary.The matrix C is always positive semi definite.
For consistency with the continuous Stokes system the matrix B should satisfy 1 ∈ null(B T ) in the case of enclosed flow (see [9,Chapter 3]).Let where D ∈ R n×n is symmetric positive definite as above, C ∈ R m×m is symmetric positive semi definite, B ∈ R m×n with m ≤ n and rank(B) = r ≤ m.
We suppose that the negative Schur complement of D, has rank p. Then under these conditions D has n positive eigenvalues, p negative eigenvalues and m -p zero eigenvalues [10].
Let the block diagonal preconditioner for D where H ∈ R n×n is a symmetric positive definite approximation to the Schur complement S. In the case where H = S and C = 0, B must be full rank for S, and hence P 1 , to be invertible, and the eigenvalues of the preconditioned system are given by [11,12] and in the case where the approximation of S (or indeed A) is inexact the preconditioner is frequently found to be extremely effective also.When the condition on C is weakened to allow the matrix to be symmetric positive semi-definite, and we have [13,14] We consider situations where the Schur complement could be singular, and we construct an inexact approximation H which is invertible.
For the Stokes equations, the approximate Schur complement is either the mass matrix associated with the pressure approximation space [9] or an approximation.Common approximations of Q are its diagonal [14,16], a lumped version [17], or a Chebyshev semi-iteration method applied to Q (see [18,40]).Instead of taking S ≈ H, we incorporate a scaling constant α > 0 such that as a potential preconditioner for D. We explain why setting a large value of á can significantly improve the performance of the iterative solver when a stabilized mixed approximation is employed.We order the eigenvalues of P α −1 D from smallest to largest, so that where We have the MINRES convergence bounds [9] where P k is the set of polynomials of at most degree k and σ(P This bound can be pessimistic, it will still provide some insight into the effect of α on preconditioned MINRES convergence.
For stable finite element discretizations Q 2 − Q 1 of the Stokes equation there exists an inf-sup constant γ , and a constant Υ resulting from the boundedness of B, such that For an unstable discretization only the upper bound holds, and a lower bound is assumed as follows : Theorem 4.1.We suppose that (4.12) holds.For the Stokes equations, we use a discretization by and, if H = diag(Q) we have Where, γ is the inf-sup constant, while φ = 1 if Dirichlet conditions are imposed on the whole boundary and φ = 2 otherwise.
Proof.Same steps of the prof of Lemma 4.1 in [38].
Theorem 4.2.We suppose that (4.13) holds.For the Stokes equations, we use a discretization by Then, the eigenvalues of Where, γ is as in (4.13) while φ = 2 if Dirichlet conditions are imposed on the whole boundary and φ = 3 otherwise.
Proof.Same steps of the prof of Theorem 4.1.

Analysis of estimators
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.
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) Hence for − → v , − → w ∈ V and q ∈ W we have (5.11) The velocity space V(P) is defined by (5.12) and the broken pressure space W (P ) is defined by Examining the previous notations reveals that 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.16) 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.23) For − → v ∈ V (P ), we have
Proof.The right-hand side of equation (5.25) vanishes en V × W . Applying theorem 5.3, we obtain (5.25).We define the linear functional R : for all ( − → w , q) ∈ V (P ) × W (P ).For ( − → w , q) ∈ V × W , we obtain (5.27) Let the lagrangian functional L : So that sup (5.30) Therefore, Using (5.31), we obtain: (5.33) We have the problems on each subdomain Suppose that the minimum exists, then the minimising element is characterized by finding − → φ K ∈ V K such that (see M. Ainsworth and J. Oden [5]) The necessary and sufficient conditions for the existence of a minimum are that the data satisfy the following equilibration condition: (5.38) When the subdomain K lies on the boundary ∂Ω the local problem (5.35) 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.36).We define Using (5.18), 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.45) For example, one might choose the piecewise bilinear pyramid functions associated with interior nodes in the partition.The relation (5.44) 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.46)where λ k KJ,A are constants to be determined.Owing the constraint (5.44), it is required that JK,A = 0, (5.47) for each A. Lemma 2. Suppose that for each X A the constants {λ Proof.The result follows immediately by using (5.46), (5.44) and (5.40).
Summarizing and incorporating the results of section 5 we have Theorem 5.5.There exists a constant C > 0 such that where We define the global error estimator η p by (5.54) 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 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 Proof.Same steps of the prof of Theorem 3.9 in [15].

Numerical simulation
Example 1 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 [31,33,34,35,40,41,42].This     Figure 1 shows the uniform streamline by MFE associated with a 64 × 64 square grid, Q 2 − Q 1 approximation.The particles in the body of the fluid move in a circular trajectory.Table 1 present the iteration numbers for MINRES solution of the leaky cavity problem using stabilized Q 1 − Q 1 elements on a uniform mesh.We observe that when α is increased, the iteration numbers clearly decrease, and there is hence a considerable benefit to applying the scaled preconditioner.This is observed for all values of mesh parameter tested.We present these results pictorially in Figure 3.The results, in Fig. 4, clearly show that increasing α reduces η, but that as we increase α beyond about 10, η decreases much more slowly.In other words, as α is increased beyond this point, we would not anticipate a further significant reduction in iteration numbers for our preconditioned solver.This therefore motivates a value of α equal to roughly 10, as this choice essentially achieves the optimal predicted convergence rate, while at the same time ensuring that the negative eigenvalues of P −1 α D are far from zero.The computational results of Figure 4 and Table 2 suggest that all two estimators seem to be able to correctly indicate the structure of the error.
Example 2 It' s a test problem with an exact solution is solved in order to compare the affectivity of two error estimation strategies: the residual estimator η r and the Poisson estimators η p .The latter approach is frequently used and is generally considered by practitioners to be one the best error estimation strategies in terms of its simplicity and reliability, especially when used as a refinement indicator in a self-adaptive refinement setting.This analytic test problem is associated with the following solution of the Stokes equation system:  To interpret the results that are presented some notation will be needed: e ωT = − → u − − → u h 2 V,ωT + P − P h 2 0,ωT .(6.4) Figure 6 shows the streamline and pressure plots, and figure 7 shows the estimated error η T associated with 64 × 64 square grid.
η G,T ≤ C( Here, the generic constant C Ω is independent of the mesh size and the exact solution but may depend on the domain and the element aspect ration.
Then the estimators η G,T is likely to be effective if it is used to drive an adaptive refinement process.
In general, if an error estimator is to be efficient then the constant on the right hand side of (6.6) should be bounded.An estimate of this constant (e.g max T ∈T h η G,T eω T ) is provided in Table 4, where we also estimate this constant for the exact error (e.g max T ∈P eT eω T ).

The local affectivity indices
η G,T eω T will be bounded above and below across the whole domain, so that elements with large errors can be singled out for local mesh refinement.This is assessed in Figure 9. Looking at the distribution of these indices it is clear that the our two estimators give a very different picture.Once again, η T is closely aligned with the exact error but η r,T is not.

Conclusion
We were interested in this work in the numeric solution for steady incompressible Stokes Equations with a new boundary condition.We use a discretization by mixed finite element methods.We use a block diagonal preconditioners and the preconditioned MINRES method for Stokes problem.We observe that when α is increased, the iteration numbers clearly decrease, and there is hence a considerable benefit to applying the scaled preconditioner.We obtain a faster convergence.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.

Remark 2. 1
If β is strictly positive constant such that β ≻≻ 1 then C β , is the Dirichlet boundary condition and if β ≺≺ 1 then the C β , is the Neumann boundary condition.For this, β is called the Dirichlet coefficient.
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 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: {y = 1, −1 ≤ x ≤ 1/u x = 1} a leaky cavity.The streamlines are computed from the velocity solution by solving the Poisson equation numerically subject to a zero Dirichlet boundary condition.

Figure 2 :
Figure 2: Velocity vectors solution by MFE (left) associated with a 64 × 64 square grid, Q 2 − Q 1 approximation and velocity vectors solution (right) computed with ADINA system.

Table 1 :
Results for the cavity problem solved with Q 1 − Q 1 finite elements, for a range of values of h and α

Figure 3 :
Figure 3: Representation of the effect of α on the MINRES iteration count for the cavity problem.

Figure 6 :
Figure 6: Fig. 6.Uniform streamline plot (left) and pressure plot (right) for the flow by MFE associated with a 64 × 64 square grid.

Table 3 :
Comparison of error estimator effectivityGride , then the Poisson problem estimator η p provides the most accurate estimate of the global error and the local estimates η T is quantitatively close to the exact error and the estimates η r is about three times larger than exact error.We see that the local error estimator η

Table 4 :
Comparison of affectivity indices T, is close to max T ∈T h eT eω T .