Augmented Stabilized and Galerkin Least Squares Formulations

We study incompressible fluid flow problems with stabilized formulations. We introduce an iterative penalty approach to satisfying the divergence free constraint in the Streamline Upwind Petrov Galerkin (SUPG) and Galerkin Least Squares (GLS) formulations, and prove the stability of the formulation. Equal order interpolations for both velocities and pressure variables are utilized for solving problems as opposed to div-stable pairs used earlier. Higher order spectral/hp approximations are utilized for solving two dimensional computational fluid dynamics (CFD) problems with the new formulations named as the Augmented SUPS (ASUPS) and Augmented Galerkin Least Squares (AGLS) formulations. Excellent conservation of mass properties are observed for the problem with open boundaries in confined enclosures. Inexact Newton Krylov methods are used as the non-linear solvers of choice for the problems studied. Faithful representations of all fields of interest are obtained for the problems tested.


Introduction
Computational fluid dynamics problems are usually solved with the finite volume and finite difference approaches.Finite element procedures have not been able to enjoy the same popularity as finite volume based methods because of ill-conditioning of the mixed weak form Galerkin formulation, satisfaction of the restrictive inf-sup condition or the LBB condition (Ladyzhenskaya-Babuska-Brezzi) condition to ensure coercivity of the functional, along with inherent central difference approximation of lower order finite element methods that become unstable at moderate to high Reynolds numbers.Various methods have been proposed to alleviate these difficulties for solving CFD problems with finite element procedures.The methods that have gained popularity are based on Penalty finite element methods (Reddy, 1982), Least Squares type approaches (Pontaza and Reddy, 2003), and Stabilized Finite element methodologies (Brooks andHughes, 1980, 1982;Hughes and Brooks, 1982).Least Squares finite element methods have issues with extensive loss of mass in contraction regions in addition to increase of the degrees of freedom to be solved with introduction of vorticity or stresses in two dimensional problems (Ranjan and Reddy, 2012).Vorticity and stresses are auxillary variables and often not of interest in flow computations.Penalty finite element formulations interpret the incompressibility as a constraint in the Navier-Stokes equations and involve solution of velocity fields as the only degrees of freedom.Penalty finite element methods have been used for solving incompressible flow however they often require a very high value of the Penalty parameter to ensure a divergence free velocity field and for open boundary problems the Penalty parameters can often lead to a very ill-conditioned global matrix which is not amenable for solution with iterative solvers.These problems prohibit the usage of such procedures for solutions of large scale problems.Stabilized finite element methods on the other hand provide a variational setting for the Navier-Stokes equations and provide a viable alternative for solving incompressible flow problems in primitive variables.Stabilization methods for solving advection-dominated problems were introduced by Hughes andBrook (1979, 1982) and Brooks andHughes (1980, 1982).The optimal upwinding schemes provide nodally exact solutions to advection diffusion problems in one dimension (Brooks and Hughes, 1982) and generalizations to multi-dimensions were also outlined.The streamline upwind finite element method added an extra term to the under-diffused Galerkin finite element formulation, which greatly improved the conditioning of the resulting stiffness matrix.
Streamline Upwind Petrov Galerkin (SUPG) finite element method was first proposed by Brooks and Hughes (1982).The SUPG formulation with additional pressure stabilization (PSPG) first appeared in Hughes et al. (1986) for the Stokes problem and was later generalized to the solution of Navier-Stokes equations.The SUPS (Streamline Upwind and Pres-sure Stabilization) (Bazilevs et al, 2012) formulation provided a variational setting for the solution of incompressible flow problems, at the same time that they admit equal order interpolations for the velocity and pressure variables.In addition to removing the requirement of satisfaction of the LBB criterion they also improve the conditioning of the discrete form.Pressure stabilization yields a better conditioned matrix, consequently facilitating the use of iterative solvers.An excellent presentation of stabilized finite element methods for solving incompressible flow equations has been provided in (Bazilevs et al, 2012) and the references therein.An earlier report highlighted results obtained from the augmented SUPS formulation for incompressible flow (Ranjan et al, 2016c) for solving both two dimensional and three dimensional problems.Further study on effective scalable pre-conditioners for solving incompressible flow problems, and introduction of inexact Newton-Krylov methods for solving spectral method based CFD problems with the augmented SUPS formulation have been provided in Ranjan et al. (2016b,c).
Lower order finite element methods require extensive meshing in areas of high gradients of the evolving fields.Spectral/hp finite element methods alleviate the extensive meshing requirements and further provide spectrally accurate results.With the advent of higher order spectral/hp finite element methods (Karniadakis and Sherwin, 1999), it was realized that very accurate numerical solutions could be obtained with coarser meshes for solving complicated problems in both structures and computational fluid dynamics areas.Similar to hp finite element methods the spectral/hp element methods use orthogonal polynomials (Legendre or Chebyshev based), which offer orthogonality and better conditioning as compared to Lagrange based expansions for higher than 4 th order polynomials.In due considerations of the advantages of spectral/hp finite element procedures we employ spectral/hp methods for solving incompressible flow problems.
In literature artificial compressibility and Penalty finite element methods have been proposed as closely related ideas which have known to provide superior convergence behavior for Navier-Stokes equations.The Penalty method which leads to the solution of a perturbed problem has been proved to provide solutions that are not too far from the original problem and the distance between the two solutions goes to zero as the penalty parameter ε goes to zero (Ern and Guermond, 2004).Even in commercial codes artificial compressibility is often utilized to obtain solutions to the stiff linear system generated from the Navier-Stokes equations when there are convergence problems during solution.A related idea to the Penalty finite element method is the iterative penalty finite element method.Iterative penalty methods have been used for solving incompressible flow problems by Codina (Codina et al, 1994;Codina Rovira et al, 1992).Discontinuous pressure and bilinear velocity div-stable pairs were utilized for solving incompressible flow problems.Recovery of the post computed pressures from a discontinuous cell based field to a continuous field for presentation purposes was obtained via a L 2 least squares projection from the center point values to the quadrature points while post processing.They eliminated the pressure from the solution and solved for the velocities only.Usage of unequal interpolation which are LBB stable requires a complex coding structure as compared to equal order interpolations if pressure is carried as a primary variable in the formulation.Further a discontinuous pressure field projected to the finite element space with a least squares projection provides a smooth pressure field at the quadrature points at a price of loss of consistency of the momentum equations.In contrast, in this work we solve for the coupled monolithic problem carrying both pressures and velocities as primary variables in the formulations.We utilize an equal order interpolation for both velocities and pressures in the formulations which is consistent with the governing differential equations being solved.We attempt to bring the advantages of stabilized finite element formulations with iterative penalization of the divergence free constraint to improve the velocity-pressure coupling in solving incompressible flow problems.We name the formulations as Augmented Stabilized SUPS (ASUPS) and Augmented Galerkin Least Squares (AGLS) formulations for solving incompressible flow equations.We provide the stability analysis of both formulations and demonstrate for benchmark problems the formulations provide better nonlinear convergence.With this improvement in the formulations we solve some benchmark CFD problems and report results.Some of the advantages of the present formulation are as follows; • It circumvents the need for artificial mass matrix enhancement with the predictor corrector step while solving incompressible flow problems.The true mass matrix is utilized for solving all problems.
• It allows the usage of standard monolithic solvers in conjugation with inexact newton krylov methods for solving problems.
• Non linear convergence is achieved with the Navier-Stokes in their natural form.This is attributed to an adequate addressing of the velocity-pressure coupling as opposed to a total loss of stability of the discrete form with SUPG formulations.
• Allows the use of alpha family of approximation for stepping the equations in time in a space time decoupled framework.
• More scales in the problem are resolved than are possible with an artificially enhanced mass matrix in effect accessing a larger bandwidth of fine scale problem.
We introduce some notation for the following discussion and proofs.
• The penalty parameter, ε is a small positive constant.
• Lower case letters (i, j, k, l, m, n) to denote positive integers (not bold).
• Utilize the lower case letter h to denote the grid-size (not bold).
• We denote by R N the Euclidean space of dimensions N = (2, 3).
) denote positive constants independent of viscosity and of any element diameter h K .
This paper is organized as follows.In section 2, we outline the incompressible Navier-Stokes equations and introduce the iterative penalty method for relaxing the divergence free condition.In section 3 we perform the stability analysis of the ASUPS formulation.In section 4 we perform the stability analysis of the AGLS formulation.Spatial discretization with the spectral/hp finite element formulation is introduced in section 5. Section 6 provides details on the non-linear solution techniques and time stepping algorithms.Section 7 provides an array of sample CFD problems solved with the ASUPS and AGLS formulations respectively.Section 8 contains numerical experiments, and Section 9 concludes the paper.

Incompressible Flow Equations and Penalty Techniques
Let us consider the flow of a Newtonian fluid with density ρ, and viscosity µ.Let Ω ⊂ R N and t ∈ [0, T ] be the spatial and temporal domains respectively, where N is the number of space dimensions.Let Γ denote the boundary of Ω.The equations of fluid motion that describe the unsteady Navier-Stokes governing incompressible flows are provided by where, u, f, and p are the velocity, body force, and pressure respectively.Re denotes the Reynolds number, ν = µ ρ is the kinematic viscosity of the fluid, p is the pressure and u is the fluid velocity.The part of the boundary at which the velocity is assumed to be specified is denoted by Penalty method relaxes the incompressibility condition as In the above equation 5 ε is a small positive constant.Penalty procedures are known to improve the convergence of the mixed finite element formulation.However, they can still provide an unstable computation for very low values of the penalty parameter.Iterative penalty methods relax this condition on the penalty parameter and allow for larger values while maintaining stability of the formulation (Dai and Cheng, 2008).Iterative penalty thus introduces a regularization step to enforce the incompressibility constraint.We employ the formulation of this method as outlined in (Dai and Cheng, 2008;Lin et al, 2004;Xiao-liang and Shaikh, 2006).The governing equations take the form where, k denotes the regularization step.The natural boundary conditions associated with eq (1) are the conditions on the stress components, and these are the conditions assumed to be imposed on the remaining part of the boundary, where Γ g and Γ h are the complementary subsets of the boundary Γ, or Γ = Γ g ∪ Γ h .In eq (8) we denote the outward normal by n, the stress tensor by σ and the traction components by h.As the initial condition, a divergence free velocity field, u 0 (x), is imposed.Error estimates for the initial value of the iterative penalty method obtained from eq.4-5 for Stokes equations (dropping the non-linearity) have been provided in (Carey and Krishnan, 1982;Oden et al, 1982;Stenberg, 1990).Similar error estimates of the iterative penalty method eq.6-7 for the Stokes problem have been provided in (Xiaoliang and Shaikh, 2006).

Stability Analysis of ASUPS Formulation
Let us study the stability of the SUPS finite element formulation resulting from the iterative regularization of the divergence free constraint.We prove the stability of the formulation based on the initial value of the algorithm (detailed below).We consider the steady state form of the equations for analysis.The linearized Navier-Stokes equations are defined on a bounded domain Ω ⊂ R N , N = 2, 3, with a polygonal or polyhedral boundary Γ.Let C h denote a partition of Ω into elements consisting of triangles or convex quadrilaterals.Let us denote L 2 (Ω) as the space of square integrable functions in Ω, and L 2 0 the space of L 2 functions with zero mean value in Ω.The L 2 inner product in Ω is (., .).Further, we also employ (., .)K the L 2 inner product and ||.|| 0,K the L 2 norm in the element domain K. Let C 0 (Ω) be the space of continuous functions in Ω. H 1 0 denotes the Sobolev space of functions with square-integrable value and derivatives in Ω with zero value on the boundary Γ.H 1 norm of a function u is Further we denote the triangular or quadrilateral decomposition of the domain by the following (Ciarlet, 2002;Johnson, 2012;Oden and Reddy, 2012) where for each integer m ≥ 0, P m and Q m have the usual meaning.
Let the finite element spaces which we wish to work with be given by where the integers k and l are greater than or equal to 1.

Consistency
Momentarily reverting back the notation introduced earlier let us denote (u, p) the solution of equation 1-2 and (u ε , p ε ) the solution to 6-7.It has been shown the following error estimates hold for the weak finite element formulation of the iterative penalty method (Xiao-liang and Shaikh, 2006) where the constant C 0 is independent of h and ε.For the above estimates we define H s (Ω) to be the standard Sobolev spaces with corresponding inner products (., .)s,Ω and norm ∥.∥ s,Ω .The formulation is consistent since the pressure converges iteratively.Further since the SUPS formulation requires addition of scaled residuals of the momentum and the continuity equations it retains consistency upon addition of these terms.

Stability Analysis
Steady-state regularized equations are written as; where a is a given velocity field, u is the unknown velocity, p is the pressure, ν the viscosity, and ε(u) is the symmetric part of the velocity gradient and f is the body force.The initial value of the algorithm follows by replacing the regularization of the incompressiblity with a one-step process as Assuming the given data a(x) the linearised velocity field and ν(x) viscosity satisfy; The finite element formulation we wish to study can be written as: find u h ∈ V h and p h ∈ P h such that B(u h , p h ; w, q) = F(w, q), (w, q) ∈ V h × P h (19) with We obtain: We note that the stabilization parameter is bounded in each domain K by a constant.Be it defined; where we define Re K (x) as and Obtaining a common bound on τ(Re K (x)), with no restrictions on Re we arrive at where Here, m k and C k are two positive constants with m k defined above.
Lemma I: For stability of the formulation we require where C ′′ K is a positive constant independent of h K .Proof: From integration by parts ((∇u)a, u) = 0, u ∈ V h (28) since u = 0 on Γ for u ∈ V h .Therefore we have: Thus we only need a measure of: Applying Schwarz inequality: Utilizing the bound on the stabilization parameter and value of m k eq.26: Inverse estimates provides us element level bounds on the derivatives as follows (Ern and Guermond, 2004) In addition: Let us denote a positive constant β such that Combining all the above: Finally we obtain: Thus we obtain the following Restrictions due to the stability analysis of the formulation require Since the regularization steps present the same problem with previous values of the pressures stability is further guaranteed across the regularization steps.

Stability Analysis of AGLS Formulation
We perform the stability analysis of the Galerkin Least Squares (GLS) formulation next.We consider the steady form of the equations.Let the weighting spaces which we wish to work with be given by (as defined above) where the integers k and l are greater than or equal to 1.

Consistency
Momentarily reverting back the notation introduced earlier let us denote (u, p) the solution of equation 1-2 and (u ε , p ε ) the solution to 6-7.It has been shown the following error estimates hold for the weak finite element formulation of the iterative penalty method (Xiao-liang and Shaikh, 2006) where the constant C 0 is independent of h and ε.For the above estimates we define H s (Ω) to be the standard Sobolev spaces with corresponding inner products (., .)s,Ω and norm ∥.∥ s,Ω .The formulation is consistent since the pressure converges iteratively.Further since the GLS formulation requires addition of scaled residuals of the momentum and the continuity equations it retains consistency upon addition of these terms.The following convergence estimates were proved in (Lin, 1997).
and the iteratively regularized pressure as From the above equation we can see we only need to assume the penalty parameter ε < 1 but not necessarily choose it very small since by iteration the error O(ε k ) can reach any accuracy that we desire.Further since we choose ε to be not very small the formulation is less stiff and more stable.

Stability Analysis
Steady-state linearised equations are written as; where a is a given velocity field, u is the unknown velocity, p is the pressure, ν the viscosity, and ε(u) is the symmetric part of the velocity gradient and f is the body force vector.Assuming the given data a(x) the linearised velocity field and ν(x) viscosity satisfy; The finite element formulation we wish to study can be written as: find u h ∈ V h and p h ∈ P h such that In the above equation 51 the parameter δ is the least squares stabilization parameter on the incompressibility constraint.We obtain: We note that the stabilization parameter is bounded in each domain K by a constant.Be it defined; Obtaining a common bound on τ(Re K (x)), with no restrictions on Re K we obtain where Lemma II: For stability of the formulation we require where C′ and C′′ are positive constants independent of h K .
Proof: From integration by parts since u = 0 on Γ for u ∈ V h .Therefore we have: The above statement proves the stability of the AGLS formulation.To derive estimates from this analysis we further examine the statement.
Classical definition of the stabilization parameter assumes a relationship of the form Eq. 22.Additional terms in the parameter dependent on the constant viscosity of the medium are traditionally defined as in (Ranjan et al, 2016a,b) τ = τ(h, ν, a).In the ensuing discussion we define the stabilization parameter to be the maximum value of this parameter obtained across all elements K and assume a constant viscosity in the domain.Let us denote this value of the stabilization parameter as τ sup .For this constant value of the parameter we obtain: Thus we only need a measure of: Applying Cauchy-Schwarz (utilizing positivity of ν): Utilizing the bound on the stabilization parameter and value of m k eq.55: Inverse estimates (on element level integrals) provides us element level bounds on the derivatives as follows (Ern and Guermond, 2004) Let us denote by the positive constant C k3 = C k1 C k2 Therefore it follows: In addition we have: Combining all the above we get: This leads to: Thus we obtain the following Further from Poincare inequality we obtain Therefore we obtain Assimilating common terms we obtain where we define the constant C k5 = 1 C k4 .For the above inequality to satisfy eq.56 we require The restrictions on viscosity from the above equations are obtained as (maximum value of viscosity) The second term from the stability estimate obtains 2ν Note, in the derivation of the above restriction on the stabilizaton parameter we did not use the definition of the parameter given by Eq. 54.This is an important note, since the estimate provided by Eq. 74 follows directly from the stability analysis for the maximum value of the stabilization parameter inside the domain.

Spectral/hp Finite Element Formulation
A class of p-type elements known as the spectral elements use the Lagrange polynomial through the zeroes of the Gauss-Lobatto polynomials (Karniadakis and Sherwin, 1999;Ranjan, 2011;Ranjan and Reddy, 2009).The spectral finite element approximation is described as follows.The primary variables are each approximated as where ψ j are the nodal expansions, which are provided by the following one-dimensional Ĉ0 spectral nodal basis as mentioned in Karniadakis et al. (1999).
where ∆ j are the nodal values due to the Kronecker delta property of the spectral basis.L n = P n (0,0) is the Legendre polynomial of order n and ξ i denote the location of the roots of (ξ The element stiffness matrices (derived in the following section) were obtained based on the Gauss-Lobatto-Legendre (GLL) rules.Gauss-Lobatto-Legendre rule includes both end points of the interval, that is, ξ i = ±1.The points and weights are respectively listed as The weights for the GLL quadrature rule are obtained as follows The zeroes of the Legendre polynomial do not have an analytical form and are commonly tabulated or evaluated separately for input in the computational framework.First and second derivatives of the Legendre polynomials are required in the formulation.The first derivative differentiation matrix is provided as The differentiation operation utilizing the differentiation matrix is obtained The computation of higher order derivatives follows the approach for the computation of the first derivative.One can compute the entries of the q th order differentiation matrix d (q) i j by evaluating the q th derivative of the Legendre polynomial at the quadrature points.Another convenient formulae for the higher order derivative is For the purposes of the current discussion we will be only concerned with the first and second derivatives.
The construction of the two dimensional spectral basis follows a tensor product with nodal expansions in either direction.Thus the construction of the tensor product follows These functions are defined on the domain (master element) which will be transformed by affine mapping to each element Ω e .

Finite Element Formulation
Let the spectral/hp spectral element spaces which we wish to work with be given by where the integers k and l are greater than or equal to 1.These spaces are defined using nodal spectral/hp basis provided by Legendre polynomials.
We consider the dimensional form of the Navier-Stokes equations for the development of the ASUPS formulation.Consider a spectral/hp element discretization of Ω into subdomains, Ω e , k = 1, 2, ...., n el where n el is the number of spectral/hp elements into which the domain is divided.Based on this discretization, for velocity and pressure, we define the trial discrete function spaces ℑ hp u and ℑ hp p , and the weighting function spaces; W hp u and P hp p defined above.These spaces are selected, by taking the Dirichlet boundary conditions into account, as subspaces of ] n el and , where is the finite-dimensional function space over Ω.We write the stabilized Galerkin formulation of eq.6-7 as follows: Find ] As can be seen from Eqn. 86, four stabilizing terms are added to the standard Galerkin formulation.In Eqn.86 the first three terms and the right hand side constitute the classical Galerkin formulation of the problem.The surface integrals on the right are obtained from the weak form development as the natural boundary conditions.The first series of elementlevel integrals comprise the SUPG and PSPG stabilization terms, which are added to the variational formulation (Brooks and Hughes, 1982;Hughes et al, 1986).Two extra terms are added to the variational formulation.The new integral adds a term (ψ i ψ j ) to the stabilized finite element formulation.A similar term is added to the forcing function which depends on last known regularization step value of the pressure.The regularization step and the non-linear iteration can be handled as a single iteration in the new formulation.Prescriptions on determinations of the stabilization parameters based on inverse estimates have been utilized by Gervasio (Gervasio and Saleri, 1998).In their paper the higher order stabilization parameter (with applications to spectral methods) was defined. where The second stabilization term on the continuity equation is known as the least squares incompressibility constraint (LSIC) term (Tezduyar and Sathe, 2003).This term is provided as In this paper we resort to the stabilization parameters defined by Eq. 87-90.In the above equations N d denotes the number of degrees of freedom in each element.The constant m k is defined above.The p norm of the velocity is defined We used the notation p to denote pressure and N to denote the space dimensions earlier.To avoid any confusion with notation we denote the p norm with pn and N d as the number of degrees of freedom in each element.
The spatial discretization of Eq. 86 leads to a system of non-linear ordinary differential equations, Here, v is the vector of unknown nodal values of v hp , and p is the vector of nodal values of p hp .The matrices N(v), K, and G are derived, respectively, from the advective, viscous, and pressure terms.The vectors F and E are due to boundary conditions.The subscripts κ 1 and κ 2 identify the SUPG and PSPG contributions, respectively.The various matrices forming the discrete finite element equations are provided.The definition of the terms in the stiffness matrices (after linearization) are outlined below; ϵq hp p k,hp dΩ (99) And the forcing functions are defined as The stabilization of the Galerkin formulation presented in this paper is a generalization of the SUPG/PSPG/LSIC formulation, employed thus far for solving incompressible flow problems.In the equations above the parameters κ are defined as With such stabilization procedures described above, it is possible to use elements that have equal order interpolation functions for the velocity and pressure, and are otherwise unstable (Hughes et al, 1986).There are various ways to linearise the non-linear terms in stabilized finite element methodology.We store the tangent stiffness matrices.The least squares term was set to zero in the present case.
The Galerkin Least Squares (GLS) formulation adds a least squares term involving the momentum equation to the Galerkin finite element formulation.The GLS stabilization is a more general stabilization approach that includes the essence of the SUPG and PSPG type stabilizations as subcomponents of its formulation.This approach has been successfully applied to Stokes flows, compressible flows (Shakib et al, 1991), and incompressible flows at finite Reynolds numbers (Tezduyar, 1992).In the GLS stabilization of incompressible flows, the stabilizing terms added are obtained by minimizing the sum of the squared residual of the momentum equation integrated over each element domain.Consequently similar to the SUPG and PSPG stabilizations, because the stabilizing terms involve the residual of the momentum equation as a factor, the GLS stabilized formulation is consistent.For time-dependent problems, a strict implementation of the GLS stabilization technique necessitates finite element discretization in both space and time, and leads to a space-time coupled finite element formulation of the problem.A slightly liberal implementation of the GLS formulation requires discretization in space only.This formulation still retains the consistency of the governing differential equations since we seek the solution of the partial differential equation minus the time dependent term in the stabilization process.Spacetime formulations have been known to generate extensive linear systems with the accompanying lethargy in solution procedures for a fully coupled simulation.In the interest of saving computational efforts in obtaining the solutions for the Navier-Stokes equations we will employ a space-time decoupled GLS formulation.
To formalize the GLS formulation for incompressible flow we introduce some notation.Let the spatial domain Ω be divided into non-overlapping subdomains Ω e , where e=1, 2, 3, ..nel where nel is total number of elements in the domain.
Based on this discretization of the velocity and pressures we define the trial function spaces Γ hp u and Γ hp p and the weighting function spaces W hp u and W hp p .These function spaces are selected by taking the Dirichlet boundary conditions into account as subsets of the finite dimensional space H 1h (Ω).The GLS stabilized finite element formulation is written as follows: Find u ∈ Γ hp u and p ∈ Γ hp p such that for w ∈ W hp u and q ∈ W hp p ] The above formulation requires the specification of the stabilization parameter.The stabilization parameter is obtained from the stability analysis of the formulation performed earlier.Classically it has been derived from a multi-dimensional generalization of the parameter derived in Shakib (1989) in one dimension.In the following developments it is also specified based on the advection and diffusion limits.The classical definition of the stabilization parameter is provided as This classical stabilization parameter is designed to provide stability to the formulation.
We employ the decoupled space-time formulation in the interest of lower compute times.We seek to obtain solutions to the incompressible Navier-Stokes with a modification of the original GLS formulation augmented with an iterative penalization of the incompressibility constraint.This augmentation is found to provide better non linear convergence histories as compared to the classical GLS formulation used to date.The idea is similar to artificial compressibility methods which have been known to improve the convergence rates for stiff problems obtained from the discretization of Navier-Stokes equations.In the iteratively penalized formulation however we maintain consistency with a regularization of the incompressibility equation.The space-time decoupled GLS stabilized finite element formulation is written as follows: u ∈ Γ hp u and p ∈ Γ hp p such that for w ∈ V hp u and q ∈ ] Two extra terms are added to the variational formulation when compared to the GLS formulation.The new integral adds a term (ψ i ψ j ) to the stabilized finite element formulation.The regularization step and the non-linear iteration can be handled as a single iteration in the augmented formulation.Augmentation is achieved with the relaxation of the divergence free condition via iterative penalization.We use higher order spectral elements where GLS and SUPS formulations are distinct.
The spatial discretization of Eq. 106 leads to a set of non-linear ordinary differential equations, Here, u is the vector of unknown nodal values of u hp , and p is the vector of nodal values of p hp .The matrices N(u), K, and G are derived, respectively, from the advective, viscous, and pressure terms.The vectors F and E are due to boundary conditions and the iterative penalization procedure.The subscripts κ 1 , κ 2 , κ 3 identify the SUPG, PSPG and GLS contributions, respectively.The various matrices forming the discrete finite element equations are provided.The definition of the terms in the stiffness matrices are outlined below; (M(u) ϵq hp p k,hp dΩ (114) And the forcing functions are defined as The stabilization of the Galerkin formulation presented in this paper is the GLS formulation for solving incompressible flow problems.In the equations above the parameters κ are defined as With such stabilization procedures described above, it is possible to use elements that have equal order interpolation functions for the velocity and pressure, and are otherwise unstable (Hughes et al, 1986).We utilize the GLS formulation for presenting results to benchmark incompressible flow problems in a following section.

Non Linear Solvers and Time Stepping Algorithms
Let us consider the discrete non-linear problem obtained from the Navier-Stokes equations by; where F(∆) : R N → R N is a non-linear mapping with the following properties: (1) There exists an (2) F is continuously differentiable in a neighborhood of ∆ * .
In the particular case of solution of incompressible flow problems, ∆ = {u, p} denotes the vector of nodal variables for the velocity field and pressures over the domain of interest.For moderate to high Reynolds numbers the non-linear convective term dominates and the solution procedures have to be able to accommodate the strong non-linearities.
We assume that F(∆) is continuously differentiable in ℜ N where N is the number of spatial dimensions.Let us denote the Jacobian matrix by J ∈ ℜ N .Classical Newton's method for solving the non-linear problem can be enunciated as (Elias et al, 2004(Elias et al, , 2006) Newton's method was designed for solving problems where the initial guess is very close to the non-linear solution ∆ * .Amongst the disadvantages of Newton's method are the requirements that the linear system be solved very accurately which can be computationally expensive.Computing an exact solution using a direct method can also be prohibitively expensive if the number of unknowns is large and may not be justified when ∆ k is far from the solution.In spectral computations where there is high computational costs involved with quadratures it is imperative to resort to efficient procedures for solving the discrete systems of equations.Thus one might prefer to compute some approximate solution leading to the following Inexact Newton Algorithm or a damped version of the final update to the new solution vector as; set ∆ k+1 = ∆ k + αs k where, α is the damping parameter.For some η k where ∥ • ∥ is the norm of choice.This new formulation automatically allows the use of an iterative linear solver method, one first chooses η k and then applies the iterative solver to the algorithm until an s k is determined for which the residual norm satisfies the convergence criterion.In this context, η k is often called the forcing term, since its role is to force the residual of the above equation to be sufficiently small.In general the non-linear forcing sequence needs to be specified to drive the solution towards the solution of the non-linear problem as; The ∥r k ∥ is defined as the residual for the non-linear problem.In our implementations we use the element by element (EBE) BiCGSTAB method to compute the s k .In the EBE-BiCGSTAB implementation we store the element stiffness matrices.The amount of storage required for the problem varied as O(e p × Here, e p is the number of elements stored in each processor, nd f is the number of degrees of freedom in each element, and p is the polynomial approximation of the spectral element.There exist several choices for the evaluation of the forcing function ∥η k ∥ as mentioned in Eisenstat and Walker (1996).Amongst the choices is the usage of a constant value of η k = 10 −4 .We stored the element tangent stiffness matrix in our computational framework.
The derivations of the convergence behavior for a generic non-linear problem have been outlined in Dembo et al. (1982).The basic requirement for the convergence of the iterative method is that F is continuously differentiable in the neighbourhood of ∆ * ∈ R N for which F(∆ * ) is non singular and that F ′ is Lipschitz continuous at x * with constant λ, i.e.
The construction of the tangent stiffness matrix required for the algorithm follows standard procedures in finite element analysis.In our implementations we store the tangent element stiffness matrices.
We employ the generalized α time integration method to solve the discretized equations.The generalized α method for Navier-Stokes equations for incompressible flows was first proposed by Jansen et al. (2000).In this method the stepping in time for the variables are approximated as In the above equation r denotes the time-step and ∆ is the vector of unknowns {u, p}.The discrete momentum and continuity equations are collocated at intermediate values of the solution at every time step.The relationship between the nodal velocity degrees-of-freedom and their time derivatives is approximated by the discrete Newmark formula Here α m , α f , and γ are the real valued parameters chosen based on the stability and accuracy considerations.It was shown by Jansen et al. (2000) that second order accuracy in time is achieved provided that while unconditional stability is achieved provided that To study the unsteady problem for the flow past a square obstruction we employ a space-time decoupled formulation.We begin by partitioning the time interval of interest [0, T ] into subintervals, or time steps.The time intervals t n and t n+1 define the end points of the time n th time step and ∆t n = t n+1 − t n is the size of the n th time step.In general the time step may vary from step to step.A one-parameter family of second-order accurate and unconditionally stable algorithms may be obtained by setting γ according to equation 125 and employing the following parametrization of the intermediate time steps

Numerical Examples
In this section we obtain solutions for some computational fluid dynamics problems with the formulations presented.First we consider the ASUPS formulation, and subsequent examples consider the AGLS formulation.In the first section we examine steady state conditions for the Kovasznay flow problem and skewed cavity problem.We provide agreement with analytical results for verification for the Kovaszany flow solutions.Steady state analysis results are presented for the flow past a square obstruction in a enclosure and unsteady state flow past a square cylinder.With the AGLS formulation we consider the performance for the driven cavity problem at aspect ratio of 2, and flow past two cylinders in tandem.

ASUPS Formulation
We consider first standard incompressible flow problems solved with the ASUPS formulation and follow with numerical examples solved with the AGLS formulation.

Kovasznay Flow
Benchmark problems in two dimensional computational fluid dynamics involve the solutions of Kovasznay flow among other problems.Analytical solutions have been provided for this flow to compare against numerical results by Kovasznay (1948).We study this standard two dimensional benchmark problem at Reynolds number of 40.
where λ = Re/2 − (Re 2 /4 − 4π 2 ) 1/2 and p 0 is the reference pressure (an arbitrary constant).Figure 1 presents the solution of the problem at Re 40.The two dimensional mesh along with contour plots of u-velocity, pressure, and v-velocity, are presented in the figure (clockwise from left).We compare the sup norm between the analytical solutions and the numerical results.The sup norm was defined as The values of this L ∞ norm was found to be 0.024170 and 0.002908 for the two components of the velocities and 0.032211 for the pressure.The non-linear tolerance for the present runs were kept at a reduction in s k to a value less than 10 −6 for generating the results in Figure 1.It should be noted that for an 'exact solve' for the problem solved with Picard method these numbers were 0.021667, 0.0028023, and 0.032027.For this comparison the p-level was set to 3. The Picard method solutions were obtained with iteration lagged linearisation.A Gauss elimination solver was used for obtaining the results for the Picard solution.Convergence study of the above formulation with increasing p-levels is provided in Figure 2. Obviously, the constant value of ϵ = 0.01 affects the convergence with p-level beyond five (5).Further studies to obtain spectral convergence of the formulation (wrt pressures) with accurate design of ϵ are desired.For the problem considered the linear convergence tolerance was set at a reduction in the residual of 10 −4 and the non linear tolerance was set at 10 −9 .
For a p-level of three (3) the inexact algorithm took from 55-95 iterations to converge for the linearised discrete systems.

Skewed Cavity
Lid driven cavity problem for orthogonal grids has been studied extensively by researchers in the past (Ghia et al, 1982).
The cavity problem for non-orthogonal grids has not been explored that extensively perhaps because of the oblique cavity orientation.In this section we test the spectral/hp Newton Krylov algorithm for solving the skewed cavity problem at angles of 30 • and 45 • respectively.We also provide the convergence history for both the SUPS and ASUPS formulations.
The length of the cavity was taken as unity, with the density of the fluid as 1, and lid velocity was taken as U L = 1 respectively.The Reynolds number was defined based on the lid velocity U L and cavity length L. The viscosity was varied from 0.01 through 0.001 which obtained Reynolds numbers 100 and 1000 respectively.The angle of skew was defined as the inclination of the side walls of the cavity.Two different values of 30 • and 45 • were used.The reason for choosing the problem included an intricate study of the difference between different types of errors when performing numerical analysis.As mentioned in Demirdžić et al. (2005) who provide benchmark results for this problem, there are three sources of error in the study of any problem namely; (a) discretization errors: the difference between the exact solution of the correctly discretized equations on a given grid and the exact solution of the differential equation, (b) convergence errors: the difference between the exact and the approximate solution of the discretized equations left over after stopping the iterations.(c) algorithmic and programming errors: errors due to incorrect discretization, erroneous implementation of boundary conditions and programming errors etc.We assume only the first two sources of errors to dominate the numerical solutions.Demirdžić (Demirdžić et al, 2005) used the control volume approach for solving the equations on a 320 × 320 mesh.
In comparison we only utilized a 40 × 40 mesh with a p level of 6 in each element.Thus, we retain the high order spatial resolution offered by the spectral/hp methods at the same time we can resort to the non-linear solvers to reach benchmark solutions.The non-linear tolerance was set to 10 −6 and the linear system solver tolerance was set to 10 −4 .We set the value of the regularization constant to ϵ = 0.0 and 0.04 respectively.Figure 3 presents the streamline and the v contour plots over the domain of the skewed cavity for Reynolds number of 100 with a skew angle of 30 • .The corresponding  Present results benchmarked with the skewed vertical and horizontal centerlines are presented in Figure 7 for Re-100 and skew angle of 45 • .The non-linear convergence history for the SUPS formulation and the ASUPS formulation have been presented in Figure 8 and Figure 9 for the velocity and pressures respectively.We can discern from the figures that the ASUPS formulation performs better for this problem and converges more rapidly to the non-linear tolerance.In general excellent agreement was found for different skew angles and also differing Reynolds numbers for this problem.Further, the length of the downstream reattachment point was found to be of larger length for the moving wall boundary condition.Figure 10 presents the development of the u and v contour plots for the stationary top wall.The development of the pressures and streamlines over the domain for this problem (BCI) are presented in Figure 11.As can be seen from Figure 11 the length of the downstream reattachment point was found to be of length 3.5 units.This is in agreement with the steady state solution of the problem solved in Laval (Laval and Quartapelle, 1990).It should be mentioned that in the above reference the upstream separation point and upstream recirculation region was not perceptible.However, due to the usage of spectral approximations we obtain more accurate resolution of the flow field in present simulation with perceptible development of a secondary circulation region upstream.The problem required from 1460 to 5452 iterations for reaching the linear convergence tolerance specified as 10 −4 for the discrete linearised system generated solved with EBE-BJCGSTAB.The inexact Newton Krylov solved required from between 3-4 iterations.In the problem subject to different boundary conditions it was realized that most of the flow fields of interest occur at the sections of x = 3 and x = 6 respectively.Figure 12 presents the U component of the velocities at the sections of x = 3 and x = 6 respectively for the stationary top wall with Reynolds number of 100. Figure 13 presents the contour plot of the U and V component of the velocity of the fluid when the top wall is set in motion with an X-component of the velocity 1.The development of the pressure and the streamline for this second boundary condition (BCII) is presented in Figure 14.From Figure 14 the length of the reattachment point downstream was found to be 4.0 units.Figure 15 presents the U component of the velocity at the section x = 3 and x = 6 respectively for BCII with the top wall in motion.The boundary conditions includes specification of the free stream velocity at the inlet of the channel of u = 1, and v = 0.At the top and bottom boundaries free stream velocity is imposed with a specification of u = 1, and v = 0. We consider the problem for a Reynolds number of Re-100 and set the regularization constant ϵ = 0.03.In general we have found the iterative regularization constant to lie in the range of 0.01-0.05.This is also in conformity with recommendations for this constant found in literature (Dai and Cheng, 2008;Lin et al, 2004;Xiao-liang and Shaikh, 2006).
The outflow boundary condition include specification of traction free surfaces t x = 0 and t y = 0.The specification of pressure p = 0 completes the problem definition by anchoring the pressure for the simulations.The problem was discretized into a macro-mesh with 50 spectral elements along the x-and y-directions respectively.The connected model Ω hp comprised of spectral elements with a uniform p level of 4 in each element.Previous numerical studies reported To address this issue we refine the mesh in the near vicinity of the square cylinder.To study the unsteady nature of the problem we employ a space-time decoupled finite element formulation.The temporal terms in the governing equations are represented by the generalized α family of approximation, which retain second order accuracy in time and allow for user controlled high-frequency damping by the single free integration parameter, 0.0 ≤ ρ hp ∞ ≤ 1.For ρ hp ∞ = 1 the method   1998) who reported a value of S t = 0.139 and also Pavlov et al. (2000) who reported a value of 0.137.For checking the predictive capabilities of the formulation we verify the conservation of mass at the exit plane.The mass conservation obtained at the exit was divided by the inlet mass to obtain the ratio m out /m in .It is well known that most finite element formulations have issues with mass conservation.From figure 20 we notice the mass conservation is also periodic in nature with a very small amplitude.The mean value of the exit mass of the fluid vs. the inlet mass was found to be 0.998.This slight difference between the input and output mass can be also attributed to the inexact solve of the non-linear problem.It is evident there is negligible errors in conservation of mass which further verifies the predictive capabilities of the ASUPS formulation.

Driven Cavity Aspect Ratio AR=2
Let us consider a benchmark problem in two dimensional incompressible flow computations, the driven cavity problem with an aspect ratio of 2. The specification of the problem is standard with the top surface of the cavity set in motion with an initial velocity of u = 1 and v = 0.No slip and no-penetration velocity boundary conditions are specified on the rest of the walls.A specification of the pressure to zero on the mid-bottom of the cavity completes the description of the problem.We consider the problem with an aspect ratio of AS=2 at a Reynolds number of 100.The dimensions of the cavity were taken as a unit length along the X direction and 2 units along Y. Reynolds number was based on the length of the cavity along the x-direction and the velocity of the driven surface.The problem was discretized into a macro-mesh of 50 × 60 elements with a p level of 3 in each element.A graded mesh was considered for analysis with a gradation towards the edges for capturing the boundary layers.The problem was subject to steady state analysis.Figure 21   Flow past two cylinders in tandem configuration is of interest to examine the effects of the cylinder interactions with each other and the accompanying flow field.Practical applications of this flow are found in off-shore structures, transmission cables with twin conductors (Mittal et al, 1997) etc.There are two configurations of interest while studying such solid fluid interaction problems.The first case is both the cylinders are placed in a tandem configuration, and the second is both cylinders placed in transverse arrangement.Different numerical and experimental studies have involved solving the flow past two cylinders in transverse configuration (Williamson, 1985).However flow past two cylinders inline with each other (in tandem) is relatively scarce.This problem has been examined in Mittal (Mittal et al, 1997) in detail.They considered the problem with a belt type boundary condition.Initial steady state flow field was perturbed with a clockwise and counter clockwise motion of the cylinders and the flow was allowed to develop.
In the following example we consider flow past two cylinders in tandem configuration without any initial perturbation of the flow field.This problem was carefully chosen as a sample problem to examine the robustness of the formulation for long time integration with the AGLS formulation with previously published results.We examine the flow interactions of the cylinders with each other and adjacent flow.The intent with an unperturbed initial condition was to be able to recover the unsteady solution obtained earlier for a similar problem.The diameters of the cylinders are denoted by D. The distance between the centers of both cylinders is denoted by S .It has been observed for cylinders in tandem configuration there exists a critical cylinder distance S between the cylinders below which there is no distinct shedding observed behind the first cylinder.Even at this separation however there is formation of a vortex street behind the second cylinder.This distance has been reported as S =3D.We present the flow metrics observed for the cylinder separation distances of 5.0D.For smaller distances lesser than 3D the flow perturbations at Reynolds number of 100 tend to reach steady state.
The problem domain was considered to be of length [32,32].The center of the first cylinder was placed at a separation of 5 units from the upstream boundary at the location [5,16].The second cylinder was placed at [10,16].The diameters of both cylinders were taken as 1.Free stream values of the velocity were specified on the top, bottom, and left faces of the computational domain.The exit face was specified traction free.The Reynolds number based on the diameter of the cylinders and the free stream value of the velocity was considered as Re-100.Figure 23 presents the contour plot of the velocity components at a specific time instant for the two cylinders in tandem.From the animation generated from the simulation the cylinders were found to be in anti-phase mode of vortex shedding.Contour plot of the pressure field has been presented in Figure 24.The obstruction in the flow downstream of the first cylinder prohibits an early development of the shedding phenomenon.For a single cylinder in free stream the Von-Karman vortex street is fully developed by the time instant t-50.However for the present case the development starts at around t-90.The values of the amplitude of the lift coefficient (Figure 25) and the Strouhal number for the first cylinder were found be close to the values for a single cylinder.The second cylinder which lies in the wake of the first cylinder was found to be significantly affected due to the presence of the first cylinder.The amplitudes of the lift coefficients for both cylinders were found be in agreement with the values presented for an initially perturbed flow condition by Mittal et al. (1997).

Conclusion
Augmented SUPS and Augmented Galerkin Least Squares formulations with iterative penalization of the incompressibility constraint were proposed for solving incompressible Navier-Stokes equations.Stability analysis was performed and better non-linear convergence history of the ASUPS formulation was demonstrated.The formulations were used for solving both steady state and unsteady state incompressible flow problems.Higher order/spectral methods provided the flexibility of obtaining exponentially accurate results with high p levels of resolutions on relatively coarse macro-meshes.We utilize inexact Newton Krylov methods for solving problems.No deterioration of results are observed for flow problems in confined enclosures with different boundary conditions.Excellent conservation of mass at the exit plane further verify the predictive capabilities of the formulations presented and the robustness of the solution algorithms.

Figure 10 .
Figure 10.Re-100 U and V contour plots with No Slip.ps

Figure 11 .
Figure 11.Re-100 Pressure and Streamlines with No Slip.ps

Figure 13 .
Figure 13.Re-100 U and V contour plots Moving Wall

Figure 16
Figure 16 presents the U-component of the velocity for the flow past a square shaped cylinder once the unsteady nature of the flow is fully developed.The development of the Von-Karman vortex street is evident.The V-component of the velocity is shown in Figure 17.The time history of the lift coefficient has been presented in Figure 18.The lift coefficient exhibits a fluctuating component with the amplitude of 0.20.The transient evolution of the U, V-velocity component with

Figure 17
Figure 17.V-contour plots flow past a square Re-100

Figure 20 .
Figure 20.Mass conservation exit plane presents the contour plot of the velocity and pressures over the domain for the problem.From the streamline plots it is evident two circulation zones develop inside the cavity.The problem was solved for two different Reynolds numbers of 100 and 400.The agreement of the centerline velocities at the mid-centerline are presented in Figure22.The agreement of the present results with benchmark results ofGupta et al. (2005) is excellent.

Figure 22 .
Figure 22. x and y component of velocities with benchmark results for Re − 100 and 400