This work explores the use of stabilized finite element formulations for the incompressible NavierStokes equations to simulate turbulent flow problems. Turbulence is a challenging problem due to its complex and dynamic nature and its simulation if further complicated by the fact that it involves fluid motions at vastly different length and time scales, requiring fine meshes and long simulation times. A solution to this issue is turbulence modeling, in which only the large scale part of the solution is retained and the effect of smaller turbulent motions is represented by a model, which is generally dissipative in nature.
In the context of finite element simulations for fluids, a second problem is the apparition of numerical instabilities. These can be avoided by the use of stabilized formulations, in which the problem is modified to ensure that it has a stable solution. Since stabilization methods typically introduce numerical dissipation, the relation between numerical and physical dissipation plays a crucial role in the accuracy of turbulent flow simulations. We investigate this issue by studying the behavior of stabilized finite element formulations based on the Variational Multiscale framework and on Finite Calculus, analyzing the results they provide for wellknown turbulent problems, with the final goal of obtaining a method that both ensures numerical stability and introduces physically correct turbulent dissipation.
Given that, even with the use of turbulence models, turbulent flow problems require significant computational resources, we also focused on programming and parallel implementation aspects of finite element codes, and in particular in ensuring that our solver can perform efficiently on distributed memory architectures and highperformance computing clusters.
Finally, we have developed an adaptive mesh refinement technique to improve the quality of unstructured tetrahedral meshes, again with the goal of enabling the simulation of large turbulent flow problems. This technique combines an error estimator based on Variational Multiscale principles with a simple refinement procedure designed to work in a distributed memory context and we have applied it to the simulation of both turbulent and nonNewtonian flow problems.
Turbulent flows are a subject of paramount interest in engineering, physics and earth sciences, since they govern many phenomena involving moving fluids like volumes of air or water. Turbulent phenomena appear in flows of all sizes, from large scale processes such as the atmospheric flow of air or oceanic currents to local problems such as wind loads over lightweight structures or the efficiency of a wind turbine. However, the complex nature of turbulent flows means that they can rarely be solved analytically and, as a result, we have to rely on experimental studies or numerical simulations for their analysis.
Experimental studies, which typically involve either campaigns of field measurements or wind tunnel tests, tend to be expensive and complex endeavors (see for example [1]) and, although they are an essential tool in providing a deeper understanding of the problem, tend to provide localized measurements such as point recordings or line averages, which are limited to the regions where sensors are placed, or indirect measurements, such as trajectories of tracers. From this point of view, numerical modeling represents a less costly alternative to experiments with the added advantage of providing a global image of the velocity and pressure fields for the entire volume of interest. Unfortunately, numerical study of turbulent problems is not without difficulty. The fundamental problem in turbulent flow simulations is the wide range of time and spatial scales involved [2], which have to be taken into account in the analysis.
Taking as an example the flow around an obstacle, the largest fluid motions, or eddies, in the flow can be assumed to have a characteristic size comparable to the size of the obstacle itself. Under a turbulent flow regime, coherent flow structures tend to break and the energy associated to these motions is transferred to successively smaller eddies until, for very localized fluctuations, viscous effects dominate and the energy is dissipated. This process can be formalized in terms of the energy associated to motions of a characteristic wavelength in what is known as the Kolmogorov energy cascade [3,4,2] and the size at which viscous dissipation dominates is defined as the Kolmogorov length .
Turbulent flows can be characterized by the Reynolds number , a dimensionless parameter indicative of the balance between inertial and viscous forces in the problem,

(1) 
where is the fluid's density and its dynamic viscosity and is a characteristic velocity of the flow. The Kolmogorov length scale can be related to the largest motions of the flow by means of the Reynolds number [2] as

(2) 
To properly represent the full range of turbulent motions in a numerical simulation, the size of the problem domain would be a multiple of the size of the obstacle on each dimension, while the grid size would have to be small enough to capture motions on the Kolmogorov scale. For a problem, this implies that the total number of elements would be on the order of . Keeping in mind that the Reynolds number in engineering applications is typically on the order of , the number of elements involved in a typical simulation would put it beyond the reach of most computers and scientific computation clusters. As a result, this approach, which is known as Direct Numerical Simulation (DNS) is only applicable to problems with a moderate Reynolds number and is not commonly used in engineering.
The only viable option in most cases is to introduce a turbulence model to reduce the computational requirements of the simulation. This typically involves neglecting the smallest motions and modifying the problem to include terms that model the effect of the neglected motions on the problem. The addition of such model increases the minimum element size required to simulate the problem, extending the range of problems that can be simulated with given computational resources. Turbulence models can be grouped into two broad families, Reynolds Averaged NavierStokes (RANS) and Large Eddy Simulation (LES).
RANS models are based on rewriting the problem in terms of timeaveraged variables (or time and spaceaveraged, if there are spatial directions of homogeneity) and introducing a model for the effect of turbulent motions on the averaged solution. This model has a dissipative effect, removing energy from the average motions, and most commonly takes the form of an added viscosity, called turbulent or eddy viscosity. The temporal and spatial distribution of the eddy viscosity is given by the specific model and is usually determined by the solution of one or more additional equations that describe the production and transport of turbulent quantities. Typical examples of RANS models are the SpalartAllmaras model [5], which involves one additional equation, or the – [6,7] and – models [8], which introduce two additional variables.
LES models are based on introducing a spatial filter to smooth out high wave number fluctuations in the solution and writing the problem equations in terms of filtered quantities. This allows the separation of larger flow features, influenced by the geometry of the fluid domain and the boundary conditions, from the smaller motions, which are assumed to have an universal behavior and thus can be replaced by local dissipation model, which is problemindependent. The most well known member of this family of models is the Smagorinsky method [9], but many variants and alternatives exist [2].
Note that, while the use of turbulence models makes the simulation of turbulent flows possible in terms of mesh resolution and computational resources, many flows of practical interest still require fine meshes and long simulation times in order to obtain statistically steady solutions and reliable statistical measurements. As a result, the simulation of complex flows at high Reynolds number still requires significant computational resources, frequently beyond what a single computer can provide, and has to be solved using High Performance Computing (HPC) clusters.
A second problem frequently encountered in the finite element simulation of incompressible flows is the apparition of numerical instabilities, which are a consequence of the incompressibility constraint and the effect of the convective term in the equations for convectiondominated flows (see for example [10] for an introduction to this topic). In the context of finite elements, one possible solution for this issue is the use of stabilization techniques, where the original Galerkin weak form of the problem is modified to obtain a stable formulation. The modified equations are characterized by the addition of new terms that are typically dissipative in nature, in a way that ensures consistency, that is, that the modified problem tends to the original equations as the simulation mesh is refined.
For turbulent flows in particular, the interaction between turbulence models and stabilization terms is an active area of research. While they have clearly different origins and motivations (turbulence models are typically motivated using physical arguments, while stabilization methods are purely numerical in nature), both types of methods have a dissipative effect on the solution, which has raised some questions regarding their interaction or the possibility of using a unified formulation to provide both stability and turbulence modeling. This subject has been studied for example for StreamlineUpwind PetrovGalerkin (SUPG) stabilization [11,12], Finite Calculus (FIC) based formulations [13,14] and stabilization techniques within the Variational Multiscale (VMS) framework.
In the case of VMS formulations, the stabilized problem is motivated by a separation of scales, differentiating a large scale part of the solution, that can be represented by the finite element mesh, from a small scale part, which remains unresolved, by means of a projection of the solution onto the finite element mesh [15,16]. This concept has intriguing parallels with the spatial filtering introduced in LES methods and can be understood as a meshinduced filtering. This has motivated multiple investigations studying the use VMS formulations in turbulence modeling [hughes2000large,hughes2001multiscale,gravemeier2006variational,codina2007time,], analyzing their relationship both from the theoretical point of view and through its application to the simulation of turbulent flows of interest.
The main topic of this monograph is the study of the applicability of stabilized finite element techniques to the solution of turbulent flow problems of practical interest in engineering, with a special focus in analyzing the behavior of different formulations as a LESlike turbulence model. The research is organized around two main lines: The first of them is centered on stabilized finite element formulations and their relationship with turbulence modeling, while the second focuses on computational techniques to facilitate the calculation of complex flows in large simulation domains.
Regarding the relationship between stabilized formulations and turbulence modeling, we choose to focus on two types of methods. The first of them is the VMS framework, where we intend to analyze the behavior of two wellknown formulations, algebraic subgrid scales and orthogonal subgrid scales, as well as dynamic subscales, while the second is a new formulation derived from the application of the FIC balance to the incompressible NavierStokes equations. We have implemented a fluid solver based on the different techniques and used it to analyze their performance.
The first topic that will be discussed in terms of the computational aspects of the present work is the use of parallel programming techniques for the simulation of finite element problems in a distributed memory context. In this sense, all developments have been made with a parallel implementation in mind, choosing algorithms and implementations that are suited to a parallel implementation in preference to alternatives that are not, and the parallel capabilities of the implemented code have been evaluated.
A related topic in relation to the calculation of large problems is the use of adaptive mesh refinement techniques to simplify the mesh generation procedure and reduce the overall number of elements required to perform the simulation. We have developed a technique based on a simple refinement technique that can work in a distributed memory environment, which has been implemented to work in combination with the fluid solver. During the course of the present work we also found the opportunity to apply these refinement techniques to nonNewtonian flow problems. The results of our investigation in this area will also be presented.
The methods presented in this document have been implemented within the open source Kratos Multiphysics finite element framework [21,22], which is based on C++ and Python and has parallel calculation capabilities, some of which have been developed or expanded as part of the present work. This implementation has been used to simulate all the examples considered and generate the presented results.
The present document is organized along the main goals outlined above. In Chapter 2 we introduce VMS methods, with a special emphasis on dynamic subscale approximations, and review the arguments that have been used in the literature to relate them to LES formulations. Our tests have driven us to propose a variant of the method characterized by a new model for the pressure subscale. We use both the standard approach and our modified formulation to simulate a turbulent channel flow benchmark problem, which allows us to compare our results to DNS data.
In Chapter 3 we present a new stabilized formulation for incompressible flows based on the FIC balance principles. This formulation is also used to simulate the channel flow problem, as well as other benchmark examples.
Chapter 4 is concerned with different aspects of the parallel implementation of the finite element flow solver that has been used to perform the simulations presented in the previous chapters, including tests of the parallel performance of the solver. In the same chapter we present the approach we used to record spatial averages and variances in a distributed memory context.
Chapter 5 presents an adaptive mesh refinement technique that combines a parallel refinement algorithm based on edge subdivision with an error indicator motivated by the VMS formulation. This technique is then used to solve both turbulent and nonNewtonian flow cases.
Finally, in Chapter 6 we summarize the work and present its main conclusions, as well as proposing future lines of research.
Variational multiscale (VMS) methods [15,16] provide a theoretical framework for the design of stabilized finite element formulations based on the separation of the solution into resolved and unresolved parts, which is achieved through the definition of large scale and small scale solution spaces. The projection of the original equations onto the large scale space gives an equivalent problem that depends on the small scale variables, while the projection of the original equations onto the small scale space is used to motivate a model for the effect of the small scale variables, which are not solved, to the large scale solution.
This methodology has interesting parallels with LES turbulence models, which use a spatial filter to introduce a separation between the resolved and unresolved parts of the solution. This fact has motivated research into the relationship between VMS stabilized formulations and LES methods, and in particular on the possible use of VMS methods as turbulence models. In the present chapter we review some current trends in VMS formulations, including the derivation of algebraic and orthogonal models for the small scales, and the arguments that have been used to relate them to turbulence modeling. We also direct our attention to dynamic subscale models [23,20], which have been shown to provide good results in turbulent flow simulations without requiring the use of an explicit turbulence model [24,25,26].
In spite of the success of the method in the mentioned tests, some aspects of the behavior of VMS methods as turbulence models are not completely understood. In particular, the solution of the problem has a degree of dependency on the precise definition of the stabilization parameters (see for example [25]) and on using or neglecting the pressure small scale [26]. We investigate this issue through the analysis of the results of turbulent channel flow simulations. Additionally, we propose a new model for the pressure small scale, which we found to provide more accurate results in our simulations.
The chapter is organized as follows: in the next pages we review the state of the art for the formulations of interest, introducing the Galerkin weak form of the incompressible NavierStokes equations Section 2.2 and using it Section 2.3 to obtain the complete stabilized equations for the different variants that will be considered in this work. We continue by presenting the arguments that have been used to justify their use as a turbulence model in Section 2.4. In Section 2.5 we present the finite element solver that we have implemented based on the described VMS formulations. In Section 2.6 we present our proposal for a new model of pressure small scale and in Section 2.7 we use both this and standard VMS methods to simulate a turbulent channel flow test case, analyzing the solution in terms of different measured statistics of the flow and their comparison with reference DNS data. Finally, we present our conclusions in Section 2.8.
The incompressible NavierStokes equations state the balance of linear momentum and mass in a fluid domain , given by

where is the fluid velocity, its density, represents the stress tensor and the external forces acting on the domain.
The problem described by Eqs. 3 and 4 must be completed with suitable initial and boundary conditions. We denote the domain boundary as and introduce its partition into Dirichlet () and Neumann () parts, verifying and . The initial and boundary conditions for the problem can be expressed as:

where is the initial velocity field, represents the imposed velocity on the Dirichlet boundary, the outer normal vector and the imposed traction acting along the Neumann boundary. Note that the initial velocity must be chosen to be divergence free to ensure that Eq. 4 is verified at all times.
For Newtonian fluids, the stress tensor can be related to the fluid velocity and pressure using the constitutive relation

(8) 
where is the second order identity tensor, the fluid's viscosity and the symmetric gradient of velocity, defined as

(9) 
To obtain a finite element formulation for the NavierStokes equations we need to rewrite them in weak form. We multiply Eq. 3 by a velocity test function , defined to be zero on the Dirichlet boundary , and Eq. 4 by a pressure test function . Integrating the resulting expressions over the fluid domain we obtain

The Neumann boundary condition can be introduced in the formulation by using the product rule on the stress term in Eq. 10 and expressing the resulting divergence as a boundary integral, which allows us to write

(12) 
Substituting Eq. 12 into Eq. 10 and introducing the definition of the stress tensor in Eq. 8 we obtain

Eqs. 14 and 15, together with the initial and Dirichlet boundary conditions, allow us to state the weak form of the problem. In addition, they also impose regularity requirements on the solution, test functions and problem data, as it must be ensured that all the integrals that appear in the equations remain bounded. This is only briefly discussed here, directing the reader to more specific literature on this topic for a rigorous formulation (see for example [27,28]). In general, for any two given functions we want to ensure that

We define the norm of a function as

and functions with bounded norm are said to be squareintegrable. The space of functions that are squareintegrable in is denoted as .
Although no proof is given here, it can be shown that, for any given time instant , it is sufficient to require that both the momentum test function , the velocity solution and their first order derivatives belong to . The space of functions verifying this property is a Hilbert space commonly denoted as in functional analysis. For the mass conservation test function and the pressure solution , it is enough to require them to be squareintegrable, as their spatial derivatives do not appear in Eqs. 14 and 15.
Considering that must verify the Dirichlet boundary condition when evaluated in the Dirichlet boundary and is zero in by definition, the solutions (for any fixed instant in time) and test functions must be contained in the spaces of functions given by

Likewise, the external forces must be such that the domain integral in the right hand side of Eq. 14 is well defined. Given that , this is equivalent to requiring the forces to belong to the dual of that space, denoted as . There is a similar requirement on the traction , as it has to be integrable when multiplied by the test function restricted to the Neumann boundary, but it will not be formally stated here.
Finally, to make sure that the dynamic problem is wellposed it is sufficient to require that and squareintegrable along the time interval of the problem. This is denoted as . In the case of the pressure it is enough to enforce that the norm is squareintegrable in time, that is, .
Before we introduce the variational form of the problem, there is an important remark to be made about the convective term in the momentum equation Eq. 3. As long as the velocity field is strongly (pointwise) divergencefree, the following three expressions are identical:

(16) 
or, in variational form,

The three expressions in Eq. 16 are respectively known as nonconservative, conservative and skewsymmetric forms of the convective term. However, in practice, the expressions in Eq. 17 are not equivalent for the discrete problem, as the velocity solution is only divergencefree in a weak (integral) sense. In fact, as shown in [29], each of them gives rise to a variational problem with different conservation properties. While the continuous NavierStokes equations enforce the balance of linear and angular momentum, as well as kinetic energy, none of the discrete variants ensures conservation of the three quantities at the same time. According to the analysis in the same reference, the variational form resulting from each expression conserves the quantities listed in Table 1.
Convective term  Linear momentum  Angular mom.  Kinetic energy 
Nonconservative  For equal  interpolations  No  No 
Conservative  Yes  Yes  No 
Skewsymmetric  For equal  interpolations  No  Yes 
In the present work we have used the skewsymmetric form, as the kinetic energy balance is an important concern in the present study and something that can be measured to to gain insight on turbulent flow simulations. In fact, we can report that using a formulation based on the skewsymmetric form resulted in a better fit to DNS data in our preliminary tests for the channel flow simulations presented in Section 2.7.
Using the skewsymmetric form for the convective term, the Galerkin weak form of the NavierStokes problem can be stated as
Find , such that, ,

Unfortunately, this problem is not straightforward to solve using finite elements, as its discrete version is not numerically stable. In fact, to ensure that that problem of finding a pair that satisfies the above weak form for all suitable choices of has a stable solution, the discrete spaces in which the solution is sought must verify the infsup or LadyzhenskayaBabuskaBrezzi (LBB) condition, given by

(20) 
where and are spaces containing the velocity and pressure solutions, respectively, and is a positive constant. In practice, satisfaction of the LBB condition implies the use of a higher order interpolation for velocity than for pressure, as is done for example in TaylorHood elements [30].
Numerical instabilities in the discrete solution of the NavierStokes may also appear in convectiondominated flows, that is, when the convective term is large in relation to the viscous term. In turbulent flows, this issue can be understood from a physical point of view by noting that viscous dissipation occurs predominantly due to high velocity gradients at small scales. If the finite element mesh is coarse, these gradients cannot be reproduced and dissipation is underestimated, which leads to energy accumulation on the larger scales and the eventual loss of convergence of the solution.
Both instabilities can be cured by the use of a stabilized formulation, which involves obtaining a modified weak form not restricted by the infsup condition of Eq. 20. This has the advantage of allowing the use of equal order interpolations for velocity and pressure. Two classical stabilized formulations for the NavierStokes equations are known as StreamlineUpwind PetrovGalerkin (SUPG) [31] and GalerkinLeast Squares (GLS) [32]. Other alternatives are methods within the VMS framework, which are the main subject of the present chapter, and the Finite Calculus (FIC) approach, presented in Chapter 3.
The Variational Multiscale (VMS) method, introduced in [15,33,16] is a theoretical framework for the development of stabilized finite element methods that has been used extensively during the last two decades in finite element formulations for fluid problems. The basic premise of this method is the separation of the problem variables into large scale and small scale, or subscale, values , which in our case corresponds to

The notation chosen to represent the large scales should be understood as a reference to the finite element mesh size . In practice, scale separation is closely related to the spatial discretization used to solve the problem. This is shown graphically in Figure 1, where the solution along a line is represented. The large scales correspond to the part of the solution that can be described using the finite element interpolation, while the small scales represent fine features of the solution that cannot be reproduced due to the limitations of the discrete interpolation.
Figure 1: Scale separation. 
Once a finite element discretization is defined, the space containing the large scale part of the solution can be identified with that of the admissible finite element functions on the discrete domain. This implies that, instead of working with test functions and solutions defined on infinitedimensional spaces of functions, we now seek a solution on a restricted, finitedimensional space of admissible solutions, expressed as

This in turn allows us to define spaces containing the small scale part of the solution , . The subscale spaces complete the corresponding large scale space, that is, , . As a result, it is clear that small scale spaces are infinitedimensional, unlike the large scale ones, and have to be approximated in order to obtain a solution. There is no single way to approximate them and, in fact, the choice of approximate space for the small scales is one of the defining features of the solution method.
We introduce the following compact notation

Note that we have modified the convective terms in to introduce a convection velocity . This notation is introduced for convenience, as it will allow us to introduce a linearized version of the operator later. Using the compact notation, the weak form of the problem can be expressed as

We introduce the scale separation of the problem variables in Eqs. 21 and 22. Testing against large scale functions, we can write

If the small scale test functions are used instead, the following expression is obtained

The large scale problem given by Eqs. 23 and 24 represents the finite element approximation of the original problem, now containing terms that describe the effect of the unresolved scales on the large scale solution. These additional terms cannot be evaluated in practice, as the small scale variables are not known. VMS methods are based on the definition of a model for the small scale variables, which is motivated by the small scale problem of Eqs. 25 and 26. This model can then be introduced in the large scale equations, closing the problem.
Before we apply this procedure to our problem, it is convenient to operate on Eqs. 23 and 24 to eliminate all spatial derivatives of small scale functions, obtaining an expression that depends on and , which will be modeled, but not on their gradients, which will remain unknown. We start by separating the terms involving the large and small scale parts of the solution

Integrating by parts within each element in the mesh, differential operators acting on and are moved to the test functions. To do so, we introduce the notation of to refer to the part of the problem domain corresponding to element and to denote its boundary.

In the present work, the boundary integrals appearing in Eqs. 30 and 31 will be neglected. This is common in VMS formulations, and can be understood as assuming that the small scale unknowns vanish on element boundaries.
In the following we express this elementbyelement integration using the notation

Now the question is to define a model for and that can be used to evaluate the domain integrals in Eqs. 30 and 31. To do this, we start from the small scale problem, given by Eqs. 25 and 26, which can be developed to read

Again, it is convenient to use integration by parts elementbyelement on some of the terms in the momentum equation, obtaining the expression

(36) 
Note that the elemental boundary integrals in Eq. 36 vanish over internal boundaries because they involve either the exact traction or the exact velocity over the boundary, which are both continuous. Similarly, they cancel out with the corresponding traction in the Neumann boundary. As a result, all boundary terms in Eq. 36 can be eliminated in the following.
Rearranging terms in Eqs. 36 and 35 to separate large and small scale unknowns we can write

Eqs. 40 and 41 can be understood as the projection onto the space of small scales of a differential equation, in the same sense as Eqs. 10 and 11 represent the projection onto of the NavierStokes equations. Moreover, observing the right hand side terms of Eq. 40, we can see that it represents the projection of the strongform linear momentum equation applied to the large scale part of the solution ,. The same can be said about the right hand side of Eq. 40, which corresponds to the mass conservation equation applied to the large scale velocity.
Denoting the projection onto the small scale spaces , with and , respectively, Eqs. 40 and 41 equation can be stated as

Since Eqs. 45 and 46 must hold for all admissible smallscale test functions and , they are equivalent to imposing that the small scale variables , verify the following problem in each element :

where and represent the residual form of the NavierStokes equations applied to the large scale variables

and and are chosen to enforce that right hand sides of Eqs. 48 and 49 belong to the corresponding small scale space.
As mentioned in the previous pages, the small scale spaces are infinitedimensional and have to be approximated before the problem given by Eqs. 48 and 49 can be solved. In practice, the definition of an approximate small scale space corresponds to a choice of projectors and . The most straightforward possibility is to use the entire residual (without projecting to a particular space), which corresponds to considering operators , equal to the identity function or, equivalently, and . This formulation gives rise to the algebraic subgrid scale (ASGS) method [34].
Another wellknown choice consists in taking a small scale space that is orthogonal to the large scale space. If and are the projection onto the large scale spaces and respectively, then the projections in Eqs. 45 and 46 are defined as , . In this case, and are chosen to subtract from the equation the part of the residuals that belongs to the finite element space, that is,

This choice leads to the orthogonal subscale (OSS) method, presented in [35,23].
Note that, due to their definition, and belong to the space of finite element functions and can be constructed from their values on mesh nodes using standard finite element interpolation functions.
A second important remark is that the calculation of and requires knowledge of the finite element solutions , and, as a result, it is coupled to the solution of the stabilized NavierStokes equations. In principle, this would double the number of nodal degrees of freedom of the problem. However, in practice, given that the NavierStokes problem is nonlinear and has to be solved iteratively anyway, the projection problem can be implemented in an staggered way, updating the projections after each nonlinear NavierStokes iteration.
The problem for the small scales, given by Eqs. 48 and 49, is not usually solved. Instead, it is approximated by an expression of the form

(55) 
where second order tensor and the scalar , known as stabilization parameters, are algorithmic quantities that have to be defined to complete the method.
A motivation for this expression can be found in [23], where the parameters in Eq. 55 are designed to ensure that the norm of the modeled subscale variables is the same as that of the exact small scale values. A different justification, based on the approximation of the Green's function of the small scale problem, is provided in [16,36].
While Eq. 55 represents the complete small scale model derived from the dynamic NavierStokes equations, many VMS formulations that can be found in the literature (see for example [34,37,36]) neglect the time variation of the velocity small scale model and define the small scales as

(56) 
This choice corresponds to the assumption that the velocity small scales adapt automatically to the large scale residual. Following the nomenclature of [34], we refer to models based on Eq. 56 as quasistatic subscales, as opposed to dynamic subscale models, based on Eq. 55.
To complete the formulation, a definition for the stabilization parameters is needed. We follow the approach of [34], where the velocity stabilization parameter is taken to be a diagonal matrix and

where is a characteristic length of the element and , are constants, which, for linear finite elements, are usually defined as , (this is the case for example in [38] or [29]). However, the studies presented in [25] for a turbulent channel flow in the low Mach number regime suggest that the choice of values for these parameters can have an impact on the solution. Based on the results presented in that reference and in [26], we have adopted , for our tests.
As pointed out in [20], the use of quasistatic subscales leaves open the possibility of instabilities appearing for small time steps once the problem is discretized in time. The same instability is studied in [39] for the Stokes problem, where it is shown that it can be neutralized if the stabilization parameter satisfies the condition

(59) 
where is the time step and is a constant. To avoid this instability, the stabilization parameter can be replaced by the modified expression

(60) 
although, as remarked in [20], this introduces a dependency of the solution on the time step, even for problems that result in a stationary solution.
If, instead of the quasistatic model of Eq. 56, we use the dynamic model of 55, the time evolution of the velocity small scale has to be taken into account. This was achieved in [23,20] by introducing a time discretization for the small scale acceleration, resulting in the timediscrete small scale model

(61) 
Taking , which corresponds to a backward Euler time scheme, we can write a closed expression for the small scale velocity, given by

(62) 
As remarked in [20], the effective stabilization parameter in Eq. 62 is

(63) 
which is precisely the expression introduced as in Eq. 60 and prevents the apparition of instabilities due to small time steps. However, unlike in quasistatic approximations, when the problem has a stationary solution, the dependency on the time step is eliminated as, in that case, and the quasistatic model of Eq. 56 is recovered.
From the point of view of its implementation, the main difference between the quasistatic model of Eq. 56 and Eq. 62 is that the latter introduces the old value of the velocity small scale in the model. This means that has to be tracked in time. In practice, this implies evaluating and storing historical values for on the integration points of the finite element mesh.
Once the small scale model is defined, it can be introduced in the large scale equations, given by Eqs. 30 and 31. The most general formulation that can be obtained from the material presented in the previous section is given by
Momentum equation

(64) 
Mass conservation

(65) 
The first two rows of Eq. 64, in combination with its right hand side, constitute the standard Galerkin discretization of the momentum equation. The terms in the third row model the effect of the velocity small scale fluctuations and the velocity small scale itself, respectively, on the large scale equations, while the first term in the last row of Eq. 64 represents the effect of pressure small scales.
Analogously, the left hand side of Eq. 65 represents the Galerkin weak form of the incompressibility equation, while its right hand side models the effect of the small scales in mass conservation.
While Eqs. 64 and 65 represent a general VMS formulation, they can be particularized to recover several wellknown stabilized methods:
For OSS formulations, the projections and are defined as the projections of the momentum and mass residuals, respectively, onto the finite element mesh. Applying this definition, they can be obtained as the solution of the projection problem

The only question that remains to close the formulation is to provide a formal definition for the auxiliary convection velocity , something that we have been deliberately avoiding up to this point. This variable was introduced to linearize the convective term, and in practice we can identify it with the last known value of velocity, which corresponds to a Picard linearization of the momentum equation. However, in a context of scale separation, should we use the large scale velocity or the full velocity ?
Both choices result in viable stabilized formulations. The classical approach is to use only the large scale velocity , as it corresponds to the finite element solution. However, from a theoretical point of view, using the full velocity has interesting implications for turbulence modeling, which will be the main focus of Section 2.4.
The concept of scale separation introduced by VMS formulations has some parallels with Large Eddy Simulation (LES) methods for the simulation of turbulent flows. LES turbulence models are also based on separating large and small motions in the flow, but in the LES approach this is is typically achieved through the introduction a filtering operation [2], defined as

(70) 
where is a filter function defined on the interval and is known as the filter width.
Applying a filter function to Eqs. 3 and 4 allows us to write the filtered NavierStokes equations, given by

where we have used the conservative form of the convective term and is the subgrid stress tensor, defined as

(73) 
The subgrid stress tensor represents the effect of the small scale motions on the large scale part of the solution and and is an unknown, since the quantity cannot be obtained from the filtered velocity . This means that Eqs. 71 and 72 do not represent a closed expression. However, if the filter width is chosen small enough that the filteredout small scale motions can be assumed to lie in the inertial subrange, Kolmogorov's hypotheses [3,4] tell us that they have an isotropic, universal (problemindependent) behavior. LES methods use this approach to motivate a model for and justify its introduction in the filtered equations, closing the formulation.
Introducing the subgrid velocity , the subgrid stress tensor can be rewritten using the Leonard decomposition as

(74) 
where each of the individual terms is defined as:

The different terms in Eqs. 75–77 represent subgrid stresses due to the interaction between resolved motions (Leonard stresses), to the interaction between large and unresolved motions (cross stresses) and to the effect of completely unresolved motions (Reynolds stresses).
Like filtering in LES methods, scale separation in VMS formulations introduces a clear division between the resolved and unresolved parts of the solution, which in this case is achieved through the projection to the finite element mesh. This projection to the mesh was introduced in writing the large scale equation, given by Eqs. 23 and 24, which is rewritten using the conservative form of the convective term to be consistent with the LES expression above, obtaining

where boundary terms and terms related to the trace of the viscous stresses have been omitted for clarity and the full velocity has been used for both arguments of the convective term.
An analogy can be established between Eqs. 71 and 72, which represent the filtered NavierStokes equations and Eqs. 79 and 80, which expresses the projection of the NavierStokes equations to the finite element mesh. This was analyzed in [17,20], where it is remarked that the convective term in Eq. 79 can be expanded as

(81) 
The last three terms in Eq. 81 can be understood as a variational version of the LES subgrid stress tensor . Ignoring the density, they can be rearranged as

(82) 
where we introduced as a variational analogue of the LES subgrid stress tensor.
Furthermore, can be decomposed into Cross and Reynolds terms as:

while an analogue of the Leonard stress, representing the contribution of the resolved velocities to the unresolved stresses, appears in the corresponding small scale equation

(85) 
Assuming that the grid size is small enough for the unresolved scales to be in the inertial subrange, we can observe that scale separation and projection to the finite element mesh play a similar role to that of filtering in classical LES methods. It is worth mentioning that the concept of meshinduced filtering has also been explored by the LES community, where it is known as implicit filtering. This is the basis of the Monotone Integrated LargeEddy Simulation (MILES) approach of Boris et al. [40], where the authors propose not using an explicit model for the subgrid stresses and relying instead on a specially designed numerical method to introduce the correct dissipation for a given mesh resolution.
A different approach to VMSbased LES, which will only be mentioned here, is based on a three–level scale separation (see for example [37,41] or the review of [19]). In such approaches, the small scales are in turn divided into resolved and unresolved small scales. The presence of two levels of small scales can be used to either introduce explicit LES model terms to represent the effect of the unresolved small scales on the problem or to calibrate the amount of dissipation that is introduced using the variational equivalent of the Germano identity [42,43].
Given that the derivation of VMS formulations is based exclusively in numerical and mathematical arguments, the physical behavior of VMS methods, in terms of reproducing the expected dissipation rates when the small scales are on the inertial subrange, has to be verified. This topic has been studied in [17,18] for the formulation that we are calling QASGS and in [24,44] for OSSbased methods.
An energy balance for the original NavierStokes problem can be obtained by taking in the Galerkin weak form of the momentum equation, given by Eq. 14, which produces

(86) 
We are interested in obtaining a balance for the kinetic energy, . Using the fact that the full velocity is incompressible and the equality , Eq. 86 can be expressed as

(87) 
Eq. 87 expresses the balance of total kinetic energy in the domain, and the individual terms represent energy storage () and convection (), the power exerted by external forces () and boundary tensions () and finally viscous dissipation ().
We can obtain an equivalent expression for the kinetic energy contained in the large scale motions, , by taking in Eq. 79, in Eq. 80 and adding the two equations, obtaining:

(88) 
Note that in Eq. 88, and for the remainder of this section, we neglected all boundary integrals to simplify the discussion. This is equivalent to considering a problem with homogeneous Dirichlet boundary conditions.
It is convenient to integrate some of the terms in Eq. 88 by parts and rearrange the convective term as follows:

where we have used the fact that the exact velocity is divergence free.
With this, Eq. 88 can be restated as

(89) 
Analogously to the complete energy balance, Eq. 89 can be understood as a balance for the kinetic energy associated to large scale motions. Terms and represent the storage and convection of large scale kinetic energy, while term represents the power exerted by the external forces on large scale motions. Term is the viscous dissipation associated to the large scale motions, which can be assumed to be negligible for high Reynolds numbers. The remaining terms represent the transfer of energy between large scale and residual motions, playing an analogous role to production terms in the filtered NavierStokes equations.
Finally, we define the residual kinetic energy as . A balance statement for can be obtained by subtracting Eq. 89 from Eq. 87, which results in

(90) 
In Eq. 90, terms and represent the storage and convection of residual kinetic energy, while term represents the power exerted by the external forces on small scale (high wavenumber) motions. The next two terms represent the difference between the total viscous dissipation () and the large scale viscous dissipation (), which was already accounted for in the large scale energy balance of Eq. 89. Again, we remark that, in practice, term is expected to be negligible in comparison to term , since viscous dissipation occurs predominantly for motions in the range of the Kolmogorov length scale, while will contain only motions on a much larger scale , lying on the inertial subrange. Finally, terms and represent the production of residual energy due to the pressure and velocity small scales, respectively, and are exactly the same (but now with opposite sign) as the production terms in Eq. 89.
Terms and can be modified by noting that, since and ,

The last step to obtain the residual kinetic energy balance is to introduce the small scale models for velocity and pressure in Eqs. 91 and 93. Noting that represents the large scale viscous dissipation and will be negligible for high Reynolds numbers, we can write

(94) 
It is clear that terms and , since the stabilization parameters and were defined as strictly positive. They play the role of energy sinks in the large scale energy balance of Eq. 89, while acting as sources on the residual energy equation Eq. 90. Term is problem dependent but, if an OSS small scale model is used, it will vanish unless the external forces have a highfrequency (small scale) component. Finally, term only exists for dynamic small scale models. This fact was used in [24] to justify that OSS formulations can account for backscatter (energy transfer from the small scales to the large ones) if a dynamic small scale model is used.
Starting from a similar reasoning, Guasch and Codina [44] use statistical and scaling arguments to show that OSS formulations extract energy from the large scale equations at the correct rate, provided that some constraints on the behavior of stabilization parameters are respected.
As a final remark on this topic, we compare the large scale energy balance to the filtered energy balance used in filterbased LES formulations. If the kinetic energy associated to the filtered velocity field is defined as , the balance for can be obtained by multiplying the filtered linear momentum equation Eq. 71 by and integrating over the fluid domain (see [2]), resulting in

(95) 
where boundary fluxes have been omitted.
We can introduce the definition of the VMS subgrid stress tensor in Eq. 88, obtaining

(96) 
where the conservative form of the convective term has been used. Rearranging some terms, Eq. 96 can be rewritten as

(97) 
where we have used the fact that is a symmetric tensor to write . Terms to have an analogous interpretation to their counterparts in Eq. 89, but here we have obtained an additional term, , which explicitly represents the contribution of the residual subgrid stresses in the energy transfer. Comparing Eq. 97 to Eq. 95, we see that both LES and VMS approaches extract energy from the large scale problem through subgrid stresses, but that the variational approach gives rise to two additional energy transfer mechanisms, represented by terms and , compared to filterbased LES.
To obtain a finite element solver based on the VMS formulation introduced in Section 2.3 we need to discretize the simulation domain, both in space and in time, and linearize the problem to obtain a system of equations that can be inverted using a linear solver. We start by introducing a finite element partition for the problem domain . Given the discrete domain , the large scale interpolation spaces and can be identified with the standard finite element interpolation functions and the large scale part of the solution, and , can be represented using a finite element interpolation as

(98) 
where represents the number of nodes in the finite element mesh, and are the nodal values of the large scale variables and respectively, represents the standard finite element basis function associated to node and its counterpart for vectorial variables, given by

Additionally, we introduce the following notation for the gradient and divergence of the finite element shape functions, which will be used to write the discrete form of the differential operators involved in the problem

(99) 
We also introduce the following operator to describe convection

(100) 
and the strain rate matrix, which will be used to write the discrete version of the viscous term

(101) 
Additionally, we define , and as the vectors of nodal values of large scale velocity , acceleration and pressure , respectively.
We introduce the finite element discretization of Eq. 98 in the variational formulation given by Eqs. 64 and 65 to obtain the matrix form of the problem. We will analyze the resulting expression for each of the variants we are considering in turn.
Note that in the present work we use linear finite elements, which can not be used to write second derivatives of the variables or test functions. As a result, terms involving in Eq. 64 or the strongform viscous term that appears in the residual introduced in Eq. 51 will be neglected in the discrete form. It must be remarked that all terms lost in this way are related to viscous stresses, which should be small in turbulent flow problems.
The quasistatic ASGS formulation is in some sense the classical VMS formulation for the NavierStokes equations. It is also relatively simple, as it does not involve dynamic terms or projections, so we will present it first.
Starting from Eqs. 64 and 65, we can neglect all terms involving the projections or . As the small scales are described using Eq. 56, terms involving or in these equations can be ignored and the momentum stabilization parameter is , given by Eq. 57. After these simplifications, we introduce the finite element interpolation of Eq. 98 to describe the problem variables , . Testing against each nodal basis function in turn we obtain a system of equations that can be expressed in matrix form as

(102) 

(103) 
Again, we have used a generic convection velocity in all terms that have a nonlinear dependence of velocity. Doing so, we leave open the possibility of using either the full velocity or the only large scale part . Note that, for linear finite elements, the latter choice is equivalent to the GalerkinLeast Squares (GLS) method [32].
The different matrices in Eqs. 102 and 103 represent the discrete version of the operators in Eqs. 64 and 65 and can be built from the assembly of elemental contributions. In general, matrix is constructed by the finite element assembly of elemental matrices of the form . For a finite element with nodes, can be defined using blocks , where and are local node indices. Using this notation, the standard Galerkin terms in the variational form of the problem give rise to the following elemental matrices

When defining the viscous matrix in Eq. 107 we introduced the constitutive matrix , which is given by

(112) 
In the same way, the stabilization terms in the QASGS give rise to additional elemental matrices, expressed here as

The next variant to be presented is the quasistatic OSS formulation. Compared to the QASGS formulation, OSS is characterized by the inclusion of the projections and , which make the small scale variables orthogonal to the large scale unknowns and should reduce the overall amount of numerical diffusion introduced in the problem. As in the previous case, we leave open the possibility of using either or for convection and the corresponding stabilization terms. The matrix form of the QOSS formulation can be expressed as

(122) 

(123) 
where the following three new terms terms, involving the projections, have been introduced:

The calculation of the projections involves the solution of an additional problem, given by Eqs. 68 and 69, which can be expressed in discrete form as

where and represent the vectors of nodal values of and respectively and the remaining matrices and vectors are given by
