A pressure-stabilized formulation of incompressible flow problems on anisotropic finite-element meshes

A pressure-stabilized formulation of incompressible flow problems on anisotropic finite element meshes Jordi Blasco ? Departament de Matem`atica Aplicada I, ETSEIB, Universitat Polit`ecnica de Catalunya, Avgda. Diagonal 647, 08038, Barcelona, Spain. Abstract We consider a pressure stabilized, finite element approximation of incompressible flow problems in primitive velocity–pressure variables, which is based on a projection of the gradient of the discrete pressure onto the space of discrete functions. Equal order interpolation for the velocity and the pressure can be employed with this formulation. The method introduced here is specially developed to be used on anisotropic finite element meshes with large element aspect ratios. Key words: Incompressible flows, Stabilized methods, Anisotropic meshes, Navier-Stokes equations 1 Introduction The numerical approximation of incompressible flow problems presents several difficulties. Apart from the unphysical oscillations in the numerical solution which may appear in convection-dominated flows, the treatment of the incompressibility constrain requires some special attention. The need to enforce the divergence-free condition on the velocity field gives the system a mixed character (see [13]). If standard approximations are employed, such as the Galerkin ?This author’s work was partially supported by the Spanish MCYT Project BFM2002-04613-C03-01 Email address: [email protected] (Jordi Blasco). URL: http://www-ma1.upc.es/blasco (Jordi Blasco). finite element method, the discrete spaces for the approximation of the velocity and the pressure are subject to the satisfaction of the well known inf-sup condition in order to get stable and convergent solutions. This condition prevents, in particular, the use of the same mesh and the same interpolation for the two variables. Equal interpolation is a very desirable choice not only in terms of ease of implementation but also from a computational viewpoint (see [35]). Several velocity-pressure finite element pairs have been developed which fulfill the inf-sup condition (see [2], [9] and [12], for instance). Techniques designed to stabilize unstable finite element pairs have also been studied (see [25], [36]). But the most successful methods have been developed under the idea of stabilized formulations, which basically modify the discrete problem in such a way that the inf-sup condition is not required any more. Some examples of pressurestabilized formulations of incompressible flows are: the Brezzi-Douglas method ([11]), the Douglas-Wang method ([24]), the Galerkin-Least-Squares or GLS method ([26], [27]), first-order least-square methods ([8]), term-by-term stabilization ([15]), (residual-free) bubble function methods ([3], [10]) and unusual stabilized methods ([4]), among several others. Some of these methods are closely related to the treatment of other sources of instability, such as convection. Another pressure stabilized, finite element formulation of steady, incompressible flows was developed and analyzed in [18], [19] and [22], and extended to the transient case in [7] and [20]. This technique was originally designed under the idea of introducing a projection of the gradient of the discrete pressure onto the space of discrete finite element functions, after which the discrete continuity equation is modified in a consistent way; this is why it was called Pressure- Gradient-Projection method (PGP). Afterwards, it was shown in [17] that in some cases this method can also be obtained within the general framework of subgrid scale models (see [29], [30]), if the so called orthogonal subscales approach is followed. Another pressure gradient projection method, based on local projections, was also considered in [6]. On the other hand, most studies of the application of the finite element method to the approximation of two and three dimensional problems, both in fluid mechanics and in other disciplines, rely on the assumption that the finite element meshes satisfy an aspect ratio condition, which essentially requires that the elements have a similar size in all the spatial directions. We will call these meshes isotropic. However, anisotropic meshes with elements having large aspect ratios are often used in many applications, such as the modelling of the hydrodynamics of the ocean in coastal regions if physical coordinates are employed (see [34]), where the horizontal and the vertical length scales may differ by 2 or 3 orders of magnitude. The design of suitable schemes to be used together with anisotropic meshes and their numerical analysis is nowadays an 2 active subject of research (see [1], [14], [16] and [32], for instance). In this paper we present a pressure stabilized, finite element formulation of incompressible flow problems which is based on a pressure gradient projection and which is specially developed to be used together with anisotropic meshes. This method is based on the consideration of different elemental stabilization parameters in each of the different spatial directions; these directional elemental stabilization parameters are computed in terms of the element size in each direction, using standard expressions for them. This way, the formulation of well adapted to mesh anisotropy. A similar anisotropically stabilized scheme was analyzed in [5] in the two dimensional, rectangular case using a stabilization technique for the Q1-Q1element. Also, an orthogonal subgrid scale method has been recently introduced in [23] which also emphasizes the use of anisotropic grids; several ways of computing the element sizes in order to obtain the stabilization parameters are compared in that reference, but these parameters are still the same in all directions within each element. The outline of the paper is as follows: Section 2 introduces the mathematical problem to approximate and some notation. In Section 3 the stabilized, pressure gradient projection method is described, first recalling the isotropic case and then introducing the anisotropic case. Some computational aspects of the discrete problem are discussed in Section 4, whereas in Section 5 numerical results are presented; these results show in particular that the computing times are significantly reduced with the anisotropic scheme compared to the isotropic one on very stretched meshes, and that better accuracy of the solution is also observed with the anisotropic method in some cases. 2 Incompressible flow equations Let us consider an open, bounded domain Ω ⊂IRd, with d= 2 or 3, which we assume to be polygonal or polyhedral. The boundary of Ω is denoted by Γ. If the region Ω is occupied by a Newtonian incompressible fluid in motion at a steady state, the velocity uand the pressure pof the fluid at each point of Ω satisfy the incompressible Navier-Stokes equations: (u·∇)u−ν∆u+∇p=fin Ω (1) ∇·u=0 in Ω (2) u=0on Γ (3) Here, and in what follows, boldface characters denote vector quantities. In the momentum equation (1), ∇and ∆ are the gradient and Laplacian operators, respectively, ν > 0 is the fluid’s kinematic viscosity (which we assume constant), and fis a given forcing term. In the continuity equation (2), ∇·uis 3 the divergence of the velocity field; this equation, therefore, enforces incompressibility. Finally, we have assumed the homogeneous Dirichlet boundary condition (3) for simplicity. In order to study problem (1)-(2)-(3), the following Hilbert spaces are usually considered: L2(Ω) denotes the space of square integrable functions on Ω, with scalar product (u, v) and norm kuk;Hm(Ω) is the space of functions of L2(Ω) with distributional derivatives up to order malso belonging to L2(Ω); H1 0(Ω) is the closed subspace of H1(Ω) of functions that vanish on Γ, a space on which k∇ukis a norm equivalent to that induced by H1(Ω); H−1(Ω) is the dual space of H1 0(Ω), the duality between the two been denoted by < , >. The corresponding vector (product) spaces are likewise denoted by boldface characters in all cases. Assuming that f∈H−1(Ω), the weak form of problem (1)-(2)-(3) consists in finding u∈V:= H1 0(Ω) and p∈Q:= L2 0(Ω) (the space of functions of L2(Ω) with zero mean on Ω) such that: ((u·∇)u,v) + ν(∇u,∇v)−(p, ∇·v) =<f,v>, ∀v∈V(4) (∇·u, q) = 0,∀q∈Q(5) If the domain Ω is Lipschitz continuous, problem (4)-(5) has at least one solution (see, for instance, [28]), which is unique under certain conditions (in particular, for sufficiently large ν). If Ω is smooth enough, or if it is a bounded convex polygon, and f∈L2(Ω), then the solution has the additional regularity u∈H2(Ω) and p∈H1(Ω) (see [28]); we assume the latter regularity for the pressure from now on. 3 Numerical approximation 3.1 Isotropic stabilized method We recall here the stabilized, pressure-gradient-projection formulation analyzed in [19]. Let Θhbe a partition of Ω (of size h > 0) into finite elements K∈Θh; we assume that each element Kis the image of a reference element ˆ Kby an affine mapping FK:ˆ K→K. We consider the following spaces of discrete functions for the numerical approximation of problem (4)-(5): Qh=nqh∈ C0(Ω) /∀K∈Θh, qh|K= ˆqh◦F−1 K,ˆqh∈Rkp(ˆ K)o Vh=nvh∈(C0(Ω))d/∀K∈Θh,vh|K=ˆ vh◦F−1 K,ˆ vh∈(Rku(ˆ K))do Vh,0=nvh∈Vh/vh|Γ=0o 4 where Rk(D) is the polynomial space of degree kon D, if ˆ Kis a simplex, or the space of polynomials of degree less than or equal to kin each variable, if ˆ Kis a square or a cube. Qhis the space of discrete pressures and Vh,0is the space of discrete velocities, whereas Vhwill be the space of the projection of the pressure gradient. Given now a set of elemental parameters {αK}K∈Θh such that αK>0∀K∈Θh, to be determined later on, we consider the following stabilized method, where for each element Kwe denote by (u, v)K the scalar product in L2(K): Find (uh, ph,rh)∈Vh,0×Qh×Vhsuch that, for all (vh, qh,sh)∈Vh,0×Qh×Vh: ((uh·∇)uh,vh) + ν(∇uh,∇vh)−(ph,∇·vh) = <f,vh>(6) (∇·uh, qh) + X K∈Θh (αK∇ph−√αKrh,∇qh)K= 0 (7) (rh,sh)−X K∈Θh √αK(∇ph,sh)K= 0 (8) It is clear from (8) that rhis the orthogonal projection in the space L2(Ω) of the gradient of the discrete pressure (∇ph), scaled within each element by √αK, onto the space of discrete vector functions Vh. This method was proved in [19] to yield stable and optimally convergent solutions, provided the approximating spaces Vh,0and Qhsatisfy a mild compatibility condition, which was proved to hold for most equal order interpolations (ku=kp). The determination of the stabilization coefficients αKis one of the key issues in the development of stabilized methods. The numerical analysis of the convergence of the scheme (6)-(7)-(8) requires that there exist two positive constants α0and α1such that, for all K∈Θh: α0h2 K≤αK≤α1h2 K where hKis the size of element K. We therefore compute these parameters by the simple expression: αK= ¯αh2 K ν,∀K∈Θh where ¯αis a fixed constant; however, more general formulations will also be considered later on. 3.2 Anisotropic stabilized method The stabilization coefficients αKof the method just described depend on the overall element size hK, irrespective of the element geometry. In several applications, very stretched meshes are often used, with very different element 5 sizes in the different space directions. In those cases, it would seem more appropriate to consider different stabilization parameters within each element in each direction. Thus, we now consider the following anisotropically stabilized scheme: Find (uh, ph,rh)∈Vh,0×Qh×Vhsuch that, for all (vh, qh,sh)∈Vh,0×Qh×Vh: ((uh·∇)uh,vh) + ν(∇uh,∇vh)−(ph,∇·vh) = <f,vh>(9) (∇·uh, qh) + X K∈ΘhX i=1,...,d (αK,i∂iph−√αK,irh,i, ∂iqh)K= 0 (10) (rh,sh)−X K∈ΘhX i=1,...,d √αK,i (∂iph, sh,i)K= 0 (11) where rh= (rh,1,...,rh,d), sh= (sh,1,...,sh,d) and now αK,i = ¯αh2 K,i ν,∀i= 1,...,d, ∀K∈Θh Here, hK,i is the size of element Kin the i-th direction (which we compute as the maximum absolute value of the difference between the i-th coordinates of pairs of nodes in K). 4 Computational aspects Both the isotropic ((6)-(7)-(8)) and the anisotropic ((9)-(10)-(11)) PGP methods can be written in matrix form as follows: K(U)U+GP =Fu(12) DU +LαP−DαR=Fp(13) −GαP+MR = 0 (14) where: •U,Pand Rrepresent the vectors of nodal values of velocity (without its prescribed values), pressure and pressure-gradient-projection, respectively. •K(U) is the stiffness matrix, accounting for convection and diffusion. •Gand Dare the gradient and divergence matrices, respectively, where the rows (respectively, columns) associated to prescribed velocity values have been omitted. •Lα,Dαand Gαare the pressure Laplacian, PGP divergence and pressure gradient matrices, respectively, which take into account the effect of the 6 weighting in each element with the stability coefficients αK(αK,i, in the anisotropic case). •Mis the mass matrix. •Fuand Fpare known vectors arising from the RHS term fin (1) and eventually from inhomogeneous boundary conditions. In Reference [22], a detailed comparison of several iterative schemes designed to solve the nonlinear, coupled problem (12)-(13)-(14) was given; the most efficient method proved to be a pair of nested loops consisting of an outer iteration to account for the nonlinearity of the problem and an inner iteration segregating the computation of the velocity and the pressure from that of the pressure gradient projection. Namely, given an initial approximation U0to the velocity solution, an iterative method is considered to solve the nonlinearity of the problem; here, we employ a Picard method, although other schemes (such as Newton-Raphson’s) were also tested. Thus, for each j≥1, we have to compute (Uj, Pj, Rj) such that: K(Uj−1)Uj+GPj=Fu(15) DUj+LαPj−DαRj=Fp(16) −GαPj+MRj= 0 (17) Then, in order to solve the linear problem (15)-(16)-(17) for each j, we take R0 j=Rj−1and obtain each new iterate (Ui j, Pi j, Ri j), for i≥1, by calculating the velocity and the pressure assuming a known value of the PGP, which is then updated at the end of the iteration, in a block-Gauss-Seidel way: K(Uj−1)Ui j+GPi j=Fu(18) DUi j+LαPi j=Fp+DαRi−1 j(19) MRi j=GαPi j(20) This way, the system matrix for the velocity-pressure problem is unaltered through the inner iteration process (and thus it only needs to be computed once every outer iteration) and the system matrix for the PGP problem is the mass matrix (which can be lumped leading to a diagonal system). 5 Numerical test cases We present in this Section some numerical results obtained with the anisotropic, PGP stabilized method (9)-(10)-(11) on some standard test cases; we include 7 also the results obtained with the isotropic scheme (6)-(7)-(8) for comparison purposes. 5.1 Kovasznay flow We first consider a problem introduced by Kovasznay (see [31]), modelling laminar flow behind a two dimensional grid, in which an analytical solution of the two-dimensional, steady, incompressible Navier-Stokes equations (1)-(2) with no forcing term is available. The velocity solution u= (u, v) is given by: u(x, y) = 1 −eλx cos(2πy) (21) v(x, y) = λ 2πeλx sin(2πy) (22) for (x, y)∈IR2, whereas the pressure is: p(x, y) = p0−1 2e2λx (23) where p0is an arbitrary constant and the parameter λis given in terms of the Reynolds Re number by: λ=Re 2−(Re2 4+ 4π2)1/2<0 We solved this flow problem in the domain Ω = [−1 2,1] ×[−1 2,1 2] for a value of the Reynolds number equal to 40, prescribing the velocity values of the exact solution (21)-(22) on the boundary of Ω. This problem was also solved in [22], where the isotropic PGP method was used on four uniform meshes and with four different types of elements of first and second order on each mesh; a comparison with the GLS method (equivalent in this case to an Algebraic Subgrid Scale approach) showed that the PGP method gives more accurate results for the pressure solution than the GLS method if first order elements are used, the velocity errors being the same for the two schemes. However, both in this and in other test cases solved it was found that the PGP formulation requires about 20% more computational time that the GLS method to reach convergence. In order to compare the performance of the isotropic and the anisotropic versions of the PGP stabilized formulation, we solved this problem on four meshes of bilinear quadrilateral elements (Q1) made up from 20 uniform subdivisions in the yvariable and 60, 300, 1500 and 3000 uniform subdivisions in the x 8 Mesh Aspect ||u−uh|| ||p−ph|| ||∇(u−uh)|| ||∇(p−ph)|| ratio 60 ×20 2 I: 7.7·10−31.0·10−24.6·10−12.7·10−2 A: 7.7·10−39.9·10−34.6·10−12.7·10−2 300 ×20 10 I: 7.7·10−31.1·10−24.6·10−12.0·10−2 A: 7.7·10−31.0·10−24.6·10−12.0·10−2 1500 ×20 50 I: 7.7·10−31.1·10−24.6·10−12.0·10−2 A: 7.7·10−31.1·10−24.6·10−12.3·10−2 3000 ×20 100 I: 7.7·10−31.0·10−24.6·10−12.2·10−2 A: 7.7·10−31.1·10−24.6·10−15.1·10−2 Table 1 Kovasznay flow: aspect ratio and errors on four increasingly stretched meshes. Mesh Isotropic Anisotropic Anisotr./Isotr. 60 ×20 3.36 2.76 82% 300 ×20 31.47 16.22 51% 1500 ×20 291.79 79.21 27% 3000 ×20 713.88 156.56 22% Table 2 Kovasznay flow: CPU time (in seconds) for convergence on each mesh for the isotropic and the anisotropic schemes, and percentage of the latter to the former. variable, respectively; the four meshes are, thus, increasingly stretched. A tolerance value of 10−3in the relative error of the velocity, the pressure and the pressure gradient projection was set for the convergence of the nonlinearity in all cases, while a value of 10−4was used for the tolerance of the inner, PGP iterations. In this numerical example, the method of nested iterations (18)- (19)-(20) was used to solve the discrete problems (6)-(7)-(8) and (9)-(10)-(11), with a direct method for the solution of the linear systems of equations for the velocity and the pressure and a consistent mass matrix for the pressure gradient projection. All the computations were performed on a PC with a AMD processor at 1.2 Gh. The aspect ratio of the elements of each mesh is given in Table (1), together with the errors with respect to the exact solution (21)-(22)-(23) for both the velocity and the pressure, in the norms of L2(Ω) and H1(Ω) and for both the isotropic (I) and the anisotropic (A) methods. Very similar accuracy is achieved with the two schemes in all cases, the error being governed by the (relatively large) vertical spacing. 9 [5] R. Becker and R. Rannacher, Finite element solution of the incompressible Navier-Stokes equations on anisotropically refined meshes, Notes Numer. Fluid Mech. 49 (1995). [6] R. Becker and M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo 38 (2001) 173-199. [7] J. Blasco and R. Codina, Space and time error estimates for a first order, pressure stabilized finite element method for the incompressible Navier-Stokes equations, Applied Numer. Math. 38 (2001) 475-497. [8] P. Bochev, Z. Cai, T.A. Manteuffel and S.F. McCormick, Analysis of velocityflux first-order system least-square principles for the Navier-Stokes equations. I, SIAM Jour. Numer. Anal. 35 (1998) 990-1009. [9] D. Boffi, Three-dimensional finite element methods for the Stokes problem, SIAM Jour. Numer. Anal. 34 (1997) 664-670. [10] F. Brezzi, M.O. Bristeau, L.P. Franca, M. Mallet and G. Rog´e, A relationship between stabilized finite element methods and the Galerkin method with bubble functions, Comp. Meth. Appl. Mech. Eng. 96 (1992) 117-129. [11] F. Brezzi and J. Douglas, Stabilized mixed methods for the Stokes problem, Numer. Math. 53 (1988) 225-235. [12] F. Brezzi and R.S. Falk, Stability of higher–order Hood–Taylor methods, SIAM Jour. Numer. Anal. 28 (1991) 581-590. [13] F. Brezzi and M. Fortin, Mixed and Hybrid finite element methods, (Springer Series in Computational Mathematics, 15. Springer–Verlag, 1991). [14] C. Canuto and A. Tabacco, An anisotropic functional setting for convection diffusion problems, East-West Jour. Numer. Math. 9 (2001) 199-231. [15] T. Chac´on, A term-by-term stabilization algorithm for the finite element solution of incompressible flow problems, Numer. Math. 79 (1998) 283-319. [16] S. Chen, Y. Zhao and D. Shi, Anisotropic interpolations with application to nonconforming elements, Applied Numer. Math. 49 (2004) 135-152. [17] R. Codina, Stabilization of incompressibility and convection through orthogonal subscales in finite element methods, Comp. Meth. Appl. Mech. Eng. 190 (2000) 1579-1599. [18] R. Codina and J. Blasco, A finite element formulation for the Stokes problem allowing equal velocity-pressure interpolation, Comp. Meth. Appl. Mech. Eng. 143, 3-4 (1997) 373-391. [19] R. Codina and J. Blasco, Analysis of a pressure–stabilized finite element approximation of the stationary Navier-Stokes equations, Numer. Math. 87 (2000) 59-81. 16 [20] R. Codina and J. Blasco, Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient projection, Comp. Meth. Appl. Mech. Eng. 182, 3-4 (2000) 277-300. [21] R. Codina and J. Blasco, Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales, Comput. Visual. Sci. 4 (2002) 167-174. [22] R. Codina, J. Blasco, G.C. Buscaglia and A. Huerta, Implementation of a stabilized finite element formulation for the incompressible Navier-Stokes equations based on a pressure gradient projection, Int. Jour. Numer. Meth. Fluids 37 (2001) 419-444. [23] R. Codina and O. Soto, Approximation of the incompressible Navier-Stokes equations using orthogonal subscale stabilization and pressure segregation on anisotropic meshes, Comp. Meth. Appl. Mech. Eng. 193 (2004), 1403-1419. [24] J. Douglas and J. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp. 52 (1989) 495-508. [25] M. Fortin and S. Boivin, Iterative stabilization of the bilinear-velocity, constantpressure element, Int. Jour. Numer. Meth. Fluids 10 (1990) 125-140. [26] L.P. Franca and S.L. Frey, Stabilized finite element methods: II. The incompressible Navier-Stokes equations, Comp. Meth. Appl. Mech. Eng. 99 (1992) 209-233. [27] L.P. Franca and T.J.R. Hughes, Convergence analyses of Galerkin least– squares methods for symmetric advective–diffusive forms of the Stokes and incompressible Navier–Stokes equations, Comp. Meth. Appl. Mech. Eng. 105 (1993) 285-298. [28] V. Girault, P.A. Raviart, Finite element methods for Navier-Stokes equations, (Springer, Berlin 1986). [29] T.J.R. Hughes, Multiscale phenomena: Green’s function, the Dirichet-to- Neumann formulation, subgrid scale models, bubbles and the origins of stabilized formulations, Comp. Meth. Appl. Mech. Eng. 127 (1995) 285-298. [30] T.J.R. Hughes, G.R. Feij´oo, L. Mazzei and J.B. Quincy, The variational multiscale method- a paradigm for computational mechanics, Comp. Meth. Appl. Mech. Eng. 166 (1998) 3-24. [31] L.I.G. Kovasznay, Laminar flow behind a two-dimensional grid, Proc. Cambridge Philos. Soc. 44 (1948). [32] G. Kunert, Robust a posteriori error estimation for a singularly perturbed reaction-diffusion equation on anisotropic tetrahedral meshes, Advances in Comp. Math. 15 (2001) 237-259. [33] B. Jiang, T.L. Lin and L. A. Povinelli, Large-scale computation of incompressible viscous flow by least-squares finite element method, Comp. Meth. Appl. Mech. Eng. 114 (1994) 213-231. 17 [34] M.A. Maidana, M. Espino and J. Blasco, Un m´etodo estabilizado en elementos finitos 3D para el estudio de las corrientes oce´anicas, Proceedings of the Congresso de M´etodos Computacionais em Engenharia, Lisboa, Portugal (2004). [35] S. Norburn and D. Silveser, Stabilized vs. stable mixed methods for incompressible flows, Comp. Meth. Appl. Mech. Eng. 166 (1998) 131-141. [36] D. Silvester and N. Kechkar, Stabilized bilinear–constant velocity–pressure finite elements for the conjugate gradient solution of the Stokes problem, Comp. Meth. Appl. Mech. Eng. 79 (1990) 71-86. 18