Photovoltaic Cells and Systems: SQU Journal for Science, 17(2) (2012) 187-199 © 2012 Sultan Qaboos University 781 Local Projection Stabilization Method for Incompressible Flow Problems A.M. Essefi* and K. Nafa** Department of Mathematics and Statistics, Sultan Qaboos University, College of Science, P.O. Box 36, Al-Khoudh 123, Muscat, Sultanate of Oman, *Email:essefi@squ.edu.om and **Email:nkamel@squ.edu.om. ABSTRACT: This paper analyzes Local Projection Stabilization (LPS) methods for the solution of Stokes problem using equal order finite elements. Their convergence, stability and accuracy properties are investigated. The resulting stabilized method is shown to lead to optimal rates of convergence for both velocity and pressure approximations. Two classes of LPS methods are distinguished: one-level and two-level methods. Numerical examples using bilinear interpolations are presented to validate the analysis and assess the accuracy of both approaches. KEYWORDS: Finite elements, Local projection stabilization, Rates of convergence, Stability. غير قابلة لالنضغاط الاستقرار طريقة االسقاط المحلي لحل مسائل السوائل اللزجة كمال نافع و ي افعفيفة الس .لعناصر محدودة متساوٌة درجة خدامستباستوكس مسألةلحل وضعًاإلسقاط الم طرقً هذا البحث ف نحلل :ملخص إلى أفضل مستقرة وتؤدي وضعًٌقة اإلسقاط المنبرهن ان طرتها. وخصائص دقتلك الطرق استقرارونتحرى تقارب مستوى واحد واثنٌن. نصف ذات اسقاطاتفٌن من الطرق: نصنمٌز بٌن نسب للتقارب سواء بالنسبة للسرعة أو للضغط. للتحقق من صحة التحلٌل وتقٌٌم مدى دقة كال النهجٌن.استقراءات شبه خطٌة أمثلة عددٌة باستخدام 1. Introduction umerical approximations of incompressible flows require a compatibility condition between the discrete velocity and pressure spaces (Girault and Raviart, 1986; Brezzi and Fortin (1991)). This condition prevents, in particular, the use of equal order interpolation spaces for the two variables, which is the most attractive choice from a computational point of view. The two-level local projection technique has been introduced by Becker and Braack (2001) and Nafa (2008) to circumvent the inf-sup condition and to allow the use of simple equal order interpolations such as 1 1 P P and 1 1 Q Q velocity-pressure approximations. This method consists in introducing the 2 L -projection of the pressure gradient as a new unknown of the problem. Hence, a third equation to enforce the projection property is added to the original discrete equations, N A.M. ESSEFI and K. NAFA 788 and a weighted difference of the pressure gradient and its local projection is introduced into the continuity equation. In this paper we analyze the LPS and distinguish two classes of LPS methods: The one-level enriched method and the two-level method. The stability and error estimate results presented here are more general than those given in Becker and Braack (2001). In fact, the proofs extend the stability and error estimate results of Becker and Braack (2001) and Nafa (2008) to the one-level enriched approximation case. Numerical results are presented to justify the order of convergence and assess the performance accuracy of both LPS methods using bilinear finite element interpolation. 2. Discrete Stokes problem Let  denote an open, bounded, connected subset of 2 , where  denotes its Lipschitz continuous boundary. Let   2 1 H      f ,   2 1 0H      V , and  20Q L  . We consider the usual Stokes problem with homogeneous Dirichlet boundary conditions. Find  , p Q u V satisfying in . 0 in on p          u f u u 0 (1) where  20L  is the set of square integrable functions with null average. The weak formulation of the above problem reads        , , . , . ,p q       u v v u f v (2) for all test functions  ,q Q v V where  .,. denotes the 2L -inner product on the region  . The compact form of the weak formulation of (1) is given by      , ; , ,B p q F q q Q   u v v v V, , (3) where        , ; - . .B p q p q     u v u v v u, , , , and    F q v f v, , 3. Local projection stabilization Let hT be a shape regular partition of the region  into quadrilateral elements K, and assume hV and hQ are finite dimensional subspaces of V and Q respectively, consisting of continuous functions. We denote by kh the local mesh size, the mesh size is then defined by max .K K Th h h   Then, the Galerkin discrete problem reads      , ; , .h h h h h h h h h hB p q F q q Q   u v v v V, , , (4) Let hM be a coarser mesh partition of the domain into quadrilateral elements M which are defined as the union of one or more neighboring elements K of the partition .hT We assume that the partition of  into macro-elements hM M is non-overlapping, shape regular and such that LOCAL PROJECTION STABILIZATION METHOD 781 , , M K hh h K M M M     where Kh and Mh denote the maximum diameter of K and M, respectively. Let  1hY H  be a scalar finite element space of continuous, piecewise polynomial functions over .hT Since we are interested in the case of equal order interpolation, we define the approximation spaces as follows 2 , .h h h hY Q Y Q V V Let  hD M be a discontinuous finite dimensional space defined on the macro-element .hM M We introduce the associated discontinuous global space  D D .h h M M h M    Figure 1. Reference macro-element (right) and macro-element (left). Then, we define the equal order approximation finite element spaces of velocity and pressure respectively by:              2 21 0 |K 1 2 0 |K / , H / , h r h h r h H Q K K T Q q L q Q K K T                   v vV (5) where, for each integer  0, rr Q K denotes standard quadrilateral finite elements by means of a reference element (as illustrated in Figure 1) as polynomials of degree less or equal to r with respect to each variable. Further, we define the space of discontinuous functions       2 22 |K 1 / , .h h h r hL Q M M Μ            D λ λ (6) We define the local 2 L -projection    2: DM hL M M  which generates the global projection  2: Dh hL   by      2|| , , .h M M hM M M L          To which, we associate the fluctuation operator    2 2:h L L    defined as ,h hid   where id stands for the identity operator in  2 .L  For simplicity we use the same notations id, ,M h and h for vector valued functions. Now, the stabilizing term is defined by     0, , .h h h M h h h h M M M h S B p q q       p , (7) Thus, the LPS discrete problem reads: Find  h h h hp Q u V, such that A.M. ESSEFI and K. NAFA 711      , ; , h h h h h h h h h h h hB p q F q q Q   u v v v V, , , (8) where      , ; , ; ,h h h h h h h h h h h hB p q B p q S B p q u v u v, , , (9a) and     .h h h h hF q F qv v, , (9b) The stabilized formulation given in (8) is written component-wise as      , , . , h h h h h h hp       u v v f v v V, (10a)      h0,, . , 0, h h M h h h h hM M M h q p p q q Q            u (10b)    h0, , 0, .M h h h h hM M M h p p          ξ ξ D (10c) In order to prove the stability and convergence of the solution of the stabilized method given in (8) we introduce the following assumptions. Assumption A1. The fluctuation operator satisfies the approximation property   ,0, , , , 0 . l l h M hl MM q Ch q q H M M M l r        (11) Assumption A2. There exists an interpolation operator  1:h hi H Y  such that    1 10 0v :h hi H Y H   with the error estimates  0, 1, , v v v v v s h K h KK K s K i h i Ch      (12) for all   v ,sH K and all , 1 1,hK T s r    where  K denotes a certain neighborhood of K. Assumption A3. Further, assume that the local inf-sup condition       1 v 0, 0, v , inf sup v h h M q D Mh h Y M h hM Mh h q q     (13) holds for all hM M with a positive constant 1 independent of the mesh size h. The following theorem defines hj and hj interpolation operators that are important for deriving error estimates for the LPS method (Matthies et al. 2007). Theorem 1. Assume that Assumptions A2 and A3 are satisfied. Then, there are interpolation operators  1:h hj H Y  and :h hj V V satisfying the following orthogonality and approximation properties    1, 0, , h h h hj D H          (14a)  0, 1, , s h M h MM M s M j h j Ch          (14b) and  , 0, , h h h h     ω j ω Φ Φ D ω V (15a) LOCAL PROJECTION STABILIZATION METHOD 717  0, 1, , . s h M h MM M s M h Ch     ω j ω ω j ω ω (15b) For all       2 , , , s s hH M H M M M   ω V and 1 1,s r   where  M denotes again a certain neighborhood of M. Next, we introduce the following mesh-dependent norm on the product space h hQV     2 2 2 2 0, 0, .h h h h M h h M M M h q q q         v v 1, , (16) 3.1 Stability Theorem 2. Assume that Assumptions A2 and A3 hold and the parameters M are such that 2 /M Mh C  for all elements .hM M Then, there is a positive constant independent of h such that             2 ,0 , , , ; , inf sup . , . , h h h h h q Qh h h h r Qh h h h h h h h B q r q r        0 v V ω V v ω v ω (17) The proof is found in (Nafa and Wathen, 2009). 3.2 Error estimates First, we introduce the following consistency error (Ganesan et al., 2008). Lemma 3. Assume that the fluctuation operator h  satisfies Assumption A1. Let       2 1r r p H Q H        u V, be the solution of Stokes problem and  h h h hp Q u V, be the solution of the LPS stabilized problem. Then, the consistency error can be estimated by     1/2 22 2 , , ; r h h h h M M h hr M M M h B p p q C h p q              u-u v v, , (18) for all   .h h h hq Q v V, Then, the error estimate of the LPS method is given by the following theorem. Theorem 4. Assume that the solution  pu, of (3) belongs to     2 1 , r r H Q H       V  h h h hp Q u V, is the solution of the local projection method (8), and M are such that 2 /M Mh C  for all elements .hM M Then, the following error estimates hold  1, 1, ,0, r h h r r p p Ch p         u u u (19) for some positive constant C independent of h. Proof. Let h hj uu and hp i p be the interpolants of the velocity and pressure respectively. By virtue of Theorem 2, there exist  h h h hq Q v V, satisfying   .h hq Cv , (20) We have A.M. ESSEFI and K. NAFA 711       2 22 1, 0, , , .h h h h h h h h h h h hp p p p S p p p p         u -u u -u (21) Hence,   22 1/2 1, 0, 1 / , .h h h h h h h hp p p p     u -u u -u (22) Thus, by equivalence of norms on finite dimensional spaces in 2  , it yields   22 1/2 1, 0, 2 / , .h h h h h h h hp p p p     u -u u -u (23) In addition, we have       2 , ; , 1 / h h h h h h h h h h h h h B p p q p p q     u -u v u -u v , , (24) i.e.           2 2 , ; , ; , 1 / 1 / . h h h h h h h h h h h h h h h h h h B p p q B p p q p p q q        u -u v u-u v u -u v v , , , , (25) Using the result of the consistency estimate with 2 ,M MCh  we obtain     , , ; rh h h h h r h h B p p q Ch p q    u-u v v , , (26) where,      , ; , ; .h h h h h h h h h h h hB p p q B p p q S B p p q    u -u v u -u v, , , The terms of  , ;h h h hB p p qu -u v , will be estimated using the approximation properties of the interpolations hj and ,hi      1, 1, 1, , r h h h h h hr Ch q           u u v u u v u v , (27) and     1, ,0, , . . r h h h h h hr p p C p p Ch p q       v v v , (28) To estimate the last term we shall use the integration by part and the orthogonality property of ,hj       . , , .h h h h h h hq q , q       u u u u u u (29) Thus,    1/2 1/2 221 0, 0, . , .h h M h M h hM M M M M Mh h q q                            u u u u (30) Hence,      1/2 1/2 2 22 1 2 2 1, 0, . , r h h M M M h hr M M M M M Mh h q C h q                               u u u     1/2 21 2 2 1, , r M M h hr M M M h C h q                 u v , (31) i.e. LOCAL PROJECTION STABILIZATION METHOD 711      1, . , , r h h h hr q C h q      u u u v (32) for some positive constant .C  Now, to estimate the stabilizing term  h h hS B p p q , we use the 2 L -stability of the fluctuation operator h and the approximation property of .hi We obtain the following:     , ,h h h M h h h h M M h S B p p q K p p K P         1/2 1/2 2 2 0,0, .M h h M h h MM M M M Mh h p p q                           (33) Thus by (11) we get       1/2 22 2 1, , r h h h M M h h hr M M M h S B p p q C h p p q                v,     1/2 22 2 , , r h h h M M h h hr M M M h S B p p q C h p p q              v, using 2 ,M MCh  we have     1/2 22 1 , , . r h h h M h hr M M M h S B p p q C h p q             v, (34) Hence     , , . r h h h h hr S B p p q Ch p q    v, (35) Using (26) and (35), we get  0,1, 1, , r h h h r r p p Ch p         u u u . □ The required error estimate then follows using the triangular inequality together with the interpolant estimates. As already observed, the existence of the interpolation operators hj and hj is fundamental, in LPS methods, for deriving the error estimates. Since, Assumption A3 may not be fulfilled, the existence of these operators is not always guaranteed. 3.3 Two-level local projection stabilization This class of stabilized schemes uses discontinuous pressure gradient projection approximations of degree 1r  which are defined on the coarser mesh 2 .h hM T Let 2 , 1h r  be the 2 L -projection onto hD then 2 , 1h h rk id   Here, for rectangular elements, the spaces hY and hD are chosen as , 1,2/ / . disc h h r h r hY D Q Q  (36) It has been proved in (Matthies et al., 2007) that the pairs given by (36) satisfy Assumptions A3. 3.4 One-level enriched equal order stabilization From the implementation point of view, a drawback of the two-level approach is the fact that the added stabilizing term produces a larger stencil which may not fit in the data structure of an available programming A.M. ESSEFI and K. NAFA 711 code. A general theory which allows the derivation of a novel class of local projection stabilization by enrichment of the approximation spaces has been established in (Matthies et al., 2007). This class of stabilized schemes uses approximation and projection spaces defined on the same mesh and leads to much more compact stencils than in the two-level approach. Let , 1, h h h rΜ T   be the 2 L -projection onto ,hD and , 1.h h rid    For rectangular elements, the spaces hY and hD are chosen as ,1 1,, / / b disc h h r hr h Y D Q Q  (37) or ,2 1,, / / b disc h h r hr h Y D Q P  (38) where      ,1 1ˆ,, ˆ ˆ ˆspan , 1, , b r r h ir h K Q K Q K b x i d     (39)      ,2 ˆ, 1,, ˆ ˆ ˆ b r h r hr h K Q K Q K b Q K  (40) and K̂ b is a biquadratic bubble function on the reference element ˆ .K Here also the pairs given by (37) and (38) satisfy Assumptions A3. Note that, the local projection method described above does not add an extra cost to the solution of Stokes algebraic system. This is due to the elimination of pressure gradient unknowns at element level (Nafa and Wathen, 2009). 4. Numerical results Numerical results for two Stokes problems are presented. The performance of the first order two -level method 1, 0,2/h hQ Q and the one-level method 1, 0, / b h hQ Q approximations are obtained for 2 0M Mh  where 0 is a constant. 4.1 Problem 1 The first problem consists of solving Stokes problem in the unit square with exact solution   1 2, ,Tu u p given by 2 2 2 2 1 2( , ) 2 (1 ) (1 ) (1 2 ), ( , ) 2 (1 )(1 2 ) (1 )u x y x x y y y y u x y x x x y y         and 2 ( , ) .p x y x x  Numerical results are obtained, using the LPS method, for 1.  Calculations were performed for 0 0.1  and 0 1  using the two-level and one-level methods. The results given in Figures 2-5 indicate that the error norms 0,h  u u and 1,h  u u converge at the predicted rates, while 0,h p p   seems to converge at a rate of about 3/2 as observed by Becker and Braack (2001) and by Nafa and Wathen (2009). LOCAL PROJECTION STABILIZATION METHOD 711 Figure 2. Velocity error norm using 1, 0,2/h hQ Q interpolation with 0 0.1  (left) and 0 1  (right). Figure 3. Pressure error norm using 1, 0,2/h hQ Q interpolation with 0 0.1  (left) and 0 1  (right). Figure 4. Velocity error norm using 1, 0,/ b h hQ Q interpolation with 0 0.1  (left) and 0 1  (right). A.M. ESSEFI and K. NAFA 711 4.2 Problem 2 To test the accuracy of the method, we consider the standard Poisseuille flow in the channel    0, 4 0,1 .   We prescribe a parabolic inflow profile and use natural outlet conditions. The Dirichlet boundary conditions are not imposed on the whole boundary and natural boundary conditions 0p     u n n are employed. The solution of the problem is given by         2 1 2, 1 4 1 2 , , 0 and , 8 2 .u x y y u x y p x y x      Figure 5. Pressure error norm using 1, 0,/ b h hQ Q interpolation with 0 0.1  (left) and 0 1  (right). Figure 6. Velocity profile using 1, 0,2/h hQ Q interpolation with 0 0.1  . LOCAL PROJECTION STABILIZATION METHOD 711 Figure 7. Pressure contours using 1, 0,2/h hQ Q interpolation with 0 0.1  . First, the problem is solved using 1, 0,2/h hQ Q approximation on a 32x16 grid. The obtained results are illustrated in Figures 6 and 7 respectively, for the velocity and pressure solutions. We observe that the approximated pressure solution is exact even near the boundaries. In addition, Figure 8 shows that the velocity profile at the outlet boundary is in perfect agreement with the exact solution there. The numerical solution of the Poiseuille problem is also performed on a 32x16 grid using the one -level enriched approximation 1, 0,/ . b h hQ Q We obtain similar results to the two-level method for the velocity profiles as shown in Figure 9. The pressure solutions obtained for 0 10,  10 2 , 10 3 , and 10 4 , respectively, are illustrated in Figure 9. We note that while varying the coefficient 2 0M Mh  over the wide range from 0.1 to 10 4 , the obtained results change considerably. These results show that the pressure solution for the enriched one -level method is more sensitive to the choice of the parameter α0 as compared to the two-level method. Figure 8. Velocity profile at the outlet of the channel using 1, 0,/ b h hQ Q interpolation with 0 0.1  (left) and 0 1  (right). A.M. ESSEFI and K. NAFA 718 5. Conclusion The convergence, stability and accuracy of LPS for the Stokes problem has been investigated using the one-level and two-level bilinear finite element interpolation. The numerical results show that the order of convergence for both approaches is as predicted theoretically for the velocity. But, the L 2 rate of convergence for the pressure is approximately 3/2 as noted by other researchers (Becker and Braack, 2001; Nafa and Wathen, 2009). Also, we observe that the computation of the Poisseuille solution produces the exact velocity independently of the parameter 0 for both methods. However, it is shown that the pressure is exact even near the boundaries for the two-level method and depends on the parameter 0 for the one-level enriched method. In fact, to obtain results similar to those obtained by the two-level method, we have to take large values for 0 . Figure 9. Pressure contours using 1, 0,/ b h hQ Q interpolation with different values of 0 . 6. References BECKER, R. and BRAACK, M. 2001. A Finite Element Pressure Gradient Stabilization for the Stokes Equations Based on Local Projections, Calcolo 38(4): 173-199. BREZZI, F. and FORTIN, M. 1991. Mixed and Hybrid Finite Element Methods, Springer Verlag, New York, USA. GANESAN, S., MATTHIES, G. and TOBISKA, L. 2008. Local Projection Stabilization of Equal Order Interpolation Applied to the Stokes Problem, Math. Comp., 77: 2039-2060. GIRAULT, V. and RAVIART, P.A. 1986. Finite Element Methods for Navier-Stokes Equations, Springer Verlag, Berlin, Germany. 0 (a) 10  2 0(b) 10  3 0( ) 10c   4 0(d) 10  LOCAL PROJECTION STABILIZATION METHOD 711 MATTHIES, G. SKRYPACZ, P., TOBISKA, L. 2007. A Unified Convergence Analysis for Local Projection Stabilisations Applied to the Oseen Problem, M2AN Math, Model. Numer. Anal., 41(4): 713-742. NAFA, K. 2008. Two-level Pressure Stabilization Method for the Generalized Stokes Problem, International Journal of Computer, Mathematics 3-4: 579-585. NAFA, K. and WATHEN, A. 2009. Local Projection Stabilized Galerkin Approximations for the Generalized Stokes Problem, Comput. Methods Appl. Mech. Engrg., 198: 877-883. Received 13 January 2011 Accepted 21 December 2011