Numerical comparisons of finite element stabilized methods for a 2D vortex dynamics simulation at high Reynolds number

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ Esta es la versión aceptada del artículo publicado en: This is an accepted manuscript of a paper published in: Computer Methods in Applied Mechanics and Engineering (2014): 13/11/2014 DOI: Copyright: El acceso a la versión publicada del artículo puede requerir la suscripción de la revista. Access to the published version may require subscription. “This is an Accepted Manuscript of an article published by Elsevier in Computer Methods in Applied Mechanics and Engineering on 13 November 2014, available at: http://dx.doi.org/10.1016/j.cma.2014.11.023” Accepted Manuscript Numerical comparisons of finite element stabilized methods for a 2D vortex dynamics simulation at high Reynolds number Naveed Ahmed, Samuele Rubino PII: S0045-7825(19)30084-2 DOI: https://doi.org/10.1016/j.cma.2019.02.013 Reference: CMA 12297 To appear in: Comput. Methods Appl. Mech. Engrg. Received date : 9 April 2018 Revised date : 29 January 2019 Accepted date : 10 February 2019 Please cite this article as: N. Ahmed and S. Rubino, Numerical comparisons of finite element stabilized methods for a 2D vortex dynamics simulation at high Reynolds number, Computer Methods in Applied Mechanics and Engineering (2019), https://doi.org/10.1016/j.cma.2019.02.013 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Numerical comparisons of finite element stabilized methods for a 2D vortex dynamics simulation at high Reynolds number Naveed Ahmeda, Samuele Rubinob aDepartment of Mathematics, School of Science and Engineering, Lahore University of Management Sciences, Opposite Sector U, DHA, Lahore Cantt., 54792, Pakistan bDepartment EDAN & IMUS, University of Seville, Avda. Reina Mercedes s/n, 41012 Seville, Spain Abstract In this paper, we consider up-to-date and classical Finite Element (FE) stabilized methods for timedependent incompressible flows. All studied methods belong to the Variational MultiScale (VMS) framework. So, different realizations of stabilized FE-VMS methods are compared using a high Reynolds number vortex dynamics simulation. In particular, a fully Residual-Based (RB)-VMS method is compared with the classical Streamline-Upwind Petrov–Galerkin (SUPG) method together with grad-div stabilization, a standard one-level Local Projection Stabilization (LPS) method, and a recently proposed LPS method by interpolation. These procedures do not make use of the statistical theory of equilibrium turbulence, and no ad-hoc eddy viscosity modeling is required for all methods. Applications to the simulation of a high Reynolds number flow with vortical structures on relatively coarse grids are showcased, by focusing on a two-dimensional plane mixing-layer flow. Both Inf-Sup Stable (ISS) and Equal Order (EO) H1-conforming FE pairs are explored using a second-order semi-implicit Backward Differentiation Formula (BDF2) in time. Based on the numerical studies conducted, it is concluded that the SUPG method using EO FE pairs performs best among all methods. Furthermore, there seems to be no reason to extend the SUPG method by the higher order terms of the RB-VMS method. Keywords: Variational multiscale methods; finite element stabilized methods; high Reynolds number incompressible flow; 2D vortex dynamics problem 1. Introduction In this paper, we consider up-to-date and classical Finite Element (FE) stabilized methods for timedependent incompressible flows fulfilling the incompressible Navier–Stokes Equations (NSE). Let Ω ∈Rd, d∈ {2,3}, be a bounded domain with Lipschitz boundary Γ and (0, T) be a bounded time interval. The incompressible NSE read as follows: Find a velocity field u: (0, T ]×Ω→Rdand a pressure field p: (0, T]×Ω→Rsuch that ∂tu−ν∆u+ (u·∇)u+∇p=fin (0, T]×Ω, ∇·u= 0 in [0, T]×Ω, u=0on [0, T]×Γ, u(0,x) = u0in Ω, (1) where νis the kinematic viscosity that is assumed to be positive and constant, fis the given body force, and u0is the given initial velocity field, assumed to be divergence-free. For simplicity of presentation, the case of homogeneous Dirichlet conditions is considered on the whole boundary. ∗Corresponding author Email addresses: [email protected] (Naveed Ahmed), [email protected] (Samuele Rubino) Preprint submitted to Elsevier January 29, 2019 *Manuscript Click here to download Manuscript: main.tex Click here to view linked References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 The main contribution of this paper is a comprehensive and thorough numerical study in the FE stabilized framework of two-scales fully Residual-Based (RB) and local projection-based Variational MultiScale (VMS) methods for a time-dependent high Reynolds number incompressible flow with a strong dynamic vortical structure. The derivation of efficient and accurate numerical schemes for the simulation of turbulent incompressible flows is a very active field of research. In particular, various realizations of VMS methods for simulating turbulent incompressible flows have been proposed in the past fifteen years (see [1] for a recent detailed review). All of these realizations obey the basic principles of VMS methods: they are based on the variational formulation of the incompressible NSE and the scale separation is defined by projections. However, apart from these common basic features, the various VMS methods look quite different. In this paper, our main goal is to focus on two-scales VMS methods, and provide a thorough numerical investigation of up- to-date and classical FE stabilized methods belonging to this category when applied to a relevant fixed setup for numerical studies such as the 2D Kelvin–Helmholtz instability problem. Indeed, even if VMS methods, despite their relatively recent development, are already well-established, and considered state-of-the-art in turbulence modeling that provides a promising and successful alternative to classical Large Eddy Simulation (LES) models, in the literature, not much has been done about a structured comparison of them in terms of numerical studies. To the best of our knowledge, the first (and only) attempt to go towards this research direction has been performed in [2], where the authors studied different realizations of VMS methods within the framework of FE in turbulent channel flow simulations. However, their main focus was on three-scales VMS models, in which the effect of the unresolved scales on the resolved ones is modeled by means of an eddy viscosity term of Smagorinsky type that only acts directly on the small resolved scales. In the present paper, it is aimed to complement and extend this research avenue, by mainly focusing on two-scales VMS methods, which use a direct modeling of the subgrid-scale flow by numerically approximating the related equations. This numerical approach, which hence relies on purely numerical artifacts without any modification of the continuous problem, was seldom followed. The MILES (Monotone Integrated LES) approach [3] being the main exception, until the RB-VMS models were introduced in the seminal papers [4, 5] and subsequently proposed as implicit LES techniques (ILES) for turbulent flows in [6]. These models provide a unified framework for the definition of spatial approximation schemes capable of preventing numerical instabilities that arise when the standard Galerkin FE method is used, and are adequate to represent the turbulence LES modeling. Thus, these models do not need any modeling of the subgrid-scales by statistical theories of turbulence, and in particular they do not include eddy viscosity. The numerical diffusion inherent to those stabilized models basically plays the role of the eddy diffusion. In this way, the present paper aims at giving a thorough numerical investigation, similar to the one performed in [2], but for two-scales VMS methods. A structured presentation is provided in this framework, with special emphasis on experience in numerical studies. After reaching almost “definitive” conclusions within this paper, a comparison of the selected “best performing” two-scales VMS method with three-scales VMS methods that use eddy viscosity (in a more or less sophisticated manner) to model the effect of subgrid-scales shall appear in a forthcoming paper. In this way, the numerical performances of different VMS methods would be assessed. The RB-VMS method was introduced in [6]. A straightforward simplification of the RB-VMS method leads to the classical Streamline-Upwind Petrov–Galerkin (SUPG) method [7, 8]. Also, another variant of the RB-VMS method, which is not fully consistent, but of optimal order with respect to the FE interpolation, is given by the so-called Local Projection Stabilization (LPS) methods [9, 29, 47]. So, different realizations of stabilized FE-VMS methods are compared in a 2D vortex dynamics simulation at high Reynolds number in this paper. In particular, the RB-VMS method [6] is compared with the classical SUPG method [7, 8] together with grad-div stabilization, a standard one-level LPS method [10], and a recently proposed LPS method by interpolation [11, 12]. To the best of our knowledge, a comparison of these methods is so far not available. However, many other alternatives exists, such as for instance the stabilized method of Franca and Frey [48] employing the negative adjoint differential operator (USFEM, Unusual Stabilized FEM) just to cite one, but have not considered here in order keep the length of the paper within reasonable limits. To keep the paper self-explanatory, a brief presentation of the cited numerical methods will be considered hereafter to provide the basic concepts. For more details on their derivation, see the up-to-date review on VMS methods for the simulation of turbulent incompressible flows [1]. To assess the different numerical methods, applications to the simulation of a high Reynolds number 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 flow with vortical structures on relatively coarse grids are presented, by focusing on a two-dimensional plane mixing-layer flow as benchmark problem, since it presents a wide range of flow scales and an interesting time evolution of the flow field. Starting from a perturbed initial condition, the transition to the development of small vortices takes place, which are then paired to larger vortices until one single eddy finally remains, rotating at a fixed position. In particular, we analyze different quantities of interest associated to this problem (i.e, temporal evolution of kinetic energy, enstrophy, palinstrophy, and vorticity thickness) in order to judge the performance of all the studied methods and draw some definitive conclusions. All the numerical results are benchmarked against a reference solution, obtained in the very recent paper [13]. There, the authors present computational studies for the long-time integration of the 2D Kelvin–Helmholtz instability problem at three different Reynolds numbers (Re = 102,103,104) with high order divergence-free H(div)- FE methods. In particular, they used a Hybrid Discontinuous Galerkin (HDG) approach for the spatial discretization and a multistep implicit-explicit (IMEX) method that combines a second-order Backward Differentiation Formula (BDF2) with a second-order accurate extrapolation in time. In our numerical simulations, we have used the results obtained with Re = 104on the finest mesh consists of 2562square as a reference solution. Note that this model problem is very sensitive to small perturbations that are almost unavoidable in numerical simulations, thus some targets, such as a mesh-independent prediction of the final pairing, are not achieved even among the very accurate simulations with the highest resolution in [13]. In this paper, we mainly focus on efficient spatial and temporal discretizations, for which both Inf-Sup Stable (ISS) and Equal Order (EO) low-order FE pairs are explored, using a second-order semi-implicit Backward Differentiation Formula (BDF2) in time, where the linearization of the fully discrete problem at each time step is done by means of temporal extrapolation. In contrast to a fully implicit scheme, this approach yields a unique linear system of equations to be solved at each time step. Altogether, performing simulations with semi-implicit schemes uses less computing time than fully implicit schemes. However, while a fully implicit approach is generally yielding a stable time discretization scheme, a semi-implicit approach may require a time step restriction due to the stability issues of the time stepping scheme. Note that semi-implicit BDF schemes for the numerical simulation of NSE with VMS turbulence modeling have already been investigated in the literature, see for instance [15], and also [16] for a stable velocity-pressure segregation version. The paper is organized as follows. Section 2 explains briefly the used VMS methods with a special focus on the derivation of two-scales VMS methods. In Section 3, a semi-implicit approach based on two-step BDF (BDF2) for the time discretization is detailed for each studied method, together with some numerical implementation aspects. Section 4 presents the numerically comparison for the simulation of two-dimensional Kelvin–Helmholtz instabilities in the high Reynolds number regime. Here, several quantities of interest are presented, evaluated and discussed. Finally, Section 5 summarizes the main conclusions of the paper and gives an outlook. 2. Variational multiscale methods As already mentioned, VMS methods are based on the variational formulation of the incompressible NSE (1). To define the variational formulation of (1), the velocity space V= [H1 0(Ω)]dand the pressure space Q=L2 0(Ω) are introduced. Let (·,·) denote the L2inner product with respect to the domain Ω. The variational formulation of (1) reads as follows: Find (u, p) : (0, T)→V×Qsuch that for all (v, q)∈V×Q d dt(u,v) + ν(∇u,∇v) + ((u·∇)u,v)−(p, ∇·v)+(∇·u, q) = hf,viin D?(0, T ),(2) with u(0,x) = u0(x) in Ω, where h·,·i denotes the duality pairing between the velocity space Vand its dual V?and D?(0, T) is the space of distribution on (0, T ). In standard conforming Finite Element (FE) formulations, the infinite-dimensional spaces (V, Q) are replaced with finite dimensional-subspaces (Vh, Qh) consisting of typically low-order piecewise polynomials with respect to a triangulation Thof Ω. In this paper, both Inf-Sup Stable (ISS, [17, 18]) and Equal Order (EO) H1-conforming FE pairs are explored, which are not exactly divergence-free, by considering in general 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 the popular Taylor–Hood FE pair Pk/Pk−1[19] and the EO FE pair Pk/Pk, respectively, with k≥2, where Pkdenotes the space of continuous functions whose restriction to each mesh cell K∈ This the Lagrange polynomial of degree less than or equal to kand Pk= [Pk]d.Thus, low-order H1-conforming FE (not exactly divergence-free) are considered in this work. Another alternative strategy is to use discontinuous approximations such as Discontinuous Galerkin (DG) methods. This has been explored for instance in [13, 36], where, in particular, an exactly divergence-free hybrid discontinuous Galerkin (HDG) method based on H(div)-FE is considered. However, this usually leads to a more expensive discretization. 2.1. Two-scales VMS methods This section discusses the basic concepts of two-scales VMS methods. A starting point of two-scales VMS methods is the separation of the flow field into resolved scales (u, p) and unresolved scales (u0, p0) such that u=u+u0and p=p+p0. Analogously, a direct-sum decomposition of velocity space V=V⊕V0 and pressure space Q=Q⊕Q0is considered. It should be emphasized that although this approach is in principle the same as in Large Eddy Simulations (LES), it is well known that the definition of the scales is different. A variational projection, either L2projection or elliptic projection, for the separation of scales and spaces is performed in VMS methods. Note that the VMS methodology allows further decompositions of the resolved scales. The most common approach of this kind is a decomposition of these scales into large resolved scales (or large scales) and small resolved scales, leading finally to a so-called three-scales VMS method. In this case, the effect of the unresolved scales on the resolved ones is modeled by means of an eddy viscosity term that only acts directly on the small resolved scales (cf. [20–22]). However, in the present paper, we just focus on the comparison between VMS methods that use a direct modeling of the subgrid-scale flow by approximating the related equations, for which no eddy viscosity is introduced to model the effect of the subgrid-scales. This is the reason why we restrict to two-scales VMS methods. For clarity of presentation, the weak formulation (2) of the NSE is expressed in a short form as follows: Given u(0,x) = u0(x), find (u, p) : (0, T)→V×Qsatisfying A(u; (u, p),(v, q)) = f(v)∀(v, q)∈V×Q. (3) Decomposing the test functions also into two scales and using the linearity with respect to the test functions, the variational formulation (3) leads to the coupled set of equations: •an equation for the resolved scales A(u; (u, p),(v, q)) + A(u; (u0, p0),(v, q)) = f(v),(4) •and an equation for the unresolved scales A(u; (u, p),(v0, q0)) + A(u; (u0, p0),(v0, q0)) = f(v0).(5) The form A(·;·,·) is decomposed into its linear part and the trilinear convective term as A(u;U,W) = Alin (U,W) + ((u·∇u),v) where the abbreviations U= (u, p)Tand W= (v, q)Tare used for simplicity. Then, the equation (5) for the unresolved scales can be written in the form AU(U0,W0) + ((u0·∇)u0,v0) = RU,W0(6) with AU(U0,W0) = Alin (U0,V0) + ((u0·∇)u,v0) + ((u·∇)u0,v0), RU,V0=f(v0)−Alin U,W0−((u·∇u),v0), 4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 where AU(U0,W0) is the Gˆateaux derivative of A(·;·,·) at Uin the direction of U0. The solution of (6) can be formally represented as U0=FURU,(7) which can be interpreted as the unresolved scales that are derived as a function of the residual of the resolved scales. Finally, inserting expression (7) in the resolved scales equations (4) leads to a single set of equations for the resolved scales. Two-scales VMS methods aim to approximate FUby models which do not rely on considerations from the physics of turbulent flow, but are derived just with mathematical arguments. In the next subsections, concrete approaches will be presented. 2.2. Residual-based VMS method The main idea in the derivation of the two-scales RB-VMS method is based on a perturbation series with respect to the norm of the residual associated with the resolved scales. It is proposed in [6] to truncate the series after the first term and to apply some modeling of this term. The resulting method can be considered as a generalization of classical stabilization methods for the NSE. A perturbation series for a potentially small quantity ε=kR(U)k(V0×Q0)∗is considered. It is assumed that the larger the space (V×Q), the better Uapproximates U, and the smaller is R(U). The perturbation series is of the form U0=εU0 1+ε2U0 2+. . . = ∞ X i=1 εiU0 i.(8) In particular, if ε= 0, i.e. R(U) = 0, then U0=FURU= 0 from (7)-(8) . Inserting the perturbation series (8) in the terms of (5) for the unresolved scales gives AU ∞ X i=1 εiU0 i,W0!= ∞ X i=1 εiAU(U0 i,W0) and ∞ X i=1 εiu0 i·∇!∞ X i=1 εiu0 i,v0!=ε2((u0 1·∇)u0 1,v0) + ε3[((u0 1·∇)u0 2,v0) + ((u0 2·∇)u0 1,v0)] + . . . = ∞ X i=2 εi  i−1 X j=1 (u0 j·∇)u0 i−j,v0 . Substituting these terms into (5) yields ∞ X i=1 εiAU(U0 i,W0) + ∞ X i=2 εi  i−1 X j=1 (u0 j·∇)u0 i−j,v0 =ε*RU RU(V0×Q0)∗ ,W0+. Collecting similar terms with respect to εleads to a system of variational problems which are coupled through the right-hand side, that is AU(U0 1,W0) = *RU RU(V0×Q0)∗ ,W0+, AU(U0 i,W0) = − i−1 X j=1 (u0 j·∇)u0 i−j,v0i≥2. 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 In the modeling of the unresolved scales, it is suggested in [6] to truncate the series (8) after the first term, and to use a linear approximation of the so-called fine-scale Green’s operator that formally represents U0 1 U0≈εU0 1=kR(U)k(V0×Q0)∗U0 1≈τRU=τRuh ph =  τm(fh−∂tuh+ν∆uh−(uh·∇)uh−∇ph) −τc(∇·uh) =  RM h RC h (9) where τis a 4×4 diagonal tensor-valued function, and the approximation of the resolved scales is computed in a standard FE space. The RB-VMS FE formulation is obtained by inserting the approximation (9) into the large scales equation (4), omitting the models of the terms (∂tu0,vh) and ν(∇u0,∇vh), neglecting the inter-element jumps of the fine scale functions, and integrating by parts the continuity equation with respect to the unresolved scale in (4), assuming that u0= 0 on Γ: Find uh: (0, T)→Vh, ph: (0, T )→Qhsatisfying (∂tuh,vh) + ν(∇uh,∇vh) + ((uh·∇)uh,vh)−(ph,∇·vh)+(∇·uh, qh) +bRM h,uh,vh+buh,RM h,vh+bRM h,RM h,vh −RC h,∇·vh−RM h,∇qh= (fh,vh) (10) for all (vh, qh)∈Vh×Qh, where bin (10) denotes the trilinear convective form given by b(u,v,w) = ((u·∇)v,w),u,v,w∈V. Concerning the actual choice of b, it is advisable from the practical point of view that one does not need to compute a derivative of the residual of the momentum equation. For this reason, it is suggested to use the following form, which is obtained from the divergence form with integration by parts: b(u,v,w)=(∇·(uvT),w) = −(uvT,∇w).(11) The two terms b(RM h,uh,vh) and b(uh,RM h,vh) are known as cross-stress terms, and bRM h,RM h,vhas the subgrid (or Reynolds-stress) term. Using (11), (uvT,∇w)=(v,(∇w)Tu) and (∇v)u= (u·∇)v, one gets for the first cross-stress term in (10): bRM h,uh,vh=−RM h(uh)T,∇vh=−uh,(∇vh)TRM h=−RM h,(∇vh)uh=−RM h,(uh·∇)vh, (12) which together with the last term in the left-hand side of (10) gives: bRM h,uh,vh−(RM h,∇qh) = −RM h,(uh·∇)vh+∇qh.(13) This term corresponds to the well known stabilization term of the Streamline-Upwind Petrov-Galerkin (SUPG) method for the convection field uh. One can also observe the contribution of the so-called grad-div stabilization term by inserting the concrete formula of the residual of the continuity equation into (10), that is: (τc∇·uh,∇·vh).(14) Similarly, using (11) and (uvT,∇w) = (v,(∇w)Tu), one obtains for the second cross-stress term and the subgrid term in (10): buh,RM h,vh=−uh(RM h)T,vh=−RM h,(∇vh)Tuh,(15) b(RM h,RM h,vh) = −RM h(RM h)T,vh=−RM h,(∇vh)TRM h.(16) Considering formulas (12) and (15) for the cross-stress terms, and formula (16) for the subgrid term, the RB-VMS method (10) can be expressed as: 6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Find uh: (0, T)→Vh, ph: (0, T )→Qhsatisfying (∂tuh,vh) + ν(∇uh,∇vh) + ((uh·∇)uh,vh)−(ph,∇vh)+(∇·uh, qh)−RM h,(uh·∇)vh+C∇qh −RM h,(∇vh)Tuh−RM h,(∇vh)TRM h+ (τc∇·uh,∇·vh)=(fh,vh),(17) for all (vh, qh)∈Vh×Qh. The formulation (17) provides the complete RB-VMS method, which retains numerical consistency in the FE equations, in the sense that the continuous solution exactly satisfies the discrete equations, whenever it is smooth enough. In this paper, both ISS and EO H1-conforming FE pairs would be explored. For this reason, we have added the constant Cin the formulation (17), so that C= 1 when using EO FE pairs, and we will drop the dependency of the pressure stabilization term from (17) when using ISS FE pairs by fixing C= 0. We recall that in (17) the terms RM h,(uh·∇)vh+∇qhand τc(∇·uh,∇·vh) are the classical stabilization terms of the SUPG and grad-div methods, respectively. In this paper, we are interested in performing numerical studies also with a simplified model arising from (17), which is the classical SUPG method together with grad-div stabilization: Find uh: (0, T)→Vh, ph: (0, T )→Qhsatisfying (∂tuh,vh) + ν(∇uh,∇vh) + ((uh·∇)uh,vh)−(ph,∇vh)+(∇·uh, qh)−RM h,(uh·∇)vh+C∇qh + (τc∇·uh,∇·vh) = (fh,vh),(18) for all (vh, qh)∈Vh×Qh, again for both ISS (C= 0) and EO (C= 1) FE pairs. 2.3. Local projection stabilization methods Local Projection Stabilization (LPS) methods are stabilization methods that provide specific stabilization of any single operator term that could be a source of instability for the numerical discretization. They were introduced in [9] and they could be viewed as simplifications of the two-scales RB-VMS method described in the previous section. Indeed, LPS methods are not fully consistent (only specific dissipative interactions are retained), but are of optimal order with respect to the FE interpolation. The fact that the stabilization enjoys the right asymptotic behavior without full consistency allows to decouple the stabilization of the pressure and the velocity without having all the residual terms coupled thus relying on a term-by-term structure. This feature could be considered an important advantage with respect to the more complex RBVMS method in view of practical implementations such as to perform the numerical analysis since it leads to a simpler and less expensive structure. Different variants of LPS methods have been investigated during the recent years for incompressible flow problems. The main common feature is that, thanks to local projection, the symmetric stabilization terms only act on the small scales of the flow, thus ensuring a higher accuracy with respect to more classical stabilization procedures, such as penalty-stabilized methods, cf. [23]. Thus, the effect of LPS is on the one hand to improve the convergence to smooth solutions and on the other hand, for rough solutions, it limit the propagation of perturbations generated in the vicinity of sharp gradients. This way these schemes remain stable and useful tools for the simulation of turbulent flows. As a single rule, the structure of LPS method is achieved by considering in the RB-VMS method (17) just the specific dissipative interactions that stabilize convection and pressure gradient, and by introducing local L2projections in the approximation of the unresolved scales, in such a way the symmetric stabilization terms only act on the small scales of the flow. This leads to a family of methods associated with the choice of the actual local L2projection. The main derivation of the LPS methods will be introduced here for the NSE (1). The stabilization effect is achieved by adding least-squares terms that give a weighted control on the fluctuations of the quantity of interest. This control is based upon a projection operation πh:L2(Ω) 7→ Dhonto a discontinuous FE space Dh(the “projection” space). This space is built on a grid Mhformed by macro-elements built from the triangulation Thof Ω. The component-wise extension of πhto vector functions is denoted by πh. The LPS approximation of the NSE reads: 7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0 50 100 150 200 250 300 350 400 0.479 0.48 0.481 0.482 EKin EO RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0.479 0.48 0.481 0.482 EKin EO SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0.477 0.478 0.479 0.48 0.481 0.482 EKin EO one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0.479 0.48 0.481 0.482 time in units EKin EO LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 1: Temporal evolution of kinetic energy with EO–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. 0 50 100 150 200 250 300 350 400 0.479 0.48 0.481 0.482 EKin ISS RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0.479 0.48 0.481 0.482 EKin ISS SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0.478 0.479 0.48 0.481 0.482 EKin ISS one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0.479 0.48 0.481 0.482 time in units EKin ISS LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 2: Temporal evolution of kinetic energy with ISS–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. monotonically decreasing kinetic energy is obtained for all the methods, which can be clearly seen in the Figures 1 and 2. For RB-VMS and SUPG methods, on all mesh levels the kinetic energy decreases very slowly, around 0.3%, as in [36]. However, results on the finest grid approaches the reference solution for the LPS methods too. 4.4. Enstrophy The temporal evolution of the enstrophy is plotted in Figures 3 and 4 for all methods respectively with EO and ISS FE on different refined meshes and the time step length ∆t= 3.125×10−3. Similar to the kinetic 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 EO RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 EO SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 EO one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 EO LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 3: Temporal evolution of enstrophy with EO–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 ISS RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 ISS SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 ISS one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 ISS LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 4: Temporal evolution of enstrophy with ISS–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. energy, the amount of the initial enstrophy is almost the same for all simulations and then behaves slightly different for different refinements. The strictly decaying behavior of the enstrophy is correctly displayed in Figures 3 and 4 for all methods. Nevertheless, only results for the EO SUPG method are almost in agreement on all grid levels with the finest solution in [36] and reference solution approximatively up to t= 200¯ t, with time unit ¯ t=δ0/U∞. However, EO RB-VMS and LPS by interpolation methods give comparable results on the finest level. 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 EO RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 EO SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 103 104 105 106 time in units 1 2k∇ωhk2 0 EO one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 EO LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 5: Temporal evolution of palinstrophy with EO–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. 4.5. Palinstrophy The palinstrophy is one of the most sensitive quantity of interest with respect to the different pairings of vortices. In Figures 5 and 6, the temporal evolution of the palinstrophy is presented for all methods using respectively EO and ISS FE on different refined meshes, and the time step length ∆t= 3.125 ×10−3. In contrast to kinetic energy and enstrophy, palinstrophy can increase. In fact local maxima are almost attained once merging processes of the vortices terminate. Reference solution indicates that the last merging process does not terminate before t= 200¯ t.As noticed in [13], we also observe here that although the magnitude of the palinstrophy is strongly mesh-dependent, the points in time where the first two pairings occur can be approximately identified independently of the particular level resolution. In terms of time intervals for local maxima, the best results for all mesh levels are obtained again by EO SUPG method, where one can see that the last merging process is still not ended at t= 200¯ tfor all grid levels. In terms of magnitude, EO LPS by interpolation method on the finest grid better approaches the reference solution. However, we reiterate that this quantity, both in magnitude and time intervals for local maxima, is highly dependent on the studied method, used mesh, refinement level, FE pair, and time step length. 4.6. Vorticity Thickness The temporal evolution of the relative vorticity thickness δ/δ0for all methods on different refinement levels is displayed in Figure 7 using EO FE and in Figure 8 using ISS FE, respectively. The computational results are obtained with the time step length ∆t= 3.125 ×10−3. The formation of succeeding peaks in the evolution of the relative vorticity thickness corresponds to the pairing process of the eddies. For the reference solution, the local maximum δ/δ0= 6.04 at t= 34¯ tindicates the pairing of two eddies from four. Comparing the relative vorticity thickness computed with all stabilization schemes on different refinement levels clearly indicates that the first pairing occurs almost at the same time. After that, the relative vorticity oscillates until the next pairing of eddies occur. The pairing into the final eddy happens somehow at a different time for different stabilization methods compared to the reference solution, where it occurs at t= 220¯ t. In comparison with the reference solution, the second pairing is closest on all grid levels for the EO SUPG method, where it occurs around t= 180¯ tfor all mesh levels. Also, EO RB-VMS and LPS by interpolation methods give comparable results on the finest level. Thus, note that the development of the extremely sensitive quantity of vorticity thickness strongly depends on the mesh refinement. In particular, the last 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 ISS RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 ISS SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 ISS one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 ISS LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 6: Temporal evolution of palinstrophy with ISS–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. pairing process, where two secondary eddies merge to become one, is very sensible with respect to how accurate the simulation is. However, the actual values of the amplitudes of the various peaks are almost identical for all refinement levels, and in agreement both with reference solution [13] and results in [36]. Similar conclusions can be drawn both for ISS and EO FE, see Figures 7 and 8. One can see that the mesh refinement have again a noticeable influence on the temporal development of the vorticity thickness, but the values of the amplitude are almost identical. Similar to the results for the palinstrophy, this quantity too is thus very sensitive with respect to vortex pairings. Altogether, the EO SUPG method is superior to all other methods, since almost approaches the very fine reference solution on relatively coarse grids. However, the EO RB-VMS method and EO/ISS LPS by interpolation methods on the finest level (Level 8) also perform quite well, being almost in agreement with the finest simulation in [36]. 4.7. Comparison of RB-VMS and SUPG methods on finer time step lengths Based on the computational studies presented in the previous sections and after having performed numerous simulations with different time step lengths, it is noticed that smaller time steps could lead to very accurate results on relatively coarse grids. Therefore, this section is devoted to the numerical comparison of the two methods (RB-VMS and SUPG) that perform better on the larger time step length of the previous sections on a rather coarse mesh (Level 6) with the finer time step lengths (∆t1= 7.8125 ×10−4and ∆t2= 5.9523 ×10−3). For the simplicity of presentation, we will use ∆t1and ∆t2as an abbreviation for the finer time step lengths. Results for EO FE are given in Figure 9 and Figure 10, respectively, with ∆t1and ∆t2. Concerning the results for ∆t1on Figure 9, observing the peaks and times at which the pairing of eddies occurs, one can clearly see that the second pairing in the case of RB-VMS occurs a bit earlier than the SUPG method, see Figure 9 (top left). However, concerning the kinetic energy, one can not observe any difference between the two methods, see Figure 9 (top right). On the other hand, for enstrophy and palinstrophy, both methods perform quite similar up to t= 180¯ ttime units. One can see that the SUPG method is more accurate with respect to the RB-VMS method when compared to the reference solution. Concerning the results for the finer time step length ∆t2on Figure 10, one can clearly observe that the results computed with the SUPG method are almost comparable to the reference data even on the rather coarse mesh for all quantities of interest, especially the relative vorticity thickness. The palinstrophy being the most sensitive to the 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0 50 100 150 200 250 300 350 400 2 4 6 8 10 δ/δ0 EO RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 2 4 6 8 10 δ/δ0 EO SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 2 4 6 8 10 δ/δ0 EO one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0 2 4 6 8 10 time in units δ/δ0 EO LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 7: Temporal evolution of vorticity thickness with EO–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. 0 50 100 150 200 250 300 350 400 2 4 6 8 10 δ/δ0 ISS RBVMS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 2 4 6 8 10 δ/δ0 ISS SUPG Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 2 4 6 8 10 δ/δ0 ISS one-level LPS Level 6 Level 7 Level 8 reference 0 50 100 150 200 250 300 350 400 0 2 4 6 8 10 time in units δ/δ0 ISS LPS by Interpolation Level 6 Level 7 Level 8 reference Figure 8: Temporal evolution of vorticity thickness with ISS–FE: RB-VMS (top left), SUPG (top right), one-level LPS (bottom left), and LPS by interpolation (bottom right), on different mesh refinement levels, ∆t= 3.125 ×10−3. numerical setup (cf. [13]) is an exception. On the other hand, compared to the time step length ∆t1, there is a slight improvement in the RB-VMS results computed with ∆t2, compare the red curves in Figures 9 and 10. Comparing these results with the corresponding results in the previous sections, it can be clearly seen that finer time step length already allowed to almost reach the reference results on a relatively coarse structured triangular grid, at least for the SUPG method with EO FE. Note that a (more than ten times) finer temporal resolution is used for the reference solution, i.e. ∆t= 3.6×10−5in [13]. In the case of ISS FE, the results are plotted in Figures 11 and 12. For both time step lengths ∆t1and ∆t2, it could be observed in all simulations that both stabilization schemes gave the same results. 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0 50 100 150 200 250 300 350 400 2 4 6 8 10 time in units δ/δ0 EO SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 0.479 0.48 0.481 0.482 time in units EKin EO SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 EO SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 EO SUPG vs RBVMS SUPG RBVMS reference Figure 9: Temporal evolution with EO–FE of: vorticity thickness (top left), kinetic energy (top right), enstrophy (bottom left), and palinstrophy (bottom right), on mesh refinement level 6, ∆t= 7.8125 ×10−4. 0 50 100 150 200 250 300 350 400 2 4 6 8 10 time in units δ/δ0 EO SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 0.48 0.481 0.482 time in units EKin EO SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 EO SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 EO SUPG vs RBVMS SUPG RBVMS reference Figure 10: Temporal evolution with EO–FE of: vorticity thickness (top left), kinetic energy (top right), enstrophy (bottom left), and palinstrophy (bottom right), on mesh refinement level 6, ∆t= 5.9523 ×10−4. 5. Summary and outlook In this paper, we compared two-scales VMS stabilized FE methods for the simulation of the incompressible NSE. These methods are widely used as one of the most promising and successful approaches that seek to simulate large-scale structures in turbulent flows. The space discretization for the studied methods using both ISS and EO FE is combined with a second-order semi-implicit time stepping scheme, based on BDF. Relatively coarse grids are chosen for the space discretization, starting from large to small time step lengths. Several variants of two-scales VMS approaches, from fully residual-based to weakly consistent, have been 19 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0 50 100 150 200 250 300 350 400 2 4 6 8 10 time in units δ/δ0 ISS SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 0.48 0.481 0.482 time in units EKin ISS SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 ISS SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 ISS SUPG vs RBVMS SUPG RBVMS reference Figure 11: Temporal evolution with ISS–FE of: vorticity thickness (top left), kinetic energy (top right), enstrophy (bottom left), and palinstrophy (bottom right), on mesh refinement level 6, ∆t= 7.8125 ×10−4. 0 50 100 150 200 250 300 350 400 2 4 6 8 10 time in units δ/δ0 ISS SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 0.48 0.481 0.482 time in units EKin ISS SUPG vs RBVMS SUPG RBVMS reference 0 50 100 150 200 250 300 350 400 20 25 30 35 time in units 1 2kωhk2 0 ISS SUPG vs RBVMS SUPG RBVMS reference 050 100 150 200 250 300 350 400 104 105 106 time in units 1 2k∇ωhk2 0 ISS SUPG vs RBVMS SUPG RBVMS reference Figure 12: Temporal evolution with ISS–FE of: vorticity thickness (top left), kinetic energy (top right), enstrophy (bottom left), and palinstrophy (bottom right), on mesh refinement level 6, ∆t= 5.9523 ×10−4. applied to the simulation of 2D Kelvin–Helmholtz instabilities, triggered by a plane mixing layer at high Reynolds number Re = 104. Section 4 presents, the detailed comparison of RB-VMS, SUPG, one-level variant of LPS and LPS by interpolation methods using both EO and ISS pair of FE on rather coarse grid levels and with different time steps, with the aim of studying their influence on the accuracy of the numerical solutions. The numerical performances of all studied methods is discussed by monitoring the relevant quantities of interest, such as relative vorticity thickness, kinetic energy, enstrophy, and palinstrophy. From the computational point of view, note that this problem is very sensitive and results strongly depend on the methods used, grid 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 refinement, and time step lengths. Through our numerical experiences, we have shown the need to consider relatively small time steps, both to prevent numerical stability issues of a less expensive semi-implicit time stepping scheme, and to guarantee not excessive numerical dissipation. Altogether, based on the presented numerical studies, it turns out that the EO SUPG method with the small time step length outperforms all other studied variants. Closest results to this best performing method are attained by RB-VMS method, for which however the extra terms seem not to provide increased accuracy for the studied problem on relatively coarse grids, and thus there seems to be no reason to extend the simpler SUPG method by the higher order terms of the more complex RB-VMS method in this case. On the other side, LPS methods, which are not fully consistent but are of optimal order with respect to the FE interpolation, despite their appealing structure both in terms of practical implementations such as to perform the numerical analysis, seems to need higher space resolutions in order to achieve the same accuracy of fully residual-based VMS stabilized methods. As a future research direction, we plan to compare the selected best performing two-scales VMS stabilized methods towards several variants of three-scales VMS methods that use turbulent eddy viscosity (in a more or less sophisticated manner) to model the effect of subgrid-scales, also on more complex problems presenting genuine 3D turbulent structure, like 3D turbulent channel flow. References [1] N. Ahmed, T. Chac´on Rebollo, V. John, S. Rubino, A review of variational multiscale methods for the simulation of turbulent incompressible flows, Archives of Computational Methods in Engineering 24 (1) (2017) 115–164. [2] V. John, A. Kindl, Numerical studies of finite element variational multiscale methods for turbulent flow simulations, Comput. Methods Appl. Mech. Engrg. 199 (13-16) (2010) 841–852. [3] J. P. Boris, F. F. Grinstein, E. S. Oran, R. L. Kolbe, New insights into large eddy simulation, Fluid Dyn. Res. 10 (4-6) (1992) 199–228. [4] T. J. R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid-scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1-4) (1995) 387–401. [5] T. J. R. Hughes, G. R. Feij´oo, L. Mazzei, J. B. Quincy, The variational multiscale method—a paradigm for computational mechanics, Comput. Methods Appl. Mech. Engrg. 166 (1-2) (1998) 3–24. [6] Y. Bazilevs, V. M. Calo, J. A. Cottrell, T. J. R. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 197 (1-4) (2007) 173–201. [7] A. N. Brooks, T. J. R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1-3) (1982) 199–259, fENOMECH ’81, Part I (Stuttgart, 1981). [8] T. J. R. Hughes, A. Brooks, A multidimensional upwind scheme with no crosswind diffusion, in: Finite element methods for convection dominated flows (Papers, Winter Ann. Meeting Amer. Soc. Mech. Engrs., New York, 1979), Vol. 34 of AMD, Amer. Soc. Mech. Engrs. (ASME), New York, 1979, pp. 19–35. [9] R. Becker, M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo 38 (4) (2001) 173–199. [10] G. Matthies, P. Skrzypacz, L. Tobiska, A unified convergence analysis for local projection stabilisations applied to the Oseen problem, M2AN Math. Model. Numer. Anal. 41 (4) (2007) 713–742. [11] T. Chac´on Rebollo, M. G´omez M´armol, V. Girault, I. S´anchez Mu˜noz, A high order term-by-term stabilization solver for incompressible flow problems, IMA J. Numer. Anal. 33 (3) (2013) 974–1007. [12] N. Ahmed, T. Chac´on Rebollo, V. John, S. Rubino, Analysis of a full space-time discretization of the Navier–Stokes equations by a local projection stabilization method, IMA J. Numer. Anal. 37 (3) (2017) 1437–1467. [13] P. W. Schroeder, V. John, P. L. Lederer, C. Lehrenfeld, G. Lube, J. Sch¨oberl, On reference solutions and the sensitivity of the 2D Kelvin–Helmholtz instability problem, Comput. Math. Appl. (2018) https://doi.org/10.1016/j.camwa.2018.10.030. [14] T. J. R. Hughes, L. P. Franca, M. Balestra, A new finite element formulation for computational fluid dynamics. V. Circumventing the Babuvska-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg. 59 (1) (1986) 85–99. [15] D. Forti, L. Ded`e, Semi-implicit BDF time discretization of the Navier-Stokes equations with VMS-LES modeling in a high performance computing framework, Comput. & Fluids 117 (2015) 168–182. [16] R. Haferssas, P. Jolivet, S. Rubino, Efficient and scalable discretization of the Navier–Stokes equations with LPS modeling, Comput. Methods Appl. Mech. Engrg. 333 (2018) 371–394. [17] I. Babuvska, Error-bounds for finite element method, Numer. Math. 16 (1970/1971) 322–333. [18] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers, Rev. Fran¸caise Automat. Informat. Recherche Op´erationnelle S´er. Rouge 8 (R-2) (1974) 129–151. 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [19] P. Hood, C. Taylor, Navier–Stokes equations using mixed interpolation, in: J. T. Oden, R. H. Gallagher, O. C. Zienkiewicz, C. Taylor (Eds.), Finite Element Methods in Flow Problems, University of Alabama in Huntsville Press, 1974, pp. 121–132. [20] V. Gravemeier, W. A. Wall, E. Ramm, A three-level finite element method for the instationary incompressible Navier- Stokes equations, Comput. Methods Appl. Mech. Engrg. 193 (15-16) (2004) 1323–1366. [21] V. Gravemeier, W. A. Wall, E. Ramm, Large eddy simulation of turbulent incompressible flows by a three-level finite element method, Internat. J. Numer. Methods Fluids 48 (10) (2005) 1067–1099. [22] V. John, S. Kaya, A. Kindl, Finite element error analysis for a projection-based variational multiscale method with nonlinear eddy viscosity, J. Math. Anal. Appl. 344 (2) (2008) 627–641. [23] T. Chac´on Rebollo, A term by term stabilization algorithm for finite element solution of incompressible flow problems, Numer. Math. 79 (2) (1998) 283–319. [24] L. He, L. Tobiska, The two-level local projection stabilization as an enriched one-level approach, Adv. Comput. Math. 36 (4) (2012) 503–523. [25] R. Becker, M. Braack, A two-level stabilization scheme for the Navier-Stokes equations, in: Numerical mathematics and advanced applications, Springer, Berlin, 2004, pp. 123–130. [26] L. R. Scott, S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (190) (1990) 483–493. [27] F. Hecht, New development in freefem++, J. Numer. Math. 20 (3-4) (2012) 251–265. [28] S. Badia, On stabilized finite element methods based on the Scott-Zhang projector. Circumventing the inf-sup condition for the Stokes problem, Comput. Methods Appl. Mech. Engrg. 247/248 (2012) 65–72. [29] M. Braack, E. Burman, Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method, SIAM J. Numer. Anal. 43 (6) (2006) 2544–2566. [30] P. Knobloch, G. Lube, Local projection stabilization for advection-diffusion-reaction problems: one-level vs. two-level approach, Appl. Numer. Math. 59 (12) (2009) 2891–2907. [31] F. E. Cellier, Continuous system modeling, Springer-Verlag, New York, 1991. [32] R. Codina, Stabilized finite element approximation of transient incompressible flows using orthogonal subscales, Comput. Methods Appl. Mech. Engrg. 191 (39-40) (2002) 4295–4321. [33] R. Codina, J. Blasco, Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales, Comput. Vis. Sci. 4 (3) (2002) 167–174. [34] R. Codina, A stabilized finite element method for generalized stationary incompressible flows, Comput. Methods Appl. Mech. Engrg. 190 (20-21) (2001) 2681–2706. [35] U. Wilbrandt, C. Bartsch, N. Ahmed, N. Alia, F. Anker, L. Blank, A. Caiazzo, S. Ganesan, S. Giere, G. Matthies, R. Meesala, A. Shamim, J. Venkatesan, V. John, Parmoon—a modernized program package based on mapped finite elements, Computers & Mathematics with Applications 74 (1) (2017) 74 – 88. [36] P. W. Schroeder, G. Lube, Divergence-free h(div)-fem for time-dependent incompressible flows with applications to high reynolds number vortex dynamics, J. Sci. Comput. 75 (2) (2018) 830–858. [37] M. Lesieur, C. Staquet, P. Le Roy, P. Comte, The mixing layer and its coherence examined from the point of view of two-dimensional turbulence, J. Fluid Mech. 192 (1988) 511–534. [38] B. J. Boersma, M. N. Kooper, F. T. M. Nieuwstadt, P. Wesseling, Local grid refinement in large-eddy simulations, J. Engrg. Math. 32 (2-3) (1997) 161–175. [39] E. Burman, Interior penalty variational multiscale method for the incompressible Navier-Stokes equation: monitoring artificial dissipation, Comput. Methods Appl. Mech. Engrg. 196 (41-44) (2007) 4045–4058. [40] V. John, An assessment of two models for the subgrid-scale tensor in the rational LES model, J. Comput. Appl. Math. 173 (1) (2005) 57–80. [41] E. Balaras, U. Piomelli, J. M. Wallace, Self-similar states in turbulent mixing layers, J. Fluid Mech. 446 (2001) 1–24. [42] M. Lesieur, O. M´etais, New trends in large-eddy simulations of turbulence, in: Annual review of fluid mechanics, Vol. 28, Annual Reviews, Palo Alto, CA, 1996, pp. 45–82. [43] K. Schneider, M. Farge, Numerical simulation of a mixing layer in an adaptive wavelet basis, Comptes Rendus de l’Acad´emie des Sciences - Series IIB - Mechanics-Physics-Astronomy 328 (3) (2000) 263 – 269. [44] E. O˜nate, A. Valls, J. Garc´ıa, Computation of turbulent flows using a finite calculus–finite element formulation, International Journal for Numerical Methods in Fluids 54 (68) (2007) 609–637. [45] C. R. Doering, J. D. Gibbon, Applied analysis of the Navier-Stokes equations, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 1995. [46] N. Ahmed, G. Matthies, Numerical study of SUPG and LPS methods combined with higher order variational time discretization schemes applied to time-dependent linear convection-diffusion-reaction equations, J. Sci. Comput. 67 (3) (2016) 988–1018. [47] G. Barrenechea, F. Valentin, Consistent local projection stabilized finite element methods, SIAM J. Numer. Anal. 48 (5) 1801–1825. [48] L. P. Franca, S. L. Frey, Stabilized Finite Element Methods: II. The incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 99 (1992) 209–233. [49] P. G. Ciarlet. The finite element method for elliptic problems. North-Holland Publishing Co., Amsterdam, 1978. Studies in Mathematics and its Applications, Vol. 4. 22