# Abstract

This work explores the use of stabilized finite element formulations for the incompressible Navier-Stokes 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 well-known 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 high-performance 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 non-Newtonian flow problems.

# 1 Introduction

## 1.1 Background and motivation

### 1.1.1 Turbulence modeling

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 ${\textstyle L}$ 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 ${\textstyle \eta }$.

Turbulent flows can be characterized by the Reynolds number ${\displaystyle {\hbox{Re}}}$, a dimensionless parameter indicative of the balance between inertial and viscous forces in the problem,

 ${\displaystyle {\hbox{Re}}={\frac {\rho \,UL}{\mu }}}$
(1.1)

where ${\textstyle \rho }$ is the fluid's density and ${\textstyle \mu }$ its dynamic viscosity and ${\textstyle U}$ 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

 ${\displaystyle {\frac {L}{\eta }}={\hbox{Re}}^{\frac {3}{4}}}$
(1.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 ${\textstyle 3D}$ problem, this implies that the total number of elements would be on the order of ${\textstyle {\hbox{Re}}^{9/4}}$. Keeping in mind that the Reynolds number in engineering applications is typically on the order of ${\textstyle 10^{4}\sim 10^{6}}$, 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 Navier-Stokes (RANS) and Large Eddy Simulation (LES).

RANS models are based on rewriting the problem in terms of time-averaged variables (or time- and space-averaged, 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 Spalart-Allmaras model [5], which involves one additional equation, or the ${\textstyle \kappa }$${\textstyle \epsilon }$ [6,7] and ${\textstyle \kappa }$${\textstyle \omega }$ 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 problem-independent. 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.

### 1.1.2 Stabilized finite element formulations for turbulent flows

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 convection-dominated 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 Streamline-Upwind Petrov-Galerkin (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 mesh-induced filtering. This has motivated multiple investigations studying the use VMS formulations in turbulence modeling [17,18,19,20,21,22,23,24,25,26,27], analyzing their relationship both from the theoretical point of view and through its application to the simulation of turbulent flows of interest.

## 1.2 Objectives and methodology

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 LES-like 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 well-known 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 Navier-Stokes 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 non-Newtonian 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 [28,29], 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.

## 1.3 Outline of this document

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 non-Newtonian flow cases.

Finally, in Chapter 6 we summarize the work and present its main conclusions, as well as proposing future lines of research.

# 2 Variational multiscale stabilization for turbulent flow problems

## 2.1 Introduction

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 [30,20], which have been shown to provide good results in turbulent flow simulations without requiring the use of an explicit turbulence model [25,31,27].

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 [31]) and on using or neglecting the pressure small scale [27]. 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 Navier-Stokes 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.

## 2.2 Variational form of the Navier-Stokes equations

### 2.2.1 Problem statement

The incompressible Navier-Stokes equations state the balance of linear momentum and mass in a fluid domain ${\textstyle \Omega }$, given by

 ${\displaystyle \rho \,\partial _{t}u+\rho \,u\cdot \nabla u-\nabla \cdot {\boldsymbol {\sigma }}={\boldsymbol {f}}{\hbox{ in }}\Omega \times \left[0,T\right)}$ (2.3) ${\displaystyle \nabla \cdot u=0{\hbox{ in }}\Omega \times \left[0,T\right)}$ (2.4)

where ${\textstyle u}$ is the fluid velocity, ${\textstyle \rho }$ its density, ${\textstyle {\boldsymbol {\sigma }}}$ represents the stress tensor and ${\textstyle {\boldsymbol {f}}}$ the external forces acting on the domain.

The problem described by Eqs. 2.3 and 2.4 must be completed with suitable initial and boundary conditions. We denote the domain boundary as ${\textstyle \partial \Omega }$ and introduce its partition into Dirichlet (${\textstyle \Gamma _{D}}$) and Neumann (${\textstyle \Gamma _{N}}$) parts, verifying ${\textstyle \partial \Omega =\Gamma _{D}\cup \Gamma _{N}}$ and ${\textstyle \Gamma _{D}\cap \Gamma _{N}=\emptyset }$. The initial and boundary conditions for the problem can be expressed as:

 ${\displaystyle u=u_{0}{\hbox{ in }}\Omega ,\,t=0}$ (2.5) ${\displaystyle u=u_{D}{\hbox{ in }}\Gamma _{D}\times \left[0,T\right)}$ (2.6) ${\displaystyle {\boldsymbol {\sigma }}\cdot {\boldsymbol {n}}={\boldsymbol {t}}{\hbox{ in }}\Gamma _{N}\times \left[0,T\right)}$ (2.7)

where ${\textstyle u_{0}}$ is the initial velocity field, ${\textstyle u_{D}}$ represents the imposed velocity on the Dirichlet boundary, ${\textstyle {\boldsymbol {n}}}$ the outer normal vector and ${\textstyle {\boldsymbol {t}}}$ the imposed traction acting along the Neumann boundary. Note that the initial velocity ${\textstyle u_{0}}$ must be chosen to be divergence free to ensure that Eq. 2.4 is verified at all times.

For Newtonian fluids, the stress tensor ${\textstyle {\boldsymbol {\sigma }}}$ can be related to the fluid velocity ${\textstyle u}$ and pressure ${\textstyle p}$ using the constitutive relation

 ${\displaystyle {\boldsymbol {\sigma }}=-p\,{\boldsymbol {I}}+2\mu \left(\nabla ^{s}u-{\frac {1}{3}}\,\left(\nabla \cdot u\right){\boldsymbol {I}}\right)}$
(2.8)

where ${\textstyle {\boldsymbol {I}}}$ is the second order identity tensor, ${\textstyle \mu }$ the fluid's viscosity and ${\textstyle \nabla ^{s}u}$ the symmetric gradient of velocity, defined as

 ${\displaystyle \nabla ^{s}u={\frac {1}{2}}\left(\nabla u+\left(\nabla u\right)^{T}\right)}$
(2.9)

To obtain a finite element formulation for the Navier-Stokes equations we need to rewrite them in weak form. We multiply Eq. 2.3 by a velocity test function ${\textstyle {\boldsymbol {w}}}$, defined to be zero on the Dirichlet boundary ${\textstyle \Gamma _{D}}$, and Eq. 2.4 by a pressure test function ${\textstyle q}$. Integrating the resulting expressions over the fluid domain we obtain

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}\cdot \left(\rho \,\partial _{t}u+\rho \,u\cdot \nabla u-\nabla \cdot {\boldsymbol {\sigma }}\right)\,{\hbox{d}}\Omega =\int _{\Omega }{\boldsymbol {w}}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega }$ (2.10) ${\displaystyle \int _{\Omega }q\,\nabla \cdot u\,{\hbox{d}}\Omega =0}$ (2.11)

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

 ${\displaystyle -\int _{\Omega }{\boldsymbol {w}}\,\nabla \cdot {\boldsymbol {\sigma }}\,{\hbox{d}}\Omega =\int _{\Omega }\nabla {\boldsymbol {w}}:{\boldsymbol {\sigma }}\,{\hbox{d}}\Omega -\int _{\Omega }\nabla \cdot \left({\boldsymbol {w}}\cdot {\boldsymbol {\sigma }}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle =\int _{\Omega }\nabla {\boldsymbol {w}}:{\boldsymbol {\sigma }}\,{\hbox{d}}\Omega +\int _{\partial \Omega }{\boldsymbol {w}}\cdot {\boldsymbol {\sigma }}\cdot {\boldsymbol {n}}\,{\hbox{d}}\Omega =\int _{\Omega }\nabla {\boldsymbol {w}}:{\boldsymbol {\sigma }}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}\cdot {\boldsymbol {t}}\,{\hbox{d}}\Omega }$
(2.12)

Substituting Eq. 2.12 into Eq. 2.10 and introducing the definition of the stress tensor in Eq. 2.8 we obtain

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}\cdot \left(\rho \,\partial _{t}u+\rho \,u\cdot \nabla u\right)\,{\hbox{d}}\Omega +\int _{\Omega }\nabla ^{s}{\boldsymbol {w}}:\,2\mu \left(\nabla ^{s}u-{\frac {1}{3}}\left(\nabla \cdot u\right)\right)\,{\hbox{d}}\Omega }$ ${\displaystyle -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}\,p\,{\hbox{d}}\Omega =\int _{\Omega }{\boldsymbol {w}}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}\cdot {\boldsymbol {t}}\,{\hbox{d}}\Omega }$ (2.13) ${\displaystyle \int _{\Omega }q\,\nabla \cdot u\,{\hbox{d}}\Omega =0}$ (2.14)

Eqs. 2.13 and 2.14, 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 [32,33]). In general, for any two given functions ${\textstyle f,g}$ we want to ensure that

 ${\displaystyle \int _{\Omega }fg\,{\hbox{d}}\Omega <\infty }$

We define the ${\textstyle L^{2}}$ norm of a function as

 ${\displaystyle \left\Vert f\right\Vert _{\Omega }=\left(\int _{\Omega }f^{2}\,{\hbox{d}}\Omega \right)^{\frac {1}{2}}}$

and functions with bounded ${\textstyle L^{2}}$ norm are said to be square-integrable. The space of functions that are square-integrable in ${\textstyle \Omega }$ is denoted as ${\textstyle L^{2}\left(\Omega \right)}$.

Although no proof is given here, it can be shown that, for any given time instant ${\textstyle t}$, it is sufficient to require that both the momentum test function ${\textstyle {\boldsymbol {w}}}$, the velocity solution ${\textstyle u}$ and their first order derivatives belong to ${\textstyle L^{2}(\Omega )}$. The space of functions verifying this property is a Hilbert space commonly denoted as ${\textstyle H^{1}(\Omega )}$ in functional analysis. For the mass conservation test function ${\textstyle q}$ and the pressure solution ${\textstyle p}$, it is enough to require them to be square-integrable, as their spatial derivatives do not appear in Eqs. 2.13 and 2.14.

Considering that ${\textstyle u}$ must verify the Dirichlet boundary condition when evaluated in the Dirichlet boundary ${\textstyle \Gamma _{D}}$ and ${\textstyle {\boldsymbol {w}}}$ is zero in ${\textstyle \Gamma _{D}}$ by definition, the solutions (for any fixed instant in time) and test functions must be contained in the spaces of functions given by

 ${\displaystyle u\in H_{D}^{1}=\left\lbrace \left.u\in H^{1}(\Omega )\,\right\vert \,u=u_{D}{\hbox{ in }}\Gamma _{D}\right\rbrace }$ ${\displaystyle {\boldsymbol {w}}\in H_{0}^{1}=\left\lbrace \left.{\boldsymbol {w}}\in H^{1}(\Omega )\,\right\vert \,{\boldsymbol {w}}={\boldsymbol {0}}{\hbox{ in }}\Gamma _{D}\right\rbrace }$ ${\displaystyle p,q\in L^{2}(\Omega )}$

Likewise, the external forces ${\textstyle {\boldsymbol {f}}}$ must be such that the domain integral in the right hand side of Eq. 2.13 is well defined. Given that ${\textstyle {\boldsymbol {w}}\in H^{1}(\Omega )}$, this is equivalent to requiring the forces to belong to the dual of that space, denoted as ${\textstyle H^{-1}(\Omega )}$. There is a similar requirement on the traction ${\textstyle {\boldsymbol {t}}}$, as it has to be integrable when multiplied by the test function ${\textstyle {\boldsymbol {w}}}$ restricted to the Neumann boundary, but it will not be formally stated here.

Finally, to make sure that the dynamic problem is well-posed it is sufficient to require that ${\textstyle \left\Vert u\right\Vert _{\Omega }}$ and ${\textstyle \left\Vert \partial u_{i}/\partial x_{j}\right\Vert _{\Omega }}$ square-integrable along the time interval of the problem. This is denoted as ${\textstyle u\in L^{2}\left(0,T,H_{D}^{1}(\Omega )\right)}$. In the case of the pressure it is enough to enforce that the ${\textstyle L^{2}}$ norm is square-integrable in time, that is, ${\textstyle p\in L^{2}\left(0,T,L^{2}(\Omega )\right)}$.

### 2.2.2 Conservation properties

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. 2.3. As long as the velocity field is strongly (point-wise) divergence-free, the following three expressions are identical:

 ${\displaystyle u\cdot \nabla u=\nabla \cdot \left(u\otimes u\right)={\frac {1}{2}}\,u\cdot \nabla u+{\frac {1}{2}}\,\nabla \cdot \left(u\otimes u\right)}$
(2.15)

or, in variational form,

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}\cdot \left(u\cdot \nabla u\right)\,{\hbox{d}}\Omega =}$ ${\displaystyle -\int _{\Omega }\nabla {\boldsymbol {w}}:\left(u\otimes u\right)\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}\left({\boldsymbol {w}}\cdot u\right)\left(u\cdot {\boldsymbol {n}}\right)\,{\hbox{d}}\Omega =}$ (2.16) ${\displaystyle {\frac {1}{2}}\int _{\Omega }{\boldsymbol {w}}\cdot \left(u\cdot \nabla u\right)\,{\hbox{d}}\Omega -{\frac {1}{2}}\int _{\Omega }\nabla {\boldsymbol {w}}:\left(u\otimes u\right)\,{\hbox{d}}\Omega +{\frac {1}{2}}\int _{\Gamma _{N}}\left({\boldsymbol {w}}\cdot u\right)\left(u\cdot {\boldsymbol {n}}\right)\,{\hbox{d}}\Omega }$ (2.17)

The three expressions in Eq. 2.15 are respectively known as non-conservative, conservative and skew-symmetric forms of the convective term. However, in practice, the expressions in Eq. 2.16 are not equivalent for the discrete problem, as the velocity solution is only divergence-free in a weak (integral) sense. In fact, as shown in [34], each of them gives rise to a variational problem with different conservation properties. While the continuous Navier-Stokes 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 Non-conservative For equal ${\displaystyle u}$-${\displaystyle p}$ interpolations No No Conservative Yes Yes No Skew-symmetric For equal ${\displaystyle u}$-${\displaystyle p}$ interpolations No Yes

In the present work we have used the skew-symmetric 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 skew-symmetric form resulted in a better fit to DNS data in our preliminary tests for the channel flow simulations presented in Section 2.7.

### 2.2.3 Galerkin weak form

Using the skew-symmetric form for the convective term, the Galerkin weak form of the Navier-Stokes problem can be stated as

Find ${\displaystyle u\in L^{2}\left(0,T,H_{D}^{1}(\Omega )\right)}$, ${\displaystyle p\in L^{2}\left(0,T,L^{2}(\Omega )\right)}$ such that, ${\displaystyle \forall {\boldsymbol {w}}\in H_{0}^{1}(\Omega ),\forall q\in L^{2}(\Omega )}$,

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}\cdot \left(\rho \,\partial _{t}u+\rho {\frac {1}{2}}u\cdot \nabla u\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla {\boldsymbol {w}}:\rho {\frac {1}{2}}\left(u\otimes u\right)\,{\hbox{d}}\Omega \qquad \qquad \qquad \qquad }$ ${\displaystyle +\int _{\Omega }\nabla ^{s}{\boldsymbol {w}}:\,2\mu \left(\nabla ^{s}u-{\frac {1}{3}}\left(\nabla \cdot u\right)\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}\,p\,{\hbox{d}}\Omega =}$ ${\displaystyle \int _{\Omega }{\boldsymbol {w}}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}\cdot \left({\boldsymbol {t}}+\rho {\frac {1}{2}}\left(u\cdot {\boldsymbol {n}}\right)u\right)\,{\hbox{d}}\Omega }$ ${\displaystyle \int _{\Omega }q\,\nabla \cdot u\,{\hbox{d}}\Omega =0}$

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 ${\textstyle u,p}$ that satisfies the above weak form for all suitable choices of ${\textstyle {\boldsymbol {w}},q}$ has a stable solution, the discrete spaces in which the solution is sought must verify the inf-sup or Ladyzhenskaya-Babuska-Brezzi (LBB) condition, given by

 ${\displaystyle \inf _{q\neq 0\;q\in Q}\sup _{{\boldsymbol {w}}\neq {\boldsymbol {0}}\;{\boldsymbol {w}}\in V}{\frac {\int _{\Omega }q\nabla \cdot {\boldsymbol {w}}\,{\hbox{d}}\Omega }{\left\Vert q\right\Vert _{\Omega }\left\Vert {\boldsymbol {w}}\right\Vert _{\Omega }}}\geq C}$
(2.18)

where ${\textstyle V}$ and ${\textstyle Q}$ are spaces containing the velocity and pressure solutions, respectively, and ${\textstyle C}$ 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 Taylor-Hood elements [35].

Numerical instabilities in the discrete solution of the Navier-Stokes may also appear in convection-dominated 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 inf-sup condition of Eq. 2.18. This has the advantage of allowing the use of equal order interpolations for velocity and pressure. Two classical stabilized formulations for the Navier-Stokes equations are known as Streamline-Upwind Petrov-Galerkin (SUPG) [36] and Galerkin-Least Squares (GLS) [37]. 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.

## 2.3 Variational multiscale stabilization

### 2.3.1 Scale separation

The Variational Multiscale (VMS) method, introduced in [15,38,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 ${\textstyle \left(\cdot \right)_{h}}$ and small scale, or subscale, values ${\textstyle \left(\cdot \right)_{s}}$, which in our case corresponds to

 ${\displaystyle u=u_{h}+u_{s}\quad p=p_{h}+p_{s}}$ ${\displaystyle {\boldsymbol {w}}={\boldsymbol {w}}_{h}+{\boldsymbol {w}}_{s}\quad q=q_{h}+q_{s}}$

The ${\textstyle \left(\cdot \right)_{h}}$ notation chosen to represent the large scales should be understood as a reference to the finite element mesh size ${\textstyle h}$. 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.

 (a) Continuous domain. (b) Discrete domain. (c) Continuous solution. (d) Large scale solution. (e) Small scale solution. 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 infinite-dimensional spaces of functions, we now seek a solution on a restricted, finite-dimensional space of admissible solutions, expressed as

 ${\displaystyle u_{h}\in W_{h}\subset W\equiv L^{2}\left(0,T,H_{D}^{1}(\Omega )\right)\qquad p_{h}\in Q_{h}\subset Q\equiv L^{2}\left(0,T,L^{2}(\Omega )\right)}$

This in turn allows us to define spaces containing the small scale part of the solution ${\textstyle {\boldsymbol {w}}_{s}\in W_{s}}$, ${\textstyle q_{s}\in Q_{s}}$. The subscale spaces complete the corresponding large scale space, that is, ${\textstyle W=W_{h}\oplus W_{s}}$, ${\textstyle Q=Q_{h}\oplus Q_{s}}$. As a result, it is clear that small scale spaces are infinite-dimensional, 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

 ${\displaystyle {\mathcal {B}}\left({\boldsymbol {w}},{\boldsymbol {a}},u\right)=\int _{\Omega }{\boldsymbol {w}}\cdot \left(\rho \,\partial _{t}u+\rho {\frac {1}{2}}{\boldsymbol {a}}\cdot \nabla u\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla {\boldsymbol {w}}:\rho {\frac {1}{2}}\left({\boldsymbol {a}}\otimes u\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Omega }\nabla ^{s}{\boldsymbol {w}}:\,2\mu \left(\nabla ^{s}u-{\frac {1}{3}}\left(\nabla \cdot u\right)\right)\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}\cdot \rho {\frac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right)u\,{\hbox{d}}\Omega }$ ${\displaystyle {\mathcal {D}}\left({\boldsymbol {w}},p\right)=\int _{\Omega }\nabla \cdot {\boldsymbol {w}}\,p\,{\hbox{d}}\Omega }$ ${\displaystyle {\mathcal {L}}\left({\boldsymbol {w}}\right)=\int _{\Omega }{\boldsymbol {w}}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}\cdot {\boldsymbol {t}}\,{\hbox{d}}\Omega }$

Note that we have modified the convective terms in ${\textstyle {\mathcal {B}}\left({\boldsymbol {w}},{\boldsymbol {a}},u\right)}$ to introduce a convection velocity ${\textstyle {\boldsymbol {a}}}$. 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

 ${\displaystyle {\mathcal {B}}\left({\boldsymbol {w}},{\boldsymbol {a}},u\right)-{\mathcal {D}}\left({\boldsymbol {w}},p\right)={\mathcal {L}}\left({\boldsymbol {w}}\right)}$ (2.19) ${\displaystyle {\mathcal {D}}\left(u,q\right)=0}$ (2.20)

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

 ${\displaystyle {\mathcal {B}}\left({\boldsymbol {w}}_{h},{\boldsymbol {a}},u_{h}+u_{s}\right)-{\mathcal {D}}\left({\boldsymbol {w}}_{h},p_{h}+p_{s}\right)={\mathcal {L}}\left({\boldsymbol {w}}_{h}\right)}$ (2.21) ${\displaystyle {\mathcal {D}}\left(u_{h}+u_{s},q_{h}\right)=0}$ (2.22)

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

 ${\displaystyle {\mathcal {B}}\left({\boldsymbol {w}}_{s},{\boldsymbol {a}},u_{h}+u_{s}\right)-{\mathcal {D}}\left({\boldsymbol {w}}_{s},p_{h}+p_{s}\right)={\mathcal {L}}\left({\boldsymbol {w}}_{s}\right)}$ (2.23) ${\displaystyle {\mathcal {D}}\left(u_{h}+u_{s},q_{s}\right)=0}$ (2.24)

The large scale problem given by Eqs. 2.21 and 2.22 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. 2.23 and 2.24. 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. 2.21 and 2.22 to eliminate all spatial derivatives of small scale functions, obtaining an expression that depends on ${\textstyle u_{s}}$ and ${\textstyle p_{s}}$, 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

 ${\displaystyle {\mathcal {B}}\left({\boldsymbol {w}}_{h},{\boldsymbol {a}},u_{h}\right)-{\mathcal {D}}\left({\boldsymbol {w}}_{h},p_{h}\right)+{\mathcal {B}}\left({\boldsymbol {w}}_{h},{\boldsymbol {a}},u_{s}\right)-{\mathcal {D}}\left({\boldsymbol {w}}_{h},p_{s}\right)={\mathcal {L}}\left({\boldsymbol {w}}_{h}\right)}$ ${\displaystyle {\mathcal {D}}\left(u_{h},q_{h}\right)+{\mathcal {D}}\left(u_{s},q_{h}\right)=0}$

Integrating by parts within each element in the mesh, differential operators acting on ${\textstyle u_{s}}$ and ${\textstyle p_{s}}$ are moved to the test functions. To do so, we introduce the notation of ${\textstyle \Omega ^{e}}$ to refer to the part of the problem domain corresponding to element ${\textstyle e}$ and ${\textstyle \Gamma ^{e}}$ to denote its boundary.

 ${\displaystyle {\mathcal {B}}\left({\boldsymbol {w}}_{h},{\boldsymbol {a}},u_{h}\right)-{\mathcal {D}}\left({\boldsymbol {w}}_{h},p_{h}\right)+\int _{\Omega }{\boldsymbol {w}}_{h}\,\rho \partial _{t}u_{s}\,{\hbox{d}}\Omega }$ ${\displaystyle -\sum _{e}\int _{\Omega ^{e}}\rho \left({\boldsymbol {a}}\cdot \nabla {\boldsymbol {w}}_{h}\right)\cdot u_{s}\,{\hbox{d}}\Omega +\sum _{e}\int _{\Gamma ^{e}}{\boldsymbol {w}}\cdot \rho {\frac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right)u_{s}\,{\hbox{d}}\Omega }$ (2.25) ${\displaystyle -\sum _{e}\int _{\Omega ^{e}}2\mu \nabla \cdot \left(\nabla ^{s}{\boldsymbol {w}}_{h}-{\frac {1}{3}}\left(\nabla \cdot {\boldsymbol {w}}_{h}\right)\right)\,u_{s}\,{\hbox{d}}\Omega -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}_{h}\,p_{s}\,{\hbox{d}}\Omega }$ (2.26) ${\displaystyle +\sum _{e}\int _{\Gamma ^{e}}2\mu {\boldsymbol {w}}_{h}\cdot \left(\nabla ^{s}u_{s}\cdot {\boldsymbol {n}}-{\frac {1}{3}}\left(\nabla \cdot u_{s}\right){\boldsymbol {n}}\right)\,{\hbox{d}}\Omega ={\mathcal {L}}\left({\boldsymbol {w}}_{h}\right)}$ (2.27) ${\displaystyle {\mathcal {D}}\left(u_{h},q_{h}\right)=\sum _{e}\int _{\Omega ^{e}}\nabla q_{h}\,u_{s}\,{\hbox{d}}\Omega -\sum _{e}\int _{\Gamma ^{e}}q_{h}u_{s}\cdot {\boldsymbol {n}}\,{\hbox{d}}\Omega }$ (2.28)

In the present work, the boundary integrals appearing in Eqs. 2.27 and 2.28 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 element-by-element integration using the notation

 ${\displaystyle \sum _{e}\int _{\Omega ^{e}}\,{\hbox{d}}\Omega =\int _{\Sigma \,\Omega ^{e}}\,{\hbox{d}}\Omega }$

### 2.3.2 Small scale equation

Now the question is to define a model for ${\textstyle u_{s}}$ and ${\textstyle p_{s}}$ that can be used to evaluate the domain integrals in Eqs. 2.27 and 2.28. To do this, we start from the small scale problem, given by Eqs. 2.23 and 2.24, which can be developed to read

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{s}\cdot \rho \left(\partial _{t}u_{h}+\partial _{t}u_{s}+{\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla \left(u_{h}+u_{s}\right)\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla {\boldsymbol {w}}_{s}:\rho {\frac {1}{2}}\left({\boldsymbol {a}}\otimes \left(u_{h}+u_{s}\right)\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Omega }\nabla ^{s}{\boldsymbol {w}}_{s}:\,2\mu \left(\nabla ^{s}\left(u_{h}+u_{s}\right)-{\frac {1}{3}}\nabla \cdot \left(u_{h}+u_{s}\right){\boldsymbol {I}}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}_{s}\,\left(p_{h}+p_{s}\right)\,{\hbox{d}}\Omega =}$ (2.29) ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{s}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}_{s}\cdot \left({\boldsymbol {t}}+\rho {\frac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right)\left(u_{h}+u_{s}\right)\right)\,{\hbox{d}}\Omega }$ (2.30) ${\displaystyle \int _{\Omega }q_{s}\,\nabla \cdot \left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega =0}$ (2.31)

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

 ${\displaystyle \int _{\Sigma \,\Omega ^{e}}{\boldsymbol {w}}_{s}\cdot \rho \left(\partial _{t}u_{h}+\partial _{t}u_{s}+{\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla \left(u_{h}+u_{s}\right)+{\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes \left(u_{h}+u_{s}\right)\right)\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Sigma \,\Omega ^{e}}{\boldsymbol {w}}_{s}\nabla \cdot \left(\left(p_{h}+p_{s}\right){\boldsymbol {I}}-2\mu \left(\nabla ^{s}\left(u_{h}+u_{s}\right)-{\frac {1}{3}}\nabla \cdot \left(u_{h}+u_{s}\right){\boldsymbol {I}}\right)\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\sum _{e}\int _{\Gamma ^{e}}{\boldsymbol {w}}_{s}\left(\left(p_{h}+p_{s}\right)\cdot {\boldsymbol {n}}+2\mu \left(\nabla ^{s}\left(u_{h}+u_{s}\right)-{\frac {1}{3}}\nabla \cdot \left(u_{h}+u_{s}\right)\right)\cdot {\boldsymbol {n}}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Gamma _{N}}{\boldsymbol {w}}_{s}\cdot \rho {\frac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right)\left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega -\sum _{e}\int _{\Gamma ^{e}}{\boldsymbol {w}}_{s}\cdot \rho {\frac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right)\left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega =}$ ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{s}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}_{s}\cdot {\boldsymbol {t}}\,{\hbox{d}}\Omega }$
(2.32)

Note that the elemental boundary integrals in Eq. 2.32 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. 2.32 can be eliminated in the following.

Rearranging terms in Eqs. 2.32 and 2.31 to separate large and small scale unknowns we can write

 ${\displaystyle \int _{\Sigma \,\Omega ^{e}}{\boldsymbol {w}}_{s}\cdot \rho \left(\partial _{t}u_{s}+{\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{s}+{\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes u_{s}\right)\right)\,{\hbox{d}}\Omega \quad \quad }$ ${\displaystyle \int _{\Sigma \,\Omega ^{e}}{\boldsymbol {w}}_{s}\cdot \left(\nabla p_{s}-2\mu \nabla \cdot \left(\nabla ^{s}u_{s}-{\frac {1}{3}}\left(\nabla \cdot u_{s}\right){\boldsymbol {I}}\right)\right)\,{\hbox{d}}\Omega =}$ (2.33) ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{s}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega -\int _{\Sigma \,\Omega ^{e}}{\boldsymbol {w}}_{s}\cdot \rho \left(\partial _{t}u_{h}+{\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{h}+{\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes u_{h}\right)\right)\,{\hbox{d}}\Omega \quad \quad }$ (2.34) ${\displaystyle -\int _{\Sigma \,\Omega ^{e}}{\boldsymbol {w}}_{s}\cdot \left(\nabla p_{h}-2\mu \nabla \cdot \left(\nabla ^{s}u_{h}-{\frac {1}{3}}\left(\nabla \cdot u_{h}\right){\boldsymbol {I}}\right)\right)\,{\hbox{d}}\Omega }$ (2.35) ${\displaystyle \int _{\Omega }q_{s}\,\nabla \cdot u_{s}\,{\hbox{d}}\Omega =-\int _{\Omega }q_{s}\,\nabla \cdot u_{h}\,{\hbox{d}}\Omega }$ (2.36)

Eqs. 2.35 and 2.36 can be understood as the ${\textstyle L^{2}}$ projection onto the space of small scales ${\textstyle W_{s}\times Q_{s}}$ of a differential equation, in the same sense as Eqs. 2.10 and 2.11 represent the ${\textstyle L^{2}}$ projection onto ${\textstyle W\times Q}$ of the Navier-Stokes equations. Moreover, observing the right hand side terms of Eq. 2.35, we can see that it represents the projection of the strong-form linear momentum equation applied to the large scale part of the solution ${\textstyle u_{h}}$,${\textstyle p_{h}}$. The same can be said about the right hand side of Eq. 2.35, which corresponds to the mass conservation equation applied to the large scale velocity.

Denoting the projection onto the small scale spaces ${\textstyle W_{s}}$, ${\textstyle Q_{s}}$ with ${\textstyle \Pi _{V_{s}}\left(\cdot \right)}$ and ${\textstyle \Pi _{Q_{s}}\left(\cdot \right)}$, respectively, Eqs. 2.35 and 2.36 equation can be stated as

 ${\displaystyle \Pi _{V_{s}}\left(\rho \partial _{t}u_{s}+\rho {\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{s}+\rho {\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes u_{s}\right)\right.}$ ${\displaystyle \left.+\nabla p_{s}-2\mu \nabla \cdot \left(\nabla ^{s}u_{s}-{\frac {1}{3}}\left(\nabla \cdot u_{s}\right){\boldsymbol {I}}\right)\right)=}$ (2.37) ${\displaystyle \Pi _{V_{s}}\left({\boldsymbol {f}}-\rho \partial _{t}u_{h}-\rho {\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{h}-\rho {\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes u_{h}\right)\right.}$ (2.38) ${\displaystyle \left.-\nabla p_{s}+2\mu \nabla \cdot \left(\nabla ^{s}u_{h}-{\frac {1}{3}}\left(\nabla \cdot u_{h}\right){\boldsymbol {I}}\right)\right)}$ (2.39) ${\displaystyle \Pi _{Q_{s}}\left(\nabla \cdot u_{s}\right)=\Pi _{Q_{s}}\left(-\nabla \cdot u_{h}\right)}$ (2.40)

Since Eqs. 2.39 and 2.40 must hold for all admissible small-scale test functions ${\textstyle {\boldsymbol {w}}_{s}}$ and ${\textstyle q_{s}}$, they are equivalent to imposing that the small scale variables ${\textstyle u_{s}}$, ${\textstyle p_{s}}$ verify the following problem in each element ${\textstyle \Omega ^{e}}$:

 ${\displaystyle \rho \partial _{t}u_{s}+\rho {\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{s}+\rho {\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes u_{s}\right)+\nabla p_{s}}$ ${\displaystyle -2\mu \nabla \cdot \left(\nabla ^{s}u_{s}-{\frac {1}{3}}\left(\nabla \cdot u_{s}\right){\boldsymbol {I}}\right)={{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}\;{\hbox{ in}}\;\Omega ^{e}\times \left[0,T\right)}$ (2.41) ${\displaystyle \nabla \cdot u_{s}={R^{c}\left(u_{h}\right)}-\delta _{h}\;{\hbox{ in}}\;\Omega ^{e}\times \left[0,T\right)}$ (2.42)

where ${\textstyle {{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}}$ and ${\textstyle {R^{c}\left(u_{h}\right)}}$ represent the residual form of the Navier-Stokes equations applied to the large scale variables

 ${\displaystyle {{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}={\boldsymbol {f}}-\rho \partial _{t}u_{h}-\rho {\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{h}-\rho {\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes u_{h}\right)-\nabla p_{h}}$ ${\displaystyle +2\mu \nabla \cdot \left(\nabla ^{s}u_{h}-{\frac {1}{3}}\left(\nabla \cdot u_{h}\right){\boldsymbol {I}}\right)}$ (2.43) ${\displaystyle {R^{c}\left(u_{h}\right)}=-\nabla \cdot u_{h}}$ (2.44)

and ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ are chosen to enforce that right hand sides of Eqs. 2.41 and 2.42 belong to the corresponding small scale space.

As mentioned in the previous pages, the small scale spaces are infinite-dimensional and have to be approximated before the problem given by Eqs. 2.41 and 2.42 can be solved. In practice, the definition of an approximate small scale space corresponds to a choice of projectors ${\textstyle \Pi _{V_{s}}\left(\cdot \right)}$ and ${\textstyle \Pi _{Q_{s}}\left(\cdot \right)}$. The most straightforward possibility is to use the entire residual (without projecting to a particular space), which corresponds to considering operators ${\textstyle \Pi _{V_{s}}\left(\cdot \right)}$, ${\textstyle \Pi _{Q_{s}}\left(\cdot \right)}$ equal to the identity function or, equivalently, ${\textstyle {\boldsymbol {\xi }}_{h}={\boldsymbol {0}}}$ and ${\textstyle \delta _{h}=0}$. This formulation gives rise to the algebraic sub-grid scale (ASGS) method [39].

Another well-known choice consists in taking a small scale space that is orthogonal to the large scale space. If ${\textstyle \Pi _{V_{h}}\left(\cdot \right)}$ and ${\textstyle \Pi _{Q_{h}}\left(\cdot \right)}$ are the ${\textstyle L^{2}}$ projection onto the large scale spaces ${\textstyle W_{h}}$ and ${\textstyle Q_{h}}$ respectively, then the projections in Eqs. 2.39 and 2.40 are defined as ${\textstyle \Pi _{V_{s}}\left(\cdot \right)\approx \Pi _{V_{h}}^{\perp }\left(\cdot \right)}$, ${\textstyle \Pi _{Q_{s}}\left(\cdot \right)\approx \Pi _{Q_{h}}^{\perp }\left(\cdot \right)}$. In this case, ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ are chosen to subtract from the equation the part of the residuals that belongs to the finite element space, that is,

 ${\displaystyle {\boldsymbol {\xi }}_{h}=\Pi _{V_{h}}\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}\right)}$ (2.45) ${\displaystyle \delta _{h}=\Pi _{Q_{h}}\left({R^{c}\left(u_{h}\right)}\right)}$ (2.46)

This choice leads to the orthogonal sub-scale (OSS) method, presented in [40,30].

Note that, due to their definition, ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ 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 ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ requires knowledge of the finite element solutions ${\textstyle u_{h}}$, ${\textstyle p_{h}}$ and, as a result, it is coupled to the solution of the stabilized Navier-Stokes equations. In principle, this would double the number of nodal degrees of freedom of the problem. However, in practice, given that the Navier-Stokes problem is non-linear and has to be solved iteratively anyway, the projection problem can be implemented in an staggered way, updating the projections after each non-linear Navier-Stokes iteration.

The problem for the small scales, given by Eqs. 2.41 and 2.42, is not usually solved. Instead, it is approximated by an expression of the form

 ${\displaystyle \rho \partial u_{s}+{\frac {1}{{\boldsymbol {\tau }}_{u}}}\,u_{s}\approx {{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}\qquad {\frac {1}{\tau _{p}}}p_{s}\approx {R^{c}\left(u_{h}\right)}-\delta _{h}}$
(2.47)

where second order tensor ${\textstyle {\boldsymbol {\tau }}_{u}}$ and the scalar ${\textstyle \tau _{p}}$, 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 [30], where the parameters in Eq. 2.47 are designed to ensure that the ${\textstyle L^{2}}$ 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,21].

### 2.3.3 Quasi-static small scale models

While Eq. 2.47 represents the complete small scale model derived from the dynamic Navier-Stokes equations, many VMS formulations that can be found in the literature (see for example [39,41,21]) neglect the time variation of the velocity small scale model and define the small scales as

 ${\displaystyle u_{s}\approx {\boldsymbol {\tau }}_{u}\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi _{h}}}\right)\qquad p_{s}\approx \tau _{p}\left({R^{c}\left(u_{h}\right)}-\delta _{h}\right)}$
(2.48)

This choice corresponds to the assumption that the velocity small scales adapt automatically to the large scale residual. Following the nomenclature of [39], we refer to models based on Eq. 2.48 as quasi-static subscales, as opposed to dynamic subscale models, based on Eq. 2.47.

To complete the formulation, a definition for the stabilization parameters is needed. We follow the approach of [39], where the velocity stabilization parameter is taken to be a diagonal matrix ${\textstyle {\boldsymbol {\tau }}_{u}=\tau _{u}{\boldsymbol {I}}}$ and

 ${\displaystyle \tau _{u}=\left({\frac {c_{1}\mu }{h^{2}}}+{\frac {c_{2}\rho \left\Vert {\boldsymbol {a}}\right\Vert }{h}}\right)^{-1}}$ (2.49) ${\displaystyle \tau _{p}={\frac {h^{2}}{c_{1}\tau _{u}}}=\mu +{\frac {c_{2}\left\Vert {\boldsymbol {a}}\right\Vert h}{c_{1}}}}$ (2.50)

where ${\textstyle h}$ is a characteristic length of the element and ${\textstyle c_{1}}$, ${\textstyle c_{2}}$ are constants, which, for linear finite elements, are usually defined as ${\textstyle c_{1}=4}$, ${\textstyle c_{2}=2}$ (this is the case for example in [42] or [34]). However, the studies presented in [31] 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 [27], we have adopted ${\textstyle c_{1}=8}$, ${\textstyle c_{2}=2}$ for our tests.

As pointed out in [20], the use of quasi-static 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 [43] for the Stokes problem, where it is shown that it can be neutralized if the stabilization parameter satisfies the condition

 ${\displaystyle \delta t\geq C\tau _{u}}$
(2.51)

where ${\textstyle \delta t}$ is the time step and ${\textstyle C}$ is a constant. To avoid this instability, the stabilization parameter ${\textstyle \tau _{u}}$ can be replaced by the modified expression

 ${\displaystyle \tau _{t}=\left({\frac {\rho }{\delta t}}+{\frac {c_{1}\mu }{h^{2}}}+{\frac {c_{2}\rho \left\Vert {\boldsymbol {a}}\right\Vert }{h}}\right)^{-1}}$
(2.52)

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.

### 2.3.4 Dynamic small-scale models

If, instead of the quasi-static model of Eq. 2.48, we use the dynamic model of 2.47, the time evolution of the velocity small scale ${\textstyle u_{s}}$ has to be taken into account. This was achieved in [30,20] by introducing a time discretization for the small scale acceleration, resulting in the time-discrete small scale model

 ${\displaystyle \rho \,{\frac {u_{s}^{n+\theta }-u_{s}^{n}}{\delta t}}+{\frac {1}{\tau _{u}}}u_{s}^{n+\theta }=\left.\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}\right)\right|_{n+\theta }}$
(2.53)

Taking ${\textstyle \theta =1}$, which corresponds to a backward Euler time scheme, we can write a closed expression for the small scale velocity, given by

 ${\displaystyle \left({\frac {\rho }{\delta t}}+{\frac {1}{\tau _{u}}}\right)u_{s}^{n+1}=\left.\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}\right)\right|_{n+1}+{\frac {\rho }{\delta t}}u_{s}^{n}}$
(2.54)

As remarked in [20], the effective stabilization parameter in Eq. 2.54 is

 ${\displaystyle \left({\frac {\rho }{\delta t}}+{\frac {1}{\tau _{u}}}\right)^{-1}=\left({\frac {\rho }{\delta t}}+{\frac {c_{1}\mu }{h^{2}}}+{\frac {c_{2}\rho \left\Vert {\boldsymbol {a}}\right\Vert }{h}}\right)^{-1}}$
(2.55)

which is precisely the expression introduced as ${\textstyle \tau _{t}}$ in Eq. 2.52 and prevents the apparition of instabilities due to small time steps. However, unlike in quasi-static approximations, when the problem has a stationary solution, the dependency on the time step is eliminated as, in that case, ${\textstyle u_{s}^{n+1}=u_{s}^{n}}$ and the quasi-static model of Eq. 2.48 is recovered.

From the point of view of its implementation, the main difference between the quasi-static model of Eq. 2.48 and Eq. 2.54 is that the latter introduces the old value of the velocity small scale ${\textstyle u_{s}^{n}}$ in the model. This means that ${\textstyle u_{s}}$ has to be tracked in time. In practice, this implies evaluating and storing historical values for ${\textstyle u_{s}}$ on the integration points of the finite element mesh.

### 2.3.5 Complete equations

Once the small scale model is defined, it can be introduced in the large scale equations, given by Eqs. 2.27 and 2.28. The most general formulation that can be obtained from the material presented in the previous section is given by

Momentum equation

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{h}\cdot \rho \left(\partial _{t}u_{h}+{\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{h}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:\rho {\frac {1}{2}}\left({\boldsymbol {a}}\otimes u_{h}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}_{h}\,p_{h}\,{\hbox{d}}\Omega +\int _{\Omega }\nabla ^{s}{\boldsymbol {w}}_{h}:\,2\mu \left(\nabla ^{s}u_{h}-{\frac {1}{3}}\left(\nabla \cdot u_{h}\right){\boldsymbol {I}}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Omega }{\boldsymbol {w}}_{h}\cdot \rho \partial _{t}u_{s}\,{\hbox{d}}\Omega -\int _{\Sigma \,\Omega ^{e}}\rho \left({\boldsymbol {a}}\cdot \nabla {\boldsymbol {w}}_{h}\right)\tau _{t}\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}+{\frac {\rho }{\delta t}}u_{s}^{n}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}_{h}\tau _{p}\left({R^{c}\left(u_{h}\right)}-\delta _{h}\right)\,{\hbox{d}}\Omega =\int _{\Omega }{\boldsymbol {w}}_{h}{\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}_{h}\left({\boldsymbol {t}}-\rho {\frac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right)u_{h}\right)\,{\hbox{d}}\Omega }$
(2.56)

Mass conservation

 ${\displaystyle \int _{\Omega }q_{h}\nabla \cdot u_{h}\,{\hbox{d}}\Omega =\int _{\Sigma \,\Omega ^{e}}\nabla q_{h}\,\tau _{t}\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}+{\frac {\rho }{\delta t}}u_{s}^{n}\right)\,{\hbox{d}}\Omega }$
(2.57)

The first two rows of Eq. 2.56, 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. 2.56 represents the effect of pressure small scales.

Analogously, the left hand side of Eq. 2.57 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. 2.56 and 2.57 represent a general VMS formulation, they can be particularized to recover several well-known stabilized methods:

• Dynamic algebraic sub-grid scales [20]. The algebraic approximation to the small scales is characterized by using an identity projector to define the small scale space, which implies that the projection terms ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ are zero.
• Quasi-static algebraic sub-grid scales [39]. A quasi-static approximation to the small scales can be recovered by neglecting all terms involving either the small scale acceleration ${\textstyle \partial _{t}u_{s}}$ or the old small-scale velocities ${\textstyle u_{s}^{n}}$ and replacing ${\textstyle \tau _{t}}$ by the static stabilization parameter ${\textstyle \tau _{u}}$. Additionally, as the small scales are algebraic, projection terms ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ are also zero.
• Dynamic orthogonal subgrid-scales [30,20]. If the small scale space is assumed to be orthogonal to the large scale space, the integral involving ${\textstyle {\boldsymbol {w}}_{h}}$ and ${\textstyle \partial _{t}u_{s}}$ is zero, as it corresponds to the ${\textstyle L^{2}}$ product of two terms belonging to orthogonal spaces.
• Quasi-static orthogonal subgrid-scales [30] can be recovered from the original expression by neglecting all terms involving ${\textstyle \partial _{t}u_{s}}$ or ${\textstyle u_{s}^{n}}$ and replacing ${\textstyle \tau _{t}}$ by ${\textstyle \tau _{u}}$.

For OSS formulations, the projections ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ are defined as the ${\textstyle L^{2}}$ 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

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{h}\cdot {\boldsymbol {\xi }}_{h}\,{\hbox{d}}\Omega =\int _{\Omega }{\boldsymbol {w}}_{h}\cdot {{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}\,{\hbox{d}}\Omega }$ ${\displaystyle =\int _{\Omega }{\boldsymbol {w}}_{h}\left({\boldsymbol {f}}-\rho \partial _{t}u_{h}-\rho {\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{h}-\rho {\frac {1}{2}}\,\nabla \cdot \left({\boldsymbol {a}}\otimes u_{h}\right)\right)\,{\hbox{d}}\Omega }$ (2.58) ${\displaystyle +\int _{\Omega }{\boldsymbol {w}}_{h}\left(2\mu \nabla \cdot \left(\nabla ^{s}u_{h}-{\frac {1}{3}}\left(\nabla \cdot u_{h}\right){\boldsymbol {I}}\right)-\nabla p_{h}\right)\,{\hbox{d}}\Omega }$ (2.59) ${\displaystyle \int _{\Omega }q_{h}\delta _{h}\,{\hbox{d}}\Omega =\int _{\Omega }q_{h}{R^{c}\left(u_{h}\right)}\,{\hbox{d}}\Omega =-\int _{\Omega }q_{h}\nabla \cdot u_{h}\,{\hbox{d}}\Omega }$ (2.60)

The only question that remains to close the formulation is to provide a formal definition for the auxiliary convection velocity ${\textstyle {\boldsymbol {a}}}$, 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 ${\textstyle u_{h}}$ or the full velocity ${\textstyle u=u_{h}+u_{s}}$?

Both choices result in viable stabilized formulations. The classical approach is to use only the large scale velocity ${\textstyle u_{h}}$, as it corresponds to the finite element solution. However, from a theoretical point of view, using the full velocity ${\textstyle u_{h}+u_{s}}$ has interesting implications for turbulence modeling, which will be the main focus of Section 2.4.

## 2.4 VMS methods and Large Eddy Simulation

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

 ${\displaystyle {\overline {\boldsymbol {u}}}\left({\boldsymbol {x}},t\right)=\int _{-\Delta /2}^{\Delta /2}\int _{-\Delta /2}^{\Delta /2}\int _{-\Delta /2}^{\Delta /2}G\left({\boldsymbol {x}}-{\boldsymbol {\chi }}\right)u\left({\boldsymbol {\chi }},t\right)\,{\hbox{d}}\Omega }$
(2.61)

where ${\textstyle G\left({\boldsymbol {x}}-{\boldsymbol {\chi }}\right)}$ is a filter function defined on the interval ${\textstyle \left[-\Delta /2,\Delta /2\right]}$ and ${\textstyle \Delta }$ is known as the filter width.

Applying a filter function to Eqs. 2.3 and 2.4 allows us to write the filtered Navier-Stokes equations, given by

 ${\displaystyle \rho \,\partial _{t}{\overline {u}}+\rho \nabla \cdot \left({\overline {u}}\otimes {\overline {u}}\right)-\rho \nabla \cdot {\boldsymbol {\tau }}^{R}-\nabla \cdot \left(2\mu \nabla ^{s}{\overline {u}}\right)+\nabla {\overline {p}}={\overline {\boldsymbol {f}}}{\hbox{ in }}\Omega \times \left[0,T\right)}$ (2.62) ${\displaystyle \nabla \cdot {\overline {u}}=0{\hbox{ in }}\Omega \times \left[0,T\right)}$ (2.63)

where we have used the conservative form of the convective term and ${\textstyle {\boldsymbol {\tau }}^{R}}$ is the subgrid stress tensor, defined as

 ${\displaystyle {\boldsymbol {\tau }}^{R}={\overline {u\otimes u}}-{\overline {u}}\otimes {\overline {u}}}$
(2.64)

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 ${\textstyle {\overline {u\otimes u}}}$ cannot be obtained from the filtered velocity ${\textstyle {\overline {u}}}$. This means that Eqs. 2.62 and 2.63 do not represent a closed expression. However, if the filter width is chosen small enough that the filtered-out 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 (problem-independent) behavior. LES methods use this approach to motivate a model for ${\textstyle {\boldsymbol {\tau }}^{R}}$ and justify its introduction in the filtered equations, closing the formulation.

Introducing the subgrid velocity ${\textstyle {\tilde {u}}=u-{\overline {u}}}$, the subgrid stress tensor ${\textstyle {\boldsymbol {\tau }}^{R}}$ can be rewritten using the Leonard decomposition as

 ${\displaystyle {\boldsymbol {\tau }}^{R}={\boldsymbol {L}}+{\boldsymbol {C}}+{\boldsymbol {R}}}$
(2.65)

where each of the individual terms is defined as:

 ${\displaystyle {\boldsymbol {L}}={\overline {u\otimes u}}-{\overline {u}}\otimes {\overline {u}}{\hbox{ Leonard stress}}}$ (2.66) ${\displaystyle {\boldsymbol {C}}={\overline {{\tilde {u}}\otimes u}}+{\overline {u\otimes {\tilde {u}}}}{\hbox{ Cross stress}}}$ (2.67) ${\displaystyle {\boldsymbol {R}}={\overline {{\tilde {u}}\otimes {\tilde {u}}}}{\hbox{ Reynolds stress}}}$ (2.68)

The different terms in Eqs. 2.662.68 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 ${\textstyle L^{2}}$ projection to the finite element mesh. This projection to the mesh was introduced in writing the large scale equation, given by Eqs. 2.21 and 2.22, which is rewritten using the conservative form of the convective term to be consistent with the LES expression above, obtaining

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{h}\cdot \rho \left(\partial _{t}u_{h}+\partial _{t}u_{s}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\rho \nabla {\boldsymbol {w}}_{h}:\left(u_{h}+u_{s}\right)\otimes \left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Omega }\nabla ^{s}{\boldsymbol {w}}_{h}:2\mu \nabla ^{s}\left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}_{h}\left(p_{h}+p_{s}\right)\,{\hbox{d}}\Omega =\int _{\Omega }{\boldsymbol {w}}_{h}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega }$ (2.69) ${\displaystyle \int _{\Omega }q_{h}\nabla \cdot \left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega =0}$ (2.70)

where boundary terms and terms related to the trace of the viscous stresses have been omitted for clarity and the full velocity ${\textstyle u_{h}+u_{s}}$ has been used for both arguments of the convective term.

An analogy can be established between Eqs. 2.62 and 2.63, which represent the filtered Navier-Stokes equations and Eqs. 2.69 and 2.70, which expresses the projection of the Navier-Stokes equations to the finite element mesh. This was analyzed in [17,20], where it is remarked that the convective term in Eq. 2.69 can be expanded as

 ${\displaystyle \int _{\Omega }\rho \nabla {\boldsymbol {w}}_{h}:\left(u_{h}+u_{s}\right)\otimes \left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega =\int _{\Omega }\rho \nabla {\boldsymbol {w}}_{h}:u_{h}\otimes u_{h}\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Omega }\rho \nabla {\boldsymbol {w}}_{h}:u_{s}\otimes u_{h}\,{\hbox{d}}\Omega +\int _{\Omega }\rho \nabla {\boldsymbol {w}}_{h}:u_{h}\otimes u_{s}\,{\hbox{d}}\Omega +\int _{\Omega }\rho \nabla {\boldsymbol {w}}_{h}:u_{s}\otimes u_{s}\,{\hbox{d}}\Omega }$
(2.71)

The last three terms in Eq. 2.71 can be understood as a variational version of the LES subgrid stress tensor ${\textstyle {\boldsymbol {\tau }}^{R}}$. Ignoring the density, they can be rearranged as

 ${\displaystyle \int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{s}\otimes u_{h}\,{\hbox{d}}\Omega +\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{h}\otimes u_{s}\,{\hbox{d}}\Omega +\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{s}\otimes u_{s}\,{\hbox{d}}\Omega }$ ${\displaystyle {}=\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:\left(u_{h}+u_{s}\right)\otimes \left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{h}\otimes u_{h}\,{\hbox{d}}\Omega }$ ${\displaystyle {}=\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:{\boldsymbol {\tau }}_{VMS}^{R}\,{\hbox{d}}\Omega }$
(2.72)

where we introduced ${\textstyle {\boldsymbol {\tau }}_{VMS}^{R}}$ as a variational analogue of the LES subgrid stress tensor.

Furthermore, ${\textstyle {\boldsymbol {\tau }}_{VMS}^{R}}$ can be decomposed into Cross and Reynolds terms as:

 ${\displaystyle \int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{s}\otimes u_{h}\,{\hbox{d}}\Omega +\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{h}\otimes u_{s}\,{\hbox{d}}\Omega {\hbox{ Cross stress}}}$ (2.73) ${\displaystyle \int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{s}\otimes u_{s}\,{\hbox{d}}\Omega {\hbox{ Reynolds stress}}}$ (2.74)

while an analogue of the Leonard stress, representing the contribution of the resolved velocities ${\textstyle u_{h}}$ to the unresolved stresses, appears in the corresponding small scale equation

 ${\displaystyle \int _{\Omega }\nabla {\boldsymbol {w}}_{h}:u_{h}\otimes u_{h}\,{\hbox{d}}\Omega -\int _{\Omega }\nabla {\boldsymbol {w}}:u_{h}\otimes u_{h}\,{\hbox{d}}\Omega =}$ ${\displaystyle -\int _{\Omega }\nabla {\boldsymbol {w}}_{s}:u_{h}\otimes u_{h}\,{\hbox{d}}\Omega {\hbox{ Leonard stress}}}$
(2.75)

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 mesh-induced filtering has also been explored by the LES community, where it is known as implicit filtering. This is the basis of the Monotone Integrated Large-Eddy Simulation (MILES) approach of Boris et al. [44], 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 VMS-based LES, which will only be mentioned here, is based on a three–level scale separation (see for example [41,45] 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 [46,47].

### 2.4.1 The VMS kinetic energy balance

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 Q-ASGS and in [25,26] for OSS-based methods.

An energy balance for the original Navier-Stokes problem can be obtained by taking ${\textstyle {\boldsymbol {w}}=u}$ in the Galerkin weak form of the momentum equation, given by Eq. 2.13, which produces

 ${\displaystyle \int _{\Omega }u\cdot \rho \partial _{t}u\,{\hbox{d}}\Omega +\int _{\Omega }u\cdot {\big (}\rho \,u\cdot \nabla u{\big )}\,{\hbox{d}}\Omega +\int _{\Omega }2\mu \nabla ^{s}u:\nabla ^{s}u\,{\hbox{d}}\Omega }$ ${\displaystyle -\int _{\Omega }\nabla \cdot up\,{\hbox{d}}\Omega =\int _{\Omega }u\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}\cdot {\boldsymbol {t}}\,{\hbox{d}}\Omega }$
(2.76)

We are interested in obtaining a balance for the kinetic energy, ${\textstyle E=\rho u\cdot u/2}$. Using the fact that the full velocity is incompressible and the equality ${\textstyle u\,\partial u/\partial x=1/2\,\partial u^{2}/\partial x}$, Eq. 2.76 can be expressed as

 ${\displaystyle \underbrace {\int _{\Omega }\partial _{t}E\,{\hbox{d}}\Omega } _{I}+\underbrace {\int _{\Omega }u\cdot \nabla E\,{\hbox{d}}\Omega } _{II}=\underbrace {\int _{\Omega }u\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega } _{III}+\underbrace {\int _{\Gamma _{N}}u\cdot {\boldsymbol {t}}\,{\hbox{d}}\Omega } _{IV}-\underbrace {\int _{\Omega }2\mu \nabla ^{s}u:\nabla ^{s}u\,{\hbox{d}}\Omega } _{V}}$
(2.77)

Eq. 2.77 expresses the balance of total kinetic energy ${\textstyle E}$ in the domain, and the individual terms represent energy storage (${\textstyle I}$) and convection (${\textstyle II}$), the power exerted by external forces (${\textstyle III}$) and boundary tensions (${\textstyle IV}$) and finally viscous dissipation (${\textstyle V}$).

We can obtain an equivalent expression for the kinetic energy contained in the large scale motions, ${\textstyle E_{h}=\rho u_{h}\cdot u_{h}/2}$, by taking ${\textstyle {\boldsymbol {w}}_{h}=u_{h}}$ in Eq. 2.69, ${\textstyle q_{h}=p_{h}}$ in Eq. 2.70 and adding the two equations, obtaining:

 ${\displaystyle \int _{\Omega }u_{h}\cdot \rho \left(\partial _{t}u_{h}+\partial _{t}u_{s}\right)\,{\hbox{d}}\Omega +\int _{\Omega }u_{h}\cdot \rho \nabla \cdot \left(\left(u_{h}+u_{s}\right)\otimes \left(u_{h}+u_{s}\right)\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Omega }\nabla ^{s}u_{h}:2\mu \nabla ^{s}\left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla \cdot u_{h}\left(p_{h}+p_{s}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Omega }p_{h}\nabla \cdot \left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega =\int _{\Omega }u_{h}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega }$
(2.78)

Note that in Eq. 2.78, 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. 2.78 by parts and rearrange the convective term as follows:

 ${\displaystyle \int _{\Omega }u_{h}\cdot \rho \nabla \cdot \left(\left(u_{h}+u_{s}\right)\otimes \left(u_{h}+u_{s}\right)\right)\,{\hbox{d}}\Omega =}$ ${\displaystyle =\int _{\Omega }u_{h}\cdot \rho \nabla \cdot \left(\left(u_{h}+u_{s}\right)\otimes u_{h}\right)\,{\hbox{d}}\Omega +\int _{\Omega }u_{h}\cdot \rho \nabla \cdot \left(\left(u_{h}+u_{s}\right)\otimes u_{s}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle =\int _{\Omega }u_{h}\cdot \rho \left(u_{h}+u_{s}\right)\cdot \nabla u_{h}\,{\hbox{d}}\Omega -\int _{\Omega }\rho \nabla u_{h}:\left(u_{h}+u_{s}\right)\otimes u_{s}\,{\hbox{d}}\Omega }$ ${\displaystyle =\int _{\Omega }\left(u_{h}+u_{s}\right)\cdot \nabla E_{h}\,{\hbox{d}}\Omega -\int _{\Omega }u_{s}\rho \left(u_{h}+u_{s}\right)\cdot \nabla u_{h}\,{\hbox{d}}\Omega }$

where we have used the fact that the exact velocity ${\textstyle u_{h}+u_{s}}$ is divergence free.

With this, Eq. 2.78 can be restated as

 ${\displaystyle \underbrace {\int _{\Omega }\partial _{t}E_{h}\,{\hbox{d}}\Omega } _{I}+\underbrace {\int _{\Omega }u\cdot \nabla E_{h}\,{\hbox{d}}\Omega } _{II}=\underbrace {\int _{\Omega }u_{h}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega } _{III}-\underbrace {\int _{\Omega }2\mu \nabla ^{s}u_{h}:\nabla ^{s}u_{h}\,{\hbox{d}}\Omega } _{IV}+\underbrace {\int _{\Omega }p_{s}\nabla \cdot u_{h}\,{\hbox{d}}\Omega } _{V}}$ ${\displaystyle \underbrace {-\int _{\Omega }u_{h}\cdot \rho \partial _{t}u_{s}\,{\hbox{d}}\Omega +\int _{\Omega }u_{s}\cdot \left(\rho u\cdot \nabla u_{h}+\nabla p_{h}+\nabla \cdot \left(2\mu \nabla ^{s}u_{h}\right)\right)\,{\hbox{d}}\Omega } _{VI}}$
(2.79)

Analogously to the complete energy balance, Eq. 2.79 can be understood as a balance for the kinetic energy associated to large scale motions. Terms ${\textstyle I}$ and ${\textstyle II}$ represent the storage and convection of large scale kinetic energy, while term ${\textstyle III}$ represents the power exerted by the external forces on large scale motions. Term ${\textstyle IV}$ 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 Navier-Stokes equations.

Finally, we define the residual kinetic energy as ${\textstyle k_{r}=E-E_{h}}$. A balance statement for ${\textstyle k_{r}}$ can be obtained by subtracting Eq. 2.79 from Eq. 2.77, which results in

 ${\displaystyle \underbrace {\int _{\Omega }\partial _{t}k_{r}\,{\hbox{d}}\Omega } _{I}+\underbrace {\int _{\Omega }u\cdot \nabla k_{r}\,{\hbox{d}}\Omega } _{II}=\underbrace {\int _{\Omega }u_{s}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega } _{III}}$ ${\displaystyle -\underbrace {\int _{\Omega }2\mu \nabla ^{s}u:\nabla ^{s}u\,{\hbox{d}}\Omega } _{IV}+\underbrace {\int _{\Omega }2\mu \nabla ^{s}u_{h}:\nabla ^{s}u_{h}\,{\hbox{d}}\Omega } _{V}-\underbrace {\int _{\Omega }p_{s}\nabla \cdot u_{h}\,{\hbox{d}}\Omega } _{VI}}$ ${\displaystyle \underbrace {+\int _{\Omega }u_{h}\cdot \rho \partial _{t}u_{s}\,{\hbox{d}}\Omega -\int _{\Omega }u_{s}\cdot \left(\rho u\cdot \nabla u_{h}+\nabla p_{h}+\nabla \cdot \left(2\mu \nabla ^{s}u_{h}\right)\right)\,{\hbox{d}}\Omega } _{VII}}$
(2.80)

In Eq. 2.80, terms ${\textstyle I}$ and ${\textstyle II}$ represent the storage and convection of residual kinetic energy, while term ${\textstyle III}$ 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 (${\textstyle IV}$) and the large scale viscous dissipation (${\textstyle V}$), which was already accounted for in the large scale energy balance of Eq. 2.79. Again, we remark that, in practice, term ${\textstyle V}$ is expected to be negligible in comparison to term ${\textstyle IV}$, since viscous dissipation occurs predominantly for motions in the range of the Kolmogorov length scale, while ${\textstyle u_{h}}$ will contain only motions on a much larger scale ${\textstyle h}$, lying on the inertial subrange. Finally, terms ${\textstyle VI}$ and ${\textstyle VII}$ 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. 2.79.

Terms ${\textstyle VI}$ and ${\textstyle VII}$ can be modified by noting that, since ${\textstyle p_{s}\in Q_{s}}$ and ${\textstyle u_{s}\in V_{s}}$,

 ${\displaystyle \int _{\Omega }p_{s}\nabla \cdot u_{h}\,{\hbox{d}}\Omega =\int _{\Omega }p_{s}\Pi _{Q_{s}}\left(\nabla \cdot u_{h}\right)\,{\hbox{d}}\Omega }$ (2.81) ${\displaystyle \int _{\Omega }u_{s}\cdot \left(\rho u\cdot \nabla u_{h}+\nabla p_{h}+\nabla \cdot \left(2\mu \nabla ^{s}u_{h}\right)\right)\,{\hbox{d}}\Omega =}$ ${\displaystyle \int _{\Omega }u_{s}\cdot \Pi _{V_{s}}{\big (}\rho u\cdot \nabla u_{h}+\nabla p_{h}+\nabla \cdot \left(2\mu \nabla ^{s}u_{h}\right){\big )}\,{\hbox{d}}\Omega }$ (2.82)

The last step to obtain the residual kinetic energy balance is to introduce the small scale models for velocity and pressure in Eqs. 2.81 and 2.82. Noting that ${\textstyle \left(2\mu \nabla ^{s}u_{h}\right)}$ represents the large scale viscous dissipation and will be negligible for high Reynolds numbers, we can write

 ${\displaystyle \int _{\Omega }u_{s}\cdot \Pi _{V_{s}}\left(\rho u\cdot \nabla u_{h}+\nabla p_{h}\right)\,{\hbox{d}}\Omega +\int _{\Omega }p_{s}\Pi _{Q_{s}}\left(\nabla \cdot u_{h}\right)\,{\hbox{d}}\Omega =}$ ${\displaystyle \underbrace {\int _{\Omega }\tau _{u}\Pi _{V_{s}}\left({\boldsymbol {f}}\right)\cdot \Pi _{V_{s}}\left(\rho u\cdot \nabla u_{h}+\nabla p_{h}\right)\,{\hbox{d}}\Omega } _{I}-\underbrace {\int _{\Omega }\tau _{u}{\big \vert }\Pi _{V_{s}}\left(\rho u\cdot \nabla u_{h}+\nabla p_{h}\right){\big \vert }^{2}\,{\hbox{d}}\Omega } _{II}}$ ${\displaystyle {}+\underbrace {\int _{\Omega }\tau _{u}\rho \partial _{t}u_{s}\cdot \Pi _{V_{s}}\left(\rho u\cdot \nabla u_{h}+\nabla p_{h}\right)\,{\hbox{d}}\Omega } _{III}-\underbrace {\int _{\Omega }\tau _{p}{\big \vert }\Pi _{Q_{s}}\left(\nabla \cdot u_{h}\right){\big \vert }^{2}\,{\hbox{d}}\Omega } _{IV}}$
(2.83)

It is clear that terms ${\textstyle II}$ and ${\textstyle IV}$, since the stabilization parameters ${\textstyle \tau _{u}}$ and ${\textstyle \tau _{p}}$ were defined as strictly positive. They play the role of energy sinks in the large scale energy balance of Eq. 2.79, while acting as sources on the residual energy equation Eq. 2.80. Term ${\textstyle I}$ is problem dependent but, if an OSS small scale model is used, it will vanish unless the external forces have a high-frequency (small scale) component. Finally, term ${\textstyle III}$ only exists for dynamic small scale models. This fact was used in [25] 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 [26] 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 filter-based LES formulations. If the kinetic energy associated to the filtered velocity field is defined as ${\textstyle {\overline {E}}=\rho {\overline {u}}\cdot {\overline {u}}/2}$, the balance for ${\textstyle {\overline {E}}}$ can be obtained by multiplying the filtered linear momentum equation Eq. 2.62 by ${\textstyle {\overline {u}}}$ and integrating over the fluid domain (see [2]), resulting in

 ${\displaystyle \int _{\Omega }\partial _{t}{\overline {E}}\,{\hbox{d}}\Omega +\int _{\Omega }{\overline {u}}\cdot \nabla {\overline {E}}\,{\hbox{d}}\Omega =\int _{\Omega }{\overline {u}}\cdot {\overline {\boldsymbol {f}}}\,{\hbox{d}}\Omega -\int _{\Omega }2\mu \nabla ^{s}{\overline {u}}:\nabla ^{s}{\overline {u}}\,{\hbox{d}}\Omega +\int _{\Omega }\nabla ^{s}{\overline {u}}:{\boldsymbol {\tau }}^{R}\,{\hbox{d}}\Omega }$
(2.84)

where boundary fluxes have been omitted.

We can introduce the definition of the VMS subgrid stress tensor ${\textstyle {\boldsymbol {\tau }}_{VMS}^{R}}$ in Eq. 2.78, obtaining

 ${\displaystyle \int _{\Omega }u_{h}\cdot \rho \left(\partial _{t}u_{h}+\partial _{t}u_{s}\right)\,{\hbox{d}}\Omega +\int _{\Omega }\rho u_{h}\nabla \cdot \left(u_{h}\otimes u_{h}\right)\,{\hbox{d}}\Omega \qquad \qquad }$ ${\displaystyle -\int _{\Omega }\rho \nabla u_{h}:{\boldsymbol {\tau }}_{VMS}^{R}\,{\hbox{d}}\Omega +\int _{\Omega }\nabla ^{s}u_{h}:2\mu \nabla ^{s}\left(u_{h}+u_{s}\right)\,{\hbox{d}}\Omega \qquad }$ ${\displaystyle -\int _{\Omega }\nabla \cdot u_{h}\left(p_{h}+p_{s}\right)+\int _{\Omega }p_{h}\nabla \cdot \left(u_{h}+u_{s}\right)=\int _{\Omega }u_{h}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega }$
(2.85)

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

 ${\displaystyle \underbrace {\int _{\Omega }\partial _{t}E_{h}\,{\hbox{d}}\Omega } _{I}+\underbrace {\int _{\Omega }u_{h}\cdot \nabla E_{h}\,{\hbox{d}}\Omega } _{II}=\underbrace {\int _{\Omega }u_{h}\cdot {\boldsymbol {f}}\,{\hbox{d}}\Omega } _{III}-\underbrace {\int _{\Omega }2\mu \nabla ^{s}u_{h}:\nabla ^{s}u_{h}\,{\hbox{d}}\Omega } _{IV}+\underbrace {\int _{\Omega }p_{s}\nabla \cdot u_{h}\,{\hbox{d}}\Omega } _{V}}$ ${\displaystyle \underbrace {-\int _{\Omega }u_{h}\cdot \rho \partial _{t}u_{s}\,{\hbox{d}}\Omega -\int _{\Omega }u_{s}\cdot \left(\nabla p_{h}-\nabla \cdot \left(2\mu \nabla ^{s}u_{h}\right)\right)\,{\hbox{d}}\Omega } _{VI}+\underbrace {\int _{\Omega }\rho \nabla ^{s}u_{h}:{\boldsymbol {\tau }}_{VMS}^{R}\,{\hbox{d}}\Omega } _{VII}}$
(2.86)

where we have used the fact that ${\textstyle {\boldsymbol {\tau }}_{VMS}^{R}}$ is a symmetric tensor to write ${\textstyle \nabla u_{h}:{\boldsymbol {\tau }}_{VMS}^{R}=\nabla ^{s}u_{h}:{\boldsymbol {\tau }}_{VMS}^{R}}$. Terms ${\textstyle I}$ to ${\textstyle VI}$ have an analogous interpretation to their counterparts in Eq. 2.79, but here we have obtained an additional term, ${\textstyle VII}$, which explicitly represents the contribution of the residual subgrid stresses in the energy transfer. Comparing Eq. 2.86 to Eq. 2.84, 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 ${\textstyle V}$ and ${\textstyle VI}$, compared to filter-based LES.

## 2.5 Discrete problem

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 ${\textstyle \Omega _{h}}$ for the problem domain ${\textstyle \Omega }$. Given the discrete domain ${\textstyle \Omega _{h}}$, the large scale interpolation spaces ${\textstyle V_{h}}$ and ${\textstyle Q_{h}}$ can be identified with the standard finite element interpolation functions and the large scale part of the solution, ${\textstyle u_{h}}$ and ${\textstyle p_{h}}$, can be represented using a finite element interpolation as

 ${\displaystyle u_{h}=\sum _{a}^{n_{n}}{\boldsymbol {N}}_{a}\left({\boldsymbol {x}}\right){\boldsymbol {u}}_{a}\qquad p_{h}=\sum _{a}^{n_{n}}N_{a}\left({\boldsymbol {x}}\right)p_{a}}$
(2.87)

where ${\textstyle n_{n}}$ represents the number of nodes in the finite element mesh, ${\textstyle {\boldsymbol {u}}_{a}}$ and ${\textstyle p_{a}}$ are the nodal values of the large scale variables ${\textstyle u_{h}}$ and ${\textstyle p_{h}}$ respectively, ${\textstyle N_{a}}$ represents the standard finite element basis function associated to node ${\textstyle a}$ and ${\textstyle {\boldsymbol {N}}_{a}}$ its counterpart for vectorial variables, given by

 ${\displaystyle {\boldsymbol {N}}_{a}={\begin{bmatrix}N_{a}&0&0\\0&N_{a}&0\\0&0&N_{a}\end{bmatrix}}}$

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

 ${\displaystyle \left(\nabla N_{a}\right)^{T}={\begin{bmatrix}{\dfrac {\partial N_{a}}{\partial x}}&{\dfrac {\partial N_{a}}{\partial y}}&{\dfrac {\partial N_{a}}{\partial z}}\end{bmatrix}}\qquad \nabla \cdot {\boldsymbol {N}}_{a}={\begin{bmatrix}{\dfrac {\partial N_{a}}{\partial x}}&{\dfrac {\partial N_{a}}{\partial y}}&{\dfrac {\partial N_{a}}{\partial z}}\end{bmatrix}}}$ ${\displaystyle \left(\nabla {\boldsymbol {N}}_{a}\right)^{T}=\left[{\begin{array}{ccccccccc}{\dfrac {\partial N_{a}}{\partial x}}&0&0&{\dfrac {\partial N_{a}}{\partial x}}&0&0&{\dfrac {\partial N_{a}}{\partial x}}&0&0\\[1.1em]0&{\dfrac {\partial N_{a}}{\partial y}}&0&0&{\dfrac {\partial N_{a}}{\partial y}}&0&0&{\dfrac {\partial N_{a}}{\partial y}}&0\\[1.1em]0&0&{\dfrac {\partial N_{a}}{\partial z}}&0&0&{\dfrac {\partial N_{a}}{\partial z}}&0&0&{\dfrac {\partial N_{a}}{\partial z}}\end{array}}\right]}$
(2.88)

We also introduce the following operator to describe convection

 ${\displaystyle {\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{a}={\begin{bmatrix}{\boldsymbol {a}}\cdot \nabla N_{a}&0&0\\0&{\boldsymbol {a}}\cdot \nabla N_{a}&0\\0&0&{\boldsymbol {a}}\cdot \nabla N_{a}\end{bmatrix}}}$
(2.89)

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

 ${\displaystyle {\boldsymbol {B}}_{a}^{T}={\begin{bmatrix}{\dfrac {\partial N_{a}}{\partial x}}&0&0&{\dfrac {\partial N_{a}}{\partial y}}&0&{\dfrac {\partial N_{a}}{\partial z}}\\[1.1em]0&{\dfrac {\partial N_{a}}{\partial y}}&0&{\dfrac {\partial N_{a}}{\partial x}}&{\dfrac {\partial N_{a}}{\partial z}}&0\\[1.1em]0&0&{\dfrac {\partial N_{a}}{\partial z}}&0&{\dfrac {\partial N_{a}}{\partial y}}&{\dfrac {\partial N_{a}}{\partial x}}\end{bmatrix}}}$
(2.90)

Additionally, we define ${\textstyle {\boldsymbol {U}}}$, ${\textstyle {\boldsymbol {\dot {U}}}}$ and ${\textstyle {\boldsymbol {P}}}$ as the vectors of nodal values of large scale velocity ${\textstyle u_{h}}$, acceleration ${\textstyle \partial _{t}u_{h}}$ and pressure ${\textstyle p_{h}}$, respectively.

We introduce the finite element discretization of Eq. 2.87 in the variational formulation given by Eqs. 2.56 and 2.57 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 ${\textstyle \nabla \cdot \nabla ^{s}{\boldsymbol {w}}_{h}}$ in Eq. 2.56 or the strong-form viscous term that appears in the residual ${\textstyle {{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}}$ introduced in Eq. 2.43 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.

### 2.5.1 Quasi-static ASGS formulation

The quasi-static ASGS formulation is in some sense the classical VMS formulation for the Navier-Stokes equations. It is also relatively simple, as it does not involve dynamic terms or projections, so we will present it first.

Starting from Eqs. 2.56 and 2.57, we can neglect all terms involving the projections ${\textstyle {\boldsymbol {\xi }}_{h}}$ or ${\textstyle \delta _{h}}$. As the small scales are described using Eq. 2.48, terms involving ${\textstyle u_{s}^{n}}$ or ${\textstyle \partial _{t}u_{s}}$ in these equations can be ignored and the momentum stabilization parameter is ${\textstyle \tau _{u}}$, given by Eq. 2.49. After these simplifications, we introduce the finite element interpolation of Eq. 2.87 to describe the problem variables ${\textstyle u_{h}}$, ${\textstyle p_{h}}$. Testing against each nodal basis function in turn we obtain a system of equations that can be expressed in matrix form as

 ${\displaystyle {\big [}{\boldsymbol {M}}+{\boldsymbol {S}}_{m}\left(\tau _{u},{\boldsymbol {a}}\right){\big ]}{\boldsymbol {\dot {U}}}+{\big [}{\boldsymbol {C}}\left({\boldsymbol {a}}\right)+{\boldsymbol {K}}+{\boldsymbol {S}}_{u}\left(\tau _{u},{\boldsymbol {a}}\right)+{\boldsymbol {H}}_{u}\left(\tau _{p}\right){\big ]}{\boldsymbol {U}}}$ ${\displaystyle +{\big [}{\boldsymbol {G}}+{\boldsymbol {S}}_{p}\left(\tau _{u},{\boldsymbol {a}}\right){\big ]}{\boldsymbol {P}}={\boldsymbol {F}}+{\boldsymbol {T}}+{\boldsymbol {S}}_{f}\left(\tau _{u},{\boldsymbol {a}}\right)}$
(2.91)

 ${\displaystyle {\boldsymbol {Q}}_{m}\left(\tau _{u}\right){\boldsymbol {\dot {U}}}+{\big [}{\boldsymbol {D}}+{\boldsymbol {Q}}_{u}\left(\tau _{u},{\boldsymbol {a}}\right){\big ]}{\boldsymbol {U}}+{\boldsymbol {Q}}_{p}\left(\tau _{u}\right){\boldsymbol {P}}={\boldsymbol {Q}}_{f}\left(\tau _{u}\right)}$
(2.92)

Again, we have used a generic convection velocity ${\textstyle {\boldsymbol {a}}}$ in all terms that have a non-linear dependence of velocity. Doing so, we leave open the possibility of using either the full velocity ${\textstyle {\boldsymbol {a}}=u_{h}+u_{s}}$ or the only large scale part ${\textstyle {\boldsymbol {a}}=u_{h}}$. Note that, for linear finite elements, the latter choice is equivalent to the Galerkin-Least Squares (GLS) method [37].

The different matrices in Eqs. 2.91 and 2.92 represent the discrete version of the operators in Eqs. 2.56 and 2.57 and can be built from the assembly of elemental contributions. In general, matrix ${\textstyle {\boldsymbol {A}}}$ is constructed by the finite element assembly of elemental matrices of the form ${\textstyle {\boldsymbol {A}}^{e}}$. For a finite element with ${\textstyle N}$ nodes, ${\textstyle {\boldsymbol {A}}^{e}}$ can be defined using ${\textstyle N\times N}$ blocks ${\textstyle {\boldsymbol {A}}_{ab}^{e}}$, where ${\textstyle a}$ and ${\textstyle b}$ 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

 ${\displaystyle {\boldsymbol {M}}_{ab}^{e}=\int _{\Omega _{e}}\rho {\boldsymbol {N}}_{a}^{T}{\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.93) ${\displaystyle {\boldsymbol {C}}\left({\boldsymbol {a}}\right)_{ab}^{e}=\int _{\Omega _{e}}\rho {\dfrac {1}{2}}\left({\boldsymbol {N}}_{a}^{T}\,{\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{b}-\left({\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{a}\right)^{T}{\boldsymbol {N}}_{b}\right)\,{\hbox{d}}\Omega }$ (2.94) ${\displaystyle {}+\int _{\Gamma _{N}}{\boldsymbol {N}}_{a}^{T}\,\rho {\dfrac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right){\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.95) ${\displaystyle {\boldsymbol {K}}_{ab}^{e}=\int _{\Omega _{e}}{\boldsymbol {B}}_{a}^{T}{\boldsymbol {C}}_{\mu }{\boldsymbol {B}}_{b}\,{\hbox{d}}\Omega }$ (2.96) ${\displaystyle {\boldsymbol {G}}_{ab}^{e}=-\int _{\Omega _{e}}\left(\nabla \cdot {\boldsymbol {N}}_{a}\right)^{T}N_{b}\,{\hbox{d}}\Omega }$ (2.97) ${\displaystyle {\boldsymbol {D}}_{ab}^{e}=\int _{\Omega _{e}}N_{a}\,\nabla \cdot {\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega =-{\big (}{\boldsymbol {G}}_{ba}^{e}{\big )}^{T}}$ (2.98) ${\displaystyle {\boldsymbol {F}}_{a}^{e}=\int _{\Omega _{e}}{\boldsymbol {N}}_{a}^{T}{\boldsymbol {f}}\,{\hbox{d}}\Omega }$ (2.99) ${\displaystyle {\boldsymbol {T}}_{a}^{e}=\int _{\Gamma _{N}}{\boldsymbol {N}}_{a}^{T}{\boldsymbol {t}}\,{\hbox{d}}\Omega }$ (2.100)

When defining the viscous matrix ${\textstyle {\boldsymbol {K}}_{ab}^{e}}$ in Eq. 2.96 we introduced the constitutive matrix ${\textstyle {\boldsymbol {C}}_{\mu }}$, which is given by

 ${\displaystyle {\boldsymbol {C}}_{\mu }={\begin{bmatrix}4\mu /3&-2\mu /3&-2\mu /3&0&0&0\\-2\mu /3&4\mu /3&-2\mu /3&0&0&0\\-2\mu /3&-2\mu /3&4\mu /3&0&0&0\\0&0&0&\mu &0&0\\0&0&0&0&\mu &0\\0&0&0&0&0&\mu \end{bmatrix}}}$
(2.101)

In the same way, the stabilization terms in the Q-ASGS give rise to additional elemental matrices, expressed here as

 ${\displaystyle {\boldsymbol {S}}_{m}\left(\tau _{u},{\boldsymbol {a}}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\rho \,{\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,{\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.102) ${\displaystyle {\boldsymbol {S}}_{u}\left(\tau _{u},{\boldsymbol {a}}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\rho \,{\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,\rho \,{\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.103) ${\displaystyle {\boldsymbol {S}}_{p}\left(\tau _{u},{\boldsymbol {a}}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\rho \,{\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,\nabla N_{b}\,{\hbox{d}}\Omega }$ (2.104) ${\displaystyle {\boldsymbol {S}}_{f}\left(\tau _{u},{\boldsymbol {a}}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\rho \,{\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,{\boldsymbol {f}}\,{\hbox{d}}\Omega }$ (2.105) ${\displaystyle {\boldsymbol {Q}}_{m}\left(\tau _{u}\right)_{ab}^{e}=\int _{\Omega }\left(\nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,{\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.106) ${\displaystyle {\boldsymbol {Q}}_{u}\left(\tau _{u},{\boldsymbol {a}}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,\rho \,{\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega ={\big (}{\boldsymbol {S}}_{p}\left(\tau _{u},{\boldsymbol {a}}\right)_{ba}^{e}{\big )}^{T}}$ (2.107) ${\displaystyle {\boldsymbol {Q}}_{p}\left(\tau _{u}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,\nabla {\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.108) ${\displaystyle {\boldsymbol {Q}}_{f}\left(\tau _{u}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\nabla {\boldsymbol {N}}_{a}\right)^{T}\tau _{u}\,{\boldsymbol {f}}\,{\hbox{d}}\Omega }$ (2.109) ${\displaystyle {\boldsymbol {H}}_{u}\left(\tau _{p}\right)_{ab}^{e}=\int _{\Omega _{e}}\left(\nabla \cdot {\boldsymbol {N}}_{a}\right)^{T}\tau _{p}\,\nabla \cdot {\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.110)

### 2.5.2 Quasi-static OSS formulation

The next variant to be presented is the quasi-static OSS formulation. Compared to the Q-ASGS formulation, OSS is characterized by the inclusion of the projections ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$, 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 ${\textstyle {\boldsymbol {a}}=u_{h}+u_{s}}$ or ${\textstyle {\boldsymbol {a}}=u}$ for convection and the corresponding stabilization terms. The matrix form of the Q-OSS formulation can be expressed as

 ${\displaystyle {\big [}{\boldsymbol {M}}+{\boldsymbol {S}}_{m}\left(\tau _{u},{\boldsymbol {a}}\right){\big ]}{\boldsymbol {\dot {U}}}+{\big [}{\boldsymbol {C}}\left({\boldsymbol {a}}\right)+{\boldsymbol {K}}+{\boldsymbol {S}}_{u}\left(\tau _{u},{\boldsymbol {a}}\right)+{\boldsymbol {H}}_{u}\left(\tau _{p}\right){\big ]}{\boldsymbol {U}}}$ ${\displaystyle +{\big [}{\boldsymbol {G}}+{\boldsymbol {S}}_{p}\left(\tau _{u},{\boldsymbol {a}}\right){\big ]}{\boldsymbol {P}}={\boldsymbol {F}}+{\boldsymbol {T}}+{\boldsymbol {S}}_{f}\left(\tau _{u},{\boldsymbol {a}}\right)-{\boldsymbol {S}}_{\Pi }\left(\tau _{u},{\boldsymbol {a}}\right)-{\boldsymbol {H}}_{\Pi }\left(\tau _{u}\right)}$
(2.111)

 ${\displaystyle {\boldsymbol {Q}}_{m}\left(\tau _{u}\right){\boldsymbol {\dot {U}}}+{\big [}{\boldsymbol {D}}+{\boldsymbol {Q}}_{u}\left(\tau _{u},{\boldsymbol {a}}\right){\big ]}{\boldsymbol {U}}+{\boldsymbol {Q}}_{p}\left(\tau _{u}\right){\boldsymbol {P}}={\boldsymbol {Q}}_{f}\left(\tau _{u}\right)-{\boldsymbol {Q}}_{\Pi }\left(\tau _{u}\right)}$
(2.112)

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

 ${\displaystyle {\boldsymbol {S}}_{\Pi }\left(\tau _{u},{\boldsymbol {a}}\right)_{a}^{e}=\int _{\Omega _{e}}\left(\rho {\boldsymbol {a}}\cdot \nabla {\boldsymbol {N}}_{a}\right)^{T}\cdot \tau _{u}\,{\boldsymbol {\xi }}_{h}\,{\hbox{d}}\Omega }$ (2.113) ${\displaystyle {\boldsymbol {Q}}_{\Pi }\left(\tau _{u}\right)_{a}^{e}=\int _{\Omega _{e}}\left(\nabla N_{a}\right)^{T}\cdot \tau _{u}\,{\boldsymbol {\xi }}_{h}\,{\hbox{d}}\Omega }$ (2.114) ${\displaystyle {\boldsymbol {H}}_{\Pi }\left(\tau _{p}\right)_{a}^{e}=\int _{\Omega _{e}}\left(\nabla \cdot {\boldsymbol {N}}_{a}\right)^{T}\tau _{p}\,{\boldsymbol {\delta }}_{h}\,{\hbox{d}}\Omega }$ (2.115)

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

 ${\displaystyle {\boldsymbol {M}}_{\xi }{\boldsymbol {\Xi }}={\boldsymbol {R}}_{\xi }}$ (2.116) ${\displaystyle {\boldsymbol {M}}_{\delta }{\boldsymbol {\Delta }}={\boldsymbol {R}}_{\delta }}$ (2.117)

where ${\textstyle {\boldsymbol {\Xi }}}$ and ${\textstyle {\boldsymbol {\Delta }}}$ represent the vectors of nodal values of ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle \delta _{h}}$ respectively and the remaining matrices and vectors are given by

 ${\displaystyle {\boldsymbol {M}}_{\xi \,ab}^{e}=\int _{\Omega _{e}}{\boldsymbol {N}}_{a}^{T}\,{\boldsymbol {N}}_{b}\,{\hbox{d}}\Omega }$ (2.118) ${\displaystyle {\boldsymbol {M}}_{\delta \,ab}^{e}=\int _{\Omega _{e}}N_{a}\,N_{b}\,{\hbox{d}}\Omega }$ (2.119) ${\displaystyle {\boldsymbol {R}}_{\xi \,a}^{e}=\int _{\Omega _{e}}{\boldsymbol {N}}_{a}^{T}\cdot {{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}\,{\hbox{d}}\Omega }$ (2.120) ${\displaystyle {\boldsymbol {R}}_{\delta \,ab}^{e}=\int _{\Omega _{e}}N_{a}\,{R^{c}\left(u_{h}\right)}\,{\hbox{d}}\Omega }$ (2.121)

Eqs. 2.116 and 2.117 represent an additional problem, coupled to Eqs 2.111 and 2.112, which effectively doubles the number of nodal unknowns in the problem. However, since the system matrices for the projection problem, defined by Eqs. 2.118 and 2.119, are effectively mass matrices, they can be replaced by the corresponding diagonal mass matrix, which allows us to obtain an approximate projection while avoiding the solution of an additional system.

Note that the matrices ${\textstyle {\boldsymbol {S}}_{m}\left(\tau _{u},{\boldsymbol {a}}\right)}$ and ${\textstyle {\boldsymbol {Q}}_{m}\left(\tau _{u}\right)}$ that appear in Eqs. 2.111 and 2.112 are not strictly necessary in OSS based formulations, since they involve the acceleration term that appears in the residual ${\textstyle {{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}}$. As noted in [39], the acceleration of the large scale lies in the large scale space ${\textstyle V_{h}}$ and therefore its projection should be itself. As a result, we could neglect both these terms if we also not take into account ${\textstyle \partial _{t}u_{h}}$ when evaluating the residual that appears in the projection right hand side vector ${\textstyle {\boldsymbol {R}}_{\xi }}$. However, as the projection is only calculated approximately using a diagonal mass matrix, we have found it convenient to take these terms into account to improve stability.

### 2.5.3 Dynamic ASGS formulation

When we use a dynamic approximation for the subscales we have keep track and update the small scale values at each integration point of the mesh. The dynamic small scale velocities are defined by the local problem given by Eq. 2.47, which was discretized in time using a Backward Euler time scheme to produce Eq. 2.53. Eq. 2.53 provides an expression for ${\textstyle u_{s}^{n+1}}$ in terms of the residual and the old subscale value ${\textstyle u_{s}^{n}}$. The same Backward Euler scheme can also be used to obtain the following time-discrete expression for the small scale acceleration ${\textstyle \partial _{t}u_{s}}$ that appears in the variational form of the problem:

 ${\displaystyle \partial _{t}u_{s}\approx {\frac {u_{s}^{n+1}-u_{s}^{n}}{\delta t}}={\frac {\tau _{t}}{\delta t}}\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}{\big |}_{n+1}-{\boldsymbol {\xi }}_{h}\right)+\rho {\frac {\tau _{t}}{\delta t^{2}}}u_{s}^{n}-{\frac {1}{\delta t}}u_{s}^{n}}$
(2.122)

Introducing Eq. 2.122 in Eqs. 2.56, we obtain a modified momentum equation with time-discrete small scales, given by

 ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{h}\cdot \rho \left(\partial _{t}u_{h}+{\frac {1}{2}}\,{\boldsymbol {a}}\cdot \nabla u_{h}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla {\boldsymbol {w}}_{h}:\rho {\frac {1}{2}}\left({\boldsymbol {a}}\otimes u_{h}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}_{h}\,p_{h}\,{\hbox{d}}\Omega +\int _{\Omega }\nabla ^{s}{\boldsymbol {w}}_{h}:\,2\mu \left(\nabla ^{s}u_{h}-{\frac {1}{3}}\left(\nabla \cdot u_{h}\right){\boldsymbol {I}}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle +\int _{\Sigma \,\Omega ^{e}}\rho {\boldsymbol {w}}_{h}\cdot \left({\frac {\tau _{t}}{\delta t}}\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}\right)+\rho {\frac {\tau _{t}}{\delta t^{2}}}u_{s}^{n}\right)\,{\hbox{d}}\Omega -\int _{\Omega }\nabla \cdot {\boldsymbol {w}}_{h}\tau _{p}\left({R^{c}\left(u_{h}\right)}-\delta _{h}\right)\,{\hbox{d}}\Omega }$ ${\displaystyle -\int _{\Sigma \,\Omega ^{e}}\rho \left({\boldsymbol {a}}\cdot \nabla {\boldsymbol {w}}_{h}\right)\tau _{t}\left({{\boldsymbol {R}}^{m}\left(u_{h},p_{h}\right)}-{\boldsymbol {\xi }}_{h}+{\frac {\rho }{\delta t}}u_{s}^{n}\right)\,{\hbox{d}}\Omega =}$ ${\displaystyle \int _{\Omega }{\boldsymbol {w}}_{h}{\boldsymbol {f}}\,{\hbox{d}}\Omega +\int _{\Gamma _{N}}{\boldsymbol {w}}_{h}\left({\boldsymbol {t}}-\rho {\frac {1}{2}}\left({\boldsymbol {a}}\cdot {\boldsymbol {n}}\right)u_{h}\right)\,{\hbox{d}}\Omega +\int _{\Sigma \,\Omega ^{e}}\rho {\boldsymbol {w}}_{h}{\frac {1}{\delta t}}u_{s}^{n}\,{\hbox{d}}\Omega }$
(2.123)

where the first term in the third row and the last term on the right hand side appear from the time discretization of the small scale acceleration. Note that all terms without an ${\textstyle n}$ index are evaluated at the current time step. Although they play no role in ASGS stabilization, we have included the projections ${\textstyle {\boldsymbol {\xi }}_{h}}$ and ${\textstyle {\boldsymbol {\delta }}_{h}}$ both in Eq. 2.122 and in Eq. 2.123, since the same expressions will be used as a starting point for the dynamic OSS method.

The matrix form of the D-ASGS formulation can be written as

 ${\displaystyle {\big [}{\boldsymbol {M}}+{\boldsymbol {\hat {S}}}_{m}\left(\tau _{t},u_{h}+u_{s}\right){\big ]}{\boldsymbol {\dot {U}}}+{\big [}{\boldsymbol {C}}\left(+\right){\boldsymbol {K}}+{\boldsymbol {\hat {S}}}_{u}\left(\tau _{t},u_{h}+u_{s}\right)+{\boldsymbol {H}}_{u}\left(\tau _{p}\right){\big ]}{\boldsymbol {U}}}$ ${\displaystyle {}+{\big [}{\boldsymbol {G}}+{\boldsymbol {\hat {S}}}_{p}\left(\tau _{t},u_{h}+u_{s}\right){\big ]}{\boldsymbol {P}}={}{\boldsymbol {F}}+{\boldsymbol {T}}+{\boldsymbol {\hat {S}}}_{f}\left(\tau _{t},u_{h}+u_{s}\right)+{\boldsymbol {S}}_{d}\left(\tau _{t},u_{s}^{n}\right)-{\boldsymbol {\hat {S}}}_{d}\left(\tau _{t}\right){u_{s}^{n}}}$
(2.124)

 ${\displaystyle {\boldsymbol {Q}}_{m}\left(\tau _{t}\right){\boldsymbol {\dot {U}}}+{\big [}{\boldsymbol {D}}+{\boldsymbol {Q}}_{u}\left(\tau _{t},u_{h}+u_{s}\right)+{\boldsymbol {Q}}_{d}\left(\tau _{t},u_{s}^{n}\right){\big ]}{\boldsymbol {U}}+{\boldsymbol {Q}}_{p}\left(\tau _{t}\right){\boldsymbol {P}}={\boldsymbol {Q}}_{f}\left(\tau _{t}\right)}$
(2.125)

where we have introduced three new terms on the right hand side involving the old subscale velocity, defined as

 ${\mathbit{S}}_{d}{\left({\tau }_{t},{u}_{s}^{n}\right)}_{a}^{e}={\int }_{{\mathrm{\Omega }}_{e}}{\left(\rho \phantom{\rule{thinmathspace}{0ex}}\mathbit{a}\cdot \mathrm{\nabla }{\mathbit{N}}_{a}\right)}^{T}\cdot \rho \frac{{\tau }_{t}}{\delta t}$