Sub-grid scale characterization and asymptotic behavior of multi-dimensional upwind schemes for the vorticity transport equations
SSub-grid scale characterization and asymptotic behavior of multi-dimensional upwindschemes for the vorticity transport equations
Daniel Foti
Department of Mechanical Engineering, University of Memphis, Memphis, TN, USA
Karthik Duraisamy
Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI, USA
We study the sub-grid scale characteristics of a vorticity-transport-based approachfor large-eddy simulations. In particular, we consider a multi-dimensional upwindscheme for the vorticity transport equations and establish its properties in the under-resolved regime. The asymptotic behavior of key turbulence statistics of velocity gra-dients, vorticity, and invariants is studied in detail. Modified equation analysis indi-cates that dissipation can be controlled locally via non-linear limiting of the gradientemployed for the vorticity reconstruction on the cell face such that low numericaldiffusion is obtained in well-resolved regimes and high numerical diffusion is realizedin under-resolved regions. The enstrophy budget highlights the remarkable ability ofthe truncation terms to mimic the true sub-grid scale dissipation and diffusion. Themodified equation also reveals diffusive terms that are similar to several commonlyemployed sub-grid scale models including tensor-gradient and hyper-viscosity mod-els. Investigations on several canonical turbulence flow cases show that large-scalefeatures are adequately represented and remain consistent in terms of spectral energyover a range of grid resolutions. Numerical dissipation in under-resolved simulationsis consistent and can be characterized by diffusion terms discovered in the modifiedequation analysis. A minimum state of scale separation necessary to obtain asymptoticbehavior is characterized using metrics such as effective Reynolds number and effectivegrid spacing. Temporally-evolving jet simulations, characterized by large-scale vorti-cal structures, demonstrate that high Reynolds number vortex-dominated flows arecaptured when criteria is met and necessitate diffusive non-linear limiting of vorticityreconstruction be employed to realize accuracy in under-resolved simulations. a r X i v : . [ phy s i c s . f l u - dyn ] F e b I. INTRODUCTION
Coherent structures are prominent in a wide range of turbulent flows including jets [1, 2], wakes [3, 4], and at-mospheric flows [5, 6]. The critical characteristics of such flows are defined to a large extent by the dynamics ofdominant coherent structures, pointing to large-eddy simulation (LES) as an ideal candidate model. In classical LES,a scale-separation is performed typically via a low-pass filter, and a sub-grid scale model is imposed to represent theimpact of the unresolved scales on the resolved scales. In practice, the sub-grid scale model provides a pathway forenergy to dissipate from the resolved scales because physical dissipation mechanisms (i.e., viscous dissipation at theKolmogorov scale) are unresolved. However, the accuracy of LES at high Reynolds numbers is significantly influencedby discretization errors [7–10], as these errors can be of a similar magnitude as the sub-grid scale model in practicalscenarios. Despite this fundamental obstacle, functional [11–18] and structural [19–25] LES sub-grid scale models areactively developed and applied to practical flows with varying degrees of success.An alternate view eschews the explicit modeling of sub-grid scales and instead focuses on tailoring numerical dissi-pation in the underlying discretization errors. Monotone integrated LES (MILES), first proposed by Boris et al. [26],utilizes functional reconstruction of the convective fluxes in a monotonicity-preserving finite volume scheme to inte-grate the effects of the sub-grid scale on the resolved scale. Results from similar numerical methods also were observedby Youngs [27] around the same time. A more broadly-defined methodology, referred to as implicit LES (ILES), re-quires a careful consideration of discretization errors. Typically, ILES employs a certain class of numerical methods,most notably monotonic upwind finite volume methods such as piece-wise parabolic [28], MPDATA [29], total vari-ation diminishing [30] and flux-corrected transport [31], which have implicit dissipation mechanisms [26, 32, 33]. Inpractice, ILES methods are employed on inertially dominated (e.g., high Reynolds number) dynamics and regularizethe under-resolved scales similar to the way in which shock-capturing, non-oscillating finite volume schemes use weaksolutions and satisfy the entropy condition.Thorough analysis of the schemes reveals the errors are often similar in form to certain sub-grid scale models [34].However, the form of the discretization errors does not necessarily has to be similar to a sub-grid scale model providedthe dissipation acts on the high wavenumber range. Drikakis et al [35] characterized the effect of numerical resolutionon the dissipation into two categories: a) When the solution is well-resolved and an adequate distinction between thestart of the inertial range and the dissipation scales exists, the numerical dissipation should not influence the largescales, and the eddies interacting with the largest scales should be reasonably resolved. The separation of scales shouldnot depend on the form of the dissipation; b) The second category involves inadequate numerical resolution, whichis often confronted in engineering applications. Given an inadequate separation between the large and dissipativescales, the numerical scheme should be designed to mimic the impact of dissipation on the large scales. Many ILESsolution fall under this category including isotropic turbulence [30, 36, 37], geophysical flows [29, 38], jets [26, 39, 40]and channel flows [40].In many vortex-dominated problems, especially those in which the vorticity distribution is compact [41–43] thevorticity-velocity formulation of the Navier–Stokes equations has the potential to be advantageous in comparisonto the pressure-velocity form [44, 45]. Despite the use of a wide array of Eulerian [46, 47], Lagrangian [48, 49]and mixed numerical implementations [50, 51] of the vorticity transport equations, modeling of unresolved dynamicsremains an outstanding issue. Further, while ILES has been prominently utilized in the context of the velocity-pressure formulation, the authors are not aware of literature investigating implicit sub-grid scale characterizations forschemes using the vorticity-velocity formulation. The goal of this work is to explore how tailored numerical schemesof vorticity-velocity formulation are suitable and can be characterized for ILES.The vorticity-velocity formulation of the incompressible vorticity transport equations (VTE), obtained by applyingthe curl operator to the momentum equations. While there is an increase to six (three vorticity and three veloc-ity) variables compared to four (three velocity and one pressure) variables in three dimensions, the formulation isadvantageous in flows with compact vorticity distributions typical of vortex-dominated flow, as the potential flowregions do not need to be part of the computational domain. In Eulerian methods, accurate boundary conditions canbe developed [43, 44, 52], allowing for efficient simulations in a compact domain. Further, the Poisson equation forpressure, normally a stiff equation, is replaced by the kinematic velocity-vorticity relationship.Numerical schemes [48, 53, 54] using the Lagrangian description of the flow field have seen considerable successand have applied the LES methodology through development of vorticity sub-grid scale models [49, 55]. Both Refs.[49, 55] introduced eddy-viscosity type sub-grid scale models for the vorticity-velocity formulation using the Lagrangiandescription. However, Lagrangian frameworks are not optimal for a wide range of turbulent flows. Furthermore, apriori tests of forced homogeneous, isotropic turbulence in Ref. [55] reveal that both vorticity convection and vortexstretching contribute to sub-grid scale dissipation, but the dissipation due to vortex stretching is inadequately capturedby the sub-grid scale model. On the other hand, relatively few Eulerian vorticity numerical schemes have beendeveloped [46, 56, 57]. Recently developed upwind finite volume methods for the VTE [45, 47, 58] have been shown toefficiently capture and preserve vortical structures with relatively coarse resolution (a few grid cells spanning vortexcores). Consistent integration of both the vorticity convection and vortex stretching terms is a prominent feature ofschemes in Refs. [45, 58]. While these schemes fit into the philosophy of ILES, an understanding of their sub-gridscale behavior has not been established in the context of turbulent flows.Margolin et al . [34] used modified equation analysis to examine leading order diffusive and dispersive error terms.Based on a Taylor series representation, the modified equation is the effective partial differential equation satisfiedby the numerical method. Their work specifically addresses finite volume approaches, and thus volume-averagedfiltering. The leading terms in the modified equation are used to elucidate the dissipative effects on the large resolvedscales. As the numerical resolution of the scheme increases, the dissipation should focus on a narrowing range ofhigh wavenumbers and not impact the large scales. Several studies [30, 36] emphasized the calculation of an effectiveviscosity or Reynolds number for a given grid resolution with resolved features. Since the filter and dissipationare based on the grid resolution, it is imperative to understand whether there is a sufficient separation of resolvedand dissipation scales. A minimum state of scale separation is necessary to reproduce high Reynolds numbers flowswith asymptotic turbulence statistics [36, 59]. In this work, we will examine the utility of a finite volume schemefor the VTE [45, 58], characterizes its ILES properties and provide a methodology to estimate effective Reynoldsnumbers. This particular formulation has a number of features pertinent for ILES: Consistent integration of vorticityconvection and vortex stretching terms, promoting stability and accuracy and the ability to represent and preservesharp gradients. Rigorous characterization of the sub-grid scale behavior in canonical turbulent flows enables the useof VTE-based schemes in complex vortex-dominated turbulent flows such as in rotorcraft and wind turbine wakes.We begin the study by introducing the governing and filtered equations in Sec. II. We provide a description of thenumerical scheme and analyze the implicit sub-grid scale model in Secs. III and IV, respectively. We detail the resultsof numerical experimentation of canonical turbulence flows in Sec. V. Finally, we conclude our work and providedetails for future studies in Sec. VI.
II. GOVERNING AND FILTERED EQUATIONS
We employ the vorticity-velocity formulation of the Navier-Stokes equations obtained by applying the curl operatorto the incompressible mass and momentum equations. In compact index notation, the unsteady, three-dimensionalincompressible VTE are as follows ( i, j = 1 , , ∂ω i ∂t + ∂∂x j ( u j ω i − u i ω j ) = f i + 1 Re ∂ ω i ∂x j ∂x j , (1)where u i is the velocity and ω i = (cid:15) ijk ∂u k /∂x j is the vorticity, the curl of the velocity where (cid:15) ijk is the Levi-Civitatensor. f i is the curl of the body force and Re = ULν is the Reynolds number defined by a characteristic velocity scale U , characteristic length L , and kinematic viscosity ν . The inviscid fluxes in Eqn. (1) can be rewritten in quasi-linearform as follows: ∂ω i ∂t + A ∂∂x ω i + A ∂∂x ω i + A ∂∂x ω i = 0 , (2)where the eigenvalues of matrix A are λ = { , u , u } , A are λ = { u , , u } , A are λ = { u , u , } . Theeigenvector matrices R i for A i are R = u u u , R = u u u , R = u u u . (3)These equations are degenerate in a hyperbolic sense. As an example, if any component of u equals zero, then theeigenvectors are linearly dependent. This presents difficulties for construction of stable upwind schemes. This issimilar to ideal magneto-hydrodynamics, for which Godunov [60] suggested the addition of a symmetrizing term.Following this idea and Ref. [58], an additional term is included to stabilize the governing equations. The equationsare thus modified in the form ∂ω i ∂t + ∂∂x j ( u j ω i − u i ω j ) + u i ∂ω j ∂x j = f i + 1 Re ∂ ω i ∂x j ∂x j , (4)Included with the VTE due to the incompressible assumption, is the solenoidality of both the velocity and vorticityfields: ∂u i ∂x i = ∂ω i ∂x i = 0 . (5)We refer to Eqn. (4) as the modified VTE because the final term on the LHS is an additional term, which isproportional to the divergence of vorticity. Analytically, the term is zero (Eqn. (5)), but numerically, this stabilizesthe hyperbolic system of equations. The modification in Eqn. (4) stabilizes the equations with eigenvalues λ i = { u i , u i , u i } and the eigenvectors of the canonical basis vector for R [58]. A similar modification is used in themagneto-hydrodynamic equations to enforce a divergence-free magnetic field [61], however, in our context we employthe modification for numerical stability.The vorticity-velocity formulation in three dimensions has six unknowns. While the three vorticity variables areobtained through the VTE, a supplemental equation is needed to obtain the velocity induced by the vorticity. Thevorticity-velocity relationship is written in the form of a Poisson equation as follows: ∂ u i ∂x j ∂x j = − (cid:15) ijk ∂ω k ∂x j . (6) II.1. The Filtered vorticity transport equations
By spatial filtering the VTE in Eqn. (4), we obtain the following filtered VTE: ∂ (cid:101) ω i ∂t + ∂∂x j ( (cid:101) u j (cid:101) ω i − (cid:101) ω j (cid:101) u i ) + (cid:101) u i ∂ (cid:101) ω j ∂x j = (cid:101) f i + 1 Re ∂ (cid:101) ω i ∂x j ∂x j − ∂τ ij ∂x j , (7)where (cid:101) · indicates the spatial filtering or a resolved quantity and the τ ij is the sub-grid scale (SGS) vorticity stressdue to the filtering operation. We note that filtering is not performed through an explicit filter function, and weassume that the filtering operation and the derivatives commute. The grid cell of the mesh at a size ∆ is an implicitphysical-space sharp cut-off filter where the velocity and vorticity fluctuations can only be resolved at a size greaterthan ∆. The SGS torque, the divergence of the SGS vorticity stress, which accounts for the unresolved velocity andvorticity fluctuations, is defined as ∂τ ij ∂x j = ∂∂x j [( (cid:103) u j ω i − (cid:101) u j (cid:101) ω i ) − ( (cid:103) u i ω j − (cid:101) u i (cid:101) ω j )] + (cid:34) (cid:94) u i ∂ω j ∂x j − (cid:101) u i ∂ (cid:101) ω j ∂x j (cid:35) . (8)The SGS torque term is composed of three different terms: 1.) the vorticity transport by unresolved velocity fluc-tuations, 2.) the unresolved vortex stretching, and 3.) the unresolved contributions for the vorticity divergencemodification in Eqn. (4). Unlike previous LES schemes for the VTE [49, 55] where there is not an additional vorticitydivergence term, the SGS vorticity stress tensor is not purely anti-symmetric. Here, we attempt to rearrange the SGSterms. The first two terms in Eqn. (8) can be combined into a single anti-symmetric tensor as follows:( τ ij ) a = ( (cid:103) u j ω i − (cid:101) u j (cid:101) ω i ) − ( (cid:103) u i ω j − (cid:101) u i (cid:101) ω j ) . (9)However, the third term introduced through the vorticity divergence modification cannot be readily decomposed intoan SGS vorticity stress tensor and as such it remains a separate term of the SGS torque: (cid:18) ∂τ ij ∂x j (cid:19) d = (cid:94) u i ∂ω j ∂x j − (cid:101) u i ∂ (cid:101) ω j ∂x j . (10)However, through manipulation, this term can be rewritten as follows: ∂ ( τ ij ) m ∂x j = (cid:18) ∂ (cid:103) ω j u i ∂x j − ∂ (cid:102) ω j (cid:101) u i ∂x j (cid:19) − (cid:0) (cid:93) ω j s ij − (cid:102) ω j (cid:102) s ij (cid:1) , (11)where s ij = ( ∂ j u i + ∂ i u j ) is the strain rate. We can obtain another two portions of the SGS terms. The firstdivergence modification contains terms that are similar to the third and fourth terms in Eqn. (9):( τ ij ) d = (cid:103) u i ω j − (cid:101) u i (cid:101) ω j (12)While the other term contains the strain rate: ∂ ( τ ij ) s ∂x j = − (cid:0) (cid:93) ω j s ij − (cid:102) ω j (cid:102) s ij (cid:1) . (13)The practical purpose of these terms in the LES methodology is to provide a pathway of energy transfer between theresolved scales to the unresolved (sub-grid) scales, which is enabled by ensuring additional dissipation. u ∆ t u ∆ t u ∆ tu ∆ t u ∆ tx i − ,j x i + ,j x i − ,j x i + ,j x i − ,j x i + ,j = + Flux F G FIG. 1. Sketch of multi-dimensional scheme using both normal F and transverse G fluxes . III. MULTI-DIMENSIONAL UPWIND FINITE VOLUME
We employ a multi-dimensional upwind finite volume approach to numerically integrate the VTE. This approach,within the philosophy of the ILES, is similar to Lax-Wendroff-like schemes where upwind differences are corrected bysecond-order transverse flux corrections to produce a solution that is second-order accurate in space and time. Ourmulti-dimensional scheme is solved in a series of dimensional sweeps. The total flux across a cell face in each directionis the sum of several fluxes accounting for both normal and transverse directional fluxes. A simplified two-dimensionalflux is demonstrated in Fig. 1. The updated vorticity is computed as shown (for clarity vectors are shown in boldand ˜ · removed from all variables): ω n +1 i,j,k = ω ni,j,k − ∆ t ∆ x (cid:104)(cid:16) ( F l ) ni + ,j,k − ( F r ) ni − ,j,k (cid:17) − (cid:16) ( G ) ni + ,j,k − ( G ) ni − ,j,k (cid:17)(cid:105) − ∆ t ∆ x (cid:104)(cid:16) ( F l ) ni,j + ,k − ( F r ) ni,j − ,k (cid:17) − (cid:16) ( G ) ni,j + ,k − ( G ) ni,j − ,k (cid:17)(cid:105) − ∆ t ∆ x (cid:104)(cid:16) ( F l ) ni,j,k + − ( F r ) ni,j,k − (cid:17) − (cid:16) ( G ) ni,j,k + − ( G ) ni,j,k − (cid:17)(cid:105) , (14)where n is the time iteration and i, j, k are the indices of the grid cell centers in three directions. The flux functions F l , F r are normal directional (with respect to the cell faces) fluxes and are obtained by solving a generalized Riemannproblem developed for the VTE [58]. In order to account for the fluxes traveling oblique to the cell faces whilesimultaneously increasing numerical stability and accuracy, transverse fluxes are included. The transverse directionalflux functions G are computed using the flux-based wave propagation approach [45]. All flux functions are stored atthe cell faces. Note that both left and right normal flux functions F li and F ri , respectively, are stored at each cellface, while a single transverse flux function G i is stored at each cell face.First, we detail the normal directional fluxes F li , F ri , which are determined by allowing the vorticity to bediscontinuous at the cell face. The sharp discontinuity-capturing scheme can be beneficial for simulating turbulencethat is dominated by large coherent structures. In flows such as these, regions that are dominated by coherentstructures can be efficiently captured by compact vorticity variables. On the other hand, in under-resolved turbulenceregions, the scheme adds additional numerical dissipation (more discussion on this in section IV).The solution to the generalized Riemann problem for vorticity, which is allowed to be discontinuous across the cellface, begins by reconstructing the vorticity stored at the cell center on the cell face. For simplicity, we will focuson the Riemann problem at cell face ( i + 1 / , j, k ) and translated to x = 0. The initial condition for the Riemannproblem is ω ( x ,
0) = (cid:26) ω L + (cid:0) x + ∆ x (cid:1) s L , if x < ω R + (cid:0) x − ∆ x (cid:1) s R , if x > , (15)where ( · ) L at i, j, k , ( · ) R is at i + 1 , j, k , ∆ x is the cell size, s is the slope of the vorticity. Achieving second order accu-racy via slope reconstruction is beneficial in maximizing the separation of scales for ILES. In areas of smooth vorticitya second-order accurate central difference can be used, however, as with most second-order accurate discontinuity-capturing schemes, an appropriate slope needs to be used, and a limiter is employed. A limiter reduces the slopecalculation to a first-order accurate difference and effectively adds numerical diffusion.Integrating the vorticity transport equations, equation (4), with the vorticity reconstruction, equation (15), overthe space (cid:2) − ∆ x , ∆ x (cid:3) × [0 , ∆ t ], where ∆ t is the time step, a left moving flux F l and right moving flux F r areobtained at the cell face as follows ( F r,l and F r,l follow similarly): F l = (cid:26) u ˆ ω L − u ˆ ω L + u s L ∆ x , if u ≥ u ˆ ω R − u ˆ ω R + u s R ∆ x + u (cid:2) (cid:0) s L − s R (cid:1) u ∆ t + (cid:0) ˜ ω L − ˜ ω R (cid:1)(cid:3) , if u < , (16) F r = (cid:26) u ˆ ω L − u ˆ ω L − u s R ∆ x − u (cid:2) (cid:0) s L − s R (cid:1) u ∆ t + (cid:0) ˜ ω L − ˜ ω R (cid:1)(cid:3) , if u ≥ u ˆ ω R − u ˆ ω R − u s R ∆ x , if u < , (17)where for conciseness, additional vorticity reconstruction variables are created:ˆ ω L = ω L + 12 s L (∆ x − u ∆ t ) , ˆ ω R = ω R − s R (∆ x + u ∆ t ) , ˜ ω L = ω L + 12 s L ∆ x , ˜ ω R = ω R − s R ∆ x . Next, we detail the computation of the transverse fluxes. In this scheme we pursue a flux-based wave decompositionintroduced in Ref. [62], where the flux difference, F l − F r , is rewritten as a linear combination of eigenvectors. Thischoice is motivated because the solution to the generalized Riemann problem is known, is linearly independent ineach direction, and can be used directly without costly manipulations. In this implementation, the fluxes, F l and F r ,only need to be computed once (per time-step), and the transverse fluxes are evaluated with that solution. The fluxdifference is decomposed into f -waves Z p , the flux wave, as follows: F l − F r = m (cid:88) p =0 β p r p ≡ m (cid:88) p =1 Z p , (18)and β = R − (cid:0) F l − F r (cid:1) , (19)where the eigenmatrix R is simply the identity matrix [45].From Eqn. (19) and the eigensystem, we obtained a simple relationship for the f -waves Z p = F l − F r exactly.The relationship allows for the multi-dimensional wave propagation corrections to be implemented equivalently to athree-dimensional advection equation described in Ref. [63] if the f -waves are normalized by the eigenvalues λ i .In a multi-dimensional problem, the fluxes propagate in the transverse x - and x -direction depending on the wavespeed given by the eigenvalues of the Jacobian, which are exactly the velocities u i at the cell face. The transverse fluxesare implemented by using the f -waves and transverse velocity at the grid cell ( i, j, k ) and updating the transverse fluxesin surrounding cells. The most important transverse fluxes are the two-dimensional fluxes given by the following:( G ) nI,J − ,k = −
12 ∆ t ∆ x u Z , ( G ) nI,j,K − = −
12 ∆ t ∆ x u Z , (20)where n is the time step iteration and the indices for the grid cell influenced by the transverse flux are given by I = (cid:26) i, if u > i − , if u < , J = (cid:26) j + 1 , if u > j, if u < , K = (cid:26) k + 1 , if u > k, if u < . Additional transverse fluxes are employed, which account for three-dimensional fluxes as well as higher-order cor-rections. For a detailed implementation see Foti and Duraisamy [45].
IV. MODIFIED EQUATION ANALYSIS
In the following, we present the modified equation analysis (MEA) for the multi-dimensional generalized Riemannproblem-based upwind finite volume scheme. MEA analysis was first proposed in Ref. [64] and subsequently usedto characterize the unresolved scales for ILES [32, 34]. The intuition is that in certain classes of numerical methods,the effects of the unresolved scales can be represented by the truncation error of the discretization. This analysismust be carefully exploited when considering high Reynolds number turbulent flows discretized with upwind finitevolume schemes. Slope limiting with a multi-dimensional scheme yields a complex, nonlinear scheme. In particular,the dissipation is proportional to the multi-dimensional interfacial wave jumps. It must also be recognized that inunder-resolved flows, the leading order terms may not necessarily be the best approximator of the truncation terms.Due to the complexity of including the exact form of the slope limiter employed in the vorticity reconstruction in thescheme, two limiting cases of the modified equation are analyzed: 1.) the vorticity at the cell interface is smooth anda second-order central difference can be employed, and 2.) the vorticity at the cell interface is discontinuous and alimiter switches the slope to a first-order upwind difference.The procedure to develop the modified equation begins with expanding the scheme with Taylor series expansions[33, 64]. For example, the scheme in Eqn. (14) contains the following, which can be substituted with a series expansion: ω n +1 i − ω ni ∆ t = ∂ω i ∂t + ∆ t ∂ ω i ∂t + ∆ t ∂ ω i ∂t + O (∆ t ) , (21)where n is the vorticity at t and n + 1 is the vorticity at t + ∆ t . Similarly, all terms in the flux functions in Eqns.(16), (17), (20) can be substituted with Taylor series expansions in terms of vorticity as a function of time t or space x . The accumulation of all expansions in the scheme can be manipulated to include terms in Eqn. (4) as follows: ∂ω i ∂t + ∂∂x j ( u j ω i − u i ω j ) + u i ∂ω j ∂x j = ν ∂ ω i ∂x j ∂x j + T ij , (22)where T ij is the term that includes all second-order accurate and higher terms from the series expansions not included inthe VTE. The remainder is manipulated to substitute temporal derivative terms with spatial terms using the governingequations. As such, the remainder T ij is a complex function that contains many high order spatial derivatives of thevorticity and velocity. In the philosophy of ILES, the leading order terms in the expansion term T ij are related to theSGS vorticity stress R ij due to truncation as follow: T ij = ∂∂x j R ij + O (∆ x ) . (23) T ij acts as an implicit sub-grid scale model, which is completely dependent on resolved variables. Note that thereis a distinction between the SGS vorticity stress R ij due to truncation and theoretical SGS vorticity stress τ ij dueto filtering. In what follows we demonstrate that by changing the form of the slope employed in the vorticity recon-struction on the cell face, the modified equation can be tailored.First, we examine terms that remain after the VTE is subtracted from the Taylor series expansion of the scheme insmooth regions and a second-order accurate central difference is employed, i.e. s l = 1 / (2∆ x ) ( ω ( x + ∆ x ) − ω ( x − ∆ x )).In turbulent flows in which we are interested, a smooth spatial region of vorticity often corresponds to a resolved (ornearly resolved) region dominated by a large-scale coherent vortical structure.In the modified equation analysis, the lowest order remaining terms are second-order accurate terms in space andare considered to have the largest influence on the scheme. Moreover, these terms can be collected to implicitlyprovide a model for the SGS torque. Eqn. (24) shows the second-order accurate terms that make up the implicit SGStorque. For clarity in presentation, only the x and x directions are shown ( j = 1 , x directionfollows the same form: ∂R j ∂x j = ∆ (cid:18) u ∂ ω ∂x − u ∂ ω ∂x − ∂ u ∂x ∂ω ∂x + 112 u ∂ ω ∂x − u ∂ ω ∂x − ∂ u ∂x ∂ω ∂x (cid:19) ,∂R j ∂x j = ∆ (cid:18) u ∂ ω ∂x − u ∂ ω ∂x − ∂ u ∂x ∂ω ∂x + 112 u ∂ ω ∂x − u ∂ ω ∂x − ∂ u ∂x ∂ω ∂x (cid:19) , (24)where ∆ i is the grid cell size in the i th direction. All terms are dispersive (containing odd-ordered derivatives ofvorticity) in nature indicating that in smooth regions where a second-order central difference is used, there are nosecond-order accurate implicit diffusive terms (even-ordered derivatives of vorticity), which can act as numerical routesfor implicit SGS dissipation. We take a closer look at the diagonal terms ( ∂ R and ∂ R ) of the SGS torque toobtain some insight into the effects of the modification of the vorticity transport equations and the effect of theanti-symmetric and non-symmetric SGS terms (due to Eqn. (9) and Eqn. (10), respectively).By rewriting terms in Eqn. (24) that are associated with the diagonal terms of R ij , an anti-symmetric part of theSGS torque associated with Eqn. (9), ∂ j R aij , can be formed and shown to be equal to zero in the following: ∂R a ∂x = 112 ∆ x u ∂ ω ∂x −
112 ∆ x u ∂ ω ∂x = 0 ,∂R a ∂x = 112 ∆ x u ∂ ω ∂x −
112 ∆ x u ∂ ω ∂x = 0 . (25)The remaining leading order terms can collected into the following: ∂R d ∂x = −
16 ∆ u ∂ ω ∂x −
18 ∆ ∂ u ∂x ∂ω ∂x ,∂R d ∂x = −
16 ∆ u ∂ ω ∂x −
18 ∆ ∂ u ∂x ∂ω ∂x , (26)where ∂ j R dij is due to the modification term added to the vorticity transport equation. Terms associated with theoff-diagonals of R ij contain dispersive terms. This modified equation of the scheme near a smooth vorticity field showsthat there is a limiting case that can be used to reduce the amount of numerical dissipation added by the numericalmethod, which can control the sub-grid scale dissipation and energy transfer.The next limiting case employs a limiter to reduce the order of the slope in regions in which vorticity gradients maybe very large. In numerical simulations of turbulence, large gradients and discontinuities arise in under-resolved regionsof the flow where gradients caused by eddies the size of the grid cell or larger cannot be smoothed by smaller scaleeddies that are physically present in the flow but are not numerically captured. Towards this end, a forward/backwarddifference is employed. For example, the slope of the vorticity is calculated as s l = 1 / ∆ x ( ω ( x ) − ω ( x − ∆ x )). Theseeddies transfer energy in a turbulent flow, and the reduction in order of the slope essentially introduces diffusion intothe solution in order to maintain non-oscillatory behavior. The lowest-order terms remaining from the modifiedequation are shown in the following (again we limit the result to two dimensions for clarity but the x -direction hasthe same form): ∂R j ∂x j = ∆ u ∂ ω ∂x + ∆ (cid:18) − ∂ u ∂x ∂ω ∂x − u ∂ ω ∂x (cid:19) +∆ u ∂ ω ∂x + ∆ (cid:18) − ∂ u ∂x ∂ω ∂x + 14 ∂u ∂x ∂ ω ∂x − ∂u ∂x ∂ ω ∂x + 13 u ∂ ω ∂x − u ∂ ω ∂x (cid:19) , (27) ∂R j ∂x j = ∆ u ∂ ω ∂x + ∆ (cid:18) − ∂ u ∂x ∂ω ∂x + 14 ∂u ∂x ∂ ω ∂x − ∂u ∂x ∂ ω ∂x + 13 u ∂ ω ∂x − u ∂ ω ∂x (cid:19) +∆ u ∂ ω ∂x + ∆ (cid:18) − ∂ u ∂x ∂ω ∂x − u ∂ ω ∂x (cid:19) , (28)which include a single first-order accurate in space term and several second-order accurate terms. We can readily seethat this case includes both diffusive and dispersive terms for vorticity. This indicates that limiting the slope of thevorticity is necessary to add diffusive terms to include dissipation to solve turbulent flows through an ILES. As withthe first limiting case, a few of the terms in Eqn. (27) can be rewritten to show a zero diagonal of an anti-symmetricportion of the SGS torque in the following: ∂R a ∂x = ∆ (cid:18) ∂u ∂x ∂ ω ∂x − ∂u ∂x ∂ ω ∂x + 13 u ∂ ω ∂x − u ∂ ω ∂x (cid:19) = 0 . (29)Similarly, terms can be rearranged in ∂ R a = 0 and ∂ R a = 0. The rest of the terms in the diagonals are accumulatedin the part of the SGS torque corresponding to the introduction of the vorticity divergence term modification, whichincludes a first-order accurate diffusive term.The analysis reveals that each component in the SGS vorticity stress tensor R ij contains diffusive terms. In thefollowing we select terms from Eqn. (27) to explore how the diffusion is present in the SGS vorticity stress tensor andits similarity to well-known SGS models. We emphasize that similarities to explicit SGS models are employed to justenhance our understanding of the ILES dissipation. In complex turbulent flows, the process of limiting the slope andmulti-dimensionality can affect the actual leading terms of the modified equation such that they may not preciselymimic the action of explicit SGS models. We show that three separate diffusive SGS vorticity stress mechanisms arepresent in the modified equation and can be summed as follows: R ij = R gij + R hij + R tij . (30)The first term is constructed as follows: ∂R g ∂x = −
18 ∆ ∂ u ∂x ∂ω ∂x + 14 ∆ ∂u ∂x ∂ ω ∂x , (31)We can manipulate the terms by introducing another dispersion term and integrating the SGS torque to obtain asingle form for the off-diagonal elements in the SGS vorticity stress tensor, which can be written in index notion asfollows: R gij = 14 ∆ ∂u j ∂x k ∂ω i ∂x k , (32)where ∆ is the grid spacing in the i th direction. The form of the Eqn. (32) can be readily seen as a term similar totensor-gradient models [15, 65, 66] for the sub-grid scale. These models have been shown to be not overly dissipativeand are often paired with an eddy viscosity model [65, 66]. They also provide both dissipation and backscatter, animportant aspect of turbulence, which is not possible in eddy viscosity-type models. However, when employed alone,the tensor-gradient model can be unstable [66].An alternate form for the numerical dissipation can be explored with another group of dispersive and diffusive termsas shown in the following: ∂R h ∂x = −
12 ∆ u ∂ ω ∂x −
14 ∆ ∂u ∂x ∂ ω ∂x , (33)where these terms again can be manipulated and integrated. The second numerical dissipation term for the SGSvorticity stress has the form : R hij = −
14 ∆ u j ∂ ω i ∂x k ∂x k . (34)This term has the form that is similar to hyper-viscosity models [49, 67], which have the form of ( − n +1 ∇ n u , where n >
1. Ref. [67] used the hyper-viscosity model for the SGS stress tensor with ∇ (cid:101) s ij , which was added to a standardSmagorinsky model and motivated by SGS dissipation of enstrophy. Hyper-viscosity models with increasing n haveless effect on dissipation and ‘bottleneck’ effects shown on the turbulence spectra [68]. However, paired with thetensor-gradient model, it may provide less dissipation than is observed with standard mixed models.With these two terms, all of the diffusive terms in the modified equation and SGS torque are accounted exceptfor the single first-order diffusive term, which appears in both the diagonal and non-diagonal elements of the SGSvorticity stress tensor. This term can be integrated to show another form of the numerical dissipation: R tij = 12 ∆ u j ∂ω i ∂x j . (35)These terms show that in the limit of a slope reconstruction with a first-order slope, there are several modes foradditional numerical diffusion that can be used to increase the energy transfer from the resolved scales to the sub-gridscales and add dissipation.The two limiting cases for the present scheme show that the dissipation of the scheme can be controlled usingwell-designed limiters for specific flows. The present scheme is developed to conserve coherent vortical structuresin flows. When regions of the flow field are dominated by large coherent structures, which are not subject to theinertial range energy transfer determined by the SGS stresses, the vorticity tends to be smooth and less implicitdissipation is included. However, in regions of small-scale vortical structures where the flow field is under-resolvedand discontinuities in the vorticity are present, the dissipation needs to be increased to account for the energy transferto the sub-grid scales. This formulation provides a plausible framework to simulate large-scale features of turbulentflows without providing dissipation at all flow scales, a common problem for Smagorinsky-like explicit methods. V. NUMERICAL EXPERIMENTATION
In this section, we discuss results obtained through numerical experimentation of several canonical flows of turbu-lence in a periodic box. We investigate a Taylor-Green vortex in Sec. V.1, forced isotropic turbulence in Sec. V.2,and a temporally evolving jet in Sec. V.3. All simulations employ the multi-dimensional wave propagation approachdescribed above.0
V.1. Taylor-Green vortex
The Taylor-Green vortex is a well-studied flow field [69–71], which encompasses large structures, transition, anddecaying turbulence, and is used to assess the ability of the scheme to capture the complexities of vortex stretchingand breakdown using under-resolved grids [36, 72, 73]. The test case is used to demonstrate how the present schemecan be employed to emulate the dominant sub-grid scales, most importantly dissipation, from the initial conditionsof large coherent structures through a transition to turbulence at high Re. The initial velocity flow field for theTaylor-Green vortex is given by u ( x , x , x ) = u cos( kx ) sin( kx ) cos( kx ) u ( x , x , x ) = − u sin( kx ) cos( kx ) cos( kx ) (36) u ( x , x , x ) = 0where u = 1 is the velocity scale and k = 1 is the length scale of the flow field. Using these scales, the Re = u /νk = 1 /ν , where ν is the dynamic viscosity specified for each simulations. Simulations are designed in a periodiccube with length 2 π with a uniform grid. We start by exploring how well the scheme captures the pertinent physicsat Re = 1600 with a non-dissipative second-order central difference (referred to as no-limiter) and a high dissipativelimiter, the minmod limiter. The first limiter is directly related with the results from the modified equation in sectionIV. On the other hand, the modified equation for the first-order backward slope detailed in section IV is only similarto the minmod limiter because the minmod limiter chooses between first-order backward and forward differences andremoves slopes where the backward and forward differences have opposite slopes.Figure 2 shows the energy spectra of simulations at several grid resolutions ( N = 64 ,
128 and 256) at two timeinstances tu k = 5, during the transition to turbulence as the large-scale structures are breaking down, and tu k =10, near the maximum dissipation. Additionally, a spectral method solver [21] is employed to provide direct numericalsimulation (DNS) comparisons. The DNS employs 256 spectral modes to fully resolve the flow. In Fig. 2(a), the lowwavenumber region of the energy spectra in all simulations is able to be captured accurately compared the DNS resultsfrom the spectral method regardless of the slope limiter employed. The decline in energy in the spectrum for eachsimulation corresponds to the grid resolution where the lowest resolution begins to fall off at lowest wavenumbers.This is expected for the under-resolved simulations. The highest resolution case, while using the same number ofdegrees of freedom as the DNS, approaches the DNS solution but has slightly less energy in inertial range modes. Testcases that employ the minmod limiter have lower energy in modes starting near a wavenumber k/k = 10 comparedto the no-limiter limiter. It becomes noticeable in the lower resolution test cases at lower wavenumbers, however,they all occur in the inertial range. This is due to the dissipative nature of the minmod limiter, which allows for moreimplicit dissipation especially in modes that have a wavenumber comparable to the inverse size of the grid cell, 1 / ∆ x .Figure 2(b) shows similar behavior for the energy spectra; however, the lowest grid resolution cases show an increasedenergy at low wavenumbers. This can be attributed to the evolution of the flow at this grid resolution because thenumber of necessary modes needed to provide an acceptable solution are higher than the number of modes resolved.The Taylor-Green vortex is simulated as several higher Re = 2000 , , ∞ , inorder to investigate the scheme in cases where the impact of viscosity is low or, in the case of Re = ∞ , nonexistent.Figures 3(a) and (b) show the energy spectra of each Re case at tu k = 5 and tu k = 10, respectively, at twodifferent grid resolutions. All simulations show similar low and high wavenumber behavior for all Re at particular gridresolution. All cases regardless of Re or grid resolution, the largest scales are similar. In fact, the Reynolds numberhas little effect on the low wavenumbers, while at the high wavenumbers, the energy in these modes increase with Re.This is due to the diminishing effects of the viscosity. However, the spectra clearly show that the scheme provides apathways for energy dissipation because there is no build up of energy in high wavenumber modes.Next, we focus on the dissipation spectra of the simulations, which is directly obtained from the vorticity flowfield. Figure 4(a) shows the dissipation spectra at tu k = 5 has several low wavenumber peaks that are well capturedby all simulations with N = 128 or 256 while there is lower dissipation in most modes for the N = 64 simulations.Overall, compared to the energy spectra, the dissipation spectra is well captured at all wavenumbers especially forthe N = 256 test cases. This is an indication that these simulations can accurately resolve most of the dissipation.Near the peak dissipation at tu k = 10, Fig 4(b) shows that grid resolution and the slope limiter start to have moreeffect on the resolved dissipation. The N = 256 no-limiter case dissipation spectrum is very comparable to the DNSdissipation spectrum at most wave numbers, while the minmod limiter test case at the same resolution reveals thatit can only capture the dissipation spectra at lower wavenumbers. A general trend shows that the no-limiter caseshave slightly higher dissipation modes. This is an indication that the no-limiter test cases can resolve slightly moredissipation compared to the minmod limiter, but we need to investigate further. While the no-limiter can captureslightly more dissipation, it also does not add any numerical dissipation on the order of ∆ x , so it may not be ableto dissipate sufficient energy. The minmod case shows slightly lower dissipation modes in the spectra, but there is1 FIG. 2. Energy spectra E ( k ) of simulations at several under-resolved grid resolutions at Re = 1600 at (a) tu k = 5 and (b) tu k = 10. Red: no-limiter, Black: minmod limiter. The magenta dotted line is the energy spectra from DNS results. numerical dissipation that is not accounted for in this metric.The total dissipation, including numerical dissipation, (cid:15) = dK/dt , where K = (cid:104) u i u i (cid:105) ( (cid:104)·(cid:105) indicates averagedover the computational domain) is the kinetic energy, is integrated over the entire domain and recast into a non-dimensional form, (cid:15) ∗ = (cid:15)/k u . Figure 5(a) shows the temporal evolution of the dissipation for the N = 256 test caseswhere both the no-limiter and minmod limiter cases provide reasonable results compared to the DNS dissipation,which is filtered using several sharp-cutoff spectral filters with widths including ∆ = 2∆ x , x , and 8∆ x . The filteredDNS dissipation (cid:101) (cid:15) is the temporal derivative of the filtered kinetic energy (cid:101) k . The filtered kinetic energy is obtainedby integrating over all wavenumbers in the sharp-cutoff spectral filtered energy spectrum with different filter widthsobtained from the DNS at each time step. While there are some discrepancies the peak dissipation is capturedreasonably well for both cases and is comparable to the DNS filtered solution at ∆ = 2∆ x . Figure 5(b), which showsthe total dissipation for the N = 128 cases, reveals that the minmod limiter case out performs the no-limiter case. Theminmod limiter case can reasonably capture the peak dissipation, which is comparable to the N = 256 cases and thefiltered DNS dissipation at ∆ = 4∆ x . Figure 5 (c) shows total dissipation for the N = 64 test cases. Neither resultsare as reasonable as the higher resolution cases indicating that very under-resolved cases may be overly dissipative incertain regimes such as near the maximum dissipation of the Taylor-Green vortex. On the other hand, the numericaldissipation in a regime of strong coherent vortex interaction at short times or decaying turbulence at long times aremore reasonably captured in very under-resolved cases. However, this low resolution case demonstrates evidence ofSGS backscatter present in numerics as the numerical dissipation becomes negative around tu k = 15. These resultsas well as the dissipation spectra indicate that when the flow is reasonably under-resolved, the numerical dissipationis necessary. The adequacy of the resolution will be explored below and is imperative to understand when employingILES in general. In our scheme, the numerical dissipation has a form that is qualitatively similar to physically devisedmodels, an aspect that offers insight in understanding the implicit SGS dissipation mechanisms.Simulations at several higher Re = 1600 , , , ∞ reveal the effects of the numerical dissipationwith diminishing viscosity for N = 128 and N = 64 in Figs. 6(a) and (b), respectively. For the N = 128 betweenthe initial condition and maximum dissipation (0 < tu k < ∞ shows the largest discrepancies between the ILES cases and the DNS and suggest that includingsome physical viscosity dissipation helps stabilize the scheme. Figure 6(b) shows the N = 64 cases with increasingRe. The cases corroborate the findings in Fig. 5(c) that a severely under-resolved case does not necessarily capturethe dissipation well for all cases. This will be further explored in further below.2 FIG. 3. Energy spectra E ( k ) of simulations at several under-resolved grid resolutions N = 128 (solid line) and N = 64 (dashedline) at Re = 1600 , , , ∞ at times (a) tu k = 5 and (b) tu k = 10.FIG. 4. Dissipation spectra E ( k ) from the resolved vorticity of simulations at several under-resolved grid resolutions at Re = 1600 at (a) tu k = 5 and (b) tu k = 10. Red: no-limiter, Black: minmod limiter. The magenta dotted line is thedissipation spectra from DNS results. Next, we provide further analysis to compare the numerical dissipation of the scheme to the form of dissipationthat is obtained through the modified equation analysis in section IV. Similar analysis [67, 74, 75] on how the resolvedvorticity field affects SGS dissipation transport has provided insights into the dynamics of the flow. We start withpresenting the transport equation for the resolved enstrophy (cid:101) E = (cid:104) (cid:101) ω i (cid:101) ω i (cid:105) that is obtained by multiplying Eqn. (4)3 FIG. 5. The filtered dissipation (cid:101) (cid:15) ∗ as a function of time for (a) N = 256, (b) N = 128, and (c) N = 64. Red dashed line:no-limiter, Black solid line: minmod limiter. The magenta line is the dissipation from DNS results spatially filtered with asharp-cutoff spectral filter with several filter widths ∆ = 2∆ x (solid line), 4∆ x (dashed line), and 8∆ x (dotted line).FIG. 6. The filtered dissipation (cid:101) (cid:15) ∗ as a function of time for (a) N = 128 and (b) N = 64 at Re = 1600 , , , ∞ . The magenta line is the dissipation from DNS results spatially filtered with a sharp-cutoff spectral filter with several filterwidths ∆ = 2∆ x (solid line), 4∆ x (dashed line), and 8∆ x (dotted line). by ω i and filtering. The resolved enstrophy transport equation is given by the following: ∂ (cid:101) E ∂t = − (cid:101) u j ∂ (cid:101) E ∂x j (cid:124) (cid:123)(cid:122) (cid:125) I + (cid:101) ω i (cid:101) ω j (cid:101) s ij (cid:124) (cid:123)(cid:122) (cid:125) II + ν ∂ (cid:101) E ∂x j ∂x j (cid:124) (cid:123)(cid:122) (cid:125) III − ν ∂ (cid:101) ω i ∂x j ∂ (cid:101) ω i ∂x j (cid:124) (cid:123)(cid:122) (cid:125) IV − (cid:101) ω i (cid:101) u i ∂ (cid:101) ω j ∂x j (cid:124) (cid:123)(cid:122) (cid:125) V − (cid:101) ω i ∂∂x j ( R ij ) (cid:124) (cid:123)(cid:122) (cid:125) VI , (37)which relates the temporal change in the resolved enstrophy to spatial changes in the vorticity, velocity, and SGSvorticity stress tensor. The equation contains several mechanisms that balance the temporal change in resolvedenstrophy: I.) convection, II.) amplification by vortex stretching, III.) diffusion by viscous effects, IV.) dissipation byviscous effects, V.) diffusion by divergence modification, and VI) SGS dissipation and SGS diffusion. Term VI can beexpanded into two terms as follows: (cid:101) ω i ∂∂x j ( R ij ) = ∂∂x j ( (cid:101) ω i R ij ) − R ij ∂ (cid:101) ω i ∂x j , (38)4 FIG. 7. (a) The temporal evolution of terms in the enstrophy transport equation for the N = 64 test case with minmodlimiter. (b) The temporal evolution of the calculated SGS dissipation and diffusion and the temporal evolution of the SGSdissipation and diffusion obtained from the modified equation for the N = 64 test case with minmod limiter. where the first term is the diffusion by SGS modes and the second term is the SGS dissipation. This expansionallows us to investigate the SGS dissipation by terms that are present in the modified equation. We discretized eachderivative (in space and time) using standard 3 point second order central differencing. Note that vorticity does notneed to be discretized, because it is the solution variable. Every term in Eqn. (37) except term VI can be directlycalculated from the instantaneous flow. However, because it is necessary to balance Eqn. (37), we can find term VI,which gives us the exact numerical SGS dissipation and SGS diffusion due to the scheme and slope limiter. Figure 7(a)shows the temporal evolution of each term in Eqn. (37) for the test case with N = 64 employing the minmod limiter.Near the initial time, the temporal change in the enstrophy is dominated by the amplification of vortex stretching.At this time, the vortical structures are large and nearly resolved by the grid resolution and there is very little SGSinteractions. However, as the flow progresses towards the transition to turbulence around tk u = 5, the SGS termbegins to increase as a balance to the amplification by vortex stretching. Around this time, there is non-negligibleviscous dissipation, however, the other terms are small compared to the vortex stretching and SGS terms. This trendcontinues towards the maximum dissipation time near tk u = 9 and after as the turbulence begins to decay. Thetrend that vortex stretching balances the SGS dissipation is consistent with experiments of high Reynolds numberturbulence [67] with multi-probe hot-wire measurements. Overall, the SGS term is balanced by the vortex stretching,which indicates that a majority of the dissipation is provided by the numerical method not the viscosity.In Fig. 7(b) the SGS dissipation and SGS diffusion term (VI) in Eqn. (37) are compared with the modeled SGSdissipation and production (i.e. modified equation sub-grid dissipation and diffusion) that uses R ij = R tij + R gij + R hij from Eqns. (35), (32), and (34), respectively. The temporal evolution of the modified equation SGS dissipation andSGS diffusion has a similar trend as the numerically calculated term while there are some discrepancies near the peakdissipation. These discrepancies arise from a) the simplified calculation of the modified equation by not incorporatingthe exact limiter definition, and b) truncation of the final modified equation to only second-order accurate terms.The modeled SGS dissipation and modeled SGS diffusion are obtained through Eqn. (38), which shows that SGSdissipation is the dominating term compared to the diffusion.The balance of enstrophy is also investigated at a higher grid resolution in Fig. 8(a) with N = 128 and the minmodlimiter. The temporal evolution of the amplification by vortex stretching term is the dominant term in the enstrophybalance. Compared the N = 64, the amount of the viscous dissipation is increased, which is expected because thehigher grid resolution resolves more viscous dissipation. The increased viscous dissipation acts as an offset for the SGSdissipation and SGS diffusion term, which is comparably less for all simulated time compared to the corresponding5 FIG. 8. (a) The temporal evolution of terms in the enstrophy transport equation for the N = 128 test case with minmodlimiter. (b) The temporal evolution of the calculated SGS dissipation and diffusion and the temporal evolution of the SGSdissipation and diffusion obtained from the modified equation for the N = 128 test case with minmod limiter. terms in the N = 64 test case. Figure 8(b) shows the comparison of the calculated SGS dissipation and diffusionterm with the modeled SGS dissipation and SGS diffusion. At this resolution, the model elucidated from the modifiedequation analysis is shown to be well represented because higher order terms in the modified equation become smallerwith higher grid resolution. As with the coarser grid resolution the SGS dissipation is dominant compared to the SGSdiffusion. Overall, the enstrophy transport equation analysis shows that the numerical results and the form of theimplicit SGS model derived from the modified equation analysis are similar and provides validation for the modifiedequation analysis for this scheme.The modified equation SGS dissipation and diffusion can be separated into three different terms based on Eqns.(32), (34) and (35) developed directly from the modified equation. This allows us to identify which mechanismdominants the dissipation in the ILES. Figure 9(a) shows the temporal evolution of the three terms for the N = 64test case with the minmod limiter. Overall the modified equation SGS dissipation is driven by terms similar toa tensor-gradient model. Similarly, for the N = 128 test case with the minmod limiter shown in Fig. 9(b), thetensor-gradient terms in Eqn. (32) are the dominant mechanism for dissipation. In general, these SGS models withtensor-gradient are not overly dissipative and normally are paired with other models. This is consistent with the SGSdissipation shown herein. The modified equation terms similar to a hyper-viscosity model in Eqn. (34) are shown tocontribute to the total SGS dissipation, which enables the total modified equation SGS dissipation to have a similartemporal evolution as the calculated SGS dissipation (term VI). The SGS dissipation contribution from Eqn. (35) isshown to relatively less significant for both resolutions.Figure 10(a) shows the enstrophy balance for the N = 64 with no-limiter case to show the temporal evolution ofthe transport of enstrophy when no limiter is employed in the scheme, which limits the numerical dissipation. Similarto the limited cases presented in Fig. 7, amplification by vortex stretching dominants the evolution of enstrophy,but is significantly larger than the limited cases. Furthermore, terms IV and V are also augmented compared tothe limited case by the lack of limiter and pathway for dissipation. This has a large influence on the balance andthe SGS dissipation/diffusion term. This term is significantly higher than than the limited case and corroboratesthe high numerical dissipation observed in Fig. 5(c). Similarly, the N = 128 case shows that the amplification ofvortex stretching is elevated by the lack of limiting. This affects the enstrophy balance by increasing the remaindercalculated in the SGS dissipation/diffusion.Figure 11 compares the SGS dissipation (cid:15) SGS with the expected SGS dissipation obtained through DNS. Thelatter is obtained by filtering the dissipation spectra Fig. 4 using a sharp-cutoff spectral filter with a length scale2 π/ ( N/ N is the grid resolution for each case. Figures 11(a), (b), and (c) show that the SGS dissipation6 FIG. 9. The temporal evolution enstrophy transport of components of the SGS dissipation and diffusion obtained frommodified equation in Eqns. (32), (34) and (35) for (a) N = 64 and (b) N = 128 test cases with minmod limiter. for the N = 256 , N = 128 and N = 64 cases with the minmod limiter ω i ∂ j R ij is included in Figs.11(b) and (c). The SGS dissipation is relatively similar to that from the filtered DNS for N = 128. It is however,under-predicted for N = 64 . This provides further evidence with the Fig. 5(c) that more analysis and clarificationis needed especially in the low resolution limit.Next, to further understand the implicit numerical dissipation and SGS model derived from the modified equationanalysis we compare several terms of the modified equation to the well-known Smagorinsky model. A Smagorinskymodel [11] was developed in Ref. [55] for the VTE and can used for closure of R ij : ∂R ij ∂x j = g i = − ∂∂x j (cid:18) ν t ∂ω i ∂x j (cid:19) − ∂ν t ∂x j ∂ω j ∂x i . (39)The eddy viscosity ν t is given by ν t = C s ∆ | (cid:101) S | , (40)where (cid:101) S ij is the filtered strain-rate tensor, C s is the Smagorinsky constant and | (cid:101) S | = (2 (cid:101) S ij (cid:101) S ij ) . While the implicitmodel is shown by modified equation analysis to consist of several terms (Eqns. (35), (32), and (34)) and provideimplicit SGS vorticity stress on all tensor element, the comparable Smagorinsky model creates an anti-symmetric SGSvorticity stress tensor.Additional LES cases are simulated using the explicit model in Eqn. (39). No limiter is used in the upwindbased scheme to ensure that the dissipation is enabled through the explicit model; however, there is still interactionbetween the scheme dissipation and explicit model dissipation. Several Smagorinsky constants are employed including C s = 0 . .
3, where the former is used in Ref. [55] with a vortex particle method discretization. Simulationsare performed at N = 64 and 128 resolution of the Taylor-Green vortex at Re = 1600 to compare to the presentresults. Figures 12(a) and (b) show the numerical filtered dissipation (cid:101) (cid:15) ∗ at N = 64 and N = 128, respectively. Thedissipation is shown to be affected significantly by the Smagorinsky constant. In fact, the constant used previousstudies is shown to be out performed by the ILES. The higher constant, especially in N = 64 performs slightly better.7 FIG. 10. The temporal evolution of terms in the enstrophy transport equation for (a) N = 64 and (b) N = 128 test caseswith no-limiter.FIG. 11. The SGS dissipation (cid:15) SGS for the Taylor-Green vortex for (a) N = 256, (b) N = 128, and (c) N = 64. The SGSdissipation obtained from the modified equation analysis is shown for the N = 128 and N = 64. The discrepancies are affected by the interactions of the modified equation of the vorticity transport scheme, evenif no limiter is employed, and the SGS model. Further insights into the performance of explicit LES compared tothe present implicit model is shown in the enstrophy transport analysis in Figs. 12(c) and (d) for the N = 64 and N = 128, respectively. Results from the explicit LES model suggest that the value of the model coefficient produces8 FIG. 12. The filtered dissipation (cid:101) (cid:15) ∗ as a function of time for (a) N = 64 and (b) N = 128 comparing test case with minmodlimiter, two explicit LES cases with C s = 0 . .
3, and DNS results from a spectral method. The temporal evolution of thecalculated SGS dissipation and diffusion obtained from the present method, the modified equation, and Smagorinsky model forexplicit LES with grid discretization of (c) N = 64 and (d) N = 128. high variability for the SGS dissipation. Moreover, a priori tests of forced homogeneous, isotropic turbulence [55]reveal that both vorticity convection and vortex stretching contribute to sub-grid scale dissipation, but the expliciteddy-viscosity model does not capture the contribution from vortex stretch adequately. While the variability maybe reduced with a dynamic model [12, 55], the overall simplification, deterministic dissipation through the modifiedequation, consistent integration both the vorticity convection and vortex stretching, and reduced computational effortof the implicit approach are desirable qualities for a scheme. Furthermore, ILES has been described by one constant—typically the effective viscosity [76]—to quantify the observed behavior of the simulated flow field.Finally, we analyze the scheme by attempting to provide methodology to characterize the flow field in terms of anumerical effective Reynolds number Re f . In the previous analysis, we demonstrated that the numerical dissipationof the present schemes and the dissipation predicted by the modified equation are similar. In what follows, we providemethodology to characterize the flow in order to determine if the simulation is sufficiently resolved to provide anaccurate solution. This is particularly important considering the N = 64 simulations where numerical dissipation isconsiderably different from DNS dissipation. If the simulation is not sufficiently resolved then we provide informationon what the simulation physically solved. In coarse-grained simulations, such as LES either explicit or implicit, wheresmall scales are not directly solved but are modeled, the Re f can be used to interpret what the simulation physicallysolved. The methodology is based on the analysis of an ILES Euler scheme of turbulence flows by Ref. [36]. Forincompressible flows with sufficiently fine grid resolution to be considered a DNS, the numerical dissipation − dK/dt FIG. 13. The effective Re, Re f , as a function of twice the enstrophy throughout the runtime of the simulations for severalgrid resolutions and viscosities. The Re f at (b) t = 5 and (c) maximum dissipation. and resolved dissipation (cid:15) E = 2 ν (cid:104) ω i ω i (cid:105) are equal. This means that the scales impacted by the numerics of scheme arein only in a narrow band of high wavenumbers dominated by the viscous scales. In the under-resolved simulationspresented herein, the ratio of the resolved enstrophy and dissipation is used to characterize the Re f as the following:Re f = u k ν f = (cid:104) ω ∗ i ω ∗ i (cid:105) (cid:15) ∗ . (41)Figure 13(a) shows the temporal evolution of Re f of several test cases including some additional test cases withdifferent Re = u /k ν than presented above in order to understand some general trends. The additional Re includedare 800, 2000, 3000, 5000, and an inviscid case where the Re = ∞ . In each case, the initial Re f is near the Re indicatingthat the initial numerical dissipation and resolved dissipation are equal. However, the Re f decreases as the enstrophyincreases towards a maximum. The grid resolution of each case has an effect on how the Re f decreases where thelargest decreases can be attributed to the lowest grid resolution cases. In Fig. 13(b), the Re f is selected at t = 5 whenthe transition to turbulence is occurring. The Re f is observed to be a function of the physical viscosity while the gridresolution has a secondary effect. However, the Re f lower than the Re, which indicates there is noticeable discrepancyin the numerical dissipation and the resolved dissipation. The Re f is selected at the location of maximum dissipationfor each case and is shown in Fig. 13(c). The results show that the grid resolution has the largest effect on Re f . TheRe f increases with the grid resolution. Regardless of the slope limiter or Re, all cases show similar behavior. Thissuggests that the numerical dissipation is highly affected by the grid resolution compared to the effects of the viscousdiffusion or limiting in the flux functions of the present scheme for the Taylor Green vortex. We observe that the N = 128 cases, which have a closer Re f to 800 than 1600, can be reasonably compared to the results for the Re = 800in some parts of the evolution of the flow. At the highest resolution case N = 256, secondary effects of the Re beginto have a larger effect.Next, we can use the ratio of the numerical dissipation and the resolved dissipation to create an effective lengthscale as follows: ∆ x f = (cid:15) / ( (cid:101) ω i (cid:101) ω i ) / . (42)The effective length scale is compared with the Kolmogorov length scale η provide by the DNS at a Re = 1600 in Fig.14. The temporal evolution of the effective length scale indicates that before the transition to turbulence, the flow0 FIG. 14. The ratio of the effective grid spacing and Kolmorgorov lengthscale (from DNS) ∆ x f /η for the Re = 1600 cases.Red: no-limiter, Black: minmod limiter. field is relatively resolved. However, as the large structures begin to break down, the grid is no longer sufficiently fineto resolve the flow and the effective length scale become larger than the Kolmogorov scale. The grid resolution has alarge effect on the evolution of the effective length scale.Further analysis into Re f of the simulations allows us to interpret the effective viscosity ν f with the Taylormicro-scale λ which is calculated as follows: λ = 13 (cid:88) i (cid:114) (cid:104) u i u i (cid:105) / (cid:104) ∂u i ∂x i ∂u i ∂x i (cid:105) . (43)The effective Taylor micro-scale Re λ = u (cid:48) λ/ν f . Many studies on high-Re turbulence have investigated the relation-ship between the energy containing eddies in flows with the dissipation especially using turbulence-resolving DNS orexperiments [77–80]. From these efforts, we are able to conclude that the time scale of the energy-containing eddies u (cid:48) /L are on the same magnitude as the time scale of the dissipation rate (cid:15)/u (cid:48) , where u (cid:48) is the characteristic velocityscale and L is the characteristic length scale. This key finding is especially beneficial for success with under-resolvedsimulations where by resolving large energy-containing eddies, there is some understanding of the magnitude of thedissipation just based on resolved characteristics. Furthermore, the DNS studies showed that with a sufficiently highRe λ , the flow reaches a minimum state [81] where turbulence seems to exhibit self-similar behavior. The minimumstate of turbulence flows where the time scale of the large scales and dissipation in terms of Re λ can be as lowas Re λ >
100 [82, 83] or Re λ >
200 [77, 80]. For the Taylor-Green vortex case with Re = 1600, the maximumRe λ = 140 [70]. This implies that in order to reach some self-similarity of small-scale turbulence the flow must havesufficient separate of scales, and the dissipation must occurs away from the large scales in the spectral sense. Ref. [36]used this criteria in order to characterize ILES schemes for velocity-pressure formulation to see if a minimum stateis reached. Furthermore, well-known relationships for Re λ to a large-scale Re L exist. Several studies have pointedtowards Re L ∼ Re λ [84, 85].We present details on the characteristics of several of the Taylor-Green simulations in Table I, which shows the effec-tive Re at transition to turbulence ( t = 5), maximum dissipation ( t ≈ t = 15 , u (cid:48) = (cid:104) ( u i u i ) / (cid:105) .The high-resolution cases with N = 256 show that the limiter has minor effects on Re f as also shown in Fig. 13.While only at the transition stage has the flow field reached a sufficiently high Re λ , this is consistent with ILESsimulations with for the Euler equations [36]. The DNS has an initial Re λ = 55 which increases to the maximumRe λ = 140 [70, 76] around t = 9. Re λ ≈
100 at larger times t >
10 as the turbulence decays. In the N = 128 testcases, the numerical dissipation provided by the minmod limiter has more effects while in the N = 64, the Re f remainsrelatively low and does not provide positive results during turbulence decay indicating that the grid resolution is nothigh enough or the dissipation is not enough for this flow.1 TABLE I. Effective Re for ILES of Taylor-Green Vortex N Limiter t Re f = u k ν f ∆ x f Re λ = u (cid:48) λν f Re L ∼ Re λ
256 No-limiter 5.00 774 2.30 × −
98 96929.03 1048 1.60 × −
58 340715.00 1624 1.61 × −
47 225120.00 1704 1.98 × −
43 1871Minmod 5.00 923 2.16 × −
122 149918.30 698 2.13 × −
52 271115.00 1135 2.08 × −
47 222620.00 1283 2.30 × −
45 2048128 No-limiter 5.00 639 2.60 × −
86 75278.49 589 2.31 × −
42 181915.00 1328 1.92 × −
45 210020.00 1596 2.04 × −
41 1725Minmod 5.00 769 2.47 × −
109 119568.29 533 2.67 × −
52 275515.00 856 2.57 × −
39 155920.00 1313 2.39 × −
47 221264 No-limiter 5.00 1356 1.84 × −
203 412968.23 329 3.28 × −
38 1503Minmod 5.00 667 2.80 × −
104 109658.33 303 3.75 × −
46 2144
V.2. Forced Isotropic Turbulence
Next, we consider the numerical experimentation of isotropic turbulence simulation employing a forcing schemedeveloped for linear forced turbulence in physical space [86]. Forced isotropic turbulence provides test cases thatallow for velocity and vorticity statistics to be collected over the simulation, which are used to assess the ILESdissipation and ability to capture turbulence statistics. A forcing term (cid:101) f = ∇ × (cid:15) / u (cid:48) u , where (cid:15) is a constantenergy injection and u (cid:48) = (cid:0) (cid:104) K (cid:105) (cid:1) / is the rms (root-mean-squared) velocity updated throughout the simulation,is included in Eqn. (4) to maintain a nearly constant energy after an initial transient time. The energy injectionconstant is (cid:15) = 5 × − . The simulation cases are initialized using the conditions for isotropic turbulence. It issimulated from an initial condition, which quickly decays into fully-developed turbulence, given in Refs. [87] and [88]where the initial energy spectrum takes the form: E ( k,
0) = 32 A k σ +1 p k σ exp (cid:32) − σ (cid:18) kk p (cid:19) (cid:33) , (44)where k p is the wavenumber at which E ( k,
0) is maximum, σ is parameter, and A is (cid:82) ∞ k σ exp( − σk / dk . A flowis initialized with σ = 4 and k p = 3. The ν = 0 . π , and the initial velocity fluctuations of Re = 6 × . Grid resolutions are employed in the 2 π -sidedperiodic box with the number of grid cells per dimension N = 256 , ,
64, and 32 with a minmod limiter for allresolutions. An additional no-limiter case is employed for the N = 256 resolution based on results above which showthat the N = 256 is nearly resolved for the present initialization.Several studies [31, 36, 39] have employed forced turbulence simulations to examine the turbulence statistics espe-cially to see sensitivities the velocity, vorticity, and strain rates have for ILES with different grid resolutions. It hasbeen observed in Ref. [36] for ILES that above a mixing transition of Re λ >
100 [83], the behavior of the turbulencestatistics in terms of velocity probability density functions (PDFs) begin to approach asymptotic Re statistics. Figure15 shows the longitudinal and transverse fluctuating velocity gradients with a comparison to DNS [89] and grid turbu-lence experimental measurements [90]. The present simulations capture the turbulence statistics, especially, the highprobability statistics well for all grid resolutions. The tails of the PDFs, which represent only a small fraction of thesolution correspond well with higher Re λ DNS and experimental cases, which suggests that the numerical dissipationin the scheme may be appropriate for simulating high Re turbulence cases in which intermittency and non-Gaussian2
FIG. 15. Longitudinal (left) and transverse (right) velocity gradient probability density functions normalized by the rmsvorticity ω (cid:48) with several grid resolutions. Markers indicate PDFs captured from DNS results from Ref. [89] are at Re λ =35 , ,
94 and 168 and Experimental measurements from Ref. [90] are at Re λ = 852. Each PDF is labeled. behavior at the tail of PDFs are observed. The negative bias of the longitudinal derivative of the fluctuating velocityis similar to results present in DNS and experimental flow field statistics. The PDF of the transverse derivative of thefluctuating velocity shown in Fig. 15 demonstrates the convergence of the velocity towards non-Gaussian behaviorobserved for this quantity. The turbulence statistics indicate that there is convergence towards high Re turbulencestatistics.Figure 16(a) show the near Gaussian behavior, similar to flow fields from DNS results, of the velocity fluctua-tions, which show asymptotic convergence as the grid resolution increases. This also suggests that the effective Reis increasing with the grid resolution as seen in the Taylor-Green vortex. The statistics for the magnitude of thevorticity fluctuations are shown in Fig. 16(b), which suggests that vorticity statistics are better represented comparedto the velocity statistics at lower grid resolutions and for this particular Re. Interestingly, this is not seen ILES ofvelocity-pressure formulations [36], which may suggest that the present scheme designed to preserve vorticity mayenhance the ability of lower grid resolution cases to capture some vorticity statistics more efficiently and DNS simula-tions from Ref. [91] have relatively low Re λ . Next, we focus on statistics related to the resolved enstrophy transportequation shown in Eqn. (37). The PDF of the vortex stretching, σ ij = ω i s ij ω j / | ω | , shows the convergence of thehigh probability portion of the PDF and tails that correspond well with different Re. Furthermore, the PDF of thestrain rate magnitude in Fig. 16(d) shows similar behavior. Overall, the velocity and vorticity statistics display thatthe turbulence characteristics have asymptotic convergence and tails of the PDF resemble higher Re statistics withincreasing grid resolution.The effect of numerics of the ILES on the correlation of turbulence statistics is shown in Fig. 17. First, Fig. 17(a)shows that the joint PDF of the vorticity magnitude and the strain rate magnitude | S | . The correlation betweenthe two quantities reveals some convergence towards in the tails of the PDF for the three lowest grid resolutions.However, because the N = 256 cases employs a central difference for the limiter, due to it being nearly resolved, thereis a difference between this case and the other three. Furthermore the correlations between the vorticity and vortexstretching and the strain rate magnitude and stretching in Figs. 17(b) and (c) respectively show similar behavior.The statistics indicate that there is high probability that the grid resolution does not have an effect on the statisticsof the flow field, while the outlying statistics have some dependence on the resolution of the simulation.In order to further assess the implicit method, invariants of the velocity gradient [92], which is related to turbulencedissipation, are analyzed. The three invariants of the velocity gradient A ij = ∂u i ∂x j are given as follows: P = − A ii = 0 R = − A ij A ji Q = − A ij A jk A ki (45)Because the first invariant P is zero for an incompressible flow, the velocity gradient is completely determined by thesecond and third invariants, R and Q , respectively. The invariants determine the local flow topology and indicatestretching and stability.Figure 18(a) and (b) shows the joint PDF for the N = 64 and N = 128 test cases, respectively. The invariants are3 FIG. 16. (a) Velocity magnitude, (b) vorticity magnitude, (c) vortex stretching, and (d) strain rate magnitude probabilitydensity function normalized by the rms vorticity ω (cid:48) of forcing case for several grid resolutions. The case N = 256 with no-limiteris shown as a solid cyan line. Marked indicate DNS results from Ref. [91] are at Re λ = 37 , , ,
142 and 168. Each PDF islabeled.FIG. 17. Turbulence statistics of the forced turbulence simulations for grid resolutions of N = 256 , ,
64 and 32. Contourlines for each plot are shown at { , − , − , − } . FIG. 18. Joint PDF of velocity gradient invariants Q and R normalized by the resolved enstrophy (cid:101) E of the forced turbulencesimulations for grid resolutions of (a) N = 64 and (b) N = 128. Contour lines for each plot are shown in log-space between 0and 10 − . determined at every grid points at t/τ = 10, where τ = u (cid:48) /(cid:15) is the eddy-turnover time. Both cases reveal a joint PDFthis consistent with previous studies of homogeneous turbulence [92]. Furthermore, the invariants calculated in eachsimulation corroborate the PDFs of the turbulence statistics in that the overall statistics are captured independent ofgrid resolution. Increased grid resolution increases the appearance of outliers and are be associated with intermittency,which should be expected by resolving smaller velocity fluctuations.Next, three new simulations are introduced with a grid resolution N = 64 but with different initial Re =8 . × , . × , and 4 . × . This is obtained by changing the kinematic viscosity. All other conditions andforcing remain consistent with the N = 64 test case at Re = 6 × . Figure 19(a) shows the longitudinal velocitygradient PDF of the four N = 64 simulations. As the Re number increases, the tails of PDF trend closer to the highRe λ results obtained from DNS and experimental measurements. Furthermore, Fig. 19(b), showing the transversevelocity gradient PDF, suggests similar, yet non-monotonic [2], convergence as the Re number increases. The higherRe demonstrate that a larger scale separation compared to lower Re is achieved at a minimum state [36] such that theturbulence statistics become asymptotic. The velocity gradient statistics as a function of Re corroborate the statisticsas a function of grid resolution (Fig. 15) to indicate that asymptotic turbulence and dissipation behavior with ILESand VTE scheme can be realized through combination of both high initial Re and sufficient grid resolution. Thisensure that the turbulence will have a sufficient separation of scales.The velocity gradient invariants of the four different Re cases at N = 64 are shown in Fig. 20. These statisticsshow that as the Re number increases, the joint PDF of the invariants resembles the canonical ‘teardrop’ shape moredefinitely. The Q-R plot at higher Re compared to Re = 6000 resemble the N = 128 case shown in Fig. 18(b) moreas the Re increases. This further suggests that once a minimum state is reached, statistics become asymptotic.The balance of terms in enstrophy transport equation introduced in Eqn. (37) is shown in Fig. 21(a) for the N = 64 case in order to assess the modified equation analysis with forced turbulence. Similar to the Taylor-Greenvortex, the amplification via vortex stretching is balanced by SGS dissipation and diffusion. Figure 21(b) showsthat the modified equation SGS dissipation and diffusion are consistent with the calculated profile throughout thesimulation. This is a further indication that the implicit model obtained through the modified equation can be usedto describe the dissipation of the scheme.The accuracy of the modified equation SGS dissipation and diffusion at N = 128 is corroborated in Fig. 22. In thiscase, the enstrophy amplification by vortex stretching is offset by both the SGS dissipation and viscous dissipation.At this resolution, there is significantly more resolved viscous dissipation than the N = 64 case. The calculated SGSdissipation and diffusion is shown to be very similar to the modified equation SGS dissipation and diffusion throughout5 FIG. 19. Longitudinal (left) and transverse (right) velocity gradient probability density functions normalized by the rmsvorticity ω (cid:48) for grid resolutions of N = 64 at several initial Re numbers. Markers indicate DNS results from Ref. [89] are atRe λ = 35 , ,
94 and 168. Experimental measurements from Ref. [90] are at Re λ = 852. Each PDF is labeled.FIG. 20. Joint PDF of velocity gradient invariants Q and R normalized by the resolved enstrophy (cid:101) E of the forced turbulencesimulations for grid resolutions of N = 64 at several Re numbers. Contour lines for each plot are shown in log-space between0 and 10 − . FIG. 21. (a) The temporal evolution of terms in the enstrophy transport equation for the N = 64 test case with minmodlimiter. (b) The temporal evolution of the calculated SGS dissipation and diffusion and the temporal evolution of the SGStransfer and diffusion obtained from the modified equation for the N = 128 test case with minmod limiter. the simulation.Finally, we examine the average energy spectra over the final half of the simulations in Fig. 23. The energyspectra indicate that there is a well-captured low wavenumber regime of all the simulations regardless of the gridresolution. A zoomed-in plot in the low wavenumber region before the inertial range show that both the N = 256and 128 spectra are nearly converged together, while the modes for the N = 64 grid resolution approaches similarenergy contents and the lowest grid resolutions have slightly higher energy contributions to the modes in this region.This suggests that very low grid resolution causes slightly more energy to be present in low wavenumber modes. Theinertial range may not be wide enough to separate the energy-containing scales from the dissipative scales. However,the higher resolutions including N = 64, the resolution becomes asymptotic. These cases exhibit the formation of aninertial range that is consistent with a − / N = 256cases suggest that the grid resolution is approaching the resolution required for a DNS.Additional inviscid simulations with the same forcing and initial conditions are performed to quantify the impactof viscosity and resolution. Figure 24 shows the original forced cases and Re = ∞ cases at grid resolutions of N = 128and 64. A forced turbulence DNS cases from Ref. [86] and explicit LES cases with N = 64 and N = 128 from Ref. [93]with a similar linear forcing in physical space are included for comparison. Note that the forcing scheme used in theDNS and explicit SGS cases are performed with the velocity-pressure formulation using the velocity variable whilethe forcing in physical space employed in the present study uses the vorticity-velocity formulation on the vorticityvariable. This may cause parameters selected in the forcing to behave with slight differences between the formulations.However, the spectra reveal that the low wavenumbers are relatively unaffected by the dissipation regardless of theReynolds number and resolution. The impact of viscosity is seen in the high wavenumber scales. Furthermore, theforced ILES simulations show similar behaviors and trends to the DNS and explicit SGS in the low wavenumberregion. The finite Reynolds number does have a larger impact on the inertial range compared to the DNS. At highRe, there is evidence of asymptotic convergence of flow statistics. V.3. Temporally Evolving Jet
A temporally evolving jet provides a more rigorous test case where both large-scale features (Kelvin-Helmholtzinstability) and small-scale fully-developed turbulence are presence. This is a flow that more closely resembles the7
FIG. 22. (a) The temporal evolution of terms in the enstrophy transport equation for the N = 128 test case with minmodlimiter. (b) The temporal evolution of the calculated SGS dissipation and diffusion and the temporal evolution of the SGStransfer and diffusion obtained from the modified equation for the N = 128 test case with minmod limiter.FIG. 23. Energy spectra of forced turbulence simulations. The thin dotted line is at a − / N = 256 withno-limiter is shown as a solid cyan line. complex flows with dominant coherent structures that the present scheme is designed to investigate.The initial velocity flow field from Ref. [74] is given as follows: u ( x ) = 12 −
12 tanh (cid:20) H θ (cid:18) − | x | H (cid:19)(cid:21) (46)where H is the initial jet thickness and θ is the initial momentum thickness. A H/θ = 35 is selected with aRe = ( U − U ) H/ν = 3200. The initial U = 1 and U = 0. Initial velocity fluctuations are prescribed from thespectrum E ( k ) ∼ k exp (cid:2) − k/k ) (cid:3) . However, for the present scheme in the vorticity-velocity formulation we impose8 FIG. 24. Energy spectra of forced turbulence simulations at Re = 6000 and Re = ∞ . The solid line is N = 128 and thedashed line is N = 64. The diamond, square, and triangle markers are forced turbulence cases 1, 2, and 3, respectively, fromRef. [86]. Circle and plus markers denote N = 64 and N = 128, respectively, and correspond to explicit LES of case 2 fromRef. [93]. the initial condition based on the transverse vorticity as follows: ω ( x ) = 14 θ | x | x cosh − (cid:20) H θ (cid:18) − | x | H (cid:19)(cid:21) . (47)The initial condition has been shown to be a good estimate of the inlet velocity profile of spatially evolving jets.However, the temporally evolving jet provides a computationally expedient solution compared to simulating thespatial jet. The temporally evolving jet employs periodic boundary conditions on all sides of a cubic domain of( L × L × L ) = (4 H × H × H ). Several grid resolution are used for comparison of the effects of the discretization: N = 64 , , and 256, where N is the number of grid cells per dimension. At N = 256, the grid resolution issufficiently fine to approach a DNS solution similar to the DNS simulation in Ref. [94]. Both a dissipative ENOslope limiter and no-limiter are employed. The ENO slope limiter is similar to the minmod slope limiter but is nottotal variation diminishing. Contours of the out-of-plane vorticity are shown for each simulation in Fig. 25 at thetime when the flow field becomes starts to become self-similar, t/T ref = 15, where T ref = H/ ∆ U , (see Fig. 26(a)).This time is equivalent to x/H = 7 . N = 64 and 128).Figure 26(a) shows profiles of the mean streamwise velocity (cid:104) u (cid:105) in x direction averaged over the x and x directions for each grid resolution and both the ENO limiter and no-limiter. The profiles are normalized by thecenterline velocity U c = (cid:104) u (cid:105) ( x = 0), while the jet half-width δ / , which is calculated with the mean streamwisevelocity (cid:104) u (cid:105) ( x = δ / ) − U ∞ = ( U c − U ∞ ) and U ∞ is the streamwise velocity far from the jet, is used to normalizethe abscissa. The profiles are chosen in the regime where the jet wake becomes self-similar. The grid resolution hasa large effect on the amount of time the temporally evolving jet requires to become self-similar and will be discussedbelow. The streamwise velocity profiles are compared with several experimental measurements of spatial evolving jetin the self-similar regime. The present simulation results compare well with experimental measurements and showconsistency between the different grid resolutions. Moreover, the spread in the results are comparable to DNS oftemporal jets performed in Refs. [74, 94]. This suggests that even the low grid resolution of N = 64 with the ENOlimiter can capture some of the large-scale features well. The mean transverse vorticity profiles are shown in Fig.9 FIG. 25. Velocity ω normalized by the jet half-width δ / and centerline velocity for the temporally evolving jet contours forthe temporally evolving jet at t/T ref = 15. N = 64 with the no-limiter case. However, the N = 64 with dissipative ENOlimiter, where the shear layer is only discretized by less than four grid cells, is captured reasonably. On the otherhand, the N = 64 and 128 with no-limiter show slightly higher mean vorticity in the shear layer suggesting thatunder-resolved simulations need additional dissipation.The second-order velocity statistics are compared with experimental measurements in Fig. 27. The rms (root-mean-square) streamwise velocity (cid:104) u (cid:105) / profiles, shown in Fig. 27(a), and rms normal velocity (cid:104) u (cid:105) / profiles,shown in Fig. 27(b), are compared with experimental measurements in the self-similar regime. While there are somedifferences between the different test cases, the overall comparison is reasonably well, especially for the lowest gridresolution case with the dissipative ENO limiter, which approximates sufficient dissipation. On the other hand, thelowest resolution with the no-limiter significantly under predicts the rms velocities. While scatter in the data canalso be attributed to the relatively low number of samples for the temporally evolving jet especially at the lowest gridresolution, the role of numerical dissipation to capture flow statistics is apparent in the jet.The grid resolution and numerical dissipation affect the temporal evolution of the jet from its initial conditionthrough the transition to turbulence. Up to this point, the result from the simulations have focused on the flowfield after the jet transitions to turbulence, the transient behavior diminishes, and a self-similar regime occurs. Intemporal jet simulations [94], the self-similar regime is obtained at times t/T ref >
20 for a similar initialization. Thisis equivalent to x/H = 10 for a spatial-evolving jet. Figure 28 shows the evolution of the mean transverse vorticityprofiles with time for each grid resolution case with the dissipative ENO limiter. The initial mean vorticity profile0
FIG. 26. Profiles of the (a) mean streamwise velocity (cid:104) u (cid:105) normalized by the centerline velocity U c compared with experimentalresults from Ref. [95] (circles) and Ref. [96] (squares) and (b) mean vorticity (cid:104) ω (cid:105) normalized by the jet half-width δ / andcenterline velocity for the temporally evolving jet. Blue - no-limiter, Black - ENO limiter.FIG. 27. Profiles of the (a) rms (root-mean-square) streamwise velocity (cid:104) u (cid:105) / and (b) rms normal velocity (cid:104) u (cid:105) / normalizedby the centerline velocity U c compared with experimental results from Ref. [95] (circles) and Ref. [96] (squares). Blue - no-limiter, Black - ENO limiter. across the shear layer is relatively high and sharp. At further instances in time, the vorticity magnitude diminishesto a self-similar solution. Each grid resolution test case due to the increase numerical dissipation and grid cell size,transitions to the self-similar solution at different times. For the most resolve case, the self-similar solution is obtainedat t/T ref >
20, while the lowest grid resolution transitions quicker, around t/T ref >
10. We are effectively solvingthe flow field at different effective Reynolds numbers through the transition to turbulence but are able to obtain aself-similar solution after the transition. The temporal impact of the numerical scheme on the flow field and theeffective Reynolds number is significant when solving temporally evolving flows.An estimate of the effective Reynolds number is determined from measurable features of the flow field. It isimportant to determine the effective Re to determine what the simulation is solving and if it has reached a minimumstate in the sense of the asymptotic turbulence statistics. Here, we assume a high Reynolds number is achieved andusing the asymptotic relationship for high Re regime of isotropic turbulence [77] where D = (cid:15)L/U ≈ similar toILES estimates for complex Richtmyer-Meshkov instabilities [36]. The velocity scale U = u (cid:48) = (cid:104) ( u i u i ) / (cid:105) and lengthscale L is chosen to the jet half-width, which is shown as a function of the time in Fig. 29(a). The temporal evolution1 FIG. 28. Temporal evolution of the vorticity profiles (cid:104) ω (cid:105) for (a) N = 256, (b) N = 128 and (c) N = 64. Each case employsthe ENO limiter.FIG. 29. The temporal evolution of the (a) jet half-width, (b) effective viscosity, and (c) effective Re. Blue - no-limiter, Black- ENO limiter. of the half-width shows that the grid resolution affects the jet spreading, which is expected based Fig. 28. Thedissipation is obtained (cid:15) = DU /δ / and the effective viscosity is obtain through ν f = (cid:15)/ω i ω i . Figure 29(b) showsthe temporal evolution of the effective viscosity. For each grid resolution and slope limiter case, the initial effectiveviscosity behaves slightly differently, however, as the jet flow field asymptotically approaches the self-similar solution,the effective viscosity for each simulation approaches a single viscosity, which is the viscosity of the simulation. Thenumerical dissipation in the lowest resolution case is able to reach a state such that the asymptotic relationships in thejet can be reached. The Re f = u (cid:48) δ / /ν f is shown in Fig. 29(c). It indicates that the lower grid resolution cases withthe increased dissipation are in fact solving a slightly different temporally evolving jet problem. However, becausethe flow comes self-similar at long times, this is not readily noticeable from profiles. Moreover, experimentation, bothnumerical and measurement, have shown that at high enough Re based on the inlet velocity and slot width (6000and Ref. [97] and 1000 in Ref. [98], respectively) the spatially evolving planar jet becomes independent of Re, whichcomparable to the same Re for the present simulations ( Re ≈ N = 256 ,
128 and 64 for both the limited and no-limitercases at t/T ref = 40. The spectra for the nearly resolved N = 256 cases are comparatively similar. Furthermore,the N = 128 case with the ENO limiter matches reasonably well in the large-scales but is slightly more dissipative2 FIG. 30. Energy spectra E ( k ) of simulations temporal simulations: (a) Re = 3200 at several discretizations, Blue - no-limiter,Black - ENO limiter, (b) N = 128 case at different Re, and (c) N = 64 case at different Re. than the N = 256 cases in the high wavenumber range as expected. The lack of a limiter in the N = 128 case showsthat a non-dissipative scheme can have a large effect on the low wavenumbers. The lowest resolved cases, N = 64show some differences in the low wavenumbers. Additional Reynolds number cases are performed to assess the lowwave number behavior and corroborate the results above that indicate the lower grid resolution are solving a slightlydifferent temporally evolving jet problem. Figure 30(b) shows the N = 128 resolution for Re = 3200 , × , × , and ∞ . The higher Re numbers cases have an asymptotic behavior in the low wavenumbers which converge to thenearly resolved N = 256 resolution. Additionally, this behavior is also found in the N = 64 cases with high Reshown in Fig. 30(c). The viscosity Re = 3200 plays a significant role in dissipation but can be mitigate by carefullyselecting the grid resolution. The spectra show that with adequate resolution and a high enough Re, the solutionconverges to correct results especially in the low wavenumber range, which is most important for performing LES. Theestimation of an effective Re becomes especially critical when attempting ILES of more complex flow fields where theflow characteristics can become independent of Re such as wind turbine wakes [59]. By selecting the grid resolutionfor the problem, an ILES with the present scheme can be employed such that the large-scale structures are capturedand small-scale turbulence are implicitly modeled to obtain physically accurate, expedient results. VI. CONCLUSIONS
The subgrid scale characteristics and effectiveness of an upwind finite volume scheme for the vorticity transportequations were investigated. The numerical scheme employs a generalized Riemann problem-based multi-dimensionalwave propagation approach. Modified equation analysis was used to characterize the dissipation and backscatter. Theanalysis revealed two limits for including dissipation implicitly through numerics: 1.) a low dissipation limit using asecond-order central difference, most appropriate in smooth areas, i.e. regions dominated by large vortical coherentstructures; and 2.) a high dissipation limit using a first order upwind difference, used when the vorticity changesrapidly across grid cells, i.e., regions of under-resolved turbulence. While the former is ideal in well-resolved areas of3the flow field, the latter is necessary in regions in which dissipation is essential to account for the transfer of energyfrom the resolved scales to the sub-grid scales in the absence of an explicit SGS model. The modified equation at thehigh dissipation limit contains many terms, some of which can be combined into forms that are similar to commonlyused explicit sub-grid scale models including the tensor-gradient models, hyper-viscosity model, and a simple gradientmodel. This serves a qualitative tool to understand the implications of using the present scheme for ILES. We arecareful to note that—for ILES in general—the terms in the modified equation do not have to be similar to explicitSGS terms. Rather, the similarities to known models provide insight in characterizing the scheme.To characterize turbulence with the present scheme for the ILES methodology, a series of turbulence-in-a-boxsimulations of the Taylor-Green vortex was performed to understand the process of energy transfer from energy-containing scales to the sub-grid scale. The Taylor-Green vortex cases revealed that the grid resolution must becarefully taken into account in order to obtain desired results. The grid resolution is the largest factor for themagnitude of the effective Re while the choice of limiter, which subtly controls the numerical dissipation, has animpact on the accuracy. However, in under-resolved simulations, a dissipative limiter is essential. The dissipationterms obtained from the modified equation analysis are shown to faithfully represent the implicit SGS torque in theTaylor-Green vortex case. In high Reynolds number flows, where there is a marked separation of scales, the methodis able to represent the high energy modes. Further numerical experiments with forced turbulence revealed thathigh-Reynolds number asymptotic turbulence statistics can be reasonably captured with the ILES methodology forthis vorticity-velocity formulation scheme.Finally, simulations of a temporally evolving jet, which contains both large-scale vortical structures and fine-scaleturbulence show that under-resolved numerics can capture asymptotic turbulence statistics and large-scale features.The method is particularly useful when the effective Reynolds number is past a threshold beyond which the flow isdependent on the Reynolds number.The simulations studied herein represent canonical flows which allow us to build our understanding of the methodto more complex flows. The simulation tests show that coarse grid resolutions provide a good estimate for largeenergy-containing modes given a large enough inertial range. This particular vorticity-velocity scheme was designedto capture and preserve large vorticity structures in flows in which fully-developed small-scale turbulence tends tolocalized and large energy-containing structures dominate the flow. Our previous work [45] showed how the schemecan be used to capture vortical structures, while this work indicates that the impact of fine-scale turbulence on theenergy-containing scales can also be reasonably represented.
ACKNOWLEDGMENTS
This work is supported through a subcontract from Continuum Dynamics, Inc. under Navy STTR Phase II contractN68335-17-C-0158 entitled “Advanced Wake Turbulence Modeling for Naval CFD Applications”, with guidance fromDr. Glen Whitehouse at Continuum Dynamics, Inc., and Drs. David Findlay and James Forsythe at the Naval AirSystems Command.This paper was prepared as an account of work sponsored by an agency of the United States (U.S.) Government.Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, express orimplied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information,apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Referenceherein to any specific commercial product, process, or service by trademark, manufacturer, or otherwise does notnecessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agencythereof. The views and opinion expressed herein do not necessarily state or reflect those of the U.S. Government orany agency thereof. [1] AJ Yule, “Large-scale structure in the mixing layer of a round jet,” Journal of Fluid Mechanics , 413–432 (1978).[2] Knut Akselvoll and Parviz Moin, “Large-eddy simulation of turbulent confined coannular jets,” Journal of Fluid Mechanics , 387–411 (1996).[3] Daniel Foti, Xiaolei Yang, Michele Guala, and Fotis Sotiropoulos, “Wake meandering statistics of a model wind turbine:Insights gained by large eddy simulations,” Physical Review Fluids , 044407 (2016).[4] Daniel Foti, Xiaolei Yang, Filippo Campagnolo, David Maniaci, and Fotis Sotiropoulos, “Wake meandering of a modelwind turbine operating in two different regimes,” Physical Review Fluids , 054607 (2018).[5] Michele Guala, M Metzger, and Beverley J McKeon, “Interactions within the turbulent boundary layer at high reynoldsnumber,” Journal of Fluid Mechanics , 573–604 (2011). [6] Jiannong Fang and Fernando Port´e-Agel, “Large-eddy simulation of very-large-scale motions in the neutrally stratifiedatmospheric boundary layer,” Boundary-Layer Meteorology , 397–416 (2015).[7] AG Kravchenko and Parviz Moin, “On the effect of numerical errors in large eddy simulations of turbulent flows,” Journalof Computational Physics , 310–322 (1997).[8] Fotini Katopodes Chow and Parviz Moin, “A further study of numerical errors in large-eddy simulations,” Journal ofComputational Physics , 366–380 (2003).[9] Martin Freitag and Markus Klein, “An improved method to assess the quality of large eddy simulations in the context ofimplicit filtering,” Journal of Turbulence , N40 (2006).[10] Siavash Toosi and Johan Larsson, “Towards systematic grid selection in LES: Identifying the optimal spatial resolution byminimizing the solution sensitivity,” Computers & Fluids , 104488 (2020).[11] Joseph Smagorinsky, “General circulation experiments with the primitive equations: I. The basic experiment,” MonthlyWeather Review , 99–164 (1963).[12] Massimo Germano, Ugo Piomelli, Parviz Moin, and William H Cabot, “A dynamic subgrid-scale eddy viscosity model,”Physics of Fluids A: Fluid Dynamics , 1760–1765 (1991).[13] Jorge Bardina, J Ferziger, and WC Reynolds, “Improved subgrid-scale models for large-eddy simulation,” in (1980) p. 1357.[14] Douglas K Lilly, “A proposed modification of the Germano subgrid-scale closure method,” Physics of Fluids A: FluidDynamics , 633–635 (1992).[15] Robert A Clark, Joel H Ferziger, and William Craig Reynolds, “Evaluation of subgrid-scale models using an accuratelysimulated turbulent flow,” Journal of Fluid Mechanics , 1–16 (1979).[16] Elie Bou-Zeid, Charles Meneveau, and Marc Parlange, “A scale-dependent Lagrangian dynamic model for large eddysimulation of complex turbulent flows,” Physics of Fluids , 025105 (2005).[17] Changping Yu, Zuoli Xiao, and Xinliang Li, “Scale-adaptive subgrid-scale modelling for large-eddy simulation of turbulentflows,” Physics of Fluids , 035101 (2017).[18] J-B Chapelier, Bono Wasistho, and Carlo Scalo, “A coherent vorticity preserving eddy-viscosity correction for large-eddysimulation,” Journal of Computational Physics , 164–182 (2018).[19] Kai Schneider, Marie Farge, Giulio Pellegrino, and Michael M Rogers, “Coherent vortex simulation of three-dimensionalturbulent mixing layers using orthogonal wavelets,” Journal of Fluid Mechanics , 39 (2005).[20] Hilde Ouvrard, Bruno Koobus, Alain Dervieux, and Maria Vittoria Salvetti, “Classical and variational multiscale LES ofthe flow around a circular cylinder on unstructured grids,” Computers & Fluids , 1083–1094 (2010).[21] Eric J Parish and Karthik Duraisamy, “A dynamic subgrid scale model for large eddy simulations based on the Mori–Zwanzig formalism,” Journal of Computational Physics , 154–175 (2017).[22] Eric J Parish and Karthik Duraisamy, “Non-Markovian closure models for large eddy simulations using the Mori-Zwanzigformalism,” Physical Review Fluids , 014604 (2017).[23] Bernard J Geurts, “Inverse modeling for large-eddy simulation,” Physics of Fluids , 3585–3587 (1997).[24] Zhideng Zhou, Shizhao Wang, and Guodong Jin, “A structural subgrid-scale model for relative dispersion in large-eddysimulation of isotropic turbulent flows by coupling kinematic simulation with approximate deconvolution method,” Physicsof Fluids , 105110 (2018).[25] Robert H Kraichnan, “Diffusion by a random velocity field,” Physics of Fluids , 22–31 (1970).[26] JP Boris, FF Grinstein, ES Oran, and RL Kolbe, “New insights into large eddy simulation,” Fluid Dynamics Research , 199–228 (1992).[27] David L Youngs, “Three-dimensional numerical simulation of turbulent mixing by Rayleigh–Taylor instability,” Physics ofFluids A: Fluid Dynamics , 1312–1320 (1991).[28] Ralf Jens Gathmann, Mohammed Si-Ameur, and Fabrice Mathey, “Numerical simulations of three-dimensional naturaltransition in the compressible confined shear layer,” Physics of Fluids A: Fluid Dynamics , 2946–2968 (1993).[29] Piotr K Smolarkiewicz and Len G Margolin, “MPDATA: A finite-difference solver for geophysical flows,” Journal ofComputational Physics , 459–480 (1998).[30] Ben Thornber, Andrew Mosedale, and Dimitris Drikakis, “On the implicit large eddy simulations of homogeneous decayingturbulence,” Journal of Computational Physics , 1902–1929 (2007).[31] Adam J Wachtor, Fernando F Grinstein, Carl R DeVore, JR Ristorcelli, and Len G Margolin, “Implicit large-eddysimulation of passive scalar mixing in statistically stationary isotropic turbulence,” Physics of Fluids , 025101 (2013).[32] Len G Margolin and William J Rider, “A rationale for implicit turbulence modelling,” International Journal for NumericalMethods in Fluids , 821–841 (2002).[33] LG Margolin and WJ Rider, “The design and construction of implicit LES models,” International Journal for NumericalMethods in Fluids , 1173–1179 (2005).[34] LG Margolin, WJ Rider, and FF Grinstein, “Modeling turbulent flow with implicit LES,” Journal of Turbulence , N15(2006).[35] Dimitris Drikakis, Marco Hahn, Andrew Mosedale, and Ben Thornber, “Large eddy simulation using high-resolution andhigh-order methods,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences , 2985–2997 (2009).[36] Ye Zhou, Fernando F Grinstein, Adam J Wachtor, and Brian M Haines, “Estimating the effective Reynolds number inimplicit large-eddy simulation,” Physical Review E , 013303 (2014).[37] Stefan Hickel, Nikolaus A Adams, and J Andrzej Domaradzki, “An adaptive local deconvolution method for implicitLES,” Journal of Computational Physics , 413–436 (2006). [38] Len G Margolin, Piotr K Smolarkiewicz, and Zbigniew Sorbjan, “Large-eddy simulations of convective boundary layersusing nonoscillatory differencing,” Physica D: Nonlinear Phenomena , 390–397 (1999).[39] Christer Fureby and Fernando F Grinstein, “Large eddy simulation of high-Reynolds-number free and wall-bounded flows,”Journal of Computational Physics , 68–97 (2002).[40] Fernando F Grinstein and Christer Fureby, “On flux-limiting-based implicit large eddy simulation,” Journal of FluidsEngineering , 1483–1492 (2007).[41] AKM Fazle Hussain, “Coherent structures and turbulence,” Journal of Fluid Mechanics , 303–356 (1986).[42] Dimitri Papamoschou and Sanjiva K Lele, “Vortex-induced disturbance field in a compressible shear layer,” Physics ofFluids A: Fluid Dynamics , 1412–1419 (1993).[43] Karthik Duraisamy and Sanjiva K Lele, “Evolution of isolated turbulent trailing vortices,” Physics of Fluids , 035102(2008).[44] Petros Koumoutsakos and A Leonard, “High-resolution simulations of the flow around an impulsively started cylinderusing vortex methods,” Journal of Fluid Mechanics , 1–38 (1995).[45] Daniel Foti and Karthik Duraisamy, “Multi-dimensional finite volume scheme for the vorticity transport equations,” Com-puters & Fluids , 17–32 (2018).[46] Richard E Brown and Andrew J Line, “Efficient high-resolution wake modeling using the vorticity transport equation,”AIAA Journal , 1434–1443 (2005).[47] Glen R Whitehouse and Alexander H Boschitsch, “Innovative grid-based vorticity–velocity solver for analysis of vorticity-dominated flows,” AIAA Journal , 1655–1669 (2014).[48] Petros Koumoutsakos, “Multiscale flow simulations using particles,” Annu. Rev. Fluid Mech. , 457–487 (2005).[49] Gr´egoire Winckelmans, TS Lund, Daniele Carati, and AA Wray, “A priori testing of subgrid-scale models for the velocity-pressure and vorticity-velocity formulations,” in Proceedings of the CTR Summer Program (1996) p. 309.[50] Georges-Henri Cottet, “A particle-grid superposition method for the Navier-Stokes equations,” Journal of ComputationalPhysics , 301–318 (1990).[51] Mohamed Lemine Ould-Salihi, Georges-Henri Cottet, and Mohammed El Hamraoui, “Blending finite-difference and vortexmethods for incompressible flow computations,” SIAM Journal on Scientific Computing , 1655–1674 (2001).[52] Steven C Rennich and Sanjiva K Lele, “Numerical method for incompressible vortical flows with two unbounded directions,”Journal of computational physics , 101–129 (1997).[53] J Zhao and C He, “A hybrid solver with combined CFD and viscous vortex particle method,” in (2011).[54] Christopher Stone, Earl Duque, Christopher Hennes, and Adrin Gharakhani, “Rotor wake modeling with a coupledEulerian and vortex particle method,” in (2010) pp. 2010–0312.[55] John R Mansfield, Omar M Knio, and Charles Meneveau, “A dynamic LES scheme for the vorticity transport equation:Formulation and a priori tests,” Journal of Computational Physics , 693–730 (1998).[56] Hubert L Meitz and Hermann F Fasel, “A compact-difference scheme for the Navier–Stokes equations in vorticity–velocityformulation,” Journal of Computational Physics , 371–403 (2000).[57] Christopher Davies and Peter W Carpenter, “A novel velocity–vorticity formulation of the Navier–Stokes equations withapplications to boundary layer disturbance evolution,” Journal of Computational Physics , 119–165 (2001).[58] Eric Parish, Karthik Duraisamy, and Praveen Chandrashekar, “Generalized Riemann problem-based upwind scheme forthe vorticity transport equations,” Computers & Fluids , 10–18 (2016).[59] Leonardo P Chamorro, REA Arndt, and F Sotiropoulos, “Reynolds number dependence of turbulence statistics in thewake of wind turbines,” Wind Energy , 733–742 (2012).[60] SK Godunov, “Symmetric form of the equations of magnetohydrodynamics,” Numerical Methods for Mechanics of Con-tinuum Medium , 26–34 (1972).[61] Kenneth G Powell, “An approximate Riemann solver for magnetohydrodynamics,” in Upwind and High-Resolution Schemes (Springer, 1997) pp. 570–583.[62] Derek S Bale, Randall J LeVeque, Sorin Mitran, and James A Rossmanith, “A wave propagation method for conservationlaws and balance laws with spatially varying flux functions,” SIAM Journal on Scientific Computing , 955–978 (2003).[63] Randall J LeVeque, “High-resolution conservative algorithms for advection in incompressible flow,” SIAM Journal onNumerical Analysis , 627–665 (1996).[64] Cyril W Hirt, “Heuristic stability theory for finite-difference equations,” Journal of Computational Physics , 339–355(1968).[65] Vadim Borue and Steven A Orszag, “Local energy flux and subgrid-scale statistics in three-dimensional turbulence,”Journal of Fluid Mechanics , 1–31 (1998).[66] Bert Vreman, Bernard Geurts, and Hans Kuerten, “Large-eddy simulation of the temporal mixing layer using the Clarkmodel,” Theoretical and Computational Fluid Dynamics , 309–324 (1996).[67] Stefano Cerutti, Charles Meneveau, and Omar M Knio, “Spectral and hyper eddy viscosity in high-Reynolds-numberturbulence,” Journal of Fluid Mechanics , 307–338 (2000).[68] AG Lamorgese, DA Caughey, and SB Pope, “Direct numerical simulation of homogeneous turbulence with hyperviscosity,”Physics of Fluids , 015106 (2005).[69] Uriel Frisch, Turbulence: the legacy of AN Kolmogorov (Cambridge University Press, 1995).[70] ME Brachet, “Direct simulation of three-dimensional turbulence in the Taylor–Green vortex,” Fluid Dynamics Research , 1–8 (1991).[71] Marc E Brachet, Daniel I Meiron, Steven A Orszag, BG Nickel, Rudolf H Morf, and Uriel Frisch, “Small-scale structure of the Taylor–Green vortex,” Journal of Fluid Mechanics , 411–452 (1983).[72] Dimitris Drikakis, Christer Fureby, Fernando F Grinstein, and David Youngs, “Simulation of transition and turbulencedecay in the Taylor–Green vortex,” Journal of Turbulence , N20 (2007).[73] FF Grinstein, AA Gowardhan, and AJ Wachtor, “Simulations of Richtmyer–Meshkov instabilities in planar shock-tubeexperiments,” Physics of Fluids , 034106 (2011).[74] Carlos B da Silva and Jos´e CF Pereira, “The effect of subgrid-scale models on the vortices computed from large-eddysimulations,” Physics of Fluids , 4506–4534 (2004).[75] Ye Zhou, Michael Groom, and Ben Thornber, “Dependence of enstrophy transport and mixed mass on dimensionalityand initial conditions in the richtmyer–meshkov instability induced flows,” Journal of Fluids Engineering (2020).[76] Ye Zhou and Ben Thornber, “A comparison of three approaches to compute the effective Reynolds number of the implicitlarge-eddy simulations,” Journal of Fluids Engineering (2016).[77] Yukio Kaneda, Takashi Ishihara, Mitsuo Yokokawa, Ken’ichi Itakura, and Atsuya Uno, “Energy dissipation rate andenergy spectrum in high resolution direct numerical simulations of turbulence in a periodic box,” Physics of Fluids ,L21–L24 (2003).[78] Gregory Falkovich, “Bottleneck phenomenon in developed turbulence,” Physics of Fluids , 1411–1414 (1994).[79] Lian-Ping Wang, Shiyi Chen, James G Brasseur, and John C Wyngaard, “Examination of hypotheses in the Kolmogorovrefined turbulence theory through high-resolution simulations. Part 1. Velocity field,” Journal of Fluid Mechanics ,113–156 (1996).[80] Katepalli R Sreenivasan, “An update on the energy dissipation rate in isotropic turbulence,” Physics of Fluids , 528–529(1998).[81] Ye Zhou, “Unification and extension of the similarity scaling criteria and mixing transition for studying astrophysics usinghigh energy density laboratory experiments or numerical simulations,” Physics of Plasmas , 082701 (2007).[82] Katepalli R Sreenivasan, “On the scaling of the turbulence energy dissipation rate,” Physics of Fluids , 1048–1051 (1984).[83] Paul E Dimotakis, “The mixing transition in turbulent flows,” Journal of Fluid Mechanics , 69–98 (2000).[84] Hendrik Tennekes and John Leask Lumley, A first course in turbulence (MIT press, 1972).[85] Pierre Sagaut and Claude Cambon,
Homogeneous turbulence dynamics , Vol. 10 (Springer, 2008).[86] Carlos Rosales and Charles Meneveau, “Linear forcing in numerical simulations of isotropic turbulence: Physical spaceimplementations and convergence properties,” Physics of Fluids , 095106 (2005).[87] NN Mansour and AA Wray, “Decay of isotropic turbulence at low Reynolds number,” Physics of Fluids , 808–814 (1994).[88] RS Rogallo, “Numerical experiments in homogeneous turbulence,” NASA STI/Recon Technical Report N (1981).[89] Javier Jim´enez, Alan A Wray, Philip G Saffman, and Robert S Rogallo, “The structure of intense vorticity in isotropicturbulence,” Journal of Fluid Mechanics , 65–90 (1993).[90] B Castaing, Y Gagne, and EJ Hopfinger, “Velocity probability density functions of high Reynolds number turbulence,”Physica D: Nonlinear Phenomena , 177–200 (1990).[91] Javier Jimenez and Alan A Wray, “On the characteristics of vortex filaments in isotropic turbulence,” Journal of FluidMechanics , 255–285 (1998).[92] Brian J Cantwell, “On the behavior of velocity gradient tensor invariants in direct numerical simulations of turbulence,”Physics of Fluids A: Fluid Dynamics , 2008–2013 (1993).[93] Ugo Piomelli, Amirreza Rouhi, and Bernard J Geurts, “A grid-independent length scale for large-eddy simulations,”Journal of fluid mechanics , 499 (2015).[94] Carlos B da Silva and Jos´e CF Pereira, “Invariants of the velocity-gradient, rate-of-strain, and rate-of-rotation tensorsacross the turbulent/nonturbulent interface in jets,” Physics of Fluids , 055101 (2008).[95] E Gutmark and I Wygnanski, “The planar turbulent jet,” Journal of Fluid Mechanics , 465–495 (1976).[96] BR Ramaprian and MS Chandrasekhara, “LDA measurements in plane turbulent jets,” Journal of Fluids Engineering ,264–271 (1985).[97] M Klein, A Sadiki, and J Janicka, “Investigation of the influence of the Reynolds number on a plane jet using directnumerical simulation,” International Journal of Heat and Fluid Flow , 785–794 (2003).[98] Ravinesh C Deo, Jianchun Mi, and Graham J Nathan, “The influence of Reynolds number on a plane jet,” Physics ofFluids20