A stable and scale-aware dynamic modeling framework for subgrid-scale parameterizations of two-dimensional turbulence

Accepted Manuscript

A stable and scale-aware dynamic modeling framework for subgrid-scale parameterizations of two-dimensional turbulence Romit Maulik, Omer San PII: DOI: Reference:

S0045-7930(16)30371-1 10.1016/j.compfluid.2016.11.015 CAF 3336

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

21 July 2016 12 October 2016 25 November 2016

Please cite this article as: Romit Maulik, Omer San, A stable and scale-aware dynamic modeling framework for subgrid-scale parameterizations of two-dimensional turbulence, Computers and Fluids (2016), doi: 10.1016/j.compfluid.2016.11.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Scale-aware dynamic eddy viscosity parametrization models are proposed for 2D turbulence. • The proposed hybrid approaches minimize error between functional and structural models. • These modular approaches include the Smagorinsky, Leith, Bardina, Layton and AD models.

CR IP T

• The robustness of these models has been tested by considering two different numerical schemes.

AC

CE

PT

ED

M

AN US

• It is shown that the proposed models could be viable tools for dynamic subgrid-scale modeling.

1

ACCEPTED MANUSCRIPT

A stable and scale-aware dynamic modeling framework for subgrid-scale parameterizations of two-dimensional turbulence Romit Maulik, Omer San∗

CR IP T

School of Mechanical and Aerospace Engineering Oklahoma State University, Stillwater, Oklahoma 74078, USA

Abstract

The physical mechanism of two-dimensional turbulence poses challenges for modeling subgrid-scale physics in largescale geophysical flows. To this end, we put forth a modular dynamic modeling approach for subgrid-scale parame-

AN US

terizations of two-dimensional turbulence. For developing a unifying dynamic modeling approach, a set of coupled subgrid-scale models are proposed by minimizing the error between the functional and structural models. This approach is fundamentally different to the test-filtering based dynamical approach, which is also included in our analysis. Our efforts include two different functional nonlinear eddy viscosity kernels: (i) Smagorinsky’s strain-rate and (ii) Leith’s vorticity-gradient based formulations. The mixing length scales associated with these eddy viscosity kernels are dynamically computed from the local flow physics by incorporating the structural models based upon the scale

M

similarity and approximate deconvolution approaches. A set of decaying turbulence experiments up to Re = 128,000 are compared against direct numerical simulations (DNSs) obtained by a resolution of 20482 . First, it is shown that

ED

less dissipative results are generally obtained using Leith’s eddy viscosity kernel due to its more scale-selective behavior. Among the proposed hybrid models, it was seen that the dynamic Bardina approach yields the least dissipative results, followed by the dynamic Layton, and the dynamic AD models. Due to its more dissipative character, the

PT

dynamic AD model seems to be an efficient approach for very large eddy simulations on coarse grid descriptions. To elucidate the effects of numerics on the subgrid-scale physics, two different high-order discretization schemes are

CE

considered, namely the fourth-order Pad´e and Arakawa schemes. Based on numerical assessments we conclude that the choice of underlying numerical discretization plays a more important role than that of the subgrid modeling in obtaining an energy spectrum that closely approximates the DNS data.

AC

Keywords: Two-dimensional turbulence; Large eddy simulations; Geophysical flows; Dynamic modeling; Smagorinsky model; Leith model

∗ Corresponding

author. Email addresses: [email protected] (Romit Maulik), [email protected] (Omer San)

Preprint submitted to Computers & Fluids

November 26, 2016

ACCEPTED MANUSCRIPT

1

1. Introduction Understanding the physics of two-dimensional turbulence is critical due to their applicability in diverse fields

3

such as in geophysical flows, astrophysical flows and plasma physics and a considerable amount of theoretical and

4

numerical attention has been provided to this concept over the last few decades [1, 2]. Two-dimensional turbulence

5

provides a suitable platform to understand the behavior of flows such as geophysical flows in the atmosphere and

6

the ocean [3, 4] or for that matter, in any particular flow where an idealized configuration can be assumed with flow

7

being restricted to two dimensions only such as with flows in rapidly rotating systems and in thin fluid films over rigid

8

objects [5].

9

CR IP T

2

An important difference is seen in the behavior of two-dimensional turbulence when compared to three-dimensional turbulence with respect to energy cascade behavior. Two-dimensional turbulence follows the Kraichnan-Batchelor-

11

Leith (KBL) theory [6–8] of energy cascade where energy is transferred from the smaller scales to the larger scales

12

(described as an inverse cascade of energy) as against the forward cascade of energy from the larger scales to the

13

smaller scales in three-dimensional turbulence through the process of vortex stretching (a phenomenon absent in two-

14

dimensional flows). In contrast, it is seen that another quadratic invariant, the enstrophy (defined as the integral of the

15

square of the vorticity) cascades from the larger scales to the smaller ones. This dual cascade (i.e. forward enstrophy

16

and inverse energy) leads to a considerable richness in the dynamics of two-dimensional turbulence as compared to

17

the three-dimensional turbulence. The KBL theory of dual cascade has been a subject of close investigation over the

18

past few years [1, 9] and has been confirmed using high fidelity computational studies for both forced (stationary)

19

turbulence [2, 10, 11] and unforced (decaying) turbulence [12–14].

M

AN US

10

In order to simulate large-scale turbulent flows with the dissipation of energy through the kinematic viscosity of

21

the fluid only, an impractically fine resolution would be required to resolve all the turbulent flow structures. Hence,

22

artificial dissipation mechanisms need to be built into the mathematical model for the purpose of computational

23

efficiency. An approximate quantification of the dissipative effects of the smaller scales can be used to simulate

24

spatially averaged flow behavior at larger scales. Reliable modelling strategies with inbuilt mechanisms for proper

25

dissipation characteristics are necessary for the modeling of turbulent flows in general circulation models for oceans

26

and the atmosphere. These models have traditionally been used along with the prescription of ad-hoc friction and eddy

27

viscosity coefficients since the horizontal scale of an ocean basin is much larger than the effective scale for molecular

29

30

PT

CE

AC

28

ED

20

diffusion caused by the kinematic viscosity of seawater. On the other hand, large eddy simulation (LES) has been proven to be a promising approach for the simulations

of 3D turbulent flows and has also been used to simulate some cases of 2D turbulence [15–18]. The basic principle

31

of an LES is to decompose the flow variables into resolved and unresolved spatial scales through the application of

32

a low pass spatial filter to the governing equations. The low-pass spatial filter is responsible for the removal of the

33

high resolution requirement of the governing equations by filtering out the high wavenumber (or finer scale) content

34

of the flow. This allows for much coarser meshes and correspondingly lower computational costs compared to direct 3

ACCEPTED MANUSCRIPT

35

numerical simulation (DNS). However, a successful LES needs to treat the closure problem where a suitable modeling

36

methodology must be chosen to quantify the effects of the finer resolutions (which have been filtered out) on the larger

37

scales (which have been preserved) [19, 20]. The choice of this model, also known as the subgrid scale model (SGS),

38

affects the quality of the LES immensely and as such has been the focus of intense research [21–28].

39

It is common to utilize user-defined heuristics to add artificial dissipation to the mathematical models of twodimensional turbulence. This approach is also known as the functional approach to artificial eddy viscosity parametriza-

41

tion where the assumption of small scale isotropy is utilized to represent the universal characteristic of the dissipation

42

of the scales in the dissipation range [29]. One of the most famous functional models based on the above assumptions

43

is the Smagorinsky model [21] which has been used for LES predictions for a variety of flows. Smagorinsky’s hypoth-

44

esis was that the modeling of the larger scales was more important since they carried a vast majority of all the energy

45

of the flow and the smaller scales could be approximated by an additional dissipative term in the spatially filtered

46

governing equations. It must be noted here that the Smagorinsky model also assumes that the rate of production of

47

energy at the largest scales equals the dissipation of energy at the smaller scales (which in itself is a cornerstone of

48

Kolmogorov’s famous theory of turbulence). The magnitude of artificial dissipation is computed from the resolved

49

strain rate and characteristic length scale which are assumed to be proportional to the filter width via a Smagorinsky

50

constant. On the other hand, the Leith model [8, 30, 31] has been applied extensively for two-dimensional geophysical

51

flows where the cascade of enstrophy from the larger to the smaller scales forms the basis of the selection of relevant

52

scales for resolution. Following a similar argument to the conservation of the turbulent energy between the larger

53

and smaller scales, it is assumed that the production of enstrophy at the larger scales is balanced by its dissipation at

54

the smallest scales. Similar to the Smagorinsky concept, a Leith coefficient is defined which controls the quantity of

55

artificial dissipation to be added to the flow. As such, the Leith formulation reflects physics more accurately since it

56

based on the resolution of the enstrophy cascade in accordance with KBL theory [32].

ED

M

AN US

CR IP T

40

Dynamic adaptations to these functional viscosities have been implemented where the coefficients which control

58

the dissipation of the model are determined adaptively by the use of a low-pass spatial test filter [24, 33]. Their

59

requirement stems from the fact that the value of the constant coefficient in the aforementioned functional models has

60

not been observed to be single-valued [34–38]. These dynamic formulations have the benefit of being applicable for

61

flows which are not perfectly isotropic and homogeneous and have thus been applied successfully to a large number

62

of engineering flows [17, 18]. Although the dynamic model allows the constant in an eddy viscosity model to vary in

64

65

CE

AC

63

PT

57

space and time depending on the local flow structures, the ensemble averaging process prohibits the true backscattering in order to limits the growth of flow instabilities in practice [39–41]. In contrast to the functional approach, structural models attempt to model the effect of the unresolved scales

66

on the larger flow features without the use of any physical assumptions or additional phenomenological arguments.

67

Approximate deconvolution (AD) is one of the most popular structural models with proven performance in large-

68

scale geophysical flows with inverse energy cascades [42–44]. Using the concept of deconvolution from the image

69

processing community, AD attempts to recreate an approximation of the subgrid field variables [45, 46] to close the 4

ACCEPTED MANUSCRIPT

LES governing equations. The deconvolution procedure is carried out by using the Van Cittert iterations which employ

71

repeated filtering operators in an iterative substitution methodology. AD has been successfully used in 3D turbulent

72

engineering flows [47–49] and atmospheric boundary layer flows [50–53] and has also been investigated thoroughly

73

from a theoretical point of view [54–61]. Similar to the concept of AD is the scale similarity models of Layton [62]

74

and Bardina [63] where the LES governing equations are closed by the process of an extrapolation from the resolved

75

scales to the subgrid scales. These models are characterized by using the filtered field variable value as the value of

76

the subgrid scale contribution and have proven accurate in a-priori testing.

CR IP T

70

Structural and functional models have often been coupled to introduce an enhanced degree of stabilization to the

78

mathematical model. This stability has been introduced in the form of a regularization to the AD method through the

79

definition of a penalty term [47] or through the coupling of a dynamic mixed scale and AD [64]. The regularization

80

may also be introduced through the use of a Smagorinsky or Leith type viscosity to arrest any leftover aliasing error

81

from the structural closure used. Another popular procedure for regularization is the use of an explicit low-pass filter

82

to account for the dissipative effect of the residual stress in LES. Infact the explicit (or relaxation) filtering approach

83

has proven successful enough to be used independently of any subgrid modeling methodology [65, 66]. This filtering

84

approach may be carried out at the end of every timestep which ensures a residual stress dissipation. We note here

85

that the efficacy of the subgrid scale modeling relies heavily on the choice of the low-pass spatial filter. A wide

86

variety of filters have been examined for their suitability to these applications and correlations between filter shape

87

and closure performance have also been elaborated [67–70]. Examples of low-pass spatial filters include differential

88

filters [71, 72], Pad´e (compact) filters [73, 74], discrete selective filters [75–77].

M

AN US

77

In this paper, we put forth a unified framework of scale-aware parameterizations by minimizing the error between

90

structural and functional models. Our hybrid models use both Smagorinsky and Leith viscosity kernels to charac-

91

terize an ‘eddy-viscosity’ based on the local strain rate deformation and vorticity gradient respectively. The free

92

modeling coefficients associated with the nonlinear viscosity kernels are determined using two methodologies. The

93

minimization of error between the structural and functional approach broadly follows the same concept used for the

94

Germano-Lilly dynamic model which minimizes the error between the low-pass filtered variable to use in the LES

95

governing equation and a low-pass test filtered variable. The standard methodology requires the prescription of fixed

96

values for the constants prior to simulation whereas the dynamic methodology involves the calculation of these coeffi-

97

cients during the course of the simulation through user-defined heuristics. We also detail the results of using structural

99

100

PT

CE

AC

98

ED

89

models which do not require the use of any phenomenological quantity during the course of the simulation. The models used in our structural framework are the Bardina, Layton and approximate deconvolution methods. A hybrid approach that combines the functional and structural approaches has been used to calculate the coefficients of the non-

101

linear viscosity kernels. Both Smagorinsky’s and Leith’s coefficients are dynamically computed within the proposed

102

hybrid approaches. These hybrid approaches are fundamentally different than the classical dynamic approaches based

103

on the test filtering procedure. The general derivation of the dynamic modeling framework proposed by Germano and

104

Lilly is also derived for the 2D turbulence case. As a secondary objective, the effect of the numerical scheme used to 5

ACCEPTED MANUSCRIPT

105

discretize the governing equations is also quantified by using both fourth-order Pad´e and energy conserving Arakawa

106

schemes. We expect to find similar trends in the behavior of the subgrid scale models but also expect a significant

107

dependence on the underlying numerical discretizations. The following paper is organized as follows: We introduce the governing equations for two-dimensional turbu-

109

lence in section 2 and detail its energy cascade process. Section 3 contains the procedure for obtain the LES governing

110

equations and the various models that will be used to close these equations. Section 4 contains the underlying numer-

111

ical methods used for the spatial discretization of the LES governing equations as well as a brief description of the

112

Poisson solver and the explicit time integration method. The low-pass spatial filtering operators used in this investi-

113

gation are also outlined here. The results obtained from this work are shown in Section 5 and concluding remarks are

114

made in Section 6.

115

2. Two-Dimensional Turbulence

117

AN US

116

CR IP T

108

On taking the curl of the 2-D Navier Stokes equations, the governing equation for incompressible flow may be expressed in the vorticity streamfunction formulation given in its following dimensionless form: ∂ω + J(ψ, ω) = ν∇2 ω, ∂t

119

with ν = 1/Re represents the inverse Reynolds number and ω = ∇ × u is the dimensionless scalar vorticity for two-dimensionless flows. The Jacobian given by J(ψ, ω) symbolizes the nonlinear interactions and can be defined as

M

118

ED

J(ψ, ω) = 120

(1)

∂ψ ∂ω ∂ψ ∂ω − . ∂y ∂x ∂x ∂y

(2)

The streamfunction ψ is related to the velocity through

PT

u=

∂ψ ∂ψ ; v=− ∂y ∂x

(3)

with u and v being the components of the two-dimensional velocity field u = (u, v). Incompressibility is enforced

122

through the divergence free condition ∇.u = 0 which imposes a kinematic relationship between the vorticity and

124

125

126

streamfunction according to the following Poisson equation

AC

123

CE

121

∇2 ψ = −ω.

(4)

The decision to use the vorticity-streamfunction formulation over the primitive variable formulation for the NavierStokes equations are due to several advantages such as an elimination of the pressure term and its associated odd-even decoupling between the pressure and velocity components. The formulation we employ also helps us avoid projection

127

inaccuracies seen in fractional step methods and automatically satisfies the divergence free condition while reducing

128

the number of equations to be solved [78]. The geometry we solve for is also assumed to have periodic boundary con-

129

ditions and a uniform Cartesian grid which eliminates errors from mesh non-uniformities and inconsistent boundary

130

schemes. 6

ACCEPTED MANUSCRIPT

131

3. Large Eddy Simulations of 2D Turbulence

132

The governing equations for LES can be obtained by low-pass filtering the two-dimensional Navier-Stokes equa-

133

tions which have been derived in their vorticity and streamfunction formulation in the previous section. On performing

134

this spatial filtering operation we have: ∂ω + J(ψ, ω) = ν∇2 ω. ∂t The above equation may further be expressed as ∂ω + J(ψ, ω) = ν∇2 ω + Π, ∂t

CR IP T

135

(5)

(6)

where the Π is considered a source term representing the effect of the subgrid scale stresses. It should be noted that

137

the primary reason for the formation of this source term is due to the fact that the spatial filtering operation does not

138

commute for nonlinear terms i.e.

AN US

136

(7)

J(ψ, ω) , J(ψ, ω). 139

Thus we can represent our source term in the following manner

(8)

Π = J(ψ, ω) − J(ψ, ω)

The modeling of this source term (which can be assumed to encompass the effects of the unresolved scales) is the primary challenge of LES closure research.

142

3.1. Functional models

143

3.1.1. Smagorinksy model

PT

145

The standard Smagorinsky model approximates the source term in the LES governing equations using a simplified eddy viscosity hypothesis. The Smagorinsky interpretation of the source term is:

CE

144

ED

141

M

140

Π = νe ∇2 ω

(9)

where ω is the filtered quantity we are attempting to model. The eddy viscosity νe is described in terms of a mixing

147

length and the absolute value of the strain tensor of the resolved flow field. This eddy viscosity controls the extra

148

149

150

AC

146

dissipation required to now account for the effect of the dissipation range (we have spatially filtered out). We have νe = l02 |S¯ |

(10)

where l0 is the mixing length usually given by l0 = C s δ in the standard Smagorinsky model and |S¯ | is based on the

second invariant of the filtered field deformation tensor given by q |S¯ | = 2S¯ i j S¯ i j . 7

(11)

ACCEPTED MANUSCRIPT

151

152

In the vorticity-streamfunction formulation used for our framework, the above equation reduces to s !2 !2 ∂2 ψ¯ ∂2 ψ¯ ∂2 ψ¯ ¯ + |S | = 4 − 2 . ∂x∂y ∂x2 ∂y

(12)

The final expression of the SM hypothesis can be represented by νe = (C s δ)2 |S¯ |,

CR IP T

(13)

where we note that the proposed model uses the respective eddy viscosity formula in the vorticity transport equation

154

and not the curl of the model in the primitive variables (i.e., see [79] for the definition of subgrid-scale torque). One

155

can observe a free-modeling parameter C s in the previous expression. This constant controls the magnitude of artificial

156

dissipation introduced into the system through the LES governing equations. However, values for this constant are

157

seen to vary in literature and therefore it is physically more appropriate to obtain this constant dynamically through

158

the process of simulation. The LES runs initiated using the Smagorinsky closure utilize the theoretically accepted

159

value of C s = 0.2 for homogeneous and isotropic turbulence.

160

3.1.2. Leith model

AN US

153

The Leith [8, 30, 31] version of an artificial viscosity is employed by focusing on resolving the forward cascade

162

of enstrophy in 2D turbulence rather than the forward cascade of energy in 3D turbulence (which forms the basic

163

assumption of Smagorinsky type models). An expression for the Leith model can be derived in a manner similar to

164

the Smagorinsky model to give

ED

where

PT

165

M

161

νe = (Cl δ)3 |∇ω|, ¯

|∇ω| =

s

∂ω ∂x

!2

!2 ∂ω + . ∂y

(14)

(15)

The runs initiated using the standard Leith closure also utilize an artificial viscosity coefficient given by Cl =

167

0.2. We note that the Leith model is similar to the well-known Baldwin-Lomax model (i.e., νe = l02 |ω|) which

169

170

171

identifies a mixing length from the filter radius used to eliminate high wavenumber content. We also note that the

AC

168

CE

166

Baldwin-Lomax model and the Smagorinsky model are both observed to be more dissipative for flows with large scale rotational structure. Our preliminary results showed that the trends in the Baldwin-Lomax (not shown in this study) and the Smagorinsky models are quite similar. Therefore, in this study, we only consider the Leith viscosity

172

since it is proportional to the vorticity gradient and l03 . It must be noted that the implementation of the functional

173

models elaborated above differ from their traditional application in the primitive variable formulation of the Navier-

174

Stokes equations in that we directly assume an eddy viscosity (and thus a source term) for the vorticity-streamfunction

175

formulation rather than applying a curl operator to an eddy viscosity obtained in the primitive approach. 8

ACCEPTED MANUSCRIPT

176

3.2. Structural models As mentioned previously, structural models differ from their functional counterparts based on the requirement of

178

phenomenological arguments for the addition of artificial dissipation. In order to correctly quantify the effect of the

179

smaller unresolved scales on the larger resolved eddies, certain approximation strategies are used to close the LES

180

governing equations.

181

3.2.1. Bardina’s scale similarity model

CR IP T

177

182

The Bardina scale similarity model [63] makes the observation that the smallest resolved scales of the LES gov-

183

erning equations are more or less similar to the largest unresolved scales of the subgrid scales. It then proposes to use

184

this information for constructing the source term. We have Π = J(ψ, ω) − J(ψ, ω) which may be approximated as

AN US

185

Π = J(ψ, ω) − J(ψ, ω) 186

(16)

(17)

The advantages of this formulation include the allowance of backscatter and proven accuracy in a-priori testing. However, it has been known to suffer from stability issues.

188

3.2.2. Layton model

M

187

Layton [62] proposed an improvement to the scale similarity model by making a simple adjustment to the Bardina

190

formulation where the LES governing equations were setup in a slightly different manner without the express use of a

191

source term. The source term of Layton’s version of the scale similarity model is given by

PT

ED

189

Π = J(ψ, ω) − J(ψ, ω)

The differences with the Bardina model is apparent and it has proven to be more stable as well.

193

3.2.3. Approximate deconvolution method

CE

192

(18)

The approximate deconvolution (AD) procedure abandons the scale similarity hypothesis (which attempts to

195

model the unresolved scales identically to the smallest resolved scales) and proceeds to recover the unresolved scales

196

197

AC

194

without the use of any physical assumption. This is obtained through the use of Van Cittert iterations, which are a popular iterative resubstitution methodology used for the deblurring of images. The basic principle of Van Cittert

198

iterations are the use of repeated filtering as shown below. If ω∗ is assumed to be the approximately recovered value

199

for the unfiltered quantity ω, we have ω∗0 = ω ¯ ω∗i = ω∗i−1 + β(ω ¯ − G ? ω∗i−1 ), 9

(19) i = 1, 2, 3, ...., N

ACCEPTED MANUSCRIPT

200

where G is our low pass spatial filtering operation (i.e., G ? ω = ω) ¯ and further iterations of the above scheme lead

201

to improved approximations of the unfiltered variable. Here, β is an over-relaxation parameter which may be used to

202

speed up convergence but must be restricted to 0 ≥ β ≥ 2 to ensure a positive and bounded transfer function for the low-pass filter. We have fixed β = 1 in our investigation. The approximately recovered field variable ω∗ now has the

204

contributions of the finer scales of the flow and may be used to determine the advective term in the LES governing

205

equations before spatial filtering. Our source term can thus be represented as Π = J(ψ, ω) − J(ψ∗ , ω∗ ).

206

3.3. Test filtering approach: The Germano-Lilly dynamic model

CR IP T

203

(20)

The functional approach to SGS modeling is limited by the fact that the constant the controls the artificial dissipa-

208

tion in the model is not observed to be single valued. Therefore it was necessary to obtain this constant dynamically

209

from the flow evolution rather than prescribing it as a fixed parameter before time integration. Germano [24] proposed

210

a dynamic estimation of the Smagorinsky coefficient from the resolved flow itself. This dynamic model was further

211

improved by Lilly [33] who proposed a least squares based calculation for C s . The traditional procedure for dynamic

212

estimation relied on an additional filtering process (also known as test filtering). The dynamic estimation process may

213

be used for both the Smagorinsky and Leith viscosity kernels as shown in the following sections.

214

3.3.1. Dynamic model using Smagorinky’s nonlinear viscosity kernel: D-S-TF

M

216

In this section, we derive the dynamic Smagorinsky test filtered (D-S-TF) model for two-dimensional flows. The LES governing equation can be written in its spatially filtered form as

ED

215

AN US

207

∂ω + J(ψ, ω) = ν∇2 ω + (C s δ)2 |S¯ |∇2 ω, ∂t

(21)

where we have used the Smagorinsky viscosity kernel to model the source term Π. The 2D turbulence governing

218

equations can also be written after a filtering operations with a low-pass spatial test filter

221

222

CE

220

∂e ω e, ω e) = ν∇2 ω e + (C se e. + J(ψ δ)2 |Se|∇2 ω ∂t

(22)

Here the filter scale is given by e δ(e δ > δ). An assumption we make here is that the difference between the test filtered fields e u and e u¯ become negligible leading to the following expression

AC

219

PT

217

e ∂ω e = ν∇2 ω e + (C se e + J(e ψ, ω) δ)2 |e S |∇2 ω. ∂t

(23)

Now we can test filter the LES governing equation in its low-pass spatially filtered form (Eq. 21) to get e ∂ω gω) = ν∇2 ω e + (C s δ)2 |S^ ¯ |∇2 ω. + J(ψ, ∂t

(24)

On subtracting equations (24) from (23) we get

gω) = (C e 2 e 2e 2 ^ e − J(ψ, ¯ 2 J(e ψ, ω) s δ) |S |∇ ω − (C s δ) |S |∇ ω. 10

(25)

ACCEPTED MANUSCRIPT

224

225

226

227

The above expression may be written in a simplified form given by H = (C s δ)2 M

(26)

gω) e − J(ψ, H = J(e ψ, ω)

(27)

where

and

e − |S^ ¯ |∇2 ω M = κ2 |e S |∇2 ω

CR IP T

223

(28)

in which κ = e δ/δ is the filter width ratio. Germano proposed the calculation of the Smagorinsky constant from Eq.

(25) but this issue could face the possible issue of a division by zero at possible locations in the fluid flow with very low vorticity gradients. Lilly proposed a procedure whereby this stability issue could be avoided by using the method

229

of least squares to minimize the average square error in the domain hE 2 i, where the error is given by E = H −(C s δ)2 M.

230

AN US

228

We utilize the latter methodology in this investigation. The averaging operator chosen here is expressed by Z lx Z ly 1 f (x, y)dx dy, hfi = l x ly 0 0

(29)

231

where l x and ly are the horizontal and vertical lengths of the periodic domain. The least squares minimization problem

232

is solved by minimizing the mean squared error with respect to our model control parameter (C s δ)2 to get

M

∂hE 2 i = −2hHMi + 2(C s δ)2 hM 2 i. ∂(C s δ)2

(30)

The left hand side of the above equation becomes zero when the error is minimized to give us the following expression

234

for the Smagorinsky constant

(C s δ)2 =

PT

235

ED

233

hHMi . hM 2 i

(31)

3.3.2. Dynamic model using Leith’s nonlinear viscosity kernel: D-L-TF A similar procedure as outlined in the previous section may be used to obtain an expression for the dynamic update

237

of the Leith constant. In order to derive the dynamic Leith test filtered (D-L-TF) model, our test filtered equation is

238

given by

240

AC

239

CE

236

∂e ω e, ω e) = ν∇2 ω e + (Cle e, + J(ψ δ)3 |∇e ω|∇2 ω ∂t

(32)

which may be considered equivalent to

e ∂ω e = ν∇2 ω e + (Cle e 2 ω. e ψ, ω) δ)3 |∇ω|∇ + J(e ∂t

(33)

The test filtering of the LES governing equation with the Leith viscosity kernel gives us e ∂ω gω) = ν∇2 ω ^ 2 ω. e + (Cl δ)3 |∇ω|∇ + J(ψ, ∂t 11

(34)

ACCEPTED MANUSCRIPT

242

243

244

In a procedure similar to the dynamic Smagorinsky procedure, we may now subtract Eq. (34) from (33) to get gω) = (C e 3 e 2e 3 ^ 2 e − J(ψ, ψ, ω) J(e l δ) |∇ω|∇ ω − (C l δ) |∇ω|∇ ω.

(35)

H = (Cl δ)3 M

(36)

The above expression can now be expressed as

with gω) e − J(ψ, ψ, ω) H = J(e

and

^ 2ω e 2ω e − |∇ω|∇ M = κ3 |∇ω|∇

CR IP T

241

(37)

(38)

where κ is the filter ratio. Using a least squares minimization procedure as shown in the dynamic Smagorinsky method

246

we get

AN US

245

(Cl δ)3 = 247

hHMi . hM 2 i

(39)

3.4. Hybrid dynamic approach: Minimizing the error between functional and structural models In this subsection we introduce a novel methodology for the closure of the LES governing equations which com-

249

bines elements of both the dynamic modeling approach (for the functional viscosity) and the structural modeling

250

approach for approximating the nonlinear Jacobian. This procedure does not require the use of a test filtering opera-

251

tor.

252

3.4.1. Dynamic hybrid model coupling Smagorinky’s nonlinear viscosity kernel with the Bardina model: D-S-B

ED

PT

The source term in the LES governing equation (Eq. 6) is written as

CE

253

M

248

Π = J(ψ, ω) − J(ψ, ω).

(40)

We focus on the reconstruction of the unresolved scales required for the second term on the right hand side of the

255

above equation. Using the Bardina scale similarity hypothesis we can write the source term as

256

257

AC

254

Π = J(ψ, ω) − J(ψ, ω).

(41)

The source term can then be equated to the Smagorinsky viscosity hypothesis to get J(ψ, ω) − J(ψ, ω) = (C s δ)2 |S |∇2 ω

(42)

L = (C s δ)2 M

(43)

which can be represented as

12

ACCEPTED MANUSCRIPT

259

where L = J(ψ, ω) − J(ψ, ω)

(44)

M = |S |∇2 ω.

(45)

and

CR IP T

258

260

We may now use an error minimization process similar to the dynamic test filtering approach to minimize an ensemble

261

averaged error given by E = L − (C s δ)2 M through the expression (C s δ)2 =

hLMi . hM 2 i

This completes the derivation of hybrid D-S-B model.

263

3.4.2. Dynamic hybrid model coupling Leith’s nonlinear viscosity kernel with the Bardina model: D-L-B As shown in the previous sections, the dynamic estimation of the source term may also be done using the Leith

264

265

AN US

262

viscosity. Using a similar procedure to the D-S-B model derived above, we obtain D-L-B model as (Cl δ)3 = where

hLMi hM 2 i

M

266

and

272

(48)

(49)

When combined with Layton’s [62] scale similarity method for the reconstruction of the finer scales we obtain the

CE

following expression for the Smagorinsky constant

AC

271

(47)

3.4.3. Dynamic hybrid model coupling Smagorinky’s nonlinear viscosity kernel with the Layton model: D-S-L

269

270

M = |∇ω|∇2 ω.

PT

268

ED

L = J(ψ, ω) − J(ψ, ω)

267

(46)

(C s δ)2 =

hLMi hM 2 i

(50)

where

L = J(ψ, ω) − J(ψ, ω)

(51)

M = |S |∇2 ω.

(52)

and

13

ACCEPTED MANUSCRIPT

3.4.4. Dynamic hybrid model coupling Leith’s nonlinear viscosity kernel with the Layton model: D-L-L For the nonlinear Leith viscosity kernel we have, when using the Layton scale similarity method for reconstruction

274

275

of unresolved scales (Cl δ)3 =

276

hLMi hM 2 i

where L = J(ψ, ω) − J(ψ, ω)

277

and M = |∇ω|∇2 ω.

(54)

(55)

3.4.5. Dynamic hybrid model coupling Smagorinky’s nonlinear viscosity kernel with the AD model: D-S-AD

AN US

278

(53)

CR IP T

273

279

Using the approximate deconvolution method to reconstruct the nonlinear advective terms gives us the approx-

280

imate values of the unfiltered variables ω∗ and ψ∗ which may now be used in our relation to compute the model

281

coefficients dynamically. We have, for the Smagorinsky viscosity, a dynamic AD (D-S-AD) model given by (C s δ)2 = where

M

282

hLMi hM 2 i

284

and

The hybrid dynamic and approximate deconvolution method when used with the Leith viscosity kernel gives us

288

289

CE

the follow expression for the dynamic update of model coefficients

AC

287

(58)

3.4.6. Dynamic hybrid model coupling Leith’s nonlinear viscosity kernel with the AD model: D-L-AD

285

286

M = |S |∇2 ω.

(57)

PT

283

ED

L = J(ψ, ω) − J(ψ∗ , ω∗ )

(56)

(Cl δ)3 =

hLMi hM 2 i

(59)

where

L = J(ψ, ω) − J(ψ∗ , ω∗ )

(60)

M = |∇ω|∇2 ω.

(61)

and

This model is referred to as D-L-AD in this paper. 14

ACCEPTED MANUSCRIPT

290

291

4. Numerical Schemes The following section outlines the various numerical schemes used in this investigation for the calculation of

292

derivatives, nonlinear terms, solution of the Poisson equation and spatial filtering operators.

293

4.1. Pad´e scheme For the purpose of the first and second order derivatives required in the governing equations of two-dimensional

295

turbulence, we use compact schemes instead of standard finite differences to restrict stencil sizes with increased

296

spectral resolution. Compact schemes also give us the benefit of reducing dissipation error at a minimal increase in

297

the computational cost required to solve a linear system of equations using a tridiagonal solver. A compact tridiagonal

298

scheme for the first derivative is given by

299

fi+2 − fi−2 fi+1 − fi−1 +b . 2h 4h

(62)

AN US

0 0 α fi−1 + fi0 + α fi+1 =a

CR IP T

294

The constants on the right hand side of the scheme are given by a = 32 (α + 2) and b = 13 (4α − 1). The subscript i,

300

represents the spatial grid index with h as its uniform grid spacing in a particular direction. Setting α = 0 recovers

301

the explicit non-compact fourth order scheme for the first derivative and the classical compact fourth-order scheme

302

popularly known as the Pad´e scheme is recovered by setting α = 1/4.

In a similar fashion, we obtain a compact scheme for the second derivative given by

M

303

00 00 α fi−1 + fi00 + α fi+1 =a

fi+2 − 2 fi + fi−2 fi+1 − 2 fi + fi−1 +b 2 h 4h2

where a = 43 (1 − α) and b = 31 (10α − 1). The classical fourth order Pad´e is recovered by setting α = 1/10.

305

4.2. Arakawa scheme

ED

304

(63)

It was proposed by Arakawa [80] that the conservation of energy, enstrophy and skew symmetry is sufficient to

307

avoid computational instabilities arising from nonlinear interactions in the Jacobian. The Jacobian in the governing

308

equations for 2D turbulence is defined as

310

CE

AC

309

PT

306

J(ψ, ω) =

∂ψ ∂ω ∂ψ ∂ω − ∂y ∂x ∂x ∂y

(64)

The second-order Arakawa scheme for the Jacobian is given by JI (ω, ψ) =

1 (J1 (ω, ψ) + J2 (ω, ψ) + J3 (ω, ψ)) 3

(65)

where the discrete parts of the Jacobian are J1 (ψ, ω) =

1 h (ωi+1, j − ωi−1, j )(ψi, j+1 − ψi, j−1 ) 4h x hy i −(ωi, j+1 − ωi, j−1 )(ψi+1, j − ψi−1, j ) 15

(66)

ACCEPTED MANUSCRIPT

J2 (ψ, ω) =

1 h ωi+1, j (ψi+1, j+1 − ψi+1, j−1 ) 4h x hy

J3 (ψ, ω) =

−ωi−1, j−1 (ψi−1, j − ψi, j−1 ) − ωi−1, j+1 (ψi, j+1 − ψi−1, j ) i +ωi+1, j−1 (ψi+1, j − ψi, j−1 )

(67)

(68)

The fourth order accurate Arakawa discretization of the Jacobian becomes JII (ψ, ω) =

312

1 h ωi+1, j+1 (ψi, j+1 − ψi+1, j ) 4h x hy

1 (J4 (ω, ψ) + J5 (ω, ψ) + J6 (ω, ψ)) 3

where

1 h (ωi+1, j+1 − ωi−1, j−1 )(ψi−1, j+1 − ψi+1, j−1 ) 8h x hy i −(ωi−1, j+1 − ωi+1, j−1 )(ψi+1, j+1 − ψi−1, j−1 )

(69)

(70)

M

J4 (ψ, ω) =

AN US

311

CR IP T

−ωi−1, j (ψi−1, j+1 − ψi−1, j−1 ) − ωi, j+1 (ψi+1, j+1 − ψi−1, j+1 ) i +ωi, j−1 (ψi+1, j−1 − ψi−1, j−1 )

1 h ωi+1, j+1 (ψi, j+2 − ψi+2, j ) 8h x hy

(71)

1 h ωi+2, j (ψi+1, j+1 − ψi+1, j−1 ) 8h x hy

−ωi−2, j (ψi−1, j+1 − ψi−1, j−1 ) − ωi, j+2 (ψi+1, j+1 − ψi−1, j+1 ) i +ωi, j−2 (ψi+1, j−1 − ψi−1, j−1 )

(72)

J(ψ, ω) = 2JI (ψ, ω) − JII (ψ, ω) + O(h4 )

(73)

ED

J5 (ψ, ω) =

313

AC

CE

PT

−ωi−1, j−1 (ψi−2, j − ψi, j−2 ) − ωi−1, j+1 (ψi, j+2 − ψi−2, j ) i +ωi+1, j−1 (ψi+2, j − ψi, j−2 ) J6 (ψ, ω) =

It is seen that JII conserves enstrophy and energy and that our Jacobian can be represented by

314

with fourth order accuracy. This scheme has successfully been used to compute two-dimensional isotropic turbulence

315

by Orlandi [81].

16

ACCEPTED MANUSCRIPT

316

317

318

4.3. Time integration scheme An optimal third-order-accurate total variation diminishing Runge-Kutta scheme (TVDRK3) is used for explicit advancement in time. In order to implement this scheme the model equations are cast in the following form du = £(u) dt

(74)

where £(u) encompasses the spatial derivatives including the nonlinear convective terms,the linear diffusive terms and

320

the source terms modeling the subgrid scale stresses. Assuming that the numerical approximation at a time level n + 1

321

is known and our time step is given by ∆t, our time stepping scheme becomes [82]u(1) = un + ∆t£(un ), 3 n u + 4 1 = un + 3

u(2) =

1 ∆t£(u(1) ), 4 2 ∆t£(u(2) ). 3

(75)

AN US

un+1

1 (1) u + 4 2 (2) u + 3

CR IP T

319

The TVDRK3 scheme has been shown to predict slightly more accurate results than some of the other third order

323

Runge-Kutta schemes for incompressible flow problems [83] and has been extensively used in hyperbolic conservation

324

laws [84–87].

325

4.4. Poisson solver

M

322

The elliptic equation obtained for the relationship between the streamfunction and vorticity is solved by a FFT-

327

based Poisson solver. The low-pass spatial filtered version of the elliptic equation given by Eq. (4) can be written in

328

the form of ∇2 u = f . The compact fourth-order discretization with nine point stencil can be written as [88]

ED

326

+ d(ui+1, j+1 + ui+1, j−1 + ui−1, j+1 + ui−1, j−1 )

(76)

= e(8 fi. j + fi+1, j + fi−1, j + fi, j+1 + fi, j−1 ),

where the coefficients are a = −10(1 + γ2 ), b = 5 − γ2 , c = 5γ2 − 1, d = (1 + γ2 )/2, e = h2x /2 with γ being the mesh

CE

329

PT

aui, j + b(ui+1, j + ui−1, j ) + c(ui, j+1 + ui, j−1 )

aspect ratio defined as the ratio of grid discretization lengths h x /hy . In our case γ = 1 for which this scheme becomes

331

the famous Mehrstellen scheme [89–91]. The presence of periodic boundary conditions in both directions suggests

332

333

334

AC

330

the use of the Fourier transform approach as against a finite sine or cosine transform. First we perform and inverse transform for the source term to get fˆ =

N x −1 N y −1 X 1 X fi, j e−i(kx xi +ky y j ) N x Ny i=0 j=0

(77)

where k x = 2πm/L x and ky = 2πn/Ly . Our elliptic equation in the Fourier space is given by (−k2x − ky2 )ˆum,n = − fˆm,n 17

(78)

ACCEPTED MANUSCRIPT

335

This equation may now be solved for the unknown uˆ k,l . Once this is achieved, a forward Fourier transform into the

336

real domain is obtained by Ny

Nx 2

ui, j =

−1 2 −1 X X

m=− N2x

uˆ m,n ei(kx xi +ky y j ) .

(79)

N n=− 2y

For the purpose of transforming to and back from the Fourier domain, we use the FFT algorithm given in Press et al

338

[92].

339

4.5. Discrete filtering operators

CR IP T

337

340

To conclude the discussion on numerical methods, we define a spatial filtering operator with a selected behavior to

341

ensure low-pass characteristics. This is to ensure the removal of high frequency content from the numerical solution

342

procedure. We use the fourth-order Pad´e compact filter family given by [73] 2 X as

AN US

α f¯j−1 + f¯j + α f¯j+1 =

s=0

343

where the constants a s are given by a0 =

1 (5 + 6α), 16

a1 =

2

( f j−s + f j+s ),

1 (1 + 2α), 2

1 a2 = − (1 − 2α). 8

(80)

(81)

In the above expression, α is a free-modeling parameter which has been kept at a value of α = 0 throughout this

345

investigation which gives us the following expression for the discrete filtering operator.

M

344

ED

 1 − f j−2 + 4 f j−1 + 10 f j + 4 f j+1 − f j+2 f¯j = 8

which is also known as the binomial smoothing filter (2,1) B [93].

347

5. Results

PT

346

(82)

This section details the performance of the various LES closures outlined in the previous sections. For the pur-

349

pose of comparison, these closures are compared to high fidelity DNS results as well as under-resolved no model

350

simulations. This process is carried out for both the Pad´e and Arakawa numerical discretizations. The 2D decaying

351

turbulence problem is specified by a square domain of side length 2π with periodic boundary conditions in both x

353

AC

352

CE

348

and y directions. The problem involves an initial energy spectrum which decays through time and we examine this decay by studying the evolution of vortices through time. In this study, we compare the energy spectra developed at

354

the end of our simulation (i.e., t = 10) in order to quantify the performance of the SGS closures and the effect of the

355

solver used. The spectra are compared to the theoretical scaling expected from a perfectly homogeneous case of 2D

356

turbulence given by KBL theory. The initial energy spectrum of the problem is given by [81]   E(k) = Ak4 exp −(k/k p )2 18

(83)

ACCEPTED MANUSCRIPT

357

where 4k−5 p

A= 358

and k = |k| =

(84)

3π

q k2x + ky2 . The maximum value of the energy spectrum is designed to occur at wavenumber k p which

is assumed to occur at k p = 10 in this study. The magnitude of vorticity Fourier coefficients related to the assumed

360

initial energy spectrum becomes |e ω(k)| =

362

363

k E(k). π

The initial vorticity distribution in Fourier space is then obtained through the introduction of a random phase r k e(k) = ω E(k)eiζ(k) π

(85)

(86)

where the phase function is given by ζ(k) = ξ(k) + η(k), where ξ(k) and η(k) are independent random values chosen

AN US

361

r

CR IP T

359

in [0, 2π] at each coordinate point in the first quadrant of the k x − ky plane. Their behavior in the other quadrants are

PT

ED

M

given by conjugate relations as shown in Fig. 1. Once the randomization process is completed, the initial vorticity

CE

Figure 1: Conjugate relations for the random phase function for the initial conditions.

364

distribution in the physical domain is obtained through an inverse FFT. It must be noted here that the randomization

366

process is identical for each run (whether LES or DNS) for the purpose of comparison. The initial condition for

367

368

AC

365

vorticity is chosen to ensure a divergence free vorticity field. In the following text, we first perform a code validation by showing the convergence of simulations with progressively finer meshes to the true DNS for both the compact and

369

Pad´e scheme. We then analyze the performance of the various prescribed SGS closures against under-resolved DNS

370

simulations for both numerical solvers. A timestep of ∆t = 5 × 10−4 is selected to ensure that the simulation remains

371

independent of errors associated with temporal discretization (i.e., negligible with respect to the spatial discretization

372

errors). 19

ACCEPTED MANUSCRIPT

373

374

In order to quantify the benefits of each SGS model combined with a particular solver, we define a statistical measure based on the energy spectrum in the wavenumber domain which is defined as ˆ t) = 1 k2 |ψ(k, ˆ t)|2 E(k, 2 and the angle averaged energy spectrum is E(k, t) =

X

1 ´ k− 21 ≤|k|≤k+ 2

´ t). ˆ k, E(

CR IP T

375

(87)

(88)

376

As mentioned previously, the result of the angle averaged energy spectrum for a run is compared to the classical KBL

377

theory scaling which approaches k−3 in the limit of infinite Reynolds number.

378

5.1. Direct numerical simulations

Figs. 2 and 3 show the time evolution of the energy spectra for both Pad´e and Arakawa discretizations of the

380

governing equations. A fully turbulent spectrum is seen to be developed by t = 2 seconds beyond which the spectrum

381

is seen to decay. We have used a resolution of N 2 = 20482 . The Pad´e scheme is seen to have a small amount of

382

aliasing error for the highest Reynolds numbers. The time evolution of the vorticity contours has been shown for the

383

highest Reynolds number (Re = 128,000) in Figs. 4 and 5 where the decay of the vorticity magnitudes is apparent.

384

The final time vorticity contours for different Reynolds number is also detailed in Figs. 6 and 7. It is seen that the

385

increase of Reynolds number corresponds to an increased magnitude of vorticity.

M

AN US

379

We also perform a code validation by carrying out a convergence study of both solvers as shown in Figs. 8 and

387

9. It is seen that both solver variants converge to the same energy spectrum scaling with increasing resolution which

388

confirms the validity of both numerical methods. A resolution of N 2 = 20482 is seen to conform to a high-fidelity DNS

389

solution. These convergence studies are also carried out for different Reynolds numbers and it is seen that the Arakawa

390

scheme is more stable than the compact scheme solver. At higher values of Reynolds numbers the compact scheme

391

fails to converge at lower resolutions. The Arakawa scheme is also more capable in arresting pileup in comparison to

392

its compact counterpart (for the same resolution).

393

5.2. Large eddy simulations

395

396

397

PT

CE

After validation of code performance using DNS runs, the generated DNS data is used (along with its low-pass

AC

394

ED

386

spatially filtered counterpart) to validate the performance of the different closures developed for this study. Both

compact and Arakawa schemes are used for the underlying discretization to ensure that the trends for the closures are universal. Three Reynolds number (Re = 32,000, Re = 64,000, Re = 128,000) are chosen for testing closure

398

performance using three coarse resolutions (N 2 = 1282 , N 2 = 2562 , N 2 = 5122 ). We note here that true DNS was

399

achieved at a minimum resolution of N 2 = 20482 and as such the coarse resolutions used here represent a large

400

reduction in computational cost. The theoretical scaling expected for 2D decaying isotropic homogeneous turbulence

401

are provided in all the figures to indicate the ideal theoretical scaling expected from DNS (and a successful LES) and 20

CR IP T

ACCEPTED MANUSCRIPT

(d) Re = 32,000

ED

M

(c) Re = 16,000

(b) Re = 8,000

AN US

(a) Re = 4,000

(e) Re = 64,000

(f) Re = 128,000

PT

Figure 2: Time evolution of energy spectrum for different Reynolds numbers using the Pad´e scheme with a resolution of N 2 = 20482 .

coarse grid resolutions with no closure modeling are included for the purpose of comparison. This section is divided

403

into two parts with the first one detailing the performance of the standard closures (i.e., without any dynamic update

404

of artificial viscosity coefficients) and the second devoted to the performance of the dynamic closures.

405

5.2.1. Standard (non-dynamic) closure performance

407

AC

406

CE

402

In this subsection, we detail the performance of the standard functional models given by Smagorinsky and Leith

viscosity kernels (without using the dynamic adaptation procedure) and the performance of the structural models given

408

by the Bardina and Layton scale similarity method along with the approximate deconvolution approach. Fig. 10 shows

409

the performance of the standard functional models for the Pad´e and Arakawa schemes respectively for the relatively

410

low Reynolds number of Re = 32,000 when we use C s = 0.2 and Cl = 0.2. Non-convergent results are not included.

411

The Smagorinsky kernel is proven to be more dissipative than the Leith kernel for all three resolutions tested. The 21

CR IP T

ACCEPTED MANUSCRIPT

(d) Re = 32,000

ED

M

(c) Re = 16,000

(b) Re = 8,000

AN US

(a) Re = 4,000

(e) Re = 64,000

(f) Re = 128,000

PT

Figure 3: Time evolution of energy spectrum for different Reynolds numbers using the Arakawa scheme with a resolution of N 2 = 20482 .

Leith kernel is seen to add no particular benefit as a closure method since it remains virtually indistinguishable from

413

no model under-resolved simulation when it converges. The standard Smagorinsky approach is thus, a more stable

414

and accurate closure approach for this test case as it is seen to converge for all the coarse resolutions shown (when

415

using the Pad´e scheme) and provide a considerable improvement in the DNS spectrum approximation.

417

418

AC

416

CE

412

Results using the Arakawa scheme also quantify the benefit of the Smagorinsky model over the Leith model where

it is obvious that the Leith approach adds very little improvement to the under-resolved DNS simulation. We remark here that the Arakawa scheme provides converged run results for all under-resolved DNS and closure approaches

419

due to its energy conserving property. A similar result is seen for the relatively higher Reynolds numbers of Re =

420

64,000 shown in Fig. 11 where the Leith model is seen to be non-convergent for all coarse resolutions using the Pad´e

421

scheme. The same behavior is also seen in the case of Re = 128,000 shown in Fig. 12. However, we emphasize that

422

the amount of dissipation added to the system depends on these modeling constants. Increasing the level of C s or Cl 22

(b) t = 2

AN US

(a) t = 0

CR IP T

ACCEPTED MANUSCRIPT

(e) t = 8

M

(d) t = 6

(c) t = 4

(f) t = 10

ED

Figure 4: Snapshots in time of the vorticity contours for Re = 128,000 using the Pad´e scheme with a resolution of N 2 = 20482 .

423

would add more dissipation to the system. Therefore, we will focus on modeling characteristics when we use dynamic

424

procedures to estimate these constants.

We then test the performance of the structural closures given by the Layton, Bardina and AD methods. The

426

common theme of these methods is the use of the approximations to determine the unresolved scales. We interject

427

here that these approximations are used to close the LES governing equations (in contrast to their use to determine the

428

Smagorinsky and Leith coefficients for the dynamic models). Fig. 13 show the performance of the aforementioned

429

closures for Re = 32,000 for both Pad´e and Arakawa schemes. The AD and Layton approach is immediately seen

430

to be more stable than the Bardina approach which does not yield converged results for the two coarsest resolutions

432

CE

AC

431

PT

425

(N 2 = 1282 and N 2 = 2562 ) when used with the Pad´e scheme. An interesting observation is the behavior of the AD and Layton methods when used with the Arakawa scheme. It is seen that the expected improvement in aliasing error

433

reduction is not pronounced when using the Arakawa scheme for these two methods. While the Bardina method is

434

shown to converge and approximate the the DNS spectrum well for N 2 = 2562 , the AD and Layton methods provide

435

a result that is more or less identical to the Pad´e approach.

436

For the case of Re = 64,000 given by Fig. 14, we see a similar behavior to the lower Reynolds number case 23

(b) t = 2

AN US

(a) t = 0

CR IP T

ACCEPTED MANUSCRIPT

(e) t = 8

M

(d) t = 6

(c) t = 4

(f) t = 10

ED

Figure 5: Snapshots in time of the vorticity contours for Re = 128,000 using the Arakawa scheme with a resolution of N 2 = 20482 .

where the Bardina method converges only for the N 2 = 5122 resolution for the Pad´e scheme. The aforementioned

438

minimal effect of the Arakawa scheme on the AD and Layton methods is still witnessed. As a positive, the AD and

439

Layton methods prove stable for all coarse resolutions for this Reynolds number as well. Fig. 15 shows the case for

440

a Reynolds number of 128,000. The cases where the Pad´e scheme is used show that the Bardina approach becomes

441

totally unstable at this high Reynolds number with no converged simulations for any of the coarse resolutions used.

442

AD and Layton, however still provide converged results with the same results of accuracy as the previous Reynolds

443

number cases. The Arakawa scheme with the Bardina closure at N 2 = 5122 provides an excellent capture of the DNS

444

spectrum as against the AD and Layton approach which show very little improvement in comparison to their Pad´e

446

CE

AC

445

PT

437

counterpart.

5.2.2. Dynamic closure performance

447

In this subsection, we document the performance of the traditional Germano-Lilly dynamic model for four differ-

448

ent eddy viscosity coefficient determination strategies. The dynamically updated coefficient controls the amount of

449

artificial dissipation added to the flow features. We first test the traditional dynamic closure, where a test filtering pro-

450

cedure is used to dynamically compute the Smagorinsky or Leith coefficients controlling the eddy viscosity (D-S-TF 24

(b) Re = 8,000

(d) Re = 32,000

(e) Re = 64,000

M

AN US

(a) Re = 4,000

CR IP T

ACCEPTED MANUSCRIPT

(c) Re = 16,000

(f) Re = 128,000

451

ED

Figure 6: Final time vorticity contours for different Reynolds numbers using the Pad´e scheme with a resolution of N 2 = 20482 .

and D-L-TF).

Fig. 16 shows the performance of the closure for Re = 32,000. At this relatively low Reynolds number, the

453

coarsest grid (N 2 = 1282 ) shows no convergence for any of the closures or the under-resolved DNS (when using

454

the Pad´e scheme). At N 2 = 2562 , we see the D-S-TF method giving us a stable but poor approximation of the

455

spectrum. At the highest resolution (N 2 = 5122 ), both the D-S-TF and D-L-TF closures give acceptable results but

456

the presence of aliasing error is seen. At this resolution, we also get an under-resolved DNS result (which is not much

457

different from the D-L-TF method). The Arakawa scheme gives us converged results for each coarse-grid resolution

458

with progressively better approximations of the DNS spectrum with increasing mesh resolution. At a resolution of

460

CE

AC

459

PT

452

N 2 = 2562 , the D-S-TF formulation with the Arakawa scheme gives us the best approximation (in terms of the coarsest

grid) for this particular Reynolds number.

461

Fig. 17 shows closure performance for the Reynolds number of 64,000. For the Pad´e scheme, the coarsest reso-

462

lutions i.e., N 2 = 1282 and N 2 = 2562 , the D-S-TF and D-L-TF methods did not converge with the exception of the

463

D-S-TF approach for N 2 = 2562 which gave very inaccurate results. At N 2 = 5122 , both D-S-TF and D-L-TF results

464

converged with the former being more dissipative than the latter. The under-resolved DNS did not converge. The 25

(b) Re = 8,000

(d) Re = 32,000

(e) Re = 64,000

M

AN US

(a) Re = 4,000

CR IP T

ACCEPTED MANUSCRIPT

(c) Re = 16,000

(f) Re = 128,000

ED

Figure 7: Final time vorticity contours for different Reynolds numbers using the Arakawa scheme with a resolution of N 2 = 20482 .

465

use of the Arakawa scheme led to converged results for both under-resolved no model simulations and closures with

466

D-S-TF being the best performing (and more dissipative) closure. Similar results were also seen for the Reynolds number of 128,000, shown in Fig. 18, where the use of the Pad´e

468

scheme resulted in no convergence for the D-L-TF approach with any of the coarse grid runs. The extra dissipative

469

characteristic of the D-S-TF approach led to converged runs for the finest resolution of N 2 = 5122 even though this

470

result was poor. In contrast, the Arakawa version of the runs saw acceptable results even at N 2 = 2562 where the

471

D-S-TF version was seen to give a remarkably good approximation of the DNS spectrum. Overall, it is seen that the

472

Smagorinsky kernel for the nonlinear viscosity proves more dissipative than the Leith kernel for the D-S-TF and D-

474

CE

AC

473

PT

467

L-TF closure. However, as against the behavior seen in the standard (non-dynamic) use of the Smagorinsky and Leith models, the Leith model used in the dynamic framework leads to a marginal improvement over no model simulations.

475

The implementation of the Smagorinsky viscosity in the dynamic framework also leads to a closer agreement with the

476

DNS spectrum in comparison to its standard implementation. On closer investigation, it is apparent that the use of the

477

underlying scheme is critical for this version of the dynamic formulation. The use of the Arakawa scheme instead of

478

the Pad´e scheme transforms a non-convergent run such as the one for N 2 = 2562 to one that closely approximates the 26

CR IP T

ACCEPTED MANUSCRIPT

(d) Re = 32,000

ED

M

(c) Re = 16,000

(b) Re = 8,000

AN US

(a) Re = 4,000

(e) Re = 64,000

(f) Re = 128,000

PT

Figure 8: DNS convergence study with varying Reynolds numbers using the Pad´e scheme.

expected scaling. However, we note that higher Reynolds numbers necessitate the use of the closure approaches even

480

while using an Arakawa scheme.

CE

479

The trends in the coefficients of the Smagorinsky and Leith kernels for the aforementioned cases (at N 2 = 5122 )

482

are given in Fig. 19 for both the Pad´e and Arakawa schemes. It is evident that the Arakawa scheme causes a smoother

483

484

485

AC

481

evolution of the eddy viscosity coefficients for all Reynolds number cases as compared to the Pad´e scheme. Maximum coefficient values for the D-L-TF cases are seen to be higher than those in the D-S-TF case. Increased Reynolds numbers are seen to add more variability in the evolution of the coefficient for the Pad´e solver and a general reduction

486

in the numerical value of the coefficient in relation to the lower Reynolds numbers. Coefficients which reach a value

487

of zero and remain so till the end of the simulation indicate non-convergence (for example for Re = 128,000 in the

488

D-L-TF Pad´e formulation). We note that the variability of the dynamic Smagorinsky coefficient is between C s = 0.15

489

and C s = 0.2 for all Reynolds numbers. However, the dynamic Leith coefficient varies from Cl = 0.15 to Cl = 0.25, 27

CR IP T

ACCEPTED MANUSCRIPT

(d) Re = 32,000

ED

M

(c) Re = 16,000

(b) Re = 8,000

AN US

(a) Re = 4,000

(e) Re = 64,000

(f) Re = 128,000

PT

Figure 9: DNS convergence study with varying Reynolds numbers using the Arakawa scheme.

resulting in smaller values for the larger Re numbers. This enhanced Re number dependence for the dynamic Leith

491

model is attributed to its more scale-selective behavior.

CE

490

The D-S-B and D-L-B closure methodologies were tested for the Reynolds number of 32,000 with details shown in

493

Fig. 20. Even at this relatively low Reynolds number, it was seen that the closure approaches did not give converged

494

495

496

AC

492

results for the coarse meshes of N 2 = 1282 and N 2 = 2562 when used with the Pad´e scheme. At N 2 = 5122 , we obtained acceptable results for the closure approaches with the D-S-B model being more dissipative than the DL-B model. The Arakawa scheme, however, was able to give converged results for all coarse grids with excellent

497

approximations of the DNS at N 2 = 2562 and N 2 = 5122 . The D-S-B approach was still seen to be more dissipative

498

than the D-L-B approach.

499

For the case of Reynolds number at 64,000 shown in Fig .(21), similar trends are seen. The Pad´e solver inhibits

500

the stability of the closure approaches which only converge for the finest case of N 2 = 5122 . We note, though, that this 28

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

PT

Figure 10: Spectral comparison of the standard functional closures (i.e., the Smagorinsky and Leith models) for Re=32,000. Non-convergent runs are not shown for the Pad´e scheme.

still represents a small improvement over the no model simulation case which fails to converge for the Pad´e scheme at

502

all of the coarse resolution cases. As expected, the Arakawa scheme allows for stable and converged cases with good

503

approximations of the DNS spectrum at N 2 = 2562 and N 2 = 5122 . For the highest Reynolds number case shown in

504

Fig. 22, the closures end up non-convergent for all coarse grid cases of the Pad´e scheme. The Arakawa scheme, gives

506

AC

505

CE

501

us stable results but the use of the closure modeling approach becomes necessary to get close to the DNS spectrum. For example, at N 2 = 5122 , the use of the closures gives us a significant improvement over the no model simulation

507

results. The Pad´e numerical scheme, however, is not a preferred underlying discretization while using the D-S-B or

508

D-L-B approaches due to very low stability.

509

Fig. 23 shows the trends for the C s and Cl parameters for the D-S-B and D-L-B approach (at N 2 = 5122 ) respec-

510

tively. A behavior similar to the test-filtering approach is seen where the Arakawa scheme causes a smoother evolution 29

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

PT

Figure 11: Spectral comparison of the standard functional closures (i.e., the Smagorinsky and Leith models) for Re=64,000. Non-convergent runs are not shown for the Pad´e scheme.

of the eddy viscosity constants with the Leith coefficients reaching larger maximum values than the Smagorinsky co-

512

efficients. Higher Reynolds numbers are also seen to evolve to lower magnitudes of the coefficients in comparison to

513

lower Reynolds numbers. The early termination of the D-S-B and D-L-B Pad´e cases for Re = 128,000 in this figure

514

show the relative instability of the closure approach. We observe, that for the highest Reynolds numbers (Re = 64,000

516

AC

515

CE

511

and Re = 128,000), using the Arakawa scheme leads to higher values of the coefficients (for both D-S-B and D-L-B) compared to their Pad´e counterparts.

517

In Fig. 24, we detail the results of the D-S-L and D-L-L methods for the Reynolds number of 32,000. It can be

518

seen that the use of the Layton scale similarity procedure for reconstruction of the subgrid contribution, leads to an

519

increased stability over the previous closure approaches with converged results even for the coarsest resolution of

520

N 2 = 1282 for both D-S-L and D-L-L methods. It is observed, that the difference between the Pad´e and Arakawa 30

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

PT

Figure 12: Spectral comparison of the standard functional closures (i.e., the Smagorinsky and Leith models) for Re=128,000. Non-convergent runs are not shown for the Pad´e scheme.

schemes for this closure combination would be restricted to the aliasing error prevention alone. Another striking

522

observation is the reversal of dissipative characteristics for the Smagorinsky and Leith formulations. We expect the

523

Smagorinsky approach to be more dissipative but the D-L-L approach obtains more dissipative results than the D-S-L

524

approach.

526

AC

525

CE

521

For Re = 64,000 shown in Fig. 25, we obtain further confirmation of the stability of the method with no non-

convergent cases for the closures in any underlying scheme. The reduction of the aliasing error though necessitates

527

the use of the Arakawa scheme. The results for the highest Reynolds number of 128,000 are shown in Fig. 26, where

528

the use of the closure methods become useful to reduce aliasing error in both Pad´e and Arakawa approaches. We can

529

state with some confidence, that the D-S-L and D-L-L closures represent a significant improvement over the no model

530

approach for all the showcased Reynolds numbers while using that Pad´e scheme. For the Arakawa scheme, at high 31

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

PT

Figure 13: Spectral comparison of the standard structural closures (i.e., the Layton, Bardina and AD models) for Re=32,000. Non-convergent runs are not shown for the Pad´e scheme.

Reynolds numbers they represent a useful tool to remove aliasing error.

CE

531

This enhanced numerical stability of the Pad´e scheme can also be observed from Fig. 27 which shows the trends

533

for the C s and Cl parameters for the D-S-L and D-L-L approach respectively. In contrast to previous time evolution

534

plots for the coefficients (refer to Figs. 19 and 23), the coefficients for the highest Reynolds numbers do not reduce

535

536

537

AC

532

to values nearing zero representing an increased stability. This is true for both the D-S-L and D-L-L approach, with the Leith coefficients reaching larger maximum values than the Smagorinsky coefficients as was seen in the previous closure models.

538

The performance of the D-S-AD and D-L-AD models is shown in Fig. 28 for the Reynolds number case of 32,000.

539

The coarsest resolution of N 2 = 1282 shows converged results for the D-S-AD formulation whereas the D-L-AD

540

approach did not converge. We can thus remark that the Smagorinsky kernel is more dissipative than the Leith kernel 32

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

PT

Figure 14: Spectral comparison of the standard structural closures (i.e., the Layton, Bardina and AD models) for Re=64,000. Non-convergent runs are not shown for the Pad´e scheme.

similar to the dynamic test-filtered and dynamic Bardina approaches. The resolutions of N 2 = 2562 and N 2 = 5122

542

show that the D-S-AD approach is not just stable but extra dissipative in comparison to other closure models. The

543

Pad´e solvers still show energy accumulation near grid cut-off which may be avoided by using the Arakawa scheme.

544

However, the downward shift of the D-S-AD and D-L-AD energy spectra (due to increased dissipation) is more

546

AC

545

CE

541

pronounced using the Arakawa scheme particularly for the coarsest resolution of N 2 = 1282 . On increasing the

resolution of the LES numerical experiment, the D-S-AD and D-L-AD approximations are closer to the DNS.

547

Fig. 29 shows the performance of these closures for Re = 64,000. Once more, the coarsest resolution of N 2 = 1282

548

combined with the Pad´e scheme shows converged results for the D-S-AD approach indicating an increased stability as

549

compared to the D-S-TF and D-S-B approaches (refer to Figs. 17 and 21). The increased dissipation of the D-S-AD

550

approach also proves promising for removing noise from converged results using both Pad´e and Arakawa schemes. A 33

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

PT

Figure 15: Spectral comparison of the standard structural closures (i.e., the Layton, Bardina and AD models) for Re=128,000. Non-convergent runs are not shown for the Pad´e scheme.

similar behavior is seen in Fig. 30 which shows the case for Re = 128,000. The coarsest resolution (N 2 = 1282 ) still

552

shows converged results for the D-S-AD approach with N 2 = 2562 onwards showing convergence for both D-S-AD

553

and D-L-AD approaches. The Arakawa scheme is efficient in preventing preventing accumulation at grid cut-off.

555

556

The use of closures adds a considerable improvement to under-resolved DNS simulations for both the Pad´e and

AC

554

CE

551

Arakawa schemes. Fig. 31 shows the trends for the C s and Cl parameters for the D-S-AD and D-L-AD approach respectively. The stability of the method at N 2 = 5122 is evident from the observation that the evolution of eddy

557

viscosity coefficients do not drop down to zero for any Reynolds number. As expected, the values of the Leith

558

coefficients are higher than the Smagorinsky coefficients.

559

To perform a quantitative analysis we define the following relationship measuring error based on the angle aver-

34

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

the Pad´e scheme.

aged energy spectrum:

561

562

AC

CE

560

PT

Figure 16: Spectral comparison of the dynamic test filtered closures (D-S-TF and D-L-TF) for Re=32,000. Non-convergent runs are not shown for

=

km P

k=k p

|E MODEL (k) − E DNS (k)| k3 km − k p

(89)

where we apply an appropriate scaling by noting that E(k) ∼ k−3 . Here, E MODEL (k) and E DNS (k) refer respectively to the energy spectrum computed by the coarse-grained LES and the high-fidelity DNS computations. The error

564

between the prescribed model and the DNS data is computed over the spectral range between the k p = 10 and √ km = N/ 2. Results are shown in tabular form in Table 1. We also include the performance for the no-model

565

simulations to clearly demonstrate the improvement of the various closures implemented. Values are presented for

566

the highest Reynolds number of 128,000. It can be seen that an increasing resolution consistently leads to a reduction

567

in error thereby validating our approach for quantification. We also note that the Leith kernel performs in a more

563

35

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

the Pad´e scheme.

PT

Figure 17: Spectral comparison of the dynamic test filtered closures (D-S-TF and D-L-TF) for Re=64,000. Non-convergent runs are not shown for

dissipative manner when combined with the dynamic formulation (D-L-L) as compared to the same formulation

569

utilizing that Smagorinsky kernel (D-S-L). It can also be observed that the Arakawa formulation provides a much

570

more accurate capture of the energy spectra through a reduced error between model and DNS as compared to the

571

compact scheme. Our numerical findings suggest that the dynamic AD models (D-S-AD and D-L-AD) yield smaller

573

AC

572

CE

568

errors (i.e., comparing against to DNS) for coarser resolutions whereas the dynamic Layton models (D-S-L and D-LL) minimize errors for finer resolutions.

574

Another measure of the dissipative characteristics of the various closures described above is the inspection of their

575

final time (t = 10 s) vorticity contours for a particular Reynolds number. For this purpose we compare the results

576

obtained from the DNS, under-resolved DNS and dynamic LES approaches at a particular resolution and Reynolds

577

number to make conclusions about the dissipation characteristics of the method employed for closure. Firstly, we use 36

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

for the Pad´e scheme.

PT

Figure 18: Spectral comparison of the dynamic test filtered closures (D-S-TF and D-L-TF) for Re=128,000. Non-convergent runs are not shown

a resolution of N 2 = 5122 for a Reynolds number of 32,000 (for both the Pad´e and Arakawa schemes) for comparison

579

as shown in Figs. 32, 33, 34 and 35. It can be seen that the fields obtained by the closures are more capable of

580

smoothing the vorticity contours in comparison to the UDNS results.

582

583

Similarly, we detail the vorticity contours obtained at the final time for simulations with the highest Reynolds

AC

581

CE

578

number (Re = 128,000) using the Arakawa scheme. We note that all the dynamic LES closures did not yield converged results for this value of the Reynolds number when combined with the Pad´e scheme even for the case with N 2 = 5122 .

584

We first show the performance of the dynamic LES closures combined with the Arakawa discretization for the coarse

585

grid given by N 2 = 1282 in Figs. 36 and 37. Similarly, coarse grid runs with N 2 = 2562 and N 2 = 5122 are shown

586

in Figs. 38, 39, 40 and 41. The smoothing ability of the D-S-AD and D-L-AD approaches is apparent as they are

587

seen to produce vorticity contours with the least amount of noise. This can be attributed to their extra dissipative 37

CR IP T

ACCEPTED MANUSCRIPT

(c) Pad´e (D-L-TF)

(b) Arakawa (D-S-TF)

AN US

(a) Pad´e (D-S-TF)

(d) Arakawa (D-L-TF)

Figure 19: Trends in the eddy viscosity coefficient for the D-S-TF (top) and D-L-TF (bottom) models using the underlying Pad´e and Arakawa schemes.

characteristics.

M

588

On a global examination of the results shown in this section, it is apparent that the choice of numerical scheme

590

is vital to the fidelity of the under-resolved no model simulation. It may even be argued that the choice of a subgrid

591

closure becomes irrelevant when the Arakawa scheme is used for the underlying numerical discretization. However,

592

the benefits of the subgrid modeling approach are apparent when using the Pad´e scheme. We highlight that all

593

LES results are also superior (i.e., showing less grid-to-grid oscillations on under-resolved grids) than the UDNS

594

results even if we use the energy conserving Arakawa scheme due to the extra stabilization coming from additional

595

dissipation mechanism. We note that the simulations with higher Reynolds numbers are generally expected to require

596

some sort of subgrid approximation when used with the Arakawa scheme if one desires a excellent representation of

597

the DNS spectrum. When such constraints of accuracy are not a concern, the requirement of a subgrid closure may

598

be superfluous. However, lower Reynolds number simulations indicate that the lone use of the energy conserving

600

601

PT

CE

AC

599

ED

589

Arakawa scheme gives us well-resolved results with coarse grids. In this regard we tend to concur with the conclusion reached in Bosshard et al. [94]. 6. Conclusion

602

To conclude, this investigation evaluated the spectral performance of various closures models including a novel

603

hybrid closure methodology based on the use of both functional and structural elements in the reconstruction of 38

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

PT

Figure 20: Spectral comparison of the dynamic Bardina scale similarity closures (D-S-B and D-L-B) for Re=32,000. Non-convergent runs are not

subgrid scale contributions to the overall flow features. Two different underlying numerical schemes were used for

605

the discretization of the nonlinear Jacobian to quantify their effect on the spectral scaling. The performance of the

606

numerical code was validated against the theoretical scaling expected for a 2D isotropic turbulence case according to

607

KBL theory. For the purpose of comparison, high fidelity DNS and filtered DNS data were used as reference for the

609

AC

608

CE

604

2D decaying turbulence problem in a square domain with periodic boundary conditions. It was seen that a majority of the cases showed an improvement when the derived LES closures were implemented as against a purely under-

610

resolved DNS case. The Arakawa scheme, however, was critical in the prevention of oscillations and aliasing error.

611

We attribute this to the energy conserving property of the Arakawa scheme.

612

Among the different closures tested, it was seen that the dynamic Bardina (D-S-B and D-L-B) approaches were

613

the least dissipative as compared to the other approaches whereas the dynamic AD approach (D-S-AD and D-L39

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

PT

Figure 21: Spectral comparison of the dynamic Bardina scale similarity closures (D-S-B and D-L-B) for Re=64,000. Non-convergent runs are not

AD) were most dissipative. The dynamic Layton methods (D-S-L and D-L-L) were seen to be the most stable LES

615

closure with converged results for all coarse resolutions and Reynolds numbers. From the point-of-view of scaling

616

law capture, the use of the Arakawa scheme can be considered to be unavoidable and its combination with the D-

617

S-L and D-L-L approaches represents a feasible approach to reduce the computational cost of the DNS approach to

619

AC

618

CE

614

two-dimensional turbulence. Results also indicate that the Leith and Smagorinsky kernels are similar to each other (when judged for their ability to capture the KBL scaling), with the slight differences in dissipative behavior when

620

combined with different closure methodologies. An important remark here is the reversal of dissipative behavior when

621

utilizing the dynamic Layton closure methodology where the Leith kernel proved to be more dissipative instead of

622

the Smagorinsky kernel (the conventional behavior of the kernels when used with other closure approaches). It must

623

be noted here that the under-dissipative nature of the Leith model when not combined in the D-L-L approach has yet 40

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

PT

Figure 22: Spectral comparison of the dynamic Bardina scale similarity closures (D-S-B and D-L-B) for Re=128,000. Non-convergent runs are not shown for the Pad´e scheme.

to be explained suitably from a physical or numerical point of view and is an important question that remains to be

625

answered.

CE

624

The use of the dynamic estimation of nonlinear viscosity coefficients through the use of structural closures (par-

627

ticularly for the Layton scale similarity approach) indicate a useful approach for subgrid scale parameterization and

628

629

AC

626

may prove useful for a wide variety of 2D turbulence applications. However, we would like to conclude by remarking that the underlying discretization plays a vital role in the feasibility of a LES in that its effect on accuracy and stability

630

far exceed that of any of the subgrid closures displayed here. We believe that the role of subgrid modeling becomes

631

superfluous for low Reynolds number simulations where an intelligent numerical discretization could reduce the grid

632

resolution required for accurate spectral scaling capture. Although each subgrid scale model utilizes a different dissi-

633

pation mechanism, we found that the choice of numerical scheme and its characteristics would seem more dominant 41

CR IP T

ACCEPTED MANUSCRIPT

(c) Pad´e (D-L-B)

(b) Arakawa (D-S-B)

AN US

(a) Pad´e (D-S-B)

(d) Arakawa (D-L-B)

Figure 23: Trends in the eddy viscosity coefficient for the D-S-B (top) and D-L-B (bottom) models using the underlying Pad´e and Arakawa schemes.

than the choice of subgrid scale model. Our findings follow the recent conclusions of Bosshard et al. [94] advocating

635

the use of improved numerical algorithms.

636

Acknowledgments

ED

The computing for this project was performed at the High Performance Computing Center at Oklahoma State

PT

637

M

634

University (OSUHPCC).

639

References

CE

638

[1] G. Boffetta, R. E. Ecke, Two-dimensional turbulence, Annu. Rev. Fluid Mech. 44 (2012) 427–451. [2] G. Boffetta, S. Musacchio, Evidence for the double cascade scenario in two-dimensional turbulence, Phys. Rev. E 82 (1) (2010) 016307.

642 643 644 645

AC

640 641

[3] R. H. Kraichnan, D. Montgomery, Two-dimensional turbulence, Rep. Prog. Phys. 43 (5) (1980) 547. [4] F. Bouchet, A. Venaille, Statistical mechanics of two-dimensional and geophysical flows, Phys. Rep. 515 (5) (2012) 227–295. [5] S. D. Danilov, D. Gurarie, Quasi-two-dimensional turbulence, Phys.-Usp. 43 (9) (2000) 863–900. [6] R. H. Kraichnan, Inertial ranges in two-dimensional turbulence, Phys. Fluids 10 (1967) 1417–1423.

646

[7] G. K. Batchelor, Computation of the Energy spectrum in Homogeneous Two-Dimensional Turbulence, Phys. Fluids 12 (12) (1969) II–233.

647

[8] C. Leith, Atmospheric predictability and two-dimensional turbulence, J. Atmos. Sci. 28 (2) (1971) 145–161.

648

[9] H. Clercx, G. van Heijst, Two-dimensional Navier-Stokes turbulence in bounded domains, Appl. Mech. Rev. 62 (2) (2009) 020802.

649

[10] S. Danilov, D. Gurarie, Forced two-dimensional turbulence in spectral and physical space, Phys. Rev. E 63 (6) (2001) 061208.

42

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

651

[11] A. Bracco, J. C. McWilliams, Reynolds-number dependency in homogeneous, stationary two-dimensional turbulence, J. Fluid Mech. 646

CE

650

PT

Figure 24: Spectral comparison of the dynamic Layton scale similarity closures (D-S-L and D-L-L) for Re=32,000. Non-convergent runs are not

(2010) 517.

[12] A. Vallgren, E. Lindborg, The enstrophy cascade in forced two-dimensional turbulence, J. Fluid Mech. 671 (2011) 168–183.

653

[13] E. Lindborg, A. Vallgren, Testing batchelors similarity hypotheses for decaying two-dimensional turbulence, Phys. Fluids 22 (9) (2010)

654 655 656 657

AC

652

091704.

[14] E. Kuznetsov, V. Naulin, A. H. Nielsen, J. J. Rasmussen, Sharp vorticity gradients in two-dimensional turbulence and the energy spectrum, Theor. Comp. Fluid Dyn. 24 (1-4) (2010) 253–258.

[15] J. P. Boris, F. F. Grinstein, E. S. Oran, R. L. Kolbe, New insights into large eddy simulation, Fluid Dyn. Res. 10 (4-6) (1992) 199–228.

658

[16] M. Lesieur, O. Metais, New trends in large-eddy simulations of turbulence, Annu. Rev. Fluid Mech. 28 (1) (1996) 45–82.

659

[17] U. Piomelli, Large-eddy simulation: achievements and challenges, Prog. Aerosp. Sci. 35 (4) (1999) 335–362.

660

[18] C. Meneveau, J. Katz, Scale-invariance and turbulence models for large-eddy simulation, Annu. Rev. Fluid Mech. 32 (1) (2000) 1–32.

661

[19] P. Sagaut, Large eddy simulation for incompressible flows: an introduction, Springer-Verlag, New York, 2006.

662

[20] L. C. Berselli, T. Iliescu, W. J. Layton, Mathematics of large eddy simulation of turbulent flows, Springer-Verlag, New York, 2006.

43

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

665 666 667 668 669 670 671 672 673

99–164.

CE

664

[21] J. Smagorinsky, General circulation experiments with the primitive equations. I. The basic experiments, Mon. Weather Rev. 91 (3) (1963) [22] A. Yakhot, S. A. Orszag, V. Yakhot, M. Israeli, Renormalization group formulation of large-eddy simulations, J. Sci. Comput. 4 (2) (1989) 139–158.

AC

663

PT

Figure 25: Spectral comparison of the dynamic Layton scale similarity closures (D-S-L and D-L-L) for Re=64,000. Non-convergent runs are not

[23] A. Yoshizawa, Subgrid-scale modeling with a variable length scale, Phys. Fluids A 1 (1989) 1293–1295. [24] M. Germano, U. Piomelli, P. Moin, W. H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids 3 (1991) 1760–1765. [25] U. Piomelli, W. H. Cabot, P. Moin, S. Lee, Subgrid-scale backscatter in turbulent and transitional flows, Phys. Fluids 3 (1991) 1766–1771. [26] T. J. R. Hughes, L. Mazzei, A. A. Oberai, A. A. Wray, The multiscale formulation of large eddy simulation: Decay of homogeneous isotropic turbulence, Phys. Fluids 13 (2001) 505–512. [27] G. S. Winckelmans, A. A. Wray, O. V. Vasilyev, H. Jeanmart, Explicit-filtering large-eddy simulation using the tensor-diffusivity model supplemented by a dynamic Smagorinsky term, Phys. Fluids 13 (2001) 1385–1403.

674

[28] B. J. Geurts, D. D. Holm, Leray and LANS-α modelling of turbulent mixing, J. Turbul. 7 (2006) 1–33.

675

[29] U. Frisch, Turbulence, Cambridge University Press, Cambridge, U.K., 1995.

44

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

PT

Figure 26: Spectral comparison of the dynamic Layton scale similarity closures (D-S-L and D-L-L) for Re=128,000. Non-convergent runs are not

[30] C. Leith, Stochastic models of chaotic systems, Physica D: Nonlinear Phenomena 98 (2) (1996) 481–491.

677

[31] C. E. Leith, Diffusion approximation for two-dimensional turbulence, Physics of Fluids 11 (1968) 671–672.

678

[32] B. Fox-Kemper, D. Menemenlis, Can large eddy simulation techniques improve mesoscale rich ocean models?, Ocean modeling in an eddying

680 681 682 683 684 685

regime (2008) 319–337.

AC

679

CE

676

[33] D. K. Lilly, A proposed modification of the Germano subgrid-scale closure method, Phys. Fluids 4 (1992) 633–635. [34] J. Smagorinsky, Some historical remarks on the use of nonlinear viscosities, in: Large eddy simulation of complex engineering and geophysical flows, eds: B. Galperin and S. A. Orszag, Cambridge University Press, New York, 1993, pp. 3–36.

[35] V. M. Canuto, Y. Cheng, Determination of the Smagorinsky-Lilly constant CS, Phys. of Fluids 9 (1997) 1368–1378. [36] A. Vorobev, O. Zikanov, Smagorinsky constant in LES modeling of anisotropic MHD turbulence, Theor. Comp. Fluid Dyn. 22 (3-4) (2008) 317–325.

686

[37] S. B. Pope, Turbulent flows, Cambridge University Press, New York, 2000.

687

[38] B. Cushman-Roisin, J.-M. Beckers, Introduction to geophysical fluid dynamics: physical and numerical aspects, Vol. 101, Academic Press,

688

2011.

45

CR IP T

ACCEPTED MANUSCRIPT

(c) Pad´e (D-L-L)

(b) Arakawa (D-S-L)

AN US

(a) Pad´e (D-S-L)

(d) Arakawa (D-L-L)

Figure 27: Trends in the eddy viscosity coefficient for the D-S-L (top) and D-L-L (bottom) models using the underlying Pad´e and Arakawa schemes.

691 692 693 694

large-eddy simulation using wrf*, Mon. Weather Rev. 140 (1) (2012) 266–284.

M

690

[39] G. Kirkil, J. Mirocha, E. Bou-Zeid, F. K. Chow, B. Kosovic, Implementation and evaluation of dynamic subfilter-scale stress models for [40] J. O’Neill, X.-M. Cai, R. Kinnersley, A generalised stochastic backscatter model: large-eddy simulation of the neutral surface layer, Q. J. Roy. Meteor. Soc. 141 (692) (2015) 2617–2629.

[41] R. Maulik, O. San, Dynamic modeling of the horizontal eddy viscosity coefficient for quasigeostrophic ocean circulation problems, J. Ocean

ED

689

Eng. Sci. (2016) http://dx.doi.org/10.1016/j.joes.2016.08.002. [42] S. Stolz, N. A. Adams, An approximate deconvolution procedure for large-eddy simulation, Phys. Fluids 11 (1999) 1699–1701.

696

[43] O. San, A. E. Staples, Z. Wang, T. Iliescu, Approximate deconvolution large eddy simulation of a barotropic ocean circulation model, Ocean

698 699

Modell. 40 (2) (2011) 120–132.

[44] O. San, A. E. Staples, T. Iliescu, Approximate deconvolution large eddy simulation of a stratified two-layer quasigeostrophic ocean model, Ocean Modell. 63 (2013) 1–20.

CE

697

PT

695

700

[45] M. Germano, A new deconvolution method for large eddy simulation, Phys. Fluids 21 (2009) 045107.

701

[46] W. J. Layton, L. Rebholz, Approximate deconvolution models of turbulence: analysis, phenomenology and numerical analysis, Springer-

703 704 705 706

Verlag, Berlin, 2012.

AC

702

[47] S. Stolz, N. A. Adams, L. Kleiser, An approximate deconvolution model for large-eddy simulation with application to incompressible wallbounded flows, Phys. Fluids 13 (2001) 997–1015.

[48] S. Stolz, N. A. Adams, L. Kleiser, The approximate deconvolution model for large-eddy simulations of compressible flows and its application to shock-turbulent-boundary-layer interaction, Phys. Fluids 13 (2001) 2985–3001.

707

[49] J. A. Domaradzki, N. A. Adams, Direct modelling of subgrid scales of turbulence in large eddy simulations, J. Turbul. 3 (24) (2002) 1–19.

708

[50] F. K. Chow, R. L. Street, M. Xue, J. H. Ferziger, Explicit filtering and reconstruction turbulence modeling for large-eddy simulation of neutral

709 710 711

boundary layer flow, J. Atmos. Sci. 62 (7) (2005) 2058–2077. [51] F. K. Chow, R. L. Street, Evaluation of turbulence closure models for large-eddy simulation over complex terrain: flow over Askervein Hill, J. Appl. Meteorol. Clim. 48 (5) (2009) 1050–1065.

46

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

715 716 717 718 719

CE

713 714

¨ okmen, Bridging the Boussinesq and primitive equations through spatio-temporal filtering, Appl. [52] J. Duan, P. Fischer, T. Iliescu, T. M. Ozg¨ Math. Lett. 23 (4) (2010) 453–456. [53] B. Zhou, F. K. Chow, Large-eddy simulation of the stable boundary layer with explicit filtering and reconstruction turbulence modeling, J. Atmos. Sci. 68 (2011) 2142–2155.

AC

712

PT

Figure 28: Spectral comparison of the dynamic AD scale similarity closures (D-S-AD and D-L-AD) for Re=32,000. Non-convergent runs are not

[54] A. A. Dunca, Y. Epshteyn, On the Stolz–Adams deconvolution model for the large-eddy simulation of turbulent flows, SIAM J. Math. Anal. 37 (6) (2006) 1890–1902.

[55] W. Layton, R. Lewandowski, Residual stress of approximate deconvolution models of turbulence, J. Turbul. 7 (2006) 1–21. [56] W. Layton, M. Neda, A similarity theory of approximate deconvolution models of turbulence, J. Math. Anal. Appl. 333 (1) (2007) 416–429.

720

[57] L. G. Rebholz, Conservation laws of turbulence models, J. Math. Anal. Appl. 326 (1) (2007) 33–45.

721

[58] I. Stanculescu, Existence theory of abstract approximate deconvolution models of turbulence, Ann. Univ. Ferrara 54 (1) (2008) 145–168.

722

[59] A. A. Dunca, On the existence of global attractors of the approximate deconvolution models of turbulence, J. Math. Anal. Appl. 389 (2)

723 724

(2011) 1128–1138. [60] L. C. Berselli, R. Lewandowski, Convergence of approximate deconvolution models to the mean Navier-Stokes equations, Ann. I. H. Poincar´e

47

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

- AN 29 (2) (2012) 171–198.

CE

725

PT

Figure 29: Spectral comparison of the dynamic AD scale similarity closures (D-S-AD and D-L-AD) for Re=64,000. Non-convergent runs are not

726

[61] A. A. Dunca, R. Lewandowski, Error estimates in approximate deconvolution models, Commun. Math. Sci. 12 (4) (2014) 757–778.

727

[62] W. Layton, R. Lewandowski, A simple and stable scale-similarity model for large eddy simulation: energy balance and existence of weak

729 730 731 732 733 734

solutions, Applied mathematics letters 16 (8) (2003) 1205–1209.

AC

728

[63] J. Bardina, J. H. Ferziger, W. Reynolds, Improved subgrid-scale models for large-eddy simulation, in: American Institute of Aeronautics and Astronautics, Fluid and Plasma Dynamics Conference, 13th, Snowmass, Colo., July 14-16, 1980, 10 p., Vol. 1, 1980.

[64] M. A. Habisreutinger, R. Bouffanais, E. Leriche, M. O. Deville, A coupled approximate deconvolution and dynamic mixed scale model for large-eddy simulation, J. Comput. Phys. 224 (1) (2007) 241–266.

[65] J. Mathew, R. Lechner, H. Foysi, J. Sesterhenn, R. Friedrich, An explicit filtering method for large eddy simulation of compressible flows, Phys. Fluids 15 (8) (2003) 2279–2289.

735

[66] J. Mathew, H. Foysi, R. Friedrich, A new approach to LES based on explicit filtering, Int. J. Heat Fluid Fl. 27 (4) (2006) 594–602.

736

[67] T. Brandt, A priori tests on numerical errors in large eddy simulation using finite differences and explicit filtering, Int. J. Numer. Meth. Fluids

737

51 (6) (2006) 635–657.

48

CR IP T

ACCEPTED MANUSCRIPT

(d) Arakawa (N 2 = 2562 )

ED

M

(c) Pad´e (N 2 = 2562 )

(b) Arakawa (N 2 = 1282 )

AN US

(a) Pad´e (N 2 = 1282 )

(e) Pad´e (N 2 = 5122 )

(f) Arakawa (N 2 = 5122 )

shown for the Pad´e scheme.

741 742 743 744 745 746

CE

739 740

[68] J. Berland, P. Lafon, F. Daude, F. Crouzet, C. Bogey, C. Bailly, Filter shape dependence and effective scale separation in large-eddy simulations based on relaxation filtering, Comput. Fluids 47 (1) (2011) 65–74. [69] O. San, A. E. Staples, T. Iliescu, A posteriori analysis of low-pass spatial filters for approximate deconvolution large eddy simulations of homogeneous incompressible flows, Int. J. Comput. Fluid D. 29 (2015) 40–66.

AC

738

PT

Figure 30: Spectral comparison of the dynamic AD scale similarity closures (D-S-AD and D-L-AD) for Re=128,000. Non-convergent runs are not

[70] O. San, Analysis of low-pass filters for approximate deconvolution closure modeling in one-dimensional decaying Burgers turbulence, Int. J. Comput. Fluid D. 30 (2016) 20–37.

[71] M. Germano, Differential filters for the large eddy numerical simulation of turbulent flows, Phys. Fluids 29 (1986) 1755–1756. [72] A. Najafi-Yazdi, M. Najafi-Yazdi, L. Mongeau, A high resolution differential filter for large eddy simulation: Toward explicit filtering on unstructured grids, J. Comput. Phys. 292 (2015) 272–286.

747

[73] S. K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1) (1992) 16–42.

748

[74] M. R. Visbal, D. V. Gaitonde, On the use of higher-order finite-difference schemes on curvilinear and deforming meshes, J. Comput. Phys.

749 750

181 (1) (2002) 155–185. [75] C. Bogey, C. Bailly, A family of low dispersive and low dissipative explicit schemes for flow and noise computations, J. Comput. Phys.

49

CR IP T

ACCEPTED MANUSCRIPT

(c) Pad´e (D-L-AD)

(b) Arakawa (D-S-AD)

AN US

(a) Pad´e (D-S-AD)

(d) Arakawa (D-L-AD)

Figure 31: Trends in the eddy viscosity coefficient for D-S-AD (top) and D-L-AD (bottom) models using the underlying Pad´e and Arakawa

756 757 758 759 760 761

dissipation, Phys. Fluids 18 (6) (2006) 065101.

[77] C. Bogey, C. Bailly, Computation of a high Reynolds number jet and its radiated noise using large eddy simulation based on explicit filtering,

ED

754 755

[76] C. Bogey, C. Bailly, Large eddy simulations of transitional round jets: influence of the Reynolds number on flow development and energy

Comput. Fluids 35 (10) (2006) 1344–1358.

[78] D. L. Brown, R. Cortez, M. L. Minion, Accurate projection methods for the incompressible Navier-Stokes equations, J. Comput. Phys. 168 (2) (2001) 464–499.

PT

752 753

194 (1) (2004) 194–214.

[79] J. R. Mansfield, O. M. Knio, C. Meneveau, A dynamic LES scheme for the vorticity transport equation: formulation and a priori tests, J. Comput. Phys. 145 (2) (1998) 693–730.

[80] A. Arakawa, Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible

CE

751

M

schemes.

flow, J. Comput. Phys. 1 (1) (1966) 119–143. [81] P. Orlandi, Fluid flow phenomena: a numerical toolkit, Vol. 55, Springer Science & Business Media, 2012.

763

[82] S. Gottlieb, C. W. Shu, Total variation diminishing Runge-Kutta schemes, Math. Comput. 67 (221) (1998) 73–85.

764 765 766 767

AC

762

[83] O. San, A. E. Staples, High-order methods for decaying two-dimensional homogeneous isotropic turbulence, Comput. Fluids 63 (2012) 105–127.

[84] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (2) (1988) 439– 471.

768

[85] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev. 43 (1) (2001) 89–112.

769

[86] X. Hu, Q. Wang, N. A. Adams, An adaptive central-upwind weighted essentially non-oscillatory scheme, J. Comput. Phys. 229 (23) (2010)

770 771 772

8952–8965. [87] K. Kara, P. Balakumar, O. A. Kandil, Effects of nose bluntness on hypersonic boundary-layer receptivity and stability over cones, AIAA J. 49 (12) (2011) 2593–2606.

50

ACCEPTED MANUSCRIPT

Table 1: Comparison of errors for different closures. DNS at N 2 = 20482 used as true solution.

256 × 256

512 × 512 Pad´e

Arakawa

Pad´e

Arakawa

Pad´e

Arakawa

UDNS

NA

90.17

NA

31.43

NA

2.38

Smagorinsky

2.19

2.24

1.24

1.347

0.61

0.69

Leith

NA

86.51

NA

31.69

NA

2.39

Layton

26.23

17.58

11.1

8.83

1.62

1.64

Bardina

NA

20.56

NA

6.99

NA

0.33

AD

214.23

9.56

10.62

4.39

0.76

0.58

D-S-TF

NA

5.12

NA

0.74

87.79

0.38

D-L-TF

NA

49.95

NA

14.11

NA

0.98

D-S-BD

NA

12.51

NA

2.28

NA

0.439

D-L-BD

NA

28.09

NA

7.54

NA

0.312

D-S-L

119.87

19.59

23.46

5.36

0.68

0.27

D-L-L

15.79

13.02

2.96

2.44

0.56

0.38

D-S-AD

23.79

2.19

1.22

1.48

0.78

0.82

D-L-AD

NA

1.36

18.87

1.31

0.61

0.67

AN US

CR IP T

Model

[88] Y. Wang, J. Zhang, Sixth order compact scheme combined with multigrid method and extrapolation technique for 2d poisson equation, J. Comput. Phys. 228 (1) (2009) 137–146.

M

773 774

128 × 128

[89] L. Collatz, The numerical treatment of differential equations, Vol. 60, Springer Science & Business Media, 2012.

776

[90] I. Tsukerman, A class of difference schemes with flexible local approximation, J. Comput. Phys. 211 (2) (2006) 659–699.

777

[91] O. San, A novel high-order accurate compact stencil Poisson solver: application to cavity flows, Int. J. Appl. Mech. 7 (01) (2015) 1550006.

778

[92] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical recipes in FORTRAN: The art of scientific computing, Cambridge University

779

Press, New York, 1992.

ED

775

[93] B. Jahne, Digital image processing: Concepts, algorithms, and scientific applications, Springer-Verlag, Berlin, Heildelberg, 1997.

781

[94] C. Bosshard, M. O. Deville, A. Dehbi, E. Leriche, UDNS or LES, that is the question, Open J. Fluid D. 5 (04) (2015) 339–352.

AC

CE

PT

780

51

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-S-TF

PT

ED

M

(a) DNS

CE

(d) D-S-B

(e) D-S-L

(f) D-S-AD

Figure 32: Vorticity contours at t = 10 s and Reynolds number of 32,000 using the Pad´e scheme. Dynamic LES closures (at N 2 = 5122 ) with the

AC

Smagorinsky kernel given by Eq. (13).

52

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-L-TF

PT

ED

M

(a) DNS

CE

(d) D-L-B

(e) D-L-L

(f) D-L-AD

Figure 33: Vorticity contours at t = 10 s and Reynolds number of 32,000 using the Pad´e scheme. Dynamic LES closures (at N 2 = 5122 ) with the

AC

Leith kernel given by Eq. (14).

53

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-S-TF

PT

ED

M

(a) DNS

CE

(d) D-S-B

(e) D-S-L

(f) D-S-AD

Figure 34: Vorticity contours at t = 10 s and Reynolds number of 32,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 5122 ) with

AC

the Smagorinsky kernel given by Eq. (13).

54

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-L-TF

PT

ED

M

(a) DNS

(e) D-L-L

CE

(d) D-L-B

(f) D-L-AD

Figure 35: Vorticity contours at t = 10 s and Reynolds number of 32,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 5122 ) with

AC

the Leith kernel given by Eq. (14).

55

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-S-TF

PT

ED

M

(a) DNS

CE

(d) D-S-B

(e) D-S-L

(f) D-S-AD

Figure 36: Vorticity contours at t = 10 s and Reynolds number of 128,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 1282 ) with

AC

the Smagorinsky kernel given by Eq. (13).

56

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-L-TF

PT

ED

M

(a) DNS

CE

(d) D-L-B

(e) D-L-L

(f) D-L-AD

Figure 37: Vorticity contours at t = 10 s and Reynolds number of 128,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 1282 ) with

AC

the Leith kernel given by Eq. (14).

57

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-S-TF

PT

ED

M

(a) DNS

CE

(d) D-S-B

(e) D-S-L

(f) D-S-AD

Figure 38: Vorticity contours at t = 10 s and Reynolds number of 128,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 2562 ) with

AC

the Smagorinsky kernel given by Eq. (13).

58

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-L-TF

PT

ED

M

(a) DNS

CE

(d) D-L-B

(e) D-L-L

(f) D-L-AD

Figure 39: Vorticity contours at t = 10 s and Reynolds number of 128,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 2562 ) with

AC

the Leith kernel given by Eq. (14).

59

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-S-TF

PT

ED

M

(a) DNS

CE

(d) D-S-B

(e) D-S-L

(f) D-S-AD

Figure 40: Vorticity contours at t = 10 s and Reynolds number of 128,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 5122 ) with

AC

the Smagorinsky kernel given by Eq. (13).

60

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b) UDNS

(c) D-L-TF

PT

ED

M

(a) DNS

CE

(d) D-L-B

(e) D-L-L

(f) D-L-AD

Figure 41: Vorticity contours at t = 10 s and Reynolds number of 128,000 using the Arakawa scheme. Dynamic LES closures (at N 2 = 5122 ) with

AC

the Leith kernel given by Eq. (14).

61