University of Southern California Dissertations and Theses

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

(USC Thesis Other)

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content LARGE EDDY SIMULATIONS OF TURBULENT FLOWS WITHOUT USE OF THE EDDY VISCOSITY CONCEPT by Guangrui Sun A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (AEROSPACE ENGINEERING) December 2019 Copyright 2019 Guangrui Sun Contents List of Tables iv List of Figures v Abstract ix 1 Introduction 1 2 Background 9 2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Large Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.1 Explicit subgrid-scale Models . . . . . . . . . . . . . . . . . 13 2.3.2 Implicit subgrid-scale Models . . . . . . . . . . . . . . . . . 15 3 Numerical Method 18 3.1 Skew-symmetric form of Navier-Stokes equation . . . . . . . . . . . 18 3.2 Channel flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.3 Temporal discretization of Navier-Stokes equation . . . . . . . . . . 20 3.3.1 Step I: nonlinear and viscous . . . . . . . . . . . . . . . . . . 20 3.3.2 Step II: pressure . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4 Quantification of numerical dissipation . . . . . . . . . . . . . . . . 23 3.5 Filtering on nonuniform mesh . . . . . . . . . . . . . . . . . . . . . 25 3.6 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4 Use filtering as an implicit LES method 30 4.1 General idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Adaptive filtering criterion . . . . . . . . . . . . . . . . . . . . . . . 36 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.3.1 Numerical dissipation introduced by explicit filters . . . . . 40 4.3.2 One-dimensional energy spectrum . . . . . . . . . . . . . . . 46 4.3.3 Different Reynolds number . . . . . . . . . . . . . . . . . . . 50 4.3.4 Performance for coarse grid resolutions . . . . . . . . . . . . 56 ii 4.3.5 Comparison of different filters . . . . . . . . . . . . . . . . . 60 4.3.6 Investigation of lower order schemes . . . . . . . . . . . . . . 65 5 Explicit interscale model 79 5.1 General idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.1.1 Scale separation by filtering . . . . . . . . . . . . . . . . . . 80 5.1.2 Derivation of the SGS stress . . . . . . . . . . . . . . . . . . 82 5.1.3 Modification to ensure Galilean invariance . . . . . . . . . . 87 5.2 Relation to other models . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2.1 Relation to Leray’s regularization . . . . . . . . . . . . . . . 92 5.2.2 Relation to Navier-Stokes-alpha regularization . . . . . . . . 97 5.2.3 Relation to Anderson and Domaradzki’s interscale energy transfer model . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.4 Relation to SFS stress and similarity model . . . . . . . . . 105 5.3 A posteriori tests for turbulent channel flow . . . . . . . . . . . . . 107 5.3.1 Different Reynolds number . . . . . . . . . . . . . . . . . . . 108 5.3.2 Different grid resolution . . . . . . . . . . . . . . . . . . . . 120 5.3.3 Energy spectra . . . . . . . . . . . . . . . . . . . . . . . . . 124 6 Conclusions 127 Reference List 131 iii List of Tables 3.1 DNSandexperimentaldatacomparedinLee&Moser(2015)(Table 1 therein). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.1 Filtered energy ratios for different energy spectrum. . . . . . . . . . 35 4.2 Current LES parameters. . . . . . . . . . . . . . . . . . . . . . . . . 41 4.3 Benchmark simulation parameters. . . . . . . . . . . . . . . . . . . 41 4.4 Relative error and numerical dissipation for Re τ = 590. . . . . . . . 41 4.5 Relative error and numerical dissipation for different filtering intervals 43 4.6 Effect of different filters for Re τ = 180, Mesh: 64× 65× 64. . . . . 60 4.7 Effect of different filters for Re τ = 590, Mesh: 96× 97× 96. . . . . 63 4.8 Different orders of accuracy for Re τ = 180, no filter case, Mesh: 64× 65× 64. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.9 Different orders of accuracy for Re τ = 180, filtered case, Mesh: 64× 65× 64. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.10 Different orders of accuracy for Re τ = 1000, no filter case, Mesh: 96× 65× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.11 Different orders of accuracy for Re τ = 1000, filtered case, Mesh: 96× 65× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.1 LES parameters for Re τ = 2000. . . . . . . . . . . . . . . . . . . . . 122 iv List of Figures 3.1 Channel geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Figure 2 from Lee & Moser (2015), mean velocity profiles. . . . . . 28 3.3 Figures 4, 5 and 6 from Lee & Moser (2015), Reynolds stresses. . . 29 4.1 Transfer functions for filters defined by equation (4.4). G trap : (red) solid line with star symbols, c = 1; G simp : (green) solid line with triangle symbols,c = 2 3 ;G sharp (blue) solid line with circle symbols, c = 1 2 ; and transfer functions for corresponding filters based on equation (4.2) for N = 3 and N = 5. . . . . . . . . . . . . . . . . . 32 4.2 Mean velocity profiles: Re τ = 590, Mesh: 96× 97× 96. . . . . . . . 42 4.3 RMS velocities and Reynolds stresses: Re τ = 590, Mesh: 96×97×96. 43 4.4 Numerical dissipation dependence on the filtering interval. . . . . . 44 4.5 Mean velocity profile and streamwise energy spectra: Re τ = 1000, Mesh: 96× 65× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.6 RMS velocities and Reynolds stresses: Re τ = 1000, Mesh: 96×65× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.7 Mean velocity profile and streamwise energy spectra: Re τ = 2000, Mesh: 256× 131× 128. . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.8 Mean velocity profiles: Re τ = 180, Mesh: 64× 65× 64. . . . . . . . 50 v 4.9 RMS velocities and Reynolds stresses: Re τ = 180, Mesh: 64×65×64. 51 4.10 Mean velocity profiles: Re τ = 1000, Mesh: 256× 97× 128. . . . . . 52 4.11 RMS velocities and Reynolds stresses: Re τ = 1000, Mesh: 256× 97× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.12 Mean velocity profiles: Re τ = 2000, Mesh: 256× 131× 128. . . . . 54 4.13 RMS velocities and Reynolds stresses: Re τ = 2000, Mesh: 256× 131× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.14 Mean velocity profiles: Re τ = 180, different resolutions. . . . . . . . 56 4.15 RMS velocities and Reynolds stresses: Re τ = 180, different resolu- tions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.16 Mean velocity profiles: Re τ = 1000, different resolutions. . . . . . . 58 4.17 RMS velocities and Reynolds stresses: Re τ = 1000, different reso- lutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.18 Mean velocity profiles: Re τ = 180, different filters. . . . . . . . . . . 61 4.19 RMS velocities and Reynolds stresses: Re τ = 180, different filters. . 62 4.20 Mean velocity profiles: Re τ = 590, different filters. . . . . . . . . . . 64 4.21 RMS velocities and Reynolds stresses: Re τ = 590, different filters. . 64 4.22 Modified wavenumber for 1st order derivatives. . . . . . . . . . . . . 66 4.23 Modified wavenumber for 2nd order derivatives. . . . . . . . . . . . 67 4.24 Mean velocity profiles: Re τ = 180, no filter, different orders of accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.25 RMS velocities and Reynolds stresses: Re τ = 180, no filter, different orders of accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.26 Mean velocity profiles: Re τ = 180, filtered with optimal period, different orders of accuracy. . . . . . . . . . . . . . . . . . . . . . . 70 vi 4.27 RMS velocities and Reynolds stresses: Re τ = 180, filtered with optimal period, different orders of accuracy. . . . . . . . . . . . . . 71 4.28 Numerical dissipation computed from the spectral method and the FDM schemes in the solver: Re τ = 1000, no filter. . . . . . . . . . . 72 4.29 Mean velocity profiles: Re τ = 1000, no filter, different orders of accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.30 RMS velocities and Reynolds stresses: Re τ = 1000, no filter, differ- ent orders of accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.31 Mean velocity profiles: Re τ = 1000, filtered with optimal period, different orders of accuracy. . . . . . . . . . . . . . . . . . . . . . . 75 4.32 RMS velocities and Reynolds stresses: Re τ = 1000, filtered with optimal period, different orders of accuracy. . . . . . . . . . . . . . 76 4.33 Onedimensionalenergyspectrathethecenterofthechannel: Re τ = 1000, comparison of no-filter and filtered (with the optimal period) results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.1 Removal of unphysical energy accumulation by scale separation. . . 81 5.2 Filtered 1D velocity profile. . . . . . . . . . . . . . . . . . . . . . . 101 5.3 Sub-grid stress terms (without taking divergence) . . . . . . . . . . 101 5.4 Mean velocity profiles: Re τ = 180, Mesh: 64× 65× 64. . . . . . . . 109 5.5 RMS velocities and Reynolds stresses: Re τ = 180, Mesh: 64×65×64.110 5.6 Mean velocity profiles: Re τ = 180, Mesh: 48× 49× 48. . . . . . . . 111 5.7 RMS velocities and Reynolds stresses: Re τ = 180, Mesh: 48×49×48.112 5.8 Mean velocity profiles: Re τ = 590, Mesh: 96× 97× 96. . . . . . . . 113 5.9 RMS velocities and Reynolds stresses: Re τ = 590, Mesh: 96×97×96.114 5.10 Mean velocity profiles: Re τ = 590, Mesh: 64× 65× 64. . . . . . . . 115 5.11 RMS velocities and Reynolds stresses: Re τ = 590, Mesh: 64×65×64.116 vii 5.12 Mean velocity profiles: Re τ = 1000, Mesh: 128× 97× 128. . . . . . 116 5.13 RMS velocities and Reynolds stresses: Re τ = 1000, Mesh: 128× 97× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.14 Mean velocity profiles: Re τ = 2000, Mesh: 256× 131× 128. . . . . 118 5.15 RMS velocities and Reynolds stresses: Re τ = 1000, Mesh: 256× 131× 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.16 Mean velocity profiles: Re τ = 1000, different resolution. . . . . . . . 120 5.17 RMS velocities and Reynolds stresses: Re τ = 1000, different reso- lution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.18 Mean velocity profiles: Re τ = 2000, different resolution. . . . . . . . 122 5.19 RMS velocities and Reynolds stresses: Re τ = 2000, different reso- lution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.20 Mean velocity profile and streamwise energy spectra: Re τ = 590, Mesh: 64× 65× 64. . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.21 Onedimensionalenergyspectrathethecenterofthechannel: Re τ = 2000, compare with corresponding no model simulations. . . . . . . 126 viii Abstract In large eddy simulations (LES) with coarse mesh, various methods can be used to model the unknown effects of small unresolved scales of turbulent flows. Due to the dominant average effect of forward energy cascade from large to small scales, unphysical energy pileup generally occurs for under-resolved LES meshes. The subgrid scale dissipation can be accounted for by explicit subgrid scale (SGS) stresses. Many widely used SGS models adopt the eddy viscosity concept, which models energy transfer in analogy to the molecular diffusion based on the Boussi- nesq hypothesis. It is well-known, however, that some eddy viscosity closures exhibit relatively weak correlation with real SGS stress, so it is desirable to seek some alternative ways that prevent the buildup of small scale energy. Implicit large eddy simulations (ILES) are increasingly popular in recent years due to ease of implementation and the relaxation of some physical assumptions, where the numerical dissipation works in a manner similar to the explicit SGS models. If spectral methods are used the numerical dissipation is negligible but it can be introduced by applying a low-pass filter in the physical space, resulting in an effective ILES. In the present work we provide a comprehensive analysis of the numericaldissipationproducedbydifferentfilteringoperationsinaturbulentchan- nel flow simulated using a non-dissipative, pseudo-spectral Navier-Stokes solver. The amount of numerical dissipation imparted by filtering can be easily adjusted ix by changing how often the filter is applied. We show that when the additional numerical dissipation is close to the subgrid-scale (SGS) dissipation of an explicit LES the overall accuracy of ILES is also comparable, indicating that periodic fil- tering can replace explicit SGS models. A new method is proposed, which does not require any prior knowledge of a flow, to determine the filtering period adaptively. Once an optimal filtering period is found, the accuracy of ILES is significantly improved at low implementation complexity and computational cost. The method is general, performing well for different Reynolds numbers, grid resolutions, and filter shapes. Another option to remove the unphysical energy buildup is through explicit SGS models based on the multi-scale method. By separating different scales with filtering, wecandirectlyremovesomeoftheenergyproductionneartheLEScutoff, whosecollectiveeffectissimilarasaddingsubgriddissipation. Thescaleseparation is facilitated by a smooth low-pass filter, which becomes increasingly more active for smaller resolved wavenumbers to reduce more production. A gradient-type modification is included to model extra dissipative effect and make the model Galilean invariant, which is found to have negligible impact to the original model. The model is compared with other nonlinear models and regularization techniques both theoretically and numerically. In a posteriori analysis of wall-bounded flow, our new model were able to provide sufficient SGS dissipation without a help of extra dissipative terms such as an eddy viscosity model. This explicit structural model makes it possible to evaluate the nonlinear energy transfer in LES more rigorously. Due to the removal of energy production in desired scales, the model consistently produces more accurate results than the similarity-type model being compared. x Chapter 1 Introduction Large Eddy Simulations (LES) is among the most widely-used technique to simulate turbulent flows. In LES, only the large, energy-containing scales are resolved, relaxing the Direct Numerical Simulations (DNS) requirements for mesh size; the effects of relatively small sub-grid scales (SGS) are modeled, leaving a new unknown in the Navier-Stokes equations. In practice, the removal of the small scales in LES is considered as the result of a low-pass convolution filter. The interaction between small SGS scales and large scale is reflected in two main phenomena: a forward energy cascade from resolved to SGS scales and an opposite backscatter. Due to the dominant average effect of forward energy cascade, there is no adequate means for energy removal in LES. Two basic ways to approximate the effects of missing scales and remove the resolution error can be classified as explicit and implicit LES approach. In explicit LES methods the motion of SGS scales is accounted for by an extra term, either reproducing the eigenvectors of the statistical correlation ten- sors (structural modeling), or just modeling their effects based on energy cascade (functional modeling). According to the review of Domaradzki & Adams (2002a), the traditional SGS modeling approaches can be classified into three main cate- gories: the eddy viscosity models, the similarity models and the mixed models. Perhaps the most widely used group among them is the eddy viscosity models, which is analogous to the Boussinesq hypothesis that models the energy transfer as similar to the molecular diffusion. The first eddy viscosity model was proposed 1 in Smagorinsky (1963) based on the large scales. Since then, many variants of eddy viscosity closure have been proposed in literatures. A detailed review can be found in Sagaut (2006). For wall-bounded flows like the turbulent channel flow, which is of interest in the present work, the eddy viscosity needs to follow the correct asymptotic behavior in the near-wall region. The simplest way is to multiply the original eddy viscosity by a damping function (Van Driest, 1956; Piomelli et al., 1990). There also exists models in which the eddy viscosity automatically vanishes near solid surfaces (Nicoud & Ducros, 1999; Vreman, 2004; Nicoud et al., 2011). Another approach is to compute the closure coefficients in SGS models through dynamic procedure. The most common choice is based on the Germano identity (Germano et al., 1991) and extended later by Lilly (1992) by minimizing the least square error. In general, as the dynamic procedure will not change the form of the model, it can be applied to the closure problem for other types SGS models as well. In this work the results from static Smagorinsky model (Smagorinsky, 1963) with the Van Driest damping function (Van Driest, 1956) and corresponding dynamic procedure (Lilly, 1992) are adopted to compare with our new methods. Alternatively, by assuming that the action of sub-grid scales on the resolved scales is purely dissipative, one can make use of the numerical dissipation as an implicit LES (ILES) model (Sagaut, 2006). The idea of implicit LES or mono- tonically integrated LES (MILES) was originally proposed by Boris et al. (1992), where the FCT scheme was used to ensure the monotonicity of the solution. A detailed overview can be found in Grinstein et al. (2007). The main objective of implicit LES is to selectively account for the effects of small unresolved scales by using the numerical dissipation, which can come from the discretization of govern- ing equations, as well as damping or filtering (Bogey & Bailly, 2006; Thornber et al., 2007; Berland et al., 2011; Kokkinakis & Drikakis, 2015). The idea can also be 2 used to form a deconvolution model to control the truncation error (Domaradzki & Adams, 2002b; Hickel et al., 2006). In general, the implementation of ILES is simpler due to lack of explicit SGS expressions, and may be more general for the situations where some physical assumptions in explicit models are suspicious. The focus of implicit LES is to manipulate the numerical dissipation to provide appropriate amount of SGS dissipation. Analytical information about the trunca- tion errors can be obtained from the modified equation analysis (Grinstein et al., 2007), which shows the relation between explicit and implicit LES. The method to quantify the numerical dissipation was first proposed by Domaradzki et al. (2003) and Domaradzki & Radhakrishnan (2005), but the application to general Navier- Stokes solvers was limited due to the required usage of an auxiliary spectral code. Schranner et al. (2015) extended the approach to the physical space, allowing an efficient estimation of the numerical dissipation rate and the numerical viscosity with the original numerical schemes in the solver. The amount of numerical dissi- pation can thus be found in a post-processing step feasible for any Navier-Stokes solver. This makes it possible to evaluate and improve the accuracy of DNS or under-resolved DNS (UDNS) during the simulation (Castiglioni & Domaradzki, 2015; Sun et al., 2017), and the comparison with viscous and explicit SGS dis- sipation can also provide guidance in the development of implicit LES methods (Cadieux et al., 2017). The design and testing of the ILES method described in the current work depends critically on the availability of this procedure to quantify the numerical dissipation in an arbitrary Navier-Stokes solver. Manyapproachescanbeusedtoremovetheaccumulatedunphysicalenergydue to insufficient grid resolution in LES. Apart from adding the numerical dissipation through the discretized convective term for specific numerical schemes, better con- trol of the numerical dissipation can be achieved either by modifying terms in the 3 governing equations, or by external processing such as filtering. Domaradzki et al. (2002) applied a periodic filtering as a simplified Truncated Navier-Stokes (TNS) approach, in which filtering and re-initialization is performed periodically with a prescribedtimeintervalbeforetheerrorofenergypileupinthesmallscalescontam- inates the large scales. Another way to add artificial dissipation near mesh cutoff is to use the hyperviscosity-modified Navier-Stokes equation (Borue & Orszag, 1998; Lamorgese et al., 2005). Lamballais et al. (2011) offered an easy way to make the use of the hyperviscosity by checking the modified wavenumbers in the viscous term such that numerical dissipation is introduced only for higher wave numbers. With a similar idea, the addition of a spectral vanishing viscosity (SVV) (Got- tlieb & Hesthaven, 2001; Karamanos & Karniadakis, 2000), can be regarded as an alternative ILES model as well. In some cases it is also beneficial to add numer- ical dissipation in other LES approaches, resulting in mixed models, for instance Pasquetti & Xu (2002) combined an approximate deconvolution model with SVV. Dairayetal.(2017)introducedSVVbythemodifiedwavenumberanalysisadjusted by an exponential kernel for the homogeneous and isotropic turbulence. A theo- retical analysis of corresponding spectral viscosity ensures the truncation of energy spectrum at the cutoff of LES as small as DNS. However, to implement the above methods, either prior knowledge of the flow or trial and error testing is necessary to achieve the best performance. As their generality cannot be assured attempts to automatically determine the amount of numerical dissipation can be made. Tantikul & Domaradzki (2010) proposed a global period filtering criterion based on the energy ratio after two different filter operations are applied. Flad et al. (2016) developed a space-dependent local filtering criterion using a similar idea for the discontinuous Galerkin method. Li & Tsubokura (2017) extend the use 4 of the energy ratio to adaptively scale the Roe upwind dissipation term for low- Mach number compressible flows. Essentially, the cited methods can automatically provide extra numerical dissipation via some specific approaches, but are still less universal because their quantitative connection with the numerical dissipation rate is unknown. The first part of the proposal is to present a more general method that imparts anappropriateamountofnumericaldissipationbasedonanadaptivefilteringcrite- rion. We correct an error in the energy ratio approach of Tantikul & Domaradzki (2010) so that the reasoning is more robust. We find that for filters that only weakly attenuate the small scales, when similar amounts of numerical dissipation are added, comparable overall accuracy can be achieved for different filters, with stronger filters applied less frequently based on the criterion. Moreover, as a global operation, computing the energy ratio every several time steps can be expensive. To make ILES as efficient as corresponding no-model UDNS, we develop a method that for a prescribed filter determines an optimal filtering period, which confines the energy ratio in a desired range so that no updates to the filtering period are necessary in the statistically steady state. To better evaluate the performance of the current approach, unlike Tantikul & Domaradzki (2010), no extra stabilization technique is used in our non-dissipative pseudo-spectral solver, with the adaptive filter found to be sufficient to stabilize the simulations. The approach is capable of improving the accuracy for lower order dissipative schemes as well. This was tested by changing the pseudo-spectral method to different finite different schemes through the modified wavenumber approach. A posteriori testing has been con- ducted for turbulent channel flow with various Reynolds numbers, grid resolutions, and filter shapes, and the predictions are compared with DNS data as well as with LES results obtained using explicit SGS models. 5 As an natural extension, the second part of this work will focus on removing the unphysicalenergypileupwithexplicitSGSmodelsbasedonthemulitscalemethod. The energy transfer between scales received much attention in both analytical analysis and numerical modeling of turbulent flows. The scale separation based on filtering is used to study the interscale energy transfer for incompressible turbulent flows (Domaradzki & Rogallo, 1990; Eyink, 1994; Domaradzki et al., 1994; Kerr et al., 1996; Eyink, 2005; Domaradzki & Carati, 2007a,b; Eyink & Aluie, 2009; Aluie & Eyink, 2009). In large eddy simulations, better understanding of energy transfer can be used to design new SGS models (Anderson & Domaradzki, 2012) or improve existing ones (Cimarelli & Angelis, 2014). In particular, many existing SGS models based on the eddy viscosity concept are purely dissipative. Piomelli et al. (1991) pointed out that accurate predictions of backscatter are required in LES modeling. Domaradzki et al. (1993) found that both forward and backward energy transfer has important impact on the resolved scales based on analyses of DNS data. Therefore, many attempts have been made to better account for the interaction between multiple scales. One strategy is to decompose velocity fields into different grid or filtering levels. For instance in the Variational Multiscale Method proposed by Hughes et al. (1998, 2000, 2001), the subgrid term is modeled by a Smagorinsky-like eddy viscosity model but allows for both small-small and large-small representations. This idea can be applied to other SGS models as well. Another approach is to regularize the Naiver-Stokes equation using smoothed or filtered velocities in order to achieve a faster decay of energy spectra at high wavenumbers, such as the family ofα models (Leray, 1934; Chen et al., 1998; Foias et al., 2001; Cheskidov et al., 2005). In other words, a DNS of the regularized equations with relaxed resolution requirements can be reinterpreted as a LES of the Navier-Stokes equations. Nevertheless, these regularizations normally focus 6 more on the mathematical properties of the equations rather than considering the actual effects of the subgrid scales in nonlinear energy transfer. Alternatively, Scale Similarity is another family of models based on physical hypothesis that the structure of subgrid scales is similar to corresponding resolved scales due to their common history in energy cascade and coherent structures. Similarity model was first proposed by Bardina et al. (1980). It showed considerable correlation with the real stress during a priori analysis (Bardina et al., 1983; Piomelli et al., 1988; Horiuti, 1989). However, the method by itself is found not sufficiently dissipative in LES, and is used as a component in mixed models (e.g. with an eddy viscosity term) (Bardina et al., 1983; Zang et al., 1993; Liu et al., 1994a). Anderson & Domaradzki (2012) proposed an interscale energy transfer model as an alternative to the scale similarity model but from a different perspective, which intends to remove the small scale production for scales near the LES cutoff to model the effect of subgrid dissipation. However, to ensure Galilean invariance, mathematically the modified interscale energy transfer model is almost equivalent to the original similarity scale model (Bardina et al., 1983) thus showing the same weaknesses (Anderer, 2015). Therefore another mathematical expression is sought inthisworktobetterrepresenttheoriginalassumptionofAnderson&Domaradzki (2012) while still ensure Galilean invariance. The idea of removing unphysical energy buildup in large wavenumbers is analogous to filtering in our ILES approach (Sun & Domaradzki, 2018), but by designing this explicit SGS model we hope to gain deeper understandings of multiscale energy transfer in under-resolved LES as well as benefit from the similarity-like structure modeling to better account for backscatter. For simulations of channel flow with the same conditions as the implicit method, the present model in general is at least as satisfactory as our 7 ILES approach and the interscale model in Anderson & Domaradzki (2012) for mean velocity profile, and the Reynolds stresses can be more accurately predicted. 8 Chapter 2 Background 2.1 Governing Equations The governing equations of fluid dynamics for incompressible flow (density ρ=const.,∇·~ u = ∂u i ∂x i = 0) are the Navier-Stokes equations ρ ∂u i ∂t +u j ∂u i ∂x j ! = ∂ ∂x j (−pδ ij +μτ ij ), i = 1, 2, 3 (2.1) where ~ u = (u 1 ,u 2 ,u 3 ) is the velocity field in the Cartesian coordinate system, p is the pressure, and μ is the dynamic viscosity. Subscript i stands for different directions in the Cartesian coordinates, for the channel flow, x 1 ,x 2 ,x 3 represent streamwise (x), vertical (y) and spanwise (z) coordinates, respectively. For a Newtonian fluid, the stress tensor is τ ij = ∂u i ∂x j + ∂u j ∂x i − ∂u k ∂x k δ ij . (2.2) Due to the nonlinear nature of the equations, there is no closed form analytical solution for the Navier-Stokes equations for general turbulent flow. In computa- tional fluid dynamics (CFD) the velocity fields are projected on a set of discrete grid points to where the governing equations are solved numerically. The accuracy of numerical simulation not only depends on the grid resolution and the quality of the mesh, but also on the numerical schemes (e.g. finite difference, finite volume, and finite element methods) used to approximate the governing equations. 9 2.2 Direct Numerical Simulation Direct Numerical Simulation (DNS) offers a complete description of the tur- bulence by solving the exact form of the governing equations without additional turbulent models. A well-resolved DNS requires that the computational grid to captures all relevant scales of flow motion. According to Kolmogorov’s hypotheis, there exists a universal equilibrium range where the statistics of the small-scale motions have a universal form uniquely determined by the viscosity ν and dissi- pation rate of kinetic energy ε. As a result the mesh size of DNS must be within the same order of magnitude as the Kolmogorov length scale η obtained from dimensional analysis η = ν 3 ε ! 1/4 , (2.3) where ν =μ/ρ is the kinematic viscosity. Typically high order numerical schemes are also used for DNS, therefore a well- resolved DNS can be regarded as the most accurate CFD tool available. However, the requirement of resolving the Kolmogrov length scaleη leads to a quick increase of a number of grid points for higher Reynolds number η L = ν 3/4 ε 1/4 L ∼Re −3/4 L , (2.4) where L is the integral length scale and Reynolds number Re L = UL ν . Therefore the number of grid points required to resolve all eddies∼Re 9/4 L for 3 dimensional homogeneous flows, so the computational cost of DNS can be prohibitively expen- sive even for moderate Reynolds numbers. Currently the finest wall-bounded DNS carried out by Lee & Moser (2015) forRe τ ≈ 5200 already used 242 billion degrees of freedom and 786K processing cores. In many circumstance a well-resolved DNS 10 is thus infeasible and it is necessary to seek other approaches for turbulence simu- lation. Note that in general a Direct Numerical Simulation can also be thought of as a numerical simulation without turbulence models. In the rest of this proposal we use the term under-resolved DNS (UDNS) to indicate no-model simulation of turbulence with a grid resolution insufficient to capture the smallest length scales, while the term DNS still stands for well-resolved simulations typically treated as benchmark. 2.3 Large Eddy Simulation A method to greatly relax DNS’s severe resolution requirement is Large Eddy Simulations (LES), which is designed based on Kolmogorov’s hypothesis that the turbulent motion of the smallest scales are dominated by viscosity and shows sta- tistically isotropic behavior. Therefore in LES, only the large, energy-containing scales are resolved, relaxing the Direct Numerical Simulations (DNS) requirements for mesh size; the relatively universal effects of small sub-grid scales (SGS) are modeled, leaving a new unknown in the Navier-Stokes equations. In explicit LES methods the motion of SGS scales is accounted for by an extra term based on the knowledge of resolved scales, either reproducing the eigenvectors of the statistical correlation tensors (structural modeling), or just modeling their effects based on energy cascade (functional modeling). Under-resolved mesh can be used in LES to reduce the degrees of freedom to a manageable number, whose effect is comparable to a low-pass filter. The effect of mesh is typically referred to as an implicit filtering, which separates larger scales computed in the filtered Navier-Stokes equations and smaller unresolved scales 11 called sub-grid scales (SGS). In spectral space this can be thought of as resolving the quantities up to a cutoff wavenumber k c = π/Δ (the Nyquist wavenumber) based on the mesh size Δ. The filtering operation with function b G is a direct multiplication in the spectral space. In physical space this can be expressed with a convolution integral with the filter kernel G (inverse Fourier transform of b G) u i (~ x,t) = ZZZ G(~ x−~ r, Δ)u i (~ r,t)d 3 ~ r, (2.5) where the velocity for the resolved scales is ¯ u i , leaving the sub-grid scale (SGS) velocity as u i − ¯ u i . The filter width Δ is generally proportional to the grid size. Accounting for the effect of sub-grid scales, the LES equation can be written as ∂¯ u i ∂t + ∂ ∂x j ¯ u i ¯ u j +τ SGS ij = ∂ ∂x j − ¯ p ρ δ ij + 2ν ¯ S ij ! , (2.6) where ν =μ/ρ is the kinematic viscosity, and τ sgs ij =u i u j − ¯ u i ¯ u j , (2.7) ¯ S ij = 1 2 ∂¯ u i ∂x j + ∂¯ u j ∂x i ! . (2.8) The subgrid-scale (SGS) stress tensor τ sgs ij is a new unknown for which an explicit SGS model or an additional equation is required to account for the effects of the small unresolved scales. In order to evaluate the performance of our new LES approaches, several com- parisons are made based on existing SGS models. 12 2.3.1 Explicit subgrid-scale Models Two classic eddy viscosity models are introduced in this subsection. We imple- mented the models in our pseudo-spectral Navier-Stokes solver to directly compare their simulation accuracy with our new methods. The Smagorinsky Model The Smagorinsky model (Smagorinsky, 1963) is among the most famous mod- els in functional modeling; it is an eddy viscosity model based on the Boussinesq approximation to relate turbulent transport to the dissipation through an addi- tional viscosity. The modeled SGS model takes the following form τ sgs ij − 1 3 τ sgs kk δ ij =−2ν sgs ( ¯ S ij − 1 3 ¯ S kk δ ij ), (2.9) where ν sgs =` 2 s | ¯ S| = (C s Δ) 2 | ¯ S|, (2.10) | ¯ S| = (2 ¯ S ij ¯ S ij ) 1 2 , (2.11) Δ = (ΔxΔyΔz) 1 3 , (2.12) C s = const. (2.13) For wall-bounded flows, additional van Driest damping function (Van Driest, 1956) is applied to ` s to achieve the correct asymptotic behavior near wall ` s =C s Δ[1− exp(−yu τ /25ν)] (2.14) 13 where for channel flow, y is the vertical (wall-normal) direction and the friction velocity u τ = q τ w /ρ, τ w is the wall shear stress. Smagorinsky Model with Germano-Lilly Dynamic Procedure To better account for the local structure of the flow, Germano identity can be used to dynamically determine the constant in some SGS models, which can be written as (Germano et al., 1991; Lilly, 1992) L ij = g ¯ u i ¯ u j − ˜ ¯ u i ˜ ¯ u j =T ij − ˜ τ ij = ( g u i u j − ˜ ¯ u i ˜ ¯ u j )− ^ (u i u j − ¯ u i ¯ u j ) (2.15) By applying a test filter denoted as ˜ with a larger characteristic filter width than Δ and assuming that stressτ ij andT ij can have the same form (e.g. the static Smagorinsky model), the coefficient C d that minimizes the error of the Germano identity using the least square approach ν sgs =C d Δ 2 | ¯ S|, (2.16) C d = min 0.2 2 , max " hL ij M ij i hM ij M ij i , 0 #! , (2.17) where M ij = 2 ^ Δ 2 | ¯ S| ¯ S ij − 2 e Δ 2 | e ¯ S| e ¯ S ij . The operatorh·i denotes a local averaging of neighboring points with equal weights. C d is restricted in the range of [0, 0.2] according to the common choices of C s (Smagorinsky constant, Eq. (2.13)). The test-filtered quantities are obtained using Simpson’s rule (filter width 2Δ), for uniform mesh e u(x)≈ 1 6 u(x− Δx) + 2 3 u(x) + 1 6 u(x + Δx). (2.18) 14 When the mesh is not uniform we first apply a local quadratic interpolation (Loh & Domaradzki, 1999) to use the same filtering kernel, the procedure is described in Section 3.5. 2.3.2 Implicit subgrid-scale Models In some cases, data from the Truncated Navier-Stokes (TNS) approach (Domaradzki et al., 2002) is also used for comparison due to its close relation to our new implicit LES method. Note that the TNS approach used for compari- son is more of a mixed model, with an explicit part based on the SGS estimation model (Domaradzki & Saiki, 1997) and an implicit part in which numerical dis- sipation is introduced by filtering. The details of the approach are presented as below. The Truncated Navier-Stokes (TNS) Approach (simplified version) TheTruncatedNavier-Stokes(TNS)approachisproposedbyDomaradzkietal. (2002). The TNS equations are just Navier-Stokes equations discretized on a mesh twice finer than the LES equations of interest, so the solution of TNS would behave asasolutionofunder-resolvedDNS(UDNS).TheinitialmethodrequiresbothLES and TNS equations to run in parallel, the estimated velocity (Domaradzki & Saiki, 1997; Domaradzki & Loh, 1999) is used as the initial condition for TNS, which is consistent which the LES velocity field at large scales. The TNS velocity is directly used to compute the SGS stress tensor for LES equation within a period T, then TNS is stopped and re-initialization is performed before error propagate toward large scales due to insufficient resolution. As TNS equation on a finer mesh already containslargescaleinformationincorrespondingLESequation, Domaradzki&Yee (2000) provided a simplified approach without running a separate LES, which is 15 extended in Domaradzki et al. (2002) to achieve more accurate simulations for turbulent channel flow. As a first step, the large scale velocity ˜ u 0 i is obtained based on an approximate deconvolution operator following Stolz & Adams (1999) based on a three-point top-hat filter approximated with trapezoidal rule: ˆ ˜ u i =G∗u i = 1 4 ˜ u i−1 + 1 2 ˜ u i + 1 4 ˜ u i−1 , (2.19) ˜ u 0 i = [I− (I−G) 2+1 ]∗ ˜ u i = 3 ˆ ˜ u i − 3 ˆ ˆ ˜ u i + ˆ ˆ ˆ ˜ u i . (2.20) The perturbation velocity ˜ u 0 i is then estimated from ˜ u 0 i ˜ u 0 i =θ(N 0 i −N 0 i ), (2.21) where the ¯ filter here is also the top-hat filter with the same width as ˆ but using Simpson’s rule of integration (Eq. (2.18)), N 0 i is the nonlinear term from which the advection effects by the large scales are removed: N 0 i =−(˜ u 0 j − ˜ u 0 j ) ∂ ∂x j ˜ u 0 i , (2.22) θ is related to the large eddy turnover time in given locations computed from θ =R ˜ u 0 i − ˜ u 0 i N 0 i −N 0 i , (2.23) in which constantR is chosen to be 1 2 for the−5/3 spectrum in the inertial range. For channel flow with constant mass flux, effects of the mean of ˜ u i are removed: ˜ u 00 i = ˜ u 0 i − ˆ ˜ u 0 i (2.24) 16 Finally thecombination of ˜ u i = ˜ u 0 i +˜ u 00 i is advanced in timein the TNS equation subject to a manually prescribed re-initialization period T. Therefore, the model only have a periodic replacement of small scale component based on the SGS estimation procedure and aimed at predicting the larges scales accurately. The method performs comparably with the parallel LES/TNS approach with standard SGS estimation model (Domaradzki & Loh, 1999). 17 Chapter 3 Numerical Method 3.1 Skew-symmetricformofNavier-Stokesequa- tion The incompressible Navier-Stokes equations are solved using a Fourier- Chebyshev pseudo-spectral solver using a time-splitting and a pressure-correction approach, whose spectral accuracy effectively eliminates numerical dissipation. The solver was previously used and validated in Cadieux & Domaradzki (2015) to simulate a laminar separation bubble flow. The details of spectral method can be found in Peyret (2002) and Canuto et al. (2007). The skew-symmetric form of Navier-Stokes equation is adopted in the current solver for a better modeling of the small scales and less aliasing (Kravchenko & Moin, 1997). For Newtonian fluids the equation can be written as ρ ∂u i ∂t + 1 2 u j ∂u i ∂x j + 1 2 ∂u i u j ∂x j ! = ∂ ∂x j " −pδ ij +μ ∂u i ∂x j + ∂u j ∂x i !# . (3.1) In general when there is a body force and a subgrid scale model ∂¯ u i ∂t =− 1 2 ¯ u j ∂¯ u i ∂x j − 1 2 ∂¯ u i ¯ u j ∂x j − 1 ρ ∂p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j + ∂¯ u j ∂x i ! +F i − ∂τ sgs ij ∂x j (3.2) =− 1 ρ ∂p ∂x i +ν ∂ 2 ∂x j ∂x j ¯ u i +N i (3.3) 18 where F i is body force, and τ sgs ij = u i u j − ¯ u i ¯ u j is the subgrid scale stress tensor, N i represent the remaining terms N i =− 1 2 ¯ u j ∂¯ u i ∂x j − ∂ ∂x j 1 2 ¯ u i ¯ u j +τ sgs ij +F i (3.4) 3.2 Channel flow Figure 3.1: Channel geometry For the channel flow (Fig.3.1), in the vertical direction (y) the Chebyshev poly- nomials are used to achieve higher near-wall resolution and Fourier polynomials in the streamwise (x) and spanwise (z) directions. No-slip boundary conditions are imposed at the walls and periodic boundary conditions in the horizontal direc- tions. In the pseudo-spectral solver all the derivatives are computed based on the Fast Fourier Transform (FFT). The evaluation of derivatives on the Gauss-Lobatto points iny-direction based on FFT can be found in Kopriva (2009). The constant mass flux in the streamwise direction (i = 1) is enforced by setting the mean pressure gradient to the instantaneous wall stress at each time step: 19 ∂¯ u i ∂t =− 1 2 ¯ u j ∂¯ u i ∂x j − 1 2 ∂¯ u i ¯ u j ∂x j − 1 ρ   ∂p ∂x i + μ 2h ∂¯ u 1 ∂y y=h y=−h δ 1i   +ν ∂ 2 ∂x j ∂x j ¯ u i − ∂τ sgs ij ∂x j (3.5) =− 1 ρ ∂p ∂x i +ν ∂ 2 ∂x j ∂x j ¯ u i +N i , (3.6) where N i =− 1 2 ¯ u j ∂¯ u i ∂x j − ν 2h ∂¯ u 1 ∂y y=h y=−h δ 1i − ∂ ∂x j 1 2 ¯ u i ¯ u j +τ sgs ij +F i . (3.7) 3.3 Temporal discretization of Navier-Stokes equation The temporal discretization is based on a combined nonlinear viscous step with pressure-correction step using a third-order backward difference scheme described in Karniadakis et al. (1991). 3.3.1 Step I: nonlinear and viscous LES equation (Eq. 3.6) without the pressure term ∂¯ u i ∂t =ν ∂ ∂x j ∂¯ u i ∂x j + ∂¯ u j ∂x i ! +N i . (3.8) A DNS simulation would follow the same framework but without the ¯ intro- duced by implicit filtering of mesh and the subgrid scale model inN i . 20 In spectral space ∂ b ¯ u i ∂t =ν ∂ 2 ∂y 2 −k 2 x −k 2 z ! b ¯ u i + c N i . (3.9) where k x and k z are wavenumbers in Fourier transform. The fully implicit scheme is applied to viscous term, a general backward differ- ence formula for nonlinear term can be written as 1 Δt   γ 0 b ¯ u ∗ i − J−1 X q=0 α q b u n−q i   −ν ∂ 2 ∂y 2 −k 2 x −k 2 z ! b ¯ u ∗ i = J−1 X q=0 β q c N n−q i , (3.10) where the superscripts represent different time steps, γ 0 , α q , β q are backward difference formula of order J coefficients available in Karniadakis et al. (1991). Rearranging the equation D 2 −k 2 − γ 0 νΔt b ¯ u ∗ i =− 1 νΔt J−1 X q=0 α q b u n−q i − 1 ν J−1 X q=0 β q c N n−q i , (3.11) whereD = ∂ ∂y andk 2 =k 2 x +k 2 z , leads to one Helmholtz equation in each direction. 3.3.2 Step II: pressure Combine the pressure part with continuity          γ 0 Δt ¯ u n+1 i − ¯ u ∗ i + 1 ρ ∂p n+1 ∂x i = 0 ∂¯ u i ∂x i = 0 (3.12) 21 Equation for Fourier modes b u i (k x ,y,k z ,t) and b p(k x ,y,k z ,t)                            γ 0 Δt (¯ u n+1 − ¯ u ∗ ) =−ik x b p ρ γ 0 Δt (¯ v n+1 − ¯ v ∗ ) =− ∂ ∂y b p ρ γ 0 Δt ( ¯ w n+1 − ¯ w ∗ ) =−ik z b p ρ ik x b u + ∂b v ∂y +ik z b w = 0 (3.13) Rearrange in spectral space, leading to one equation for v (vertical velocity), which is equivalent to a Poisson equation in physical space D 2 −k 2 b ¯ v n+1 =−D ik x b ¯ u ∗ +ik z b ¯ w ∗ −k 2 b ¯ v ∗ , (3.14) where no-slip Dirichlet boundary condition is applied b v| y=−h,h = 0. (3.15) Once b ¯ v n+1 is computed, the pressure and streamwise and spanwise velocities can be obtained as follows based on previous equations b p ρ =− D b ¯ v n+1 +ik x b ¯ u ∗ +ik z b ¯ w ∗ / k 2 Δt/γ 0 , (3.16) b ¯ u n+1 = b ¯ u ∗ −ik x Δt γ 0 b p ρ , (3.17) b ¯ w n+1 = b ¯ w ∗ −ik y Δt γ 0 b p ρ . (3.18) 22 3.4 Quantification of numerical dissipation The transport equation for the kinetic energy (e kin = 1 2 u i u i ) can be obtained by multiplying the Navier-Stokes equation (2.1) by u i ∂ρe kin ∂t + ∂ (ρe kin u i ) ∂x i =−u i ∂p ∂x i +μu j ∂τ ij ∂x i . (3.19) The second term in the right hand side can be re-written as μu j ∂τ ij ∂x i =μ ∂ (u j τ ij ) ∂x i −μτ ij ∂u j ∂x i . (3.20) For incompressible flow and Newtonian fluids the kinetic energy equation becomes ∂e kin ∂t + ∂ (e kin u i ) ∂x i + 1 ρ ∂ (pu i ) ∂x i − ∂ ∂x i ν ∂e kin ∂x i ! +ν ∂u j ∂x i ∂u j ∂x i = 0. (3.21) Integrating the kinetic energy transport equation over an arbitrary control vol- ume leads to ∂ ∂t E kin +F e kin +F ac −F ν +E ν = 0, (3.22) where each term corresponds to the respective terms in equation (3.21) integrated over the same sub-domain. E kin is the kinetic energy in a sub-domain, F e kin , F ac and F ν are contributions from the advective, pressure, and viscous fluxes, respectively, andE ν is the integral form of the viscous (or physical) dissipation. Based on the results of Schranner et al. (2015), when the above equation is discretized using finite difference or finite volume methods, the discretized terms (denoted as (d)) do not add to zero but leave a residual: 23 −E n = ∂ ∂t E (d) kin +F (d) e kin +F (d) ac −F (d) ν +E (d) ν . (3.23) The residualE n can be regarded as the numerical dissipation rate as shown by (Schranner et al., 2015; Castiglioni & Domaradzki, 2015; Sun et al., 2017; Cadieux et al., 2017). In analyzing the effects of the numerical dissipation it is also useful to consider the ratio of numerical and physical dissipationR n =E n /E ν . The kinetic energy equation for LES simulations of turbulent channel flow with aconstantmassfluxcanbeobtainedinthesamewaybymultiplyingequation(3.6) (without additional body force in current simulations) by ¯ u i ∂¯ e kin ∂t + ∂ (¯ e kin ¯ u i ) ∂x i +¯ u i ∂τ sgs ij ∂x j +   1 ρ ∂ (¯ p¯ u i ) ∂x i + ν 2h ¯ u 1 ∂¯ u 1 ∂y y=h y=−h   − ∂ ∂x i ν ∂¯ e kin ∂x i ! +ν ∂¯ u j ∂x i ∂¯ u j ∂x i = 0. (3.24) In the present work, all numerical dissipation terms are obtained by integrating the above equation over a half channel −E n = Z Lx 0 Z h 0 Z Lz 0   ∂¯ e kin ∂t + ∂¯ e kin ¯ u i ∂x i + ¯ u i ∂τ sgs ij ∂x j + ¯ u i ρ ∂ ¯ p ∂x i + ν 2h ¯ u 1 ∂¯ u 1 ∂y y=h y=−h −ν ∂ 2 ¯ e kin ∂x 2 i +ν ∂¯ u j ∂x i ∂¯ u j ∂x i # dxdydz (3.25) = Z Lx 0 Z h 0 Z Lz 0   ∂¯ e kin ∂t +ν ∂¯ u i ∂x j ∂¯ u i ∂x j −τ sgs ij ∂¯ u i ∂x j +   ¯ u i ρ ∂ ¯ p ∂x i + ν 2h ¯ u 1 ∂¯ u 1 ∂y y=h y=−h     dxdydz + ZZ A n j " ¯ e kin ¯ u j −ν ∂e kin ∂x j + ¯ u i τ sgs ij # dA (3.26) = ∂ ¯ E (d) kin ∂t + ¯ E (d) ν − ¯ E (d) sgs + ¯ W (d) ac + ¯ F (d) e kin − ¯ F (d) ν + ¯ F (d) sgs , (3.27) where n j with j ∈ {1, 2, 3} are the components of the outward pointing unit vector normal to each surface A bounding the half channel. The pressure flux ¯ F ac in equation (3.23) can no longer be written as a flux term after adding the mean 24 pressure gradient, so the combined pressure term is denoted as pressure work ¯ W ac in equation (3.27). Note that the SGS dissipation term ¯ E sgs =τ sgs ij ∂¯ u i ∂x j can also be written using a more common definition τ sgs ij ¯ S ij , given that τ sgs ij is symmetric. For the current analysis, the unsteady term in the kinetic energy equation is computed by central differences ∂ ¯ E kin ∂t = ¯ E n+1 − ¯ E n−1 2Δt . (3.28) The numerical dissipation can be regarded as the effect of an extra numerical viscosity ν n = E n E ν ν =R n ν. (3.29) To compare the numerical dissipation with the dissipation of an explicit SGS modelE sgs , the effective viscosity associated with the SGS dissipation is computed in the same way ν sgs = E sgs E ν ν. (3.30) ν sgs can be regarded as an effective SGS viscosity used to quantify the SGS dissipa- tioncomparingwiththeviscousdissipation, whichissimilartothesubgrid-activity parameter s defined in Geurts & Fröhlich (2002). 3.5 Filtering on nonuniform mesh When filtering is applied for nonuniform mesh, it can be implemented after a local quadratic interpolation. For the three-point top-hat filter used in the current work, we define Δ − = y n −y n−1 and Δ + = y n+1 −y n . Then it is pos- sible to interpolate velocities to uniform mesh (denoted by superscript ‘*’) using 25 min{Δ − , Δ + }. The result when Δ + > Δ − is as follows, which simply replacesy n+1 with y =y n + Δ − based on the closer neighboring point U n =aU ∗ n−1 +bU ∗ n +aU ∗ n+1 =c 0 U n−1 +c 1 U n +c 2 U n+1 . (3.31) where the filtering coefficients c 0 =a (y−y n )(y−y n+1 ) (y n−1 −y n )(y n−1 −y n+1 ) +a =a Δ − (Δ − − Δ + ) Δ − (Δ − + Δ + ) +a = 2aΔ − Δ + + Δ − , (3.32) c 1 =a (y−y n−1 )(y−y n+1 ) (y n −y n−1 )(y n −y n+1 ) +b = a2Δ − (Δ − − Δ + ) −Δ − Δ + +b = 2a 1− Δ − Δ + ! +b, (3.33) c 2 =a (y−y n−1 )(y−y n ) (y n+1 −y n )(y n+1 −y n−1 ) = 2aΔ − Δ − Δ + (Δ + + Δ − ) . (3.34) Then the filtered velocity U = 2a Δ − Δ + + Δ − U n−1 + " 2a 1− Δ − Δ + ! +b # U n + 2a Δ − Δ − Δ + (Δ + + Δ − ) U n+1 . (3.35) Similarly, when Δ + < Δ − U = 2a Δ + Δ + Δ − (Δ − + Δ + ) U n+1 + " 2a 1− Δ + Δ − ! +b # U n + 2a Δ + Δ − + Δ + U n−1 . (3.36) According to Loh & Domaradzki (1999), for incompressible turbulent channel flow, the commutation errors due to this local interpolation are insignificant for plan averaged quantities. The procedure is used in the vertical direction of the channel flow in order to apply the same filter to all three directions. 26 3.6 Validation For turbulent channel flow results are generally presented in a dimensionless form based on the friction velocity u τ = q τ w /ρ, where τ w is the wall shear stress. Dimensional analysis leads to the normalized velocity U + = U/u τ , length y + =yu τ /ν, and Reynolds stresses R + ij =R ij /u 2 τ =u 0 i u 0 j /u 2 τ , in which U(y) is the streamwise mean velocity averaged in homogeneous directions x andz, andu 0 i are the fluctuating velocities for i = 1, 2, 3 in each Cartesian coordinate x,y,z. In order to quantify the relative error in comparison with DNS benchmarks, we follow approach of Bose et al. (2010) and Toosi & Larsson (2017) with the overall error defined as Err DNS = 1 5 R b a |U + −U + DNS |dy + R b a U + DNS dy + + X i,j 1 5 R b a |R + ij −R + ij,DNS |dy + R b a 1 2 R + kk,DNS dy + . (3.37) In the formula above only Reynolds stressesR uu ,R vv ,R ww ,R uv are included inR ij , and R kk =R uu +R vv +R ww . All terms are integrated over half channel width. These lower order statistics are compared with benchmark results from Direct Numerical Simulation (DNS), which are available for the canonical incompressible turbulence channel flow used to test our models. The DNS data are obtained from Moser et al. (1999) and Lee & Moser (2015) for different frictional Reynolds number Re τ = u τ h/ν. For relatively low Reynolds numbers, DNS data is more commonly used to evaluate simulation results since it can provided very detailed information at every grid point. In Lee & Moser (2015), statistics for turbulent channel flow for several DNS results are compared with experimental data from laser Doppler velocimetry measurements in Schultz & Flack (2013). The lower order statistics we interested matches closely with each other. All their test cases 27 are summarized in Table 3.1, some of the comparisons are show in figures 3.2 and 3.3. These comparisons confirm that the DNS statistics almost collapses on results from experiments. Therefore, we only compare with DNS data in the current work. Table 3.1: DNS and experimental data compared in Lee & Moser (2015) (Table 1 therein). Figure 3.2: Figure 2 from Lee & Moser (2015), mean velocity profiles. Note that in equation (3.37), all terms are treated as if they had equal weights. However, itmustbeunderstoodthatthemeanvelocityprofileshouldbeconsidered as a priority because the Reynolds stresses such as R + uu hardly could be trusted 28 Figure 3.3: Figures 4, 5 and 6 from Lee & Moser (2015), Reynolds stresses. when the mean velocity U + is off. As a reference, the solver used in the present work produced an overall error (Err DNS ) of roughly 0.7% for Re τ = 180 for the same resolution as used in Moser et al. (1999). 29 Chapter 4 Use filtering as an implicit LES method 4.1 General idea In implicit large eddy simulations (ILES) numerical dissipation is necessary to remove the energy accumulation in the small scales near the cutoff before they can contaminate large scales of interest. If spectral methods are used the numerical dissipation is negligible but it can be introduced by applying a low-pass filter, resulting in an effective ILES. Itwasshownpreviously(Domaradzkietal.,1993)thatitissufficienttoconsider the energy pileup between wavenumber 1 2 k c and k c , where k c =π/Δ is the cutoff wavenumber based on the mesh size Δ. In order to construct filters that mostly affect this region, given any primary filter G, a series of secondary filters obtained from the deconvolution method is used in the current work. For the so-called soft deconvolution problem (Sagaut, 2006), if one assumes thatG is an invertible filter, then theoretically one may recover the full velocity field from its filtered version using a filter inverse u =G −1 ∗ ¯ u (e.g., the inverse modeling in Geurts (1997) and the Approximate Deconvolution Model (ADM) of Stolz & Adams (1999); Stolz et al. (2001)). In practice, deconvolution is frequently used in image processing. In general, the inversion of a continuous filter is ill-conditioned and singular problem may occur, so the original filter cannot be accurately recovered, but it is possible 30 to generate a weaker filter through approximate deconvolution. Here we choose the discrete van Cittert deconvolution operator. The filtering and inverse filtering is facilitated by various iterative methods based on polynomial approximations of the filter kernel G −1 ≈Q N = N X ν=0 (I−G) ν = I− (I−G) N+1 G , (4.1) Q N G = N X ν=0 (I−G) ν G =I− (I−G) N+1 . (4.2) A secondary filter is constructed as a product of Q N and G: ¯ U = (Q N G)∗U = " N X ν=0 (I−G) ν # G∗U. (4.3) It retains most of the large scale information and mainly damps the small scales. In terms of the primary filter G, several three-point filters based on different quadrature rules were tested for a uniform gird in 1D G∗U i = c 4 U i−1 + 1− c 2 U i + c 4 U i+1 . (4.4) In particular c = 1 leads to a filter based on the trapezoidal rule with coefficients ( 1 4 , 1 2 , 1 4 ) denoted as G trap , c = 2 3 leads to a filter based on the Simpson’s rule with coefficients ( 1 6 , 2 3 , 1 6 ) denoted as G simp , c = 1 2 leads to a even sharper filter with coefficients ( 1 8 , 3 4 , 1 8 ) denoted as G sharp . When the mesh is not uniform we apply locally a quadratic interpolation to a uniformly spaced points before applying the above formulas (Loh & Domaradzki, 1999), the procedure is described in Section 3.5. It was tested that an approximate deconvolution based inverse of orderN = 5 (denoted as Q 5 ) is sufficient to recover almost all the information in the range 31 k < 1 2 k c . The shapes of both primary and secondary filters using inverse filters of order 3 and 5 are shown in Fig. 4.1. Figure 4.1: Transfer functions for filters defined by equation (4.4). G trap : (red) solid line with star symbols,c = 1;G simp : (green) solid line with triangle symbols, c = 2 3 ; G sharp (blue) solid line with circle symbols, c = 1 2 ; and transfer functions for corresponding filters based on equation (4.2) for N = 3 and N = 5. In Fourier space, these three-point primary filters in 1-D have a general form for uniform grid ˆ G(k) = 1−c sin 2 k Δ 2 ! (4.5) Therefore all the primary filters in Fig. 4.1 are second order operators, and based on equation (4.2) the corresponding secondary filters are of 2(N +1)th order [ Q N G(k) = 1− " c sin 2 k Δ 2 !# N+1 . (4.6) These filters only weakly affect the large wavenumbers k for increasing the decon- volution operator order N. While all the secondary filters shown in Fig. 4.1 can be used to remove the energy between 1 2 k c and k c , in actual simulations it is not necessary to apply 32 them at each time step. This is because the energy accumulation near k c does not contaminate large scales immediately but gradually. In the TNS simulations of Domaradzki et al. (2002) (Section 2.3.2) the filtering period was determined manually for a given filter. Tantikul & Domaradzki (2010) developed an automatic filtering criterion based on the ratio of removed kinetic energy for two different filters. For a given energy spectrum E(k), after filtering by filters G 1 and G 2 E 1 (k) = ˆ G 2 1 E(k), (4.7) E 2 (k) = ˆ G 2 2 E(k), (4.8) where ˆ G is the transfer function of a physical space filter G. The energy loss due to filtering is I 1 = Z kc 0 (E−E 1 )dk, (4.9) I 2 = Z kc 0 (E−E 2 )dk. (4.10) TheratioI 1 /I 2 canbecomputedforagivenenergyspectrumE(k)andtwoselected filters ˆ G 1 and ˆ G 2 . For numerical simulations on under-resolved grids the energy will gradually accumulate in the small scales near the cutoff k c , changing the energy spectrum in that range and consequently the value of a computed ratio I 1 /I 2 . Therefore, the ratio I 1 /I 2 can be used as a diagnostic quantity indicating the energy accumulation near cutoff. Tantikul & Domaradzki (2010) computed this quantity using three theoretical spectral forms corresponding to the inertial, the Batchelor, and the dissipation range and found that for the particular filters chosen the ratio was confined to a narrow range for all spectral forms. The ratio can be monitored in actual simulations. When we choose G 2 to be a stronger 33 filter than G 1 thus I 1 /I 2 < 1, if the value exceeds values from the theoretically determined range, indicating an unphysical energy build-up, an explicit filtering is applied to attenuate small resolved scales. In the work of Tantikul & Domaradzki (2010), the energy ratio is based on kernel Q 5 G trap with filter width Δ (mesh size) and 2Δ, in which the filtering operation of width 2Δ is approximated as Q 5 G sharp on mesh size Δ assuming U n±1 = 1 2 (U n±2 +U n ). Theoretically, any two filters that have relatively low effects on wave numbers k < 1 2 k c can be used to compute the energy ratio. However, in order to make the energy ratio independent of the shape of the energy spectrum, we should avoid using too weak filters such as Q 5 G sharp . To see this, consider a Taylor’s expansion of 1D secondary filter [ Q N G [ Q N G(k) = 1− [c(kΔ/2) 2 ] N+1 +.... (4.11) which quickly converges to 1 for largeN whenc< 1 and makes the filtered energy spectrum close to the original E(k). Then the approximate energy ratio is I 1 I 2 = R {2[c 1 (kΔ/2) 2 ] N+1 E(k) +....}dk R {2[c 2 (kΔ/2) 2 ] M+1 E(k) +...}dk , (4.12) which should vary weakly on small changes of power k in E(k) given order N and M large enough to make the higher order terms are negligible. Thus the ratio can be regarded as a constant purely related to the filters. However, it was tested that when at least one of c 1 or c 2 were smaller than 1, the energy ratio was affected by computational errors and higher order terms. This is consistent with the observation that, in real simulations, the subtle energy difference after one application of a weak filter is more difficult to be computed accurately, making the result less robust. 34 In the present work we decided to use only the trapezoidal filter as our pri- mary filter so that c = 1. For 3D case, the primary filter G can be implemented sequentially in each direction as a product in the spectral space before applying a deconvolution operator Q N . Then the filtered transfer function is integrated over a sphere to be able to compare the filter strength with typical 1-D filters [ Q N G 3D (k) = 1 2π Z π 0 Z π 0 {1− (1− [1− sin 2 (k Δ 2 sin(φ) cos(θ))][1− sin 2 (k Δ 2 sin(φ) sin(θ))] [1− sin 2 (k Δ 2 cos(φ))]) N+1 } sin(φ)dφdθ. (4.13) To develop a filtering criterion, in this work we choose the same primary filter G trap with the deconvolution order N = 5 and M = 3, respectively, i.e. G 1 = Q 5 G trap andG 2 =Q 3 G trap . The energy ratiosI 1 /I 2 based on equations (4.7 - 4.10) for the dissipation (∼ k −5/3 exp(−5ηk)), the inertial (∼ k −5/3 ) and the Batchelor spectrum (∼k −1 ) are summarized in Table. 4.1. Table 4.1: Filtered energy ratios for different energy spectrum. energy spectrum Dissipation Inertial Batchelor I 1 0.00043738 0.00579159 0.29930728 I 2 0.00055926 0.00809234 0.40573954 I 1 /I 2 0.7820691628 0.7156879222 0.7376832931 In practice, an easy way to implement the filter (Jeanmart & Winckelmans, 2007) is to compute the coefficients of the geometric sequence in equation (4.2) Q N G =I− (I−G) N+1 =a 1 G +a 2 G 2 +... +a N+1 G N+1 . (4.14) Coefficients a i can be calculated recursively, thus the filtering operations can be implemented easily and efficiently also in the physical space. 35 In actual simulations, the energy loss integrals I 1 and I 2 are computed using the filtered velocities in the physical space ¯ u i,1 =G 1 ∗ ¯ u i and ¯ u i,2 =G 2 ∗ ¯ u i , giving the energy ratio I 1 I 2 = RRR P 3 i=1 (¯ u i − ¯ u i,1 )(¯ u i − ¯ u i,1 ) RRR P 3 i=1 (¯ u i − ¯ u i,2 )(¯ u i − ¯ u i,2 ) . (4.15) While developing the filtering criterion for this investigation we repeated calcu- lations for the ratio I 1 /I 2 performed in Tantikul & Domaradzki (2010) and found an error. The filter kernel used by Tantikul & Domaradzki (2010) in a form equa- tion (4.13), i.e., the product of 1-D filters averaged over angles φ and θ, had a square root applied to the equivalent integral on the r.h.s. of the equation. This leads to an erroneous value of the ratioI 1 /I 2 = 0.007− 0.010 reported in Tantikul & Domaradzki (2010) rather than the value of around 0.1 computed with the for- mula without a square root for the inertial range spectrum. As a result the filtering period in that work should have been longer. Nevertheless no major changes in the results are expected because as observed in Domaradzki et al. (2002) changing the filtering period even by a factor of two led to negligible changes in simulation results. 4.2 Adaptive filtering criterion For the current pair of filters in Table. 4.1 we choose energy ratio 0.7 as a theoretical benchmark to develop a filtering criterion. However, if we regard 0.7 as an upper limit and filter whenever the energy ratio exceeds this value like the criterion used in Tantikul & Domaradzki (2010), it would be necessary to compute the energy ratio every few time steps, with the corresponding interval determined empirically. As the filtering and volume integral in energy ratio calculation can be computationally expensive, to reduce the computation time, we decided to seek 36 an optimal fixed filtering period when the statistically steady state is reached in simulations. The idea is motivated by the fact that, the simulation accuracy in general is not very sensitive to the filtering interval (Domaradzki et al., 2002). Besides, even if we monitor the energy ratio at every time step, the filtering period varies very little in the statistically steady state. So rather than monitoring the energy ratio every one or several time steps to determine if filtering should be applied, we can instead start with a prescribed filtering period and allow a fixed relative error, e.g., 3% of the benchmark value 0.7. The corresponding pseudo-code is as follows er_theoritical = 0.7 tolerance = 0.03 call er = average(energy_ratio) after filtered n times ! averaged energy ratio if er < (1 − tolerance )∗er_theoritical filter_interval += 1 else if er > (1 + tolerance )∗er_theoritical filter_interval −= 1 else optimal_filter_interval = filter_interval The process is then repeated until an optimal filtering interval is found, which takes at most hundreds or thousands of time steps for the present cases, to keep that energy ratio in a narrow range when filtering is applied. The method can be easily applied to a statistically steady UDNS run, either during the simulation, or separately based on the output velocity fields as a diagnostic function, to signifi- cantly improve accuracy with very low computational cost. For the former case, 37 it is observed that once this optimal period is found, it never changes in the stati- cally steady state. Therefore, we can continue simulations using the fixed filtering period without needing to compute the energy ratio. For simulations with a filter- ing interval of 10 time steps, the computational cost of the current ILES approach is only 1.04 times the cost of UDNS with the same mesh, and much cheaper than all the explicit SGS models we tested. This is because the computational cost of periodic filtering is much less than performing more complex calculations of SGS stresses at each time step in traditional LES. In principle, any filter that only affects large wavenumbers near the cutoff k c can be adopted to add numerical dissipation in ILES. However, because the energy pileup does not occur exclusively in a small scales range, smoother filters are more favorable than the spectral cutoff filters. Different filtering interval can be deter- mined for different filters, but the resulting time-averaged numerical dissipation should be comparable, with stronger filters applied less frequently. Therefore, for the same adaptive filtering criterion, ideally different filters should produce com- parable results under the same filtering criterion. However, we found that this is true only to some extent, in particular, if a filter is too strong, it may remove excessive information and deteriorate the results. In the present work, most of the results are from a relatively weak filterQ 5 G sharp , which is treated as a benchmark. The investigation of different filters Q 5 G sharp , Q 5 G simp and Q 5 G trap (see Fig. 4.1) is reported later in Section 4.3.5. Note that the hyperparameters in the above filtering criterion are designed for the weakest filter we tested (Q 5 G sharp ), leading to an optimal filtering interval around10timestepsformostofoursimulations. Whenastrongerfilterisadopted, it is advisable (but not necessary) to reduce the tolerance, the step size and the 38 number of samples for averaging to find a narrower range of a filtering interval in a shorter time. 4.3 Results The Fourier-Chebyshev pseudo-spectral solver introduced in Chapter 3 is adopted to solve the incompressible Navier-Stokes equations, whose spectral accu- racy in general leads to negligible numerical dissipation, despite possible effects from aliasing errors when grid resolution is relatively low. In order to to quantify the numerical dissipation added by purely explicit filtering based on the criterion in Section 4.2, no other sources of numerical dissipation were applied. Therefore, the current adaptive filter not only adds the numerical dissipation as a surrogate for explicit SGS models, but also serves to stabilize the simulations and reduce aliasing errors. The amount of the numerical dissipation can be easily controlled by changing the filtering period. In Section 4.3.1 we show that the numerical dissipation introduced by filtering is linearly dependent on the filtering interval. Moreover, when the additional numerical dissipation is close to the sub-grid scale (SGS) dissipation in an explicit LES, the overall accuracy is also comparable, indi- cating that simple filtering can potentially be used to replace explicit SGS models. A posteriori tests of the adaptive filter are shown in Section 4.3.3 for different Reynolds numbers by comparing with DNS as well as classic explicit LES models (static Smagorinksy (Smagorinsky, 1963) with the van Driest near wall damping function (Van Driest, 1956) for C s = 0.1, dynamic Smagorinsky (Germano et al., 1991; Lilly, 1992) and Truncated Navier-Stokes (TNS) model (Domaradzki et al., 2002), see Section 2.3). The approach works well even when a grid resolution is 39 further reduced, as tested in Section 4.3.4. As discussed earlier, the adaptive fil- tering criterion for different filters should produce similar results as long as the filter is not too strong and removes energy mostly near the cutoff wavenumber. This is confirmed by comparison of results obtained for different filters discussed in Section 4.3.5. On the other hand, unlike spectral method, lower order methods such as finite differences already have a relatively large amount of the numerical dissipation and its effects on the energy transfer cannot be neglected. To be able to address dissipative properties of low order numerical schemes in the pseudo- spectral environment we use the modified wavenumber concept (Lele, 1992). The numerical dissipation analysis of simulations performed with lower order schemes can be found in Section 4.3.6. In the plots in the following sections the label "Re*" is associated with the current ILES approach. The parameters for all simulations are summarized in Table. 4.2 and parameters for reference DNS/TNS data are listed in Table. 4.3. The benchmark DNS data is obtained from Moser et al. (1999) for Re τ = 180, 590 and Lee & Moser (2015) for Re τ = 1000, 2000, and the TNS results are from Domaradzki et al. (2002) where the filtering interval was set manually. Lower orderstatisticsforthemeanvelocityprofileU + , rmsvelocitiesu + rms ,v + rms ,w + rms and the Reynolds stress R + uv are compared with the reference data following approach described in Section 3.6. The rms velocities are computed as the square root of the diagonal terms of the Reynolds stress tensor R + ii for i = 1, 2, 3, respectively. 4.3.1 Numerical dissipation introduced by explicit filters To check if the secondary filters in Section 4.1 can effectively add numerical dissipation for different conditions, a numerical dissipation analysis of a series of UDNS/LES simulations was performed. To achieve more generality, the approach 40 Table 4.2: Current LES parameters. Case Domain Resolution Δx + × Δy + w × Δz + N/N DNS (%) Re180h 4π× 2× 1.5π 64× 65× 64 36.0× 0.22× 13.5 12.6 Re180 4π× 2× 1.5π 48× 49× 48 46.0× 0.38× 17.3 5.23 Re180l 2π× 2×π 32× 33× 32 36.8× 0.90× 18.4 1.6 Re590 2π× 2×π 96× 97× 96 39.9× 0.33× 19.9 2.36 Re590l 2π× 2×π 64× 65× 64 59.8× 0.74× 29.8 0.702 Re1000h 2.5π× 2× 0.5π 256× 97× 128 30.9× 0.54× 16.5 0.132 Re1000 2.5π× 2× 0.5π 128× 97× 128 60.9× 0.53× 16.3 0.066 Re1000l 2.5π× 2× 0.5π 96× 65× 128 77.3× 1.10× 22.8 0.033 Re2000 2.5π× 2× 0.5π 256× 131× 128 61.2× 0.58× 24.1 0.044 Table 4.3: Benchmark simulation parameters. Case Re τ Domain Resolution DNS180 180 4π× 2× 1.5π 128× 129× 128 DNS590 590 2π× 2×π 384× 257× 384 DNS1000 1000 8π× 2× 3π 2304× 512× 2048 DNS2000 2000 8π× 2× 3π 4096× 768× 3072 TNSZ1 180 4π× 2× 2π 64× 65× 64 TNSZ2 1000 2.5π× 2× 0.5π 96× 65× 128 Table 4.4: Relative error and numerical dissipation for Re τ = 590. Error (comparison with DNS) no model static Smagorinsky C s =0.1 dynamic Smagorinsky U + (%) 5.817805 1.094381 1.193594 R + uu (%) 15.40269 6.201147 5.82024 R + vv (%) 1.966532 3.419877 2.298193 R + ww (%) 5.341539 2.887433 2.101117 R + uv (%) 1.340966 1.329854 0.6068324 Average (%) 5.973906 2.986538 2.403995 Numerical Dissipation E ν 3.697E-003 2.476E-003 2.673E-003 E SGS 0.000E+000 4.971E-004 4.400E-004 ν SGS 0.000E+000 1.807E-005 1.482E-005 E n -3.143E-005 3.310E-005 1.636E-005 ν n 2.524E-007 1.201E-006 5.498E-007 RatioE n /E ν (%) 0.280480 1.334613 0.6108837 was first tested for a moderate Reynolds numberRe τ = 590 with an adequate grid resolution 96× 97× 96. The corresponding DNS resolution (Moser et al., 1999) 41 Figure 4.2: Mean velocity profiles: Re τ = 590, Mesh: 96× 97× 96. is 384× 257× 384, so the current resolution is roughly 2% of the DNS bench- mark. Results for simulations without explicit filtering are listed in Table. 4.4. The relative errors and the numerical dissipation analysis are based on terms in equations (3.37) and (3.27) and corresponding plots are shown in Fig. 4.2 and 4.3. For eddy viscosity models like the Smagorinsky models, only a term for the explicit SGS tensor τ ij is added, thus we are solving the Navier-Stokes equations with a very similar form as for a no-model case. Therefore, the value of the numerical dissipation is similar among all cases, though contributions from each term can be slightlydifferent. Wecanconcludethatforasufficientlyhighresolutionthenumer- ical dissipation is mostly determined by grid resolution irrespective of models used. Note that all terms are time averaged, including the ratioE n /E ν . Thus when the numerical dissipation is negligibly small and oscillates around zero, it is possible forE n andE n /E ν to have different signs. In the current cases, results are relatively accurate when explicit SGS models are used, whereas the no-model UDNS is much less satisfactory. Since the numerical dissipation is negligible compared with the viscous dissipation, either an explicit or implicit LES method would be beneficial. 42 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.3: RMS velocities and Reynolds stresses: Re τ = 590, Mesh: 96× 97× 96. Table 4.5: Relative error and numerical dissipation for different filtering intervals Filtering interval (time steps) 2 4 6 8 10 12 Filtering period (t + ) 6.521E-004 1.304E-003 1.956E-003 2.608E-003 3.260E-003 3.912E-003 Error (comparison with DNS) U + (%) 2.906691 0.5320195 0.414884 1.311051 1.638773 2.544442 R + uu (%) 13.03592 7.800148 7.130713 7.087964 7.70797 8.715973 R + vv (%) 1.685311 1.527238 1.124426 1.30913 1.058724 0.7714084 R + ww (%) 1.207151 1.746765 1.941009 2.131751 2.310026 2.523582 R + uv (%) 0.4383451 0.3256939 0.5087641 0.5962698 0.5486774 0.7617005 Average (%) 3.854684 2.386373 2.223959 2.487233 2.652834 3.063421 Numerical Dissipation E ν 2.409E-003 2.522E-003 2.620E-003 2.710E-003 2.786E-003 2.851E-003 E n 4.715E-004 4.491E-004 4.294E-004 4.100E-004 3.887E-004 3.703E-004 ν n 1.761E-005 1.602E-005 1.474E-005 1.360E-005 1.254E-005 1.167E-005 RatioE n /E ν (%) 19.56663 17.79539 16.37606 15.10800 13.93206 12.96244 43 Figure 4.4: Numerical dissipation dependence on the filtering interval. In Table. 4.5 we present results for the same conditions for simulations with filtering. A relatively weak secondary filter Q 5 G sharp was used and a filtering interval was set manually to several time steps (a corresponding filtering period is nondimensionalized in wall time units t + = tu τ /h). It can be observed that the amount of additional numerical dissipation introduced by the same filter is almost linearly dependent on the filtering period. In the present case, when the interval decreasesby2timesteps, thenumericaldissipationincreasesbyaconstantamount around 2.0× 10 −5 (see Fig. 4.4). This makes the current filter a good candidate for introducing artificial dissipation. This filter has similar strengths as the weak filter used in Cadieux & Domaradzki (2016). The overall error reaches a minimum for the filtering interval of 6, and it is slightly less than the errors for the explicit SGS models in Table. 4.4. At this filtering interval, as the numerical viscosity ν n provided by the filtering is close to the effective SGS viscosity ν sgs for the dynamic Smagorinsky model their overall simulation accuracy is also comparable. This further confirms that the numerical dissipation added by our explicit filtering 44 approach can serve as a reasonable surrogate for explicit SGS models. For a non- dissipative solver there should exist a global minimum of the error for some value of the filtering period, which can be regarded as the optimal filtering period for the current setup in terms of the sub-grid scale dissipation provided by the filter. The optimal filtering interval found using the criterion introduced in Section 4.2 is 7. The corresponding results are labeled as Re590 in Fig. 4.2 and 4.3 and are very accurate. Therefore we can conclude that the adaptive filtering approach did a decent job in finding the optimal period. 45 4.3.2 One-dimensional energy spectrum The goal of the periodic filtering is to remove the unphysical energy pileup near the LES cutoff, which can be verified by checking the energy spectrum. In this section, the one-dimensional energy spectrum for relatively high Reynolds number and low grid resolution test cases are presented, so that the cutoff wavenumber is located at the inertial range. All energy spectra in this section is streamwise spectra (x-direction) at the center of the channel, but results at other locations show similar trends. Fig. 4.5 shows results forRe τ = 1000 with very low grid resolution (corresponds to case Re1000l in Table 4.2). The relatively weak filter Q 5 G sharp is still adopted in the present section, and “F*” in the label of the plots indicates filtering interval. Apparentlytheno-modelUDNSandthesimulationwiththedynamicSmagorinsky model suffers from excessive energy near the LES cutoff, which contaminates the energy spectra for moderate wavenumbers as well. As a result, even the mean velocity profile cannot be accurately predicted. On the other hand, the application of this weak filter at low frequency were already able to effectively remove the energy pileup. However, if the frequency is too high, for instance the “F2” cases, energy at moderated wavenumbers might be affected. Recall that the coarser grid used in LES is interpreted as a filtering operation. Therefore, if low-pass filtering is applied almost every time step, it is similar to using a lower resolution grid, which reduces the accuracy as well. For the optimal filtering interval determined by the energy spectrum in the inertial range (“F13”), the energy spectra indeed follows closely to the -5/3 law, so the mean velocity profile show excellent agreement with the DNS benchmark. For even larger filtering interval (“F25”), the results are satisfactory as well, indicating that in the present case the accuracy is not very sensitive to the filtering period. This also justifies the use of a range of energy 46 ratio, instead of necessarily try to find the best filtering interval at every time step. If we further check the the Reynolds stresses in Fig. 4.6, the results of the optimal filtering interval determined by our adaptive filtering approach is still the most accurate. (a) U + (b) E uu (c) E vv (d) E ww Figure 4.5: Mean velocity profile and streamwise energy spectra: Re τ = 1000, Mesh: 96× 65× 128. We can draw similar conclusions based on the results for Re τ = 2000 case in Fig. 4.7. The no model simulation (UDNS) and overly frequent filtering (“F10”) still leads to energy pileup and excessive energy removal, respectively. The dynam- ics Smagorinsky model in this case were able to alleviate the energy accumulation a little bit, in comparison with UDNS, so its accuracy improved. Note that since 47 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure4.6: RMSvelocitiesandReynoldsstresses: Re τ = 1000, Mesh: 96×65×128. the optimal filtering interval 34 is already relatively large, further increasing the filtering period does not have very significant impact to the energy spectra, and the optimal filtering periodic still leads to more accurate results. Therefore, we believe thattheadaptivefilteringapproachisabletocorrectlyreproducethedesiredshape of energy spectrum, the unphysical energy accumulation due to under-resolved mesh is removed with little influence to lower wavenumbers, so the predication of lower order statistics can be significantly improved. 48 (a) U + (b) E uu (c) E vv (d) E ww Figure 4.7: Mean velocity profile and streamwise energy spectra: Re τ = 2000, Mesh: 256× 131× 128. 49 4.3.3 Different Reynolds number In this section, results for the friction Reynolds number Re τ = 180, 1000, and 2000 are compared with the benchmarks in Table. 4.3. Results for the dynamic Smagorinsky model (sometimes the static Smagorinsky model with C s = 0.1 as well) obtained by the same code are also used to evaluate the performance of our new method. The Re τ = 180 case is simulated on a fine grid (Fig. 4.8 and 4.9, denoted as Re180h in Table. 4.2), so results with either explicit or implicit LES meth- ods are rather satisfactory and in general offer visible improvements to UDNS. A closer observation confirms that, when the grid resolution is adequate, the filtering approach in implicit LES (Re180h, TNSZ1) gives satisfactory results, outperform- ing the classic dynamic Smagorinsky model. Figure 4.8: Mean velocity profiles: Re τ = 180, Mesh: 64× 65× 64. For a higher Reynolds number Re τ = 1000 (Fig. 4.10 and 4.11, Re1000h in Table. 4.2), the ratio of LES and benchmark DNS resolution decreases. The rel- atively lower resolution provides more severe conditions for comparison among different methods. Again, the implicit LES Re1000h and TNSZ2 more faithfully 50 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.9: RMS velocities and Reynolds stresses: Re τ = 180, Mesh: 64× 65× 64. reproduced the mean velocity profile than other methods. The Reynolds stress results from our new approach are also the most consistent with DNS benchmarks, except for a small under-prediction of w + rms near the wall. On the other hand, in this relatively low resolution case, the improvements we get using the explicit SGS models are less significant. Therefore, the filtering seems to offer a better capability in lower resolution LES. We then further increase Re τ to 2000 but only slightly increase the number of grid points in the vertical direction, (case Re2000 in Table. 4.2, plots in Fig. 4.12 and 4.13). In this case the overall accuracy of our approach and the dynamic 51 Figure 4.10: Mean velocity profiles: Re τ = 1000, Mesh: 256× 97× 128. Smagorinsky model is similar. Both methods predicts the mean velocity profile U + (y + ) accurately as opposed to the poor performance of UDNS. The adaptive filter in the current case seems to damp the fluctuation velocity u 0 too strongly. This makes the profile ofu + rms andR + uw to deviate from the DNS benchmark away from the wall, while the dynamic Smagorinsky model shows good agreement with these two quantities. However, it under-predicts v + rms and w + rms . Recall that our filtering operation is global, theoretically it only affects very small scales which are expected to be almost homogeneous. However, when we have more inhomogeneous conditions in the vertical direction for high Reynolds number with low resolution, the homogeneity assumption may not hold. To improve results a further mesh refinement would be beneficial for this inhomogeneous direction. In general for simulations at various Reynolds numbers, the present method agreeswellwithDNSdataandoftenoutperformsthedynamicSmagorinskymodel. The robustness of the method on coarser meshes is investigated in the next section. 52 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.11: RMS velocities and Reynolds stresses: Re τ = 1000, Mesh: 256× 97× 128. 53 Figure 4.12: Mean velocity profiles: Re τ = 2000, Mesh: 256× 131× 128. 54 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.13: RMS velocities and Reynolds stresses: Re τ = 2000, Mesh: 256× 131× 128. 55 4.3.4 Performance for coarse grid resolutions The major objective of LES is to accurately predict the low order statistics of turbulent flows with the smallest number of grid points possible. To further evaluate the performance of the adaptive filtering we decreased the resolution for Re τ = 180 and 1000 cases. Figure 4.14: Mean velocity profiles: Re τ = 180, different resolutions. Simulations cases Re180h, Re180 and Re180l in Table. 4.2 at Re τ = 180 were performed with decreasing grid resolutions. Note that while we use the same computational domain as DNS reference (Moser et al., 1999) for the first two, in the last case the computational domain is somewhat smaller. For this Re180l case, although our no model simulation becomes unstable due to the lack of any stabilization techniques, we were able to get reasonable results by turning on the adaptive filtering, so the filter by itself is capable of stabilizing the simulation. For the mean velocity profiles on Fig. 4.14, all results are almost on the top of each other, indicating the adaptive filter was able to adjust properly to the different resolutions. For the Reynolds stresses in Fig. 4.15, it can be observed that the 56 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.15: RMS velocities and Reynolds stresses: Re τ = 180, different resolu- tions. highest resolution case Re180h in general still captured best the DNS profiles, but the method performed reasonably well also at lower resolutions. When Re τ increases to 1000, all the previous observations still hold. In gen- eral all the results are still satisfactory, simulations on coarser meshes Re1000 and Re1000l can capture the mean velocity profile almost as well as simulations on the finer mesh Re1000h (Fig. 4.16), but larger differences can be seen for Reynolds stresses (Fig. 4.17). This is expected as the corresponding UDNS becomes less accurate when a grid resolution decreases and more under-resolved scales are 57 Figure 4.16: Mean velocity profiles: Re τ = 1000, different resolutions. ignored. The filtering operations can only remove their negative effects but will not recover them. 58 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.17: RMS velocities and Reynolds stresses: Re τ = 1000, different resolu- tions. 59 4.3.5 Comparison of different filters According to the analysis in Section 4.1 the energy ratio can be computed in principle for any two filters. Therefore, it would be straightforward to build a new filtering criterion based on the energy ratio of any reasonable pair of filters. Moreover, one may expect that the filtering criterion should work in the same manner for different pairs of filters, as long as these filters substantially damp the small scales without noticeable impact on the large scales. In other words, the filtering approach described here should be largely independent of the actual filters employed. In the previous sections we only used a weak filter Q 5 G sharp . In this section we investigate the capability of two other filters,Q 5 G simp andQ 5 G trap , shown in Fig. 4.1. In all cases the filters were applied to the same no model simulation with the identical, fixed time step, to see if the filtering criterion can find the optimal filtering interval based on the strength of different filters. Table 4.6: Effect of different filters for Re τ = 180, Mesh: 64× 65× 64. Filter Q 5 G trap Q 5 G simp Q 5 G sharp Filtering interval (time steps) 57 21 8 Filtering period (t + ) 5.806E-002 2.139E-002 8.149E-003 Error (comparison with DNS) U + (%) 0.9951294 0.5417327 1.201423 R + uu (%) 11.71444 5.108443 3.698822 R + vv (%) 0.8341516 1.076409 0.3762055 R + ww (%) 1.14583 1.208954 1.400624 R + uv (%) 0.6271194 0.3880678 0.5955131 Average (%) 3.063335 1.664721 1.454518 Numerical Dissipation E ν 3.989E-003 4.113E-003 4.237E-003 E n 3.259E-004 3.120E-004 2.896E-004 ν n 2.947E-005 2.728E-005 2.458E-005 RatioE n /E ν (%) 8.185376 7.57844 6.828093 60 Figure 4.18: Mean velocity profiles: Re τ = 180, different filters. Results for the case Re180h are shown in Table. 4.6 and Fig. 4.18, and 4.19. Apparently, with the filter strength increasing from Q 5 G sharp ,Q 5 G simp toQ 5 G trap (see Fig. 4.1), i.e. more energy removed by a single application of a filter, cor- responding optimal filtering interval becomes larger, as expected. In principle, filtering less frequently is more desirable to slightly improve the efficiency, but with the strongest filter Q 5 G trap we tested, which is also the one used in Tantikul &Domaradzki(2010), resultsarelesssatisfactoryincomparisonwiththeothertwo filters, while simulations using filters Q 5 G sharp and Q 5 G simp have a similar overall accuracy. Therefore in the present cases, a relatively weak filter with more frequent filtering would be a safer choice. When the optimal filtering interval approaches unity, it would be equivalent to an LES model with the SGS term replaced by a filtering operation at each time step. On the other hand, although the filter such as Q 5 G trap might be better theoretically since it vanishes at k c , if we completely remove all the energy near the cutoff, it affects the energy cascade process more significantly and may take a long time for flows to recover. 61 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.19: RMS velocities and Reynolds stresses: Re τ = 180, different filters. In terms of the numerical dissipation, all filters generate comparable amount of numerical dissipation averaged over time, with slightly larger values for stronger filters despite lower filtering frequency, which is caused by the excessive energy they can potentially remove. The same analysis was performed for the case Re590 (Table. 4.7 and Fig. 4.20 and 4.21). As the relative grid resolution of the current case is smaller than the previous Re180h case (see Table. 4.2), for the same filter, the filtering period decreases compared with results in Table. 4.6, implying that the filter needs to be applied more frequently to remove larger amount of energy accumulation near 62 Table 4.7: Effect of different filters for Re τ = 590, Mesh: 96× 97× 96. Filter Q 5 G trap Q 5 G simp Q 5 G sharp Filtering interval 63 28 7 Filtering period (t + ) 2.054E-002 9.129E-003 2.282E-003 Error (comparison with DNS) U + (%) 3.172722 0.2720504 1.16982 R + uu (%) 13.49817 4.727628 5.948215 R + vv (%) 1.550555 1.137231 0.8470457 R + ww (%) 1.582716 2.064209 2.29037 R + uv (%) 0.4500272 0.2977598 0.6291394 Average (%) 4.050838 1.699776 2.176918 Numerical Dissipation E ν 2.463E-003 2.646E-003 2.672E-003 E n 4.277E-004 4.201E-004 4.195E-004 ν n 1.567E-005 1.423E-005 1.411E-005 RatioE n /E ν (%) 17.4059 15.81157 15.68028 k c . For Re τ = 590, the strongest filter Q 5 G trap is still the worst, but results from Q 5 G simp is actually more accurate than Q 5 G sharp , with very similar numerical dissipation. Recall that for the Re180h case, Q 5 G simp also better predicted the mean velocity profile. We may thus conclude that Q 5 G simp can overall be as good as the weaker filterQ 5 G sharp , especially when the amount of numerical dissipation added by filtering is comparable. Of course filters may perform differently under other flow conditions, but in general one may assume relatively weak filters like Q 5 G sharp and Q 5 G simp to have comparable accuracy, but when a filter becomes stronger like Q 5 G trap , results might deteriorate. 63 Figure 4.20: Mean velocity profiles: Re τ = 590, different filters. (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.21: RMS velocities and Reynolds stresses: Re τ = 590, different filters. 64 4.3.6 Investigation of lower order schemes All previous results were obtained using the pseudo-spectral solver, whose spec- tral accuracy effectively eliminates the numerical dissipation, so the numerical dis- sipation is primarily due to the filters. However, as the amount of native numerical dissipation in a CFD code is also related to the order of the solver, it is of interest to check how the method performs for lower order schemes where the numeri- cal dissipation is already not negligible. For (compact) finite difference method (FDM), the error of the schemes can be evaluated using a modified wavenumber analysis in the Fourier space, in which the modified wavenumbers k 0 and k 00 are compared with wavenumbers in the spectral method. On the other hand, for a pseudo-spectral code, if we replace the wavenumbers k by the modified wavenum- bers, it will be equivalent to applying FDM schemes corresponding to a given modified wavenumber. The relation is exact, which offers an easy way to evaluate different low order schemes within the framework of a spectral scheme. Here we changed the numerical scheme in the horizontal directions (x and z) to common finite difference schemes with different orders of accuracy by replacing horizontal wavenumbers k x and k z in the code by modified wavenumbers. The discretiza- tion in the vertical direction (y) is based on Chebyshev polynomials and remains unchanged. Pure central difference schemes are tested in this section, which is feasible for the periodic boundary conditions in horizontal directions. For the 1st order derivatives, the 5 term central difference equation for compact difference schemes up to 10th order is generalized as follows (Lele, 1992) β(f 0 i−2 +f 0 i+2 ) +α(f 0 i−1 +f 0 i+1 ) +f 0 i =c f i+3 −f i−3 6Δ +b f i+2 −f i−2 4Δ +a f i+1 −f i−1 2Δ (4.16) in which Δ =L/(N− 1) is uniform mesh size 65 Figure 4.22: Modified wavenumber for 1st order derivatives. In this section, order 2, 4, 6, 8 and 10 based on equation (4.16) are tested. For each order, we use the smallest number of non-zero coefficients to form the most common schemes (order 2 and 4 will become purely central difference schemes). The corresponding coefficients can be found in or derived from Lele (1992), and the modified wavenumber k 0 is given by the formula: k 0 Δ = a sin(kΔ) + (b/2) sin(2kΔ) + (c/3) sin(3kΔ) 1 + 2α cos(kΔ) + 2β cos(2kΔ) , (4.17) in which the wavenumber k = 0, 2π/L, 4π/L,..., and kΔ∈ [0,π]. The modified wavenumbers for the schemes tested in this paper are shown in Fig. 4.22 Similarly, for the 2nd order derivatives β(f 00 i−2 +f 00 i+2 )+α(f 00 i−1 +f 00 i+1 )+f 00 i =c f i−3 − 2f i +f i+3 9Δ 2 +b f i−2 − 2f i +f i+2 4Δ 2 +a f i−1 − 2f i +f i+1 Δ 2 . (4.18) and the modified wavenumber k 00 is given by the formula: k 00 Δ 2 = 2a(1− cos(kΔ)) + (b/2)(1− cos(2kΔ)) + (2c/9)(1− cos(3kΔ)) 1 + 2α cos(kΔ) + 2β cos(2kΔ) , (4.19) 66 Figure 4.23: Modified wavenumber for 2nd order derivatives. and is shown Fig. 4.23. The approach is first investigated for the case Re180h and results from the modified solvers with different orders of accuracy are shown in Table. 4.8 and Fig. 4.24 and 4.25. The numerical dissipation rate is computed using the same numerical method as used by the Navier-Stokes solver itself (i.e., with modified wavenumbers) for self-consistency. As seen in Table. 4.8 the numerical dissipation increasesmonotonicallyastheorderofaschemedecreases. However,theadditional numerical dissipation for the lower order schemes in general does not improve the overall accuracy. Table 4.8: Different orders of accuracy for Re τ = 180, no filter case, Mesh: 64× 65× 64. order or accuracy spectral o10 o8 o6 o4 o2 Error (comparison with DNS) U + (%) 4.781803 3.846254 4.657559 4.637896 3.800464 1.209274 R + uu (%) 18.09592 15.93359 18.09063 16.1243 12.12179 18.93923 R + vv (%) 1.035923 0.6763259 0.8043295 0.5818885 1.07008 2.71913 R + ww (%) 4.043432 3.013137 3.381043 2.311485 2.504524 4.335509 R + uv (%) 1.226828 1.28007 1.622934 1.728694 1.174829 0.3975939 Average (%) 5.83678 4.949876 5.711299 5.076852 4.134338 5.520148 Numerical Dissipation RatioE n /E ν (%) 1.267855 2.597434 2.811173 3.382733 6.144547 8.751022 67 Figure 4.24: Mean velocity profiles: Re τ = 180, no filter, different orders of accu- racy. In Table. 4.9 and Fig. 4.26 and 4.27 we show results of implicit LES obtained using the adaptive filtering approach for all cases discussed above. Although the current approach in general can still make some improvements to corresponding UDNS simulations, the results are not as accurate as those obtained using the spectral code. For finite difference schemes, especially for lower orders 2 and 4, the inherent numerical dissipation is already larger than the one found in the spectral method with the adaptive filtering, so adding more numerical dissipation by filter- ingcannolongerbeexpectedtobeeffective. Asallthelowerorderschemesalready have large numerical dissipation, the explicit filtering is triggered less frequently, thus less dissipation is added through filtering. However, in analogy to the rela- tively strong filters we tested in Section 4.3.5, the 1st order derivatives for FDM all vanish atk =k c , therefore if we consider their effects as equivalent filters, they may still be too strong and lead to the loss of information in high wavenumbers. On the other hand, the numerical dissipation itself does not provide information about the distribution of dissipative effects in different range of wavenumbers. However, for 68 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.25: RMS velocities and Reynolds stresses: Re τ = 180, no filter, different orders of accuracy. FDM schemes, the energy spectrum would gradually go to 0 rather than exhibit a sharp drop at k c , making it unlikely to have a severe energy accumulation close tok c . Therefore, for numerical schemes with low order of accuracy, the filtering is expected to mainly smooth the profile to improve accuracy. The approach was subsequently tested on the case Re1000l, which possesses the lowest grid resolution relative to the DNS benchmark. For no filter simula- tions (Table. 4.10 and Fig. 4.29,4.30) most of the observations are consistent with the previous high resolution case. A notable issue is the numerical dissipation for order 2, which became negative when we use the same numerical scheme as 69 Table 4.9: Different orders of accuracy for Re τ = 180, filtered case, Mesh: 64× 65× 64. order or accuracy spectral o10 o8 o6 o4 o2 Filtering interval 8 23 35 52 53 54 Error (comparison with DNS) U + (%) 1.201423 2.565452 3.091512 3.782421 3.546077 1.014171 R + uu (%) 3.698822 10.7382 12.35304 16.5383 14.55257 16.85912 R + vv (%) 0.3762055 0.424793 0.4387869 0.4817185 1.348159 2.695042 R + ww (%) 1.400624 2.002394 2.369451 1.549874 2.52258 4.496868 R + uv (%) 0.5955131 0.9052288 1.158041 1.196133 0.8972027 0.3227613 Average (%) 1.454518 3.327214 3.882166 4.70969 4.573318 5.077591 Numerical Dissipation RatioE n /E ν (%) 6.828093 5.270255 4.794535 4.603238 6.817811 9.247348 Figure 4.26: Mean velocity profiles: Re τ = 180, filtered with optimal period, different orders of accuracy. in the Navier-Stokes solver to compute it. A negative time averaged numerical dissipation theoretically indicates the inverse energy transfer, which is inconsistent with the physics of turbulence and can make the simulations unstable. Based on the Reynolds stresses plot in Fig. 4.30, we can conclude that order 2 is too low for the current Reynolds number and grid resolution, and since it cannot produce reliable results the numerical dissipation calculation should not be trusted either. Previous analysis of Schranner et al. (2015) shows that the numerical dissipation analysis for a reasonable simulation is independent of the order of corresponding 70 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.27: RMS velocities and Reynolds stresses: Re τ = 180, filtered with opti- mal period, different orders of accuracy. numerical schemes, however, when the simulation itself is suspicious, it’s not safe to apply the same schemes to compute it. Therefore, if the solver is of a low order and results are not accurate, it is advisable to use an extra post-processing code based on higher order schemes for the numerical dissipation analysis. In Fig. 4.28, the numerical dissipation ratioE n /E ν computed with the original pseudo-spectral method is compared with values obtained in a self-contained calculations. Though the values for high orders cases were in a good agreement with the spectral calcu- lations, the dramatic difference in sign is observed for the 2nd order case. 71 Table 4.10: Different orders of accuracy for Re τ = 1000, no filter case, Mesh: 96× 65× 128. order or accuracy spectral o10 o8 o6 o4 o2 Error (comparison with DNS) U + (%) 3.523027 1.292342 2.488192 3.319122 6.446266 4.581886 R + uu (%) 10.49554 10.99817 15.38092 14.92271 26.06426 38.3471 R + vv (%) 4.019252 5.238255 3.589799 5.553783 6.689112 14.85287 R + ww (%) 6.221551 6.877904 4.180001 8.271229 10.15309 18.22385 R + uv (%) 0.7948734 1.286577 0.8907049 1.125829 1.074343 0.4613919 Average (%) 5.01085 5.13865 5.305924 6.638534 10.08541 15.29342 Numerical Dissipation (by FDM) RatioE n /E ν (%) 7.344825 11.63679 11.85087 12.02124 15.10917 -8.44594 Figure 4.28: Numerical dissipation computed from the spectral method and the FDM schemes in the solver: Re τ = 1000, no filter. When the adaptive filter is applied (Table. 4.11 and Fig. 4.31, 4.32), we observe the increase of the filtering interval and reduced accuracy for lower order schemes, as expected, but differently than in the previous cases, the numerical dissipation of FDM schemes decreases when the explicit filter is applied. If explicit filters can indeed reduce the intrinsic numerical dissipation of a scheme if it is large, it could potentially be useful. However, because information is lost due to the sharp cutoffs for the modified wave number approach, it is almost impossible to make FDM as accurate as spectral method for the same low resolution with the current 72 Figure 4.29: Mean velocity profiles: Re τ = 1000, no filter, different orders of accuracy. approach because filtering only removes information. For the present cases, our filtering method could improve the overall accuracy compared with corresponding no-model cases by about 30%. Table 4.11: Different orders of accuracy for Re τ = 1000, filtered case, Mesh: 96× 65× 128. order or accuracy spectral o10 o8 o6 o4 o2 Filtering interval 13 29 40 64 145 255 Error (comparison with DNS) U + (%) 0.8465734 2.70001 3.222855 3.895813 3.384768 3.621044 R + uu (%) 9.698326 12.15836 11.35095 16.54547 22.22139 39.35841 R + vv (%) 1.753032 1.21703 1.227448 1.222356 2.330911 3.605784 R + ww (%) 3.0109 4.134004 4.75554 4.737139 5.349884 7.439873 R + uv (%) 0.6222812 0.5585319 0.6474082 0.9562947 0.4037544 0.5686227 Average (%) 3.186222 4.153587 4.24084 5.471415 6.738142 10.91875 Numerical Dissipation (by FDM) RatioE n /E ν (%) 16.21293 4.566073 2.174174 -0.4368408 -0.3372216 3.564017 The deficiency for the lower order FDM schemes to capture the resolved scales can be explained by the energy spectra, similar as what we did in Section 4.3.2. The results for case Re1000l are shown in Fig. 4.33. As we can see in Fig. 4.22, the modified wavenumber for the 1st order derivatives decreases to zero, unlike the spectral method or the weak filters we use. Therefore, the energy spectra for 73 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.30: RMS velocities and Reynolds stresses: Re τ = 1000, no filter, different orders of accuracy. no-filter FDM simulations also decreased to a much lower value than the spectra from the spectral method. As a result, no energy pileup is observed. The only exception is order 2. The application of this low order scheme and coarse resolution is apparently insufficient to capture the physics of the flow, so the energy spectra is suspicious, which explains why we got the inaccurate results showed before. For other FDM schemes, higher order scheme removes less energy near the LES cutoff, as expected. However, all these schemes observe some oscillations at moderate wavenumbers, indicating that despite the energy accumulation did not occur near the cutoff, the excessive loss of information from moderate to high wavenumbers 74 Figure 4.31: Mean velocity profiles: Re τ = 1000, filtered with optimal period, different orders of accuracy. still contaminate the results. On the other hand, when filtering is applied, the unphysical energy accumulation at the smallest resolved scales for spectral and 2nd order FDM methods is effectively removed, so much better accuracy is achieved. For other FDM schemes, the agreement with DNS at moderate wavenumbers also improved, though the loss caused by the application of the modified wavenumbers cannot be recovered by our additional filtering. For this reason, our filtering were able to improve the results for different FDM schemes as well, though in order to get comparable results as the ones from the pseudo-spectral method, finer grid resolution is necessary to provide more spectral support. Note that the spectra for highorderFDMsisclosetothespectralmethoduptorelativelyhighwavenumbers, after which some small energy accumulations are observed, which also explains why their results are less satisfactory. Nevertheless, it is still good to see that the combination of our adaptive with non-spectral methods also leads to improved results. The similarity between the spectra of weakly filtered FDM and spectral method is also encouraging. This indicates that despite canonical compact FDM 75 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 4.32: RMS velocities and Reynolds stresses: Re τ = 1000, filtered with optimal period, different orders of accuracy. schemes (with the least number of nonzero entries) is purely dispersive, a slight addition of dissipation can effectively regularize the results, and we should be able to include it in the scheme itself. This would be analogous to the additional dissipation in WENO schemes (Liu et al., 1994b; Jiang & Shu, 1996), in which extra dissipation is added when the stencils are not smooth, or the approach in Lamballais et al. (2011), where hyper-viscosities are added through the 2nd order derivatives. An investigation of whether similar idea can be used to find optimal schemes to automatically reduce the energy accumulation and maintain the correct 76 shape of energy spectra, without the help of external filtering, can be a further work. 77 (a) E uu (b) E uu (filtered) (c) E vv (d) E vv (filtered) (e) E ww (f) E ww (filtered) Figure 4.33: One dimensional energy spectra the the center of the channel: Re τ = 1000, comparison of no-filter and filtered (with the optimal period) results. 78 Chapter 5 Explicit interscale model 5.1 General idea This section will focus on designing an explicit structural SGS model that removes excessive energy production by the energy cascade due to lack of small dissipative scales in LES. The approach is facilitated by the evaluation of energy transfer among multi-level scales. The idea follows closely the interscale energy transfer model proposed by Anderson & Domaradzki (2012), which adopts a form close to the similarity model (Bardina et al., 1980) to maintain the Galilean invari- ance. It is well recognized that the similarity-type models suffer from lack of dissi- pation and reported to have stability issue in many circumstances; this weakness is unfortunately maintained in the interscale model. Here we seek different math- ematical representations to compute the energy accumulation near the LES mesh cutoff more accurately, while still satisfying the necessary condition of Galilean invariance. Unlike most of the eddy-viscosity models, the model accounts for the backscatter as well, which is a crucial part of the actual SGS energy transfer. We show that theoretically the new model is more stable than the previous model in Anderson & Domaradzki (2012), which is subsequently confirmed byaprosteriori tests of turbulent channel flow. The effective removal of the energy generated by large scales near the LES cutoff enhances stability and leads to better prediction of low order statistics. 79 5.1.1 Scale separation by filtering Velocityfieldsofturbulentflowscanbeanalyzedasmulti-leveldata. Ingeneral, if we useN filtersG 1 ,...,G N with corresponding cutoff lengths ¯ Δ 1 ≤...≤ ¯ Δ N , the velocity fields of each level can be written as u 1 i = (G 1 ∗G 2 ∗...∗G N )∗u i (5.1) u m i = (G 1 ∗...∗G m )∗u i − (G 1 ∗...∗G m ∗G m+1 )∗u i , for 2≤m<N (5.2) In practice, normally only a limited number of filtering levels are considered, and a three-level decomposition for filtered velocity ¯ u i in LES equation (equation (2.6)) can be written as u 1 i = b ¯ u i , (5.3) u 2 i = ¯ u i − b ¯ u i , (5.4) u 3 i =u i − ¯ u i , (5.5) where ¯ stands for implicit filtering due to the mesh, whileb is explicit filtering. In this way, the decomposed velocity fieldu 1 i ,u 2 i andu 3 i stands for the resolved filtered field, the resolved scales near the grid cutoff, and the unresolved field, respectively. The full velocity field is u i =u 1 i +u 2 i +u 3 i , therefore the nonlinear term in the Navier-Stokes equation (equation (2.1)) can be written as ∂ ∂x j (u i u j ) = ∂ ∂x j ( 3 X p=1 u p i 3 X q=1 u q j ). (5.6) 80 Figure 5.1: Removal of unphysical energy accumulation by scale separation. Following Anderson & Domaradzki (2012), to remove energy mostly in the scales near the cutoff (region 2), we can further decompose the above nonlinear term in the same way ∂ ∂x j (u i u j ) = ∂ ∂x j 3 X m=1 ( 3 X p=1 u p i 3 X q=1 u q j ) m . (5.7) There are 27 possible combinations ofp,q andm in the above equation. How- ever, as the unresolved scales are missing, only 2 3 = 8 terms are known in actual LES runs. The effects of remaining unresolved terms needs to be modeled prop- erly. The net result of nonlinear interactions is to extract the energy from large scales to small scales and dissipate energy in under-resolved scales. Therefore, we aim at reducing some of the production of energy in small scales to model this overall dissipative effect. For simplicity, the nonlinear term without the divergence operator can be written as N pqm = (u p i u q j ) m . (5.8) 81 Specifically, the 8 nonlinear terms available in LES runs are written explicitly as N 111 = d b ¯ u i b ¯ u j , N 112 = b ¯ u i b ¯ u j − d b ¯ u i b ¯ u j , N 121 = d b ¯ u i ¯ u 0 j , N 122 = b ¯ u i ¯ u 0 j − d b ¯ u i ¯ u 0 j , N 211 = d ¯ u 0 i b ¯ u j , N 212 = ¯ u 0 i b ¯ u j − d ¯ u 0 i b ¯ u j , N 221 = d ¯ u 0 i ¯ u 0 j , N 222 = ¯ u 0 i ¯ u 0 j − d ¯ u 0 i ¯ u 0 j . (5.9) Following similar idea as the implicit LES method in the previous section, it is sufficient to consider scale interaction between neighboring bands, so here we only consider the production term in region m = 2 to reduce the unphysical energy pileup (see Fig. 5.1), i.e. N 112 and N 122 . Then the desired stress tensor takes the following form, τ ori ij =−N 112 −N 122 =−( b ¯ u i b ¯ u j − d b ¯ u i b ¯ u j )− ( b ¯ u i ¯ u 0 j − d b ¯ u i ¯ u 0 j ) = d b ¯ u i ¯ u j − b ¯ u i ¯ u j . (5.10) Unfortunately, the above tensor is not invariant to Galilean transformation, so we have to construct an alternative mathematical expression for τ ij while still maintaining the same idea. 5.1.2 Derivation of the SGS stress Intuitively, the removal of energy accumulation due to the interactions between different turbulent scales can have similar effects as an addition of a diffusive subgrid-scale stress. This can be theoretically justified by relating the resolved and the subgrid scales using the assumption of statistically stationary state 82 (Domaradzki & Liu, 1995), i.e. in each region (denoted byn) the energy is approx- imately constant, ∂E n ∂t ≈ 0, (5.11) where n = 1, 2, 3, following the notation from the previous section for different scales. We start from the Navier-Stokes equation for incompressible flow (equation (2.1)), which is written here again ∂u i ∂t +u j ∂u i ∂x j =− ∂ ∂x i p ρ ! +ν ∂ 2 u i ∂x j ∂x j i = 1, 2, 3 (5.12) Corresponding kinetic energy equation (E = 1 2 u i u i ) ∂E ∂t +u i ∂u i u j ∂x j =−u i ∂ ∂x i p ρ ! +νu i ∂ 2 u i ∂x j ∂x j , (5.13) which can be filtered to different regions to evaluate separately. Applying the statistical stationarity assumption (equation (5.11)), the kinetic energy equation in region 2 (m = 2) is approximately zero ¯ u i ∂ ∂x j 3 X p=1 3 X q=1 N pq2 + ¯ u i " 1 ρ ∂ (¯ p) ∂x i −ν ∂¯ u i ∂x j ∂x j # m=2 ≈ 0. (5.14) where ¯ u i is the grid filtered velocity that is supported by the LES mesh (region 1 and 2). For incompressible flow−ν¯ u i ∂¯ u i ∂x j ∂x j =− ∂ ∂x j ν ∂¯ e kin ∂x j +ν ∂¯ u i ∂x j ∂¯ u i ∂x j , where the last term is normally considered as the viscous dissipation. The above equation is equivalent to the spectral space derivation in Domaradzki & Liu (1995), except that the effect of mean profile (wavenumber k = 0) is included in the nonlinear term. 83 Note that for channel flow the mean pressure gradient will not contribute to the above equation, since region 2 does not contain the zeroth mode. For incom- pressible flow, in spectral space the pressure term exactly balances the parallel component of the nonlinear term, which is equivalent to solving the Poisson equa- tion in the physical space. Therefore, in actual implementations of the SGS stress, it is not necessary to account for the pressure term explicitly. The SGS energy transfer can be written as all the interactions involving region 3 affecting regions 1 and 2 T sgs = ¯ u i    3 X p=1 3 X q=1 2 X m=1 N pqm i 3∈p [ q    , (5.15) whereN i =− ∂u i u j ∂x j , which contains 10 terms. Since the LES cutoff is between region 2 and 3, subgrid scales normally have small impact on region 1. Based on the spectral analysis in Domaradzki et al. (1993), it is sufficient to consider the energy pileup between wavenumber 1 2 k c and k c , which in general holds for the very weak filters used in the present model. By neglecting the nonlinear interactions in region 1, T sgs ≈ ¯ u i    3 X p=1 3 X q=1 N pq2 i 3∈p [ q    . (5.16) 84 In this way, the unknown terms containing region 3 in equation (5.14) can be canceled out, leading to T sgs based only on computable terms involving LES velocity T sgs ≈−¯ u i (N 112 i +N 122 i +N 212 i +N 222 i )− ¯ u i " − 1 ρ ∂ (¯ p) ∂x i +ν ∂¯ u i ∂x j ∂x j # m=2 (5.17) = ¯ u i ∂ ∂x j (N 112 +N 122 +N 212 +N 222 ) + ¯ u i (I−G)∗ " 1 ρ ∂ (¯ p) ∂x i − 2ν ∂ ¯ S ij ∂x j # (5.18) =−¯ u i ∂τ sgs ij ∂x j , (5.19) where G is a filtering operator in physical space. Assuming that the application of taking derivatives and filtering commutes, the current model can be written explicitly as τ sgs ij ≈ d ¯ u i ¯ u j −¯ u i ¯ u j −(I−G)∗ " ¯ p ρ δ ij − 2ν ¯ S ij # = d ¯ u i ¯ u j −¯ u i ¯ u j − " ¯ p ρ δ ij − 2ν ¯ S ij # 0 . (5.20) The modification of nonlinear term can be considered as replacing the stress ¯ u i ¯ u j with d ¯ u i ¯ u j in the LES equation so that the production in region 2 will be smaller. As discussed earlier, the pressure term does not have to be included in the SGS stress. On the other hand, the above stress also reduces dissipation (replace ¯ S ij by b ¯ S ij ) thus makes simulations less stable. Recall that the statistical stationarity assumption does not necessarily hold instantaneously at each grid point where the SGS stress tensor is added. The current derivation is based on DNS so that the viscous dissipation in region 2 is removed to balance the effect of SGS scales. However, in LES due to insufficient spectral support from coarser mesh, extra dissipation at the smallest resolved scales is necessary. In addition, theapriori analysis in Domaradzki et al. (1993, 1994); Domaradzki & Liu (1995) 85 was based only relatively low Reynolds number. For high Reynolds numbers, the effect of the filtered viscous term is smaller. Therefore, in actual implementations this small anti-dissipative term is removed. On the other hand, the SGS energy transfer in region 1 is also neglected and should be in general dissipative. It is expected that due to more significant nonlocal interactions for higher Reynolds numbers, a further compensation in region 1 might be necessary. Therefore, the removal of the viscous term can also be regarded as adding extra dissipation to account for the interaction in region 1. Note that apart from the additional pressure and viscous terms, the nonlinear term also changed a little in comparison with the model in the previous section. The model in equation (5.10) removes only the production terms (N 112 andN 122 ), while the above equation contains all the terms in region 2. An additional filtering of the convective velocity in term N 212 = ¯ u 0 i b ¯ u j − d ¯ u 0 i b ¯ u j does not lead to reduced production, which appears due to the derivation based on the Navier-Stokes equa- tions that guarantees the stress to be symmetric. Thus the removal of this term will make the model asymmetric. Apparently,N 212 is asymptotically smaller than N 112 . According to Anderson & Domaradzki (2012), the removal of this term causes negligible difference in (τ ij −τ ji ) in their model. Besides, when the SGS dissipation is computed by multiplying a symmetric mean stain rate tensor S ij , the non-symmetry in off-diagonal terms is only going to impact the redistribution of energy rather than causing a net effect on energy transfer. The term is also not Galilean invariant, therefore we also choose to ignoreN 212 . The interaction among scales in region 2 only (N 222 ) is asymptotically much smaller than those involving 86 the well-resolved region 1, thus it is neglected. As a result, the model can still be reduced to the original interscale model (equation (5.10)) in the previous section τ ori ij =−N 112 −N 122 = d b ¯ u i ¯ u j − b ¯ u i ¯ u j . (5.21) 5.1.3 Modification to ensure Galilean invariance Galilean invariance is a property satisfied by the governing equations of fluid dynamics. It guarantees that the underlying physics in the governing equations remain the same under a transformation of inertial reference frame. Since the Navier-Stokes equations are frame-invariant, it is reasonable to have the same requirementforareducedsysteminLES.Thiscanbecheckedbyaddingaconstant velocity V i , the divergence of the SGS stress tensor should remain unchanged to ensure Galilean invariance. The original interscale energy transfer model (equation (5.10)) is non-Galilean invariant thus cannot be directly used as an SGS model. Here we proposed two general idea to ensure the Galilean invariance with minimal change to the original model. A relatively straightforward way to eliminate the effect ofV is to simply remove the mean velocity fieldhui. Then the resulting u−hui =u +V−h(u +V )i will not change as V =hVi. The same idea is also tested in Thiry & Winckelmans (2016). The modified SGS stress then become τ m0 ij = \ (¯ u i −h¯ u i i)(¯ u j −h¯ u j i) V − \ (¯ u i −h¯ u i i)(¯ u j −h¯ u j i), (5.22) inwhichhirepresentsthevolumeaverageoperator(foranonuniformgrid,integrate for the entire domain and then divide by its volume), sohu i i is the mean velocity 87 field (or bulk velocity) in each direction. Checking the Galilean invariance can be done by replacing u i with u i +V i and using the fact that V i =hV i i =const. τ m0 ij = (¯ u i +V i −h¯ u i +V i i) V (¯ u j +V j −h¯ u j +V j i) V − (¯ u i +V i −h¯ u i +V i i) V (¯ u j +V j −h¯ u j +V j i) (5.23) = (¯ u i −h¯ u i i +V i −hV i i) V (¯ u j −h¯ u j i +V j −hV j i) V − (¯ u i −h¯ u i i +V i −hV i i) V (¯ u j −h¯ u j i +V j −hV j i) (5.24) = \ (¯ u i −h¯ u i i)(¯ u j −h¯ u j i) V − \ (¯ u i −h¯ u i i)(¯ u j −h¯ u j i). (5.25) Thus the SGS stress in equation (5.22) can be recovered. As the mean velocity fieldhui can be regarded as a constant in each step, it is possible for the difference to be less significant after taking derivatives. However, to remove the mean veloci- ties for non-equilibrium or transient flows, it is necessary to compute the quantity every time step. Besides, it is speculated that the mean velocity gradient would be necessary to correctly represent the energy transfer. In a posteriori analyses for turbulent channel flow, the above model was able to yield similar results as the original model. But for relatively high Reynolds numbers and low grid resolu- tions, it seems that the current SGS stress is not sufficiently dissipative, so further modifications to control the backscatter is likely necessary. Because of these defi- ciencies, a different modification without a need to remove the mean velocities is sought. The violation of Galilean invariance for our nonlinear model is apparently due to filtering. Recall that a divergence operator will be applied to the SGS stress tensor. For incompressible flow assuming that filtering commutes with derivative operator, we basically need to bring the the solution u i in different terms to the same filtering level, so that they can be canceled. While for advection velocityu j it 88 does not make much sense to modify it, and it would not matter since a divergence operation will be taken anyway. The most straightforward way is to replace b ¯ u i in the last term with b b ¯ u i τ m1 ij =−N 112 −N 122 ≈ b ¯ u i ¯ u j V − b b ¯ u i ¯ u j . (5.26) Then τ m1 ij = (¯ u i +V i ) V (¯ u j +V j ) V − (¯ u i +V i ) V V (¯ u j +V j ) (5.27) = b ¯ u i ¯ u j +V i ¯ u j +V j b ¯ u i +V i V j V − ( b b ¯ u i ¯ u j +V i ¯ u j +V j b b ¯ u i +V i V j ) (5.28) = ( b ¯ u i ¯ u j V − b b ¯ u i ¯ u j ) +V i ( b ¯ u j − ¯ u j ). (5.29) The last term vanishes after taking the divergence operator ∂ ∂x j . The modification can be regarded as adding b ¯ u 0 i ¯ u j (where b ¯ u 0 i = b ¯ u i − b b ¯ u i ). We regard this as an attempt to model the SGS energy transfer effective in region 1 (equation (5.15)), which was neglected in the following derivations but should be reconsidered due to the more significant nonlocal interactions when Reynolds number becomes high. Basically, the additional term can be interpreted as the interaction between the subfiltered velocity in region 2 and resolved scales in region 1, acting on regions 1 and 2, which should act similarly as the interaction between region 3 (subgrid) and 1 based on energy cascade and similarity assumptions. In general ¯ u 0 i is proportional to higher order derivatives of ¯ u i . For common smooth filters such as 3-point top-hat filters or differential filters, the leading order is the 2nd order, and it leads to a dissipative behavior after integration by parts. The original interscale model τ ori ij , despite being non-Galilean invariant, is unstable in actual simulations when relatively high Reynolds number and low grid resolution is used. This is likely due to the absence of nonlocal energy transfer which we are 89 trying to add. Therefore, our modification can further stabilize simulations. To see this, let’s still use the 3-point top-hat filter in equation (4.4) as an example. Based on Taylor’s expansion we can write the filtered velocity as b b ¯ u i = b ¯ u i + c 2 " 3 X k=1 (Δx k ) 2 2! ∂ 2 b ¯ u i ∂x 2 k +O(Δx 4 k ) # . (5.30) Then the additional term b ¯ u 0 i ¯ u j =−¯ u j c 2 " 3 X k=1 (Δx k ) 2 2! ∂ 2 b ¯ u i ∂x 2 k +O(Δx 4 k ) # (5.31) = c 2 3 X k=1 (Δx k ) 2 2! " ∂¯ u j ∂x k ∂ b ¯ u i ∂x k − ∂ ∂x k ¯ u j ∂ b ¯ u i ∂x k ! +O(Δx 2 k ) # (5.32) = c 2 3 X k=1 (Δx k ) 2 2! " ∂¯ u j ∂x k ∂¯ u i ∂x k − ∂ ∂x k ¯ u j ∂¯ u i ∂x k ! +O(Δx 2 k ) # . (5.33) The first term can be considered as an extra dissipative term, whose strength is characterized by coefficientc depending on the choice of filter. The second term is a flux term, which should be relatively small for incompressible flow. It is note- worthy that the same leading order terms can also be obtained from other 2nd order discrete filters, for instance the differential filter (see equation (5.35) in the discussion of the α model, in which the coefficient in the leading order term is replaced byα 2 and high order terms are also different). Recall that the removal of the subfiltered energy in the original interscale model does not change the mean energy thus can hardly be sufficiently dissipative, so a stabilizing term would be necessary. This additional term may look like the gradient model, which was first developedbyLeonard(1975)andextendedinClarketal.(1979)(canalsobeclassi- fied as a general tensor diffusivity model or Clark model). A generalized expansion for different filters can be found in Carati et al. (1999). It is reported that the gradient model also encounters a blowup of the kinetic energy due to improperly 90 accounting for the energy transfer (Iliescu et al., 2003; Chumakov, 2006). There- fore, it is usually implemented with an eddy viscosity term, resulting in a mixed model (Winckelmans et al., 2001), which is what our model is trying to overcome. Since coefficientc is decided based on filtering in the interscale model and different series of higher order terms also exist, it should not be regarded as an coefficient of the gradient model. Like its original form, the modified model (equation (5.26)) still maintains a neat expression with only filtered velocity available from LES. In the rest of this chapter, the mathematical properties and validations through turbulent channel flow with the same pseudo-spectral code described in Section 3 will be presented. 5.2 Relation to other models Before proceeding to evaluate the model through actual simulations, in this section we first compare it with several other models based on filtered nonlinear terms. Since the present model is derived using properties of SGS energy transfer, it is fundamentally different from the models being compared. Nevertheless, the effectofsubgridscalesonresolvedscalesismostlyenergeticaction, theSGSmodels modify the spectral distribution of energy. We believe that for any SGS model to work on under-resolved LES grid, a removal or faster decay of energy at the smallest resolved scales is necessary to model the collective effect of subgrid scales and stabilize simulations. We show that common nonlinear models can also be reinterpreted as a removal of nonlinear energy transfer using the notation from Section 5.1.1. Our analysis is qualitatively consistent with the spectral analysis and simulation results in previous works. Since the models being compared did not take nonlinear interaction explicitly into account, they all suffers the same issue 91 of unreasonable energy removal. This motivates our model to remove unphysical energy production mostly near the LES cutoff, with smaller impact on larger scales through the application of smooth filters. 5.2.1 Relation to Leray’s regularization In Leray’s regularization (Leray, 1934) the solution u i is convected by a smoothed (filtered) velocity ˜ u j , and the governing equations (for incompressible flow and Newtonian fluids) become ∂u i ∂t + ˜ u j ∂u i ∂x j =− ∂p ∂x i +ν ∂ ∂x j ∂u i ∂x j ! ; ∂˜ u i ∂x i = 0, (5.34) where p =P/ρ, P is pressure and ρ is density. Note that even though u i is unfiltered, it cannot be regarded as full DNS solution from the Navier-Stokes equations (up to the smallest Kolmogorov length scale) due to the modification of the convective term. If we regard the equation as an alternative to LES, it is not appropriate to treat the filtered velocity ˜ u j as grid-filtered velocity either. In fact, historically the filtering in Leray regularization is defined through a Helmholtz operator (1−α 2 ∂ j ∂ j )˜ u i =u i . (5.35) The parameter α is used to characterize the so-called α-models (e.g. Leray-α, LANS-α, etc.). In general, solving extra Helmholtz equations in 3D is computa- tionally expensive, and it is possible to either use other types of filters as well (e.g. a top-hat filter used in Geurts & Holm (2003, 2006) and a Gaussian filter is used in Liu et al. (2008)), or apply it with only one or two Jacobi iterations (Verstappen, 92 2008). Therefore, there are two ways to interpret Leray’s equation as a surrogate LES model: (1). Unfiltered velocity u i in equation (5.34) is considered as velocity from DNS of the Leray equation. The translational energy spectrum for the filtered velocity decays at a rate of k −13/3 for kα 1 (Cheskidov et al., 2005) (based on Kraichnan’s theory of cascade Kraichnan (1967)), faster than the energy spectra of the Navier-Stokes equation, thus the grid resolution demands can be reduced. Therefore, a DNS of the Leray’s regularized equation can be considered as a LES for the Navier-Stokes equation. The grid-filtered velocity (denoted by bar) ¯ u i is obtained from LES runs. Then we can directly write the LES equation as ∂¯ u i ∂t + e ¯ u j ∂¯ u i ∂x j =− ∂ ¯ p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j ! ; ∂ e ¯ u i ∂x i = 0. (5.36) Then the SGS stress τ Leray1 ij = ¯ u i e ¯ u j − ¯ u i ¯ u j =−¯ u i ¯ u 0 j , where ¯ u 0 j = ¯ u j − e ¯ u j , while the true SGS stress for the Navier-Stokes equation is τ ij =u i u j − ¯ u i ¯ u j . Note that this requires the divergence of explicitly filtered velocity ˜ ¯ u i to be divergence free, which can cause issue for nonuniform grid and filters, see e.g. van Reeuwijk et al. (2006). (2). Filtered velocity ˜ u j in equation (5.34) is considered as the grid-filtered veloc- ity ¯ u j obtained from LES runs. Then following Geurts & Holm (2006), denoting the explicit filtering as ¯ u i =G∗(u i ), assuming that filtering commutes with deriva- tives, the equation (5.34) can be written as G −1 ∗ ∂¯ u i ∂t + ∂ ¯ p ∂x i −ν ∂ ∂x j ∂¯ u i ∂x j !! =−¯ u j ∂G −1 ∗ (¯ u i ) ∂x j (5.37) ⇒G −1 ∗ ∂¯ u i ∂t + ¯ u j ∂¯ u i ∂x j + ∂ ¯ p ∂x i −ν ∂ ∂x j ∂¯ u i ∂x j !! =G −1 ∗ ¯ u j ∂¯ u i ∂x j ! − ¯ u j ∂G −1 ∗ (¯ u i ) ∂x j (5.38) 93 Rewritten in LES form ∂¯ u i ∂t + ¯ u j ∂¯ u i ∂x j =− ∂ ¯ p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j ! − ∂ ∂x j (u i ¯ u j − ¯ u i ¯ u j ); ∂¯ u i ∂x i = 0. (5.39) Then the SGS stressτ Leray2 ij =u i ¯ u j − ¯ u i ¯ u j =τ ij −u i u 0 j , which is invariant with respect to to Galilean transformations as well as to rotating reference frame. Apparently the two interpretations are equivalent since we are still solving exactly the same equation. Note that even if the inverse filtering G −1 is exact and commutes with the derivatives, the two SGS stresses are not identical due to the nonlinear term. The difference can be considered as a subfilter stress (SFS), which can be computed exactly (see e.g. Berselli et al. (2006)). In practice, due to lack of spectral support for under-resolved LES mesh, it will not be possible to recover the unfiltered velocityu i (can be considered as DNS velocity) and subgrid velocity u 0 i , so one needs to either seek an approximate inversion ofG −1 (Geurts & Holm, 2003) (impossible if we regard the bar filter as a projector, such as sharp spectral cutoff filter for coarse mesh), or further modify the stress. The modifications depend on actual implementation and lead to additional uncertainties, so only the SGS stress from equation (5.34) is compared here. The idea of Leray’s model is to regularize on convective velocity u j while the interscale method focuses on the production effect ofu i , so they are fundamentally different. In fact, the equivalent SGS stress in equation (5.36) can be rewritten using the notation in equation (5.9) as τ Leray1 ij =−¯ u i ¯ u 0 j =−( e ¯ u i + ¯ u 0 i )¯ u 0 j =−N 121 −N 122 −N 221 −N 222 . (5.40) Similar as the interscale model in equation (5.10), Leray’s regulation also attempts to reduce some nonlinear interaction. The removal ofN 122 is the same as 94 in our model. However, instead of reducing production term N 112 for the smallest resolved wavenumbers, Leray’s regulation removed some nonlinear effect in region 1 (N 121 and N 221 ), which are in general dissipative interactions; removing these terms can be considered as adding backscatter. Besides, we think reducing the production effect in the smallest resolved scales is more imperative, as opposed to regularizing the nonlinear interaction in region 1. Based on the theoretical analysis from Leslie & Quarini (1979) and verifications through numerical data (e.g. Domaradzki et al. (1987); Domaradzki & Rogallo (1990); Domaradzki et al. (1993)), in LES the majority of energy is transfered locally from the intermediate wavenumber before the cutoff toward unresolved wavenumber just after it. The interaction between the small wavenumber and subgrid scales through nonlocal triad contribute less to the energy transfer. On the other hand, the application of filtering for our implicit LES method in Section 4 would lead to a similar effect, in the following sections we will show that the two approaches can cause similar decay of the energy spectra. As shown in Section 4.3.5, when stronger filtering is applied, which is analogous to removing more energy production including region 1, the results deteriorate. In other words, a goal of LES is to reduce the resolution requirement while still maintaining the large scale properties as much as possible, so we believe it is better to focus more on the energy transfer in region 2. The con- sideration in region 1, however, can be facilitated by using a smooth filter rather than a sharp spectral cutoff filter, which allows for some energy removal at large scales. For this reason, only low-pass smooth filters that become increasingly more active at higher wavenumbers are used for the present model. TheanalysisofenergytransferisconsistentwiththespectralanalysisinCheski- dovetal.(2005). InLeray’sregularization, therangeoffasterdecayofenergyspec- tra depends on the coefficientα, which can also be considered as a partial removal 95 of energy depending on the filtering strength. Thek −13/3 decay of the energy spec- trum of the filtered velocity forkα 1 (Cheskidov et al., 2005) indicates that the nominal kinetic energy∼k −1/3 in that region. This is not always observed though, if a coarse LES resolution is used, but it limits the choice ofα to be not much larger than the Kolmogorov dissipation scalek η . Nevertheless, the Leray’s regularization only offers a redistribution of energy rather than adding dissipation. The addi- tional SGS stress in the kinetic energy equation equals ¯ u i ∂ j (−¯ u i ¯ u 0 j ) =−∂ j ( 1 2 ¯ u 2 i ¯ u 0 j ), which shows the kinetic energy is still conserved after the regularization. As a result, insufficient removal of energy at the smallest resolved scales would gener- ate energy accumulation and gradually contaminate large scales due to the energy redistribution. Specifically, the removal of energy production from terms−N 122 and−N 222 is compensated by the inverse transfer from terms−N 121 and−N 221 . This is consistent with the simulations from Geurts & Holm (2006), in which a small over-prediction of large scales is accompanied by a slight under-prediction of small scales. However, improper energy distribution at smaller wavenumbers has larger impact on the lower order statistics interest in LES. In another investiga- tion from van Reeuwijk et al. (2009) with a 2nd order central difference scheme in space, Leray’s regularization is found unable to improve results over correspond- ing under-resolved DNS for channel flow, and the results deteriorate even more when α increases. Therefore, further modifications of the approach would be nec- essary to get accurate results, for instance the modified-Leray-α model (Ilyin et al., 2006) and the Leray deconvolution model (Dunca & Epshteyn, 2006; Layton & Lewandowski, 2008). 96 5.2.2 Relation to Navier-Stokes-alpha regularization The NS-α model is also known as the viscous Camassa-Holm equations (Holm et al., 1998a,b). The equations can be obtained through a nonlinear dispersive modification by filtering the velocity of the fluid loop in Kelvin’s circulation theo- rem. Chen et al. (1998, 1999,a,b) introduced the Navier-Stokes model by adding a Newtonian viscous term. The model can be written as ∂u i ∂t + ˜ u j ∂u i ∂x j +u k ∂˜ u k ∂x i =− ∂p ∂x i +ν ∂ ∂x j ∂u i ∂x j ! ; ∂˜ u i ∂x i = 0. (5.41) The NS-α model can be considered as the Leray’s regularization (Leray, 1934) with an additional term u k ∂˜ u k ∂x i . According to Guermond et al. (2003), the Leray regularization is not frame-indifferent, and the NS-α model can be considered as a frame-indifferent perturbation of the Leray regularization. The equation can also be derived through some other ways (Montgomery & Pouquet, 2002; Marsden & Shkoller, 2003). The filtered velocity ˜ u is still obtained from the Helmholtz operator in equation (5.35) characterized by the coefficent α. In this modified equation, it can be shown that a translational kinetic energy E α is bounded (see e.g. (Rebholz, 2007)). E α = 1 2 Z u i ˜ u i dx i = 1 2 Z [|˜ u i | 2 +α 2 |∂ j ∂ j ˜ u i | 2 ]dx i . (5.42) According to the spectral analysis of Foias et al. (2001), it decays at a rate of k −3 whenkα 1, faster than Komogolov’sk −5/3 law for the kinetic energy in the NS equations. Therefore, the required degrees of freedom to conduct DNS for the NS-α equation can be reduced, and the application of lower resolution grid can be in return regarded as LES. Similarly as Leray’s regularization in the previous section, two interpretations of the filtered velocity in the NS-α model can lead to 97 different SGS stresses. The results are available in Domaradzki & Holm (2001) and repeated here. (1). Unfiltered velocity u i in equation (5.41) is considered as velocity from DNS of the NS-α equation. Then the grid-filtered velocity (denoted by bar) ¯ u i obtained from LES runs, so we can directly rewrite the equation as ∂¯ u i ∂t + e ¯ u j ∂¯ u i ∂x j + ¯ u k ∂ ˜ ¯ u k ∂x i =− ∂ ¯ p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j ! ; ∂ e ¯ u i ∂x i = 0. (5.43) Then the SGS stress τ α1 ij = ¯ u i e ¯ u j − ¯ u i ¯ u j −α 2∂ ˜ ¯ u k ∂x i ∂ ˜ ¯ u k ∂x j =−¯ u i ¯ u 0 j −α 2∂ ˜ ¯ u k ∂x i ∂ ˜ ¯ u k ∂x j . Note that in the above equation, 1 2 ˜ ¯ u k ˜ ¯ u k has been included in the pressure term. (2). Filtered velocity ˜ u j in equation (5.41) is considered as the grid-filtered velocity ¯ u j obtained from LES runs. Then following the same procedure as in the previous section, we can get ∂¯ u i ∂t + ¯ u j ∂¯ u i ∂x j =− ∂ ¯ p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j ! − ∂ ∂x j u i u j − ¯ u i ¯ u j −u i ¯ u 0 j −α 2 ∂¯ u k ∂x i ∂¯ u k ∂x j ! . (5.44) Then the equivalent SGS stress τ α2 ij = u i u j − ¯ u i ¯ u j − u i u 0 j − α 2∂¯ u k ∂x i ∂¯ u k ∂x j = ¯ u i ¯ u j − ¯ u i ¯ u j +u 0 i u j −α 2∂¯ u k ∂x i ∂¯ u k ∂x j , where the u 0 i can be rewritten as LES velocities via the Helmholtz operator (equation (5.35)). Similarly as what we indicated in the previous section for Leray’s regularization, the difference is still simply a sub- filtered scale stress. To be consistent, let’s still analyzeτ α1 ij , which is Leray’s model minusα 2∂ ˜ ¯ u k ∂x i ∂ ˜ ¯ u k ∂x j , which looks like a gradient model (Leonard, 1975) (see a brief discussion in Section 5.1.3). The additional term made it possible for the equation to conserve a weaker translational energy in equation (5.42), which rolls off faster∼ k −3 when kα > 1 than the -5/3 law for the Navier-Stokes equations in the inertial range, making the 98 equation more computable. Nevertheless, this is equivalent to a k −1 decay of the nominal kinetic energy, indicating the approach is still not capable of removing the energy production near the LES cutoff if regarded as a SGS closure of the Navier- Stokes equation. In Graham et al. (2007, 2008), the spectrum is observed to be ∼k +1 for small scales (kα 1), which was attributed to the frozen-in-turbulence closure employed in deriving the model. Since the approach can be considered as the Leray’s regularization with subtracted gradient-type model, it is reasonable for the NS-α model to have comparable level of energy removal as Leray’s model and be less dissipative near the LES cutoff. According to Gibbon & Holm (2006) and Graham et al. (2007), the NS-α equation cannot asymptomatically reduce the number of degrees of freedom (DOF). Due to the slower decay (or even growing) of the energy spectrum at small scales, the choice ofα cannot be much larger than the dissipation scale either, similar to previous discussion of the Leray’s regularization. Note that in certain applications when the spectrum is not of a great concern, it is still possible to use relatively highα, making the computational savings a function of the Reynolds number. In actual simulations (e.g. Rebholz et al. (2017) for turbulent channel flow), the optimizedα (by comparing with the DNS) for the NS- α model and its modified version are generally also found to be small. Therefore, the model should be more precisely regarded as a quasi-DNS as opposed to a traditional LES. 5.2.3 Relation to Anderson and Domaradzki’s interscale energy transfer model The previous interscale energy transfer model proposed in Anderson & Domaradzki (2012) also originates from equation (5.10) that removes energy only 99 in the region near the LES cutoff. However, to ensure Galilean invariance, they adopted a modified form of N 112 −N 112∗ = d b ¯ u i b ¯ u j − b b ¯ u i b b ¯ u j (5.45) ≈ d b ¯ u i b ¯ u j − b ¯ u i b ¯ u j =−N 112 , (5.46) while N 122 remains the same, which leads to τ int ij =−N 112∗ −N 122 =−( b b ¯ u i b b ¯ u j − d b ¯ u i b ¯ u j )−( b ¯ u i ¯ u 0 j − d b ¯ u i ¯ u 0 j ) = d b ¯ u i ¯ u j − b b ¯ u i b b ¯ u j − b ¯ u i ¯ u 0 j . (5.47) However, b ¯ u = b b ¯ u holds only for projectors, the most common choice being the sharp spectral cutoff filters (more examples can be found in Buzzicotti et al. (2018)), whichisnotapplicabletocommonnon-spectralsolversthuswasnotadoptedinthe previous work. If a smooth filter is used (e.g. the three-point filter equation (5.48) used in Anderson & Domaradzki (2012)), double filtering in general is equivalent to stronger filter. To compare the modifications more clearly,τ ij is computed based on two simple 1D profiles (Gaussian, square wave) filtered by equation (5.48): b u n = 1 12 (u n−1 +u n+1 ) + 5 6 u n , (5.48) which is an optimized filtered used in Anderson & Domaradzki (2012). Our new modification in equation (5.26) is also compared. Original and filtered velocity profiles are shown in Fig. 5.2, and the 1D subgrid scale stress can be found in Fig. 5.3. Apparently our new modification agrees much better with the original model, especially for the smooth Gaussian profile. For interactions from the same profiles 100 (a) Gaussian (b) Square wave Figure 5.2: Filtered 1D velocity profile. (a) Gaussian (b) Square wave (c) interaction of Gaussian and square wave Figure 5.3: Sub-grid stress terms (without taking divergence) 101 (first 2 subfigures in Fig. 5.3), the modification also automatically reduces the negative region while leaving the positive region mostly unchanged. This can be explainedbytheadditionaldissipativeterminequation(5.33), whoseleadingorder is the square of a derivative of u for diagonal terms. In practice, this behavior is favorable as a means to automatically reduce some backscatter, while leaving the forward energy transfer mostly unchanged. On the other hand, the expression proposed in Anderson & Domaradzki (2012) leads to either completely positive or negative results. It is hard to expect the model to work in a similar way as the original model for more complex real turbulent flows, even after an application of the divergence operator. Therefore, we believe that our new model can more faithfully represent the original idea of energy removal. In fact, rather than an approximation of N 112 , it makes more sense to regard −N 112∗ as a pseudo-Leonard term (original Leonard stress ¯ u i ¯ u j − ¯ u i ¯ u j ). In this way, we may draw analogy between−N 112∗ and the very first similarity model from Bardina et al. (1983) τ ij = g e u i e u j − e e u i ˜ ˜ u j , (5.49) wheretheeffectofthemeshfilterisneglected(thefullvelocityfieldisapproximated with the LES resolved components, i.e. u≈ ¯ u in the previous notation), and ˜ here represents explicit filtering operations (Gaussian filter in Bardina et al. (1983)) as well. Therefore, the only difference between the original similarity model (Bardina et al., 1983) and the−N 112∗ term (Anderson & Domaradzki, 2012) is the filter shape. Liu et al. (1994a) suggests the application of a different filter in the second level with a larger filter width, which was tested with various filters τ ij =C( g ¯ u i ¯ u j − e ¯ u i e ¯ u j ), (5.50) 102 where C is a modeling coefficient, which can either be evaluated with spectral analysis by requiring the average of the trace of the model to be exact (Cook, 1997), or obtained with a dynamic procedure based on Germano identity (similar to the dynamic Smagorinsky model in Section 2.3.1). Apparently,−N 112∗ can simply be regarded as a similarity model using a test filtered velocity b ¯ u i , which can be used in the dynamic procedure. Therefore, the additional filtering likely does not make much change to the model, so the−N 112∗ part in Anderson & Domaradzki (2012)’s model can still be considered as a similarity model and share the same weaknesses. The second term,−N 122 , is asymptotically smaller thus relatively less important. Similarly to the analysis in the previous sections, the model can also be rein- terpreted as a combination of nonlinear terms in different regions. The double filtered velocity b b ¯ u i is denoted as region 0 for smooth filters, since it is essentially a filtered quantity from region 1. To better distinguish the two, here we use region 0 ∗ to indicate in general even larger filtered scales, with the new region 1 ∗ repre- senting the intermediate range for b ¯ u 0 i = b ¯ u i − b b ¯ u i . In other words, the combination of region 0 ∗ and 1 ∗ is the original region 1. Then the model can be rewritten with the previous notation N pqm with p,q∈{0 ∗ , 1 ∗ , 2} τ int ij =−N 112∗ −N 122 =−( b b ¯ u i b b ¯ u j − d b ¯ u i b ¯ u j )− ( b ¯ u i ¯ u 0 j − d b ¯ u i ¯ u 0 j ) (5.51) = d b b ¯ u i b ¯ u 0 j + d b ¯ u 0 i b ¯ u 0 j + d b ¯ u 0 i b b ¯ u j − ( b b ¯ u i b b ¯ u j − d b b ¯ u i b b ¯ u j )−N 122 (5.52) =N 0 ∗ 1 ∗ 1 +N 1 ∗ 1 ∗ 1 +N 1 ∗ 0 ∗ 1 −N 0 ∗ 0 ∗ 2 −N 0 ∗ 22 −N 1 ∗ 22 . (5.53) Apparently, the model attempts to removeN 0 ∗ 0 ∗ 2 andN 0 ∗ 22 . Since region 0 ∗ is just region 1 with extra filtering, these terms should have similar effect asN 112 and N 122 , as indicated in Anderson & Domaradzki (2012). However, duplicating the 103 nonlinear interaction at larger scales in region 1 is hardly desirable. In particular, the additional production term N 11 ∗ 1 (the combination of the first two terms in the above equation) is counterproductive, which is likely responsible for the energy blowup in actual simulations. On the other hand, no such terms were added in our refined model in equation (5.26), so it make sense to expect the model to be more stable under the same conditions. More precisely, we can compare it directly with our new model in Section 5.1.3 τ int ij −τ new ij =− b b ¯ u i b b ¯ u j − b ¯ u i ¯ u 0 j + b b ¯ u i ¯ u j , (5.54) = b b ¯ u i b ¯ u 0 j − b ¯ u 0 i ¯ u 0 j . (5.55) The first term is asymptotically larger than the second one, which in gen- eral causes an unnecessary enhancement of production. Therefore, the results of the new modification is much closer to the original non-Galilean invariant form (equation (5.10)), which we attempt to reproduce. For the same filtering strength (related to filtering coefficient c in equation (5.30)), our new model was able to reduce the unphysical energy accumulation near the cutoff wavenumber more effec- tively. As a consequence, the optimal filtering strength for the present model is smaller than that found in Anderson & Domaradzki (2012). More importantly, for the same Fourier-Chebyshev pseudo-spectral solver, our model was able to keep the simulations stable at low grid resolution using relatively weak filters, when the model in Anderson & Domaradzki (2012) sometimes fails. More detailed com- parisons between the two models will be presented in the a posteriori analysis in Section 5.3. 104 5.2.4 Relation to SFS stress and similarity model As discussed in the previous section, the Anderson & Domaradzki (2012)’s interscale model can be regarded as a variation of Bardina’s similarity model (Bar- dina et al., 1983) using a test-filtered velocity b ¯ u i . The similarity model (equation (5.49)) itself is mathematically the same as a subfilter-scale (SFS) stress. How- ever, they are applied to different equations. For SFS stress, the LES equation is further filtered by a filter denoted bye, where e u = G∗u can be any invertible filter. Assume that the filtering operation commutes with the derivatives ∂ e ¯ u i ∂t + e ¯ u j ∂ e ¯ u i ∂x j =− ∂ e ¯ p ∂x i +ν ∂ ∂x j ∂ e ¯ u i ∂x j ! − ∂ ∂x j (τ sfs ij + g τ sgs ij ), ∂ e ¯ u i ∂x i = 0, (5.56) where τ sfs ij = g ¯ u i ¯ u j − e ¯ u i e ¯ u j = τ sim ij . Since equation (5.56) is simply the LES equa- tion with an application of the tilde filter, there is no change to the equation. Apparently, if the SGS stress tensorτ SGS ij is neglected and the filtered e ¯ u i is used to compute other terms, it is still equivalent to solving for the filtered velocity with- out any modeling. Therefore, despite looking the same as the similarity model, the SFS stress itself should not be interpreted as a means of modeling. AnotherwaytoshowwhytheSFSmodeldoesnotofferanymodelingcapability is to rewrite the equation back to the LES form by taking the inverse of filtering G −1 ∗ () ∂¯ u i ∂t +G −1 ∗ e ¯ u j ∂ e ¯ u i ∂x j ! =− ∂ ¯ p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j ! − ∂ ∂x j G −1 ∗τ sfs ij +τ sgs ij . (5.57) ⇒ ∂¯ u i ∂t + ¯ u j ∂¯ u i ∂x j =− ∂ ¯ p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j ! − ∂ ∂x j h G −1 ∗ τ sfs ij + e ¯ u i e ¯ u j − ¯ u i ¯ u j +τ sgs ij i . (5.58) 105 In some works, e.g. Bull & Jameson (2016), only τ sfs ij (or G −1 ∗τ sfs ij ) is used, by assuming the effect of SGS stress on filtered velocity e ¯ u i is small. However, equation (5.58) is still just the LES equation since we applied G∗ and G −1 ∗ to rewrite the equation. The last term in the equation remains to be τ sgs ij since the other 3 terms simply canceled each other, but it does not make sense to assume τ sgs ij ≈ 0. If G −1 ∗τ sfs ij is the only term retained, then it is equivalent to saying τ sgs ij =G −1 ∗τ sfs ij = ¯ u i ¯ u j −G −1 ∗ e ¯ u i e ¯ u j . Therefore, the SFS stress in LES equation can be considered as a deconvoluted similarity model, due to the application to different filtering level. If the commutation error between filtering and derivative operators are negligible and no approximate deconvolution is used (which normally adds viscous or hyper-viscous terms that have a stabilizing effect), the results should be the same as under-resolved DNS (UDNS). On the other hand, if the similarity model is applied to the LES level (¯ u i ) with the assumption of true SGS stress τ sgs ij = ¯ u i ¯ u j − ¯ u i ¯ u j ≈ g ¯ u i ¯ u j − e ¯ u i e ¯ u j . In this way, the similarity model is able to offer some modeling capability of the SGS stress. However, despite achieving very high correlations in a priori analysis, in practice the similarity model is widely reported to cause stability issues. For this reason, normally the modifications of Bardina’s similarity model are used as a SGS model. TheconversionofBardina’ssimilaritymodeltoournonlinearnotationisalmost the same as equation (5.53), which is also outlined in Domaradzki (2010) and repeated as follows τ sim ij = g ¯ u i ¯ u j − e ¯ u i e ¯ u j (5.59) = g e ¯ u i ¯ u 0 j + g ¯ u 0 i e ¯ u j + g ¯ u 0 i ¯ u 0 j − ( e ¯ u i e ¯ u j − g e ¯ u i e ¯ u j ) (5.60) =N 121 +N 211 +N 221 −N 112 . (5.61) 106 The similarity model also removes N 112 , which is the most direct reduction of the unphysical energy buildup for the smallest resolved wavenumbers. However, the model adds nonlinear interaction in region 1 in a similar way as the interscale model in Anderson & Domaradzki (2012). This explains why similarity models is widely reported to be insufficiently dissipative and needs to be combined with an eddy viscosity model. 5.3 A posteriori tests for turbulent channel flow In this section, the results using different SGS models are compared. The simulations are based on the pseudo-spectral solver described in Section 3. The test cases are similar to those for our implicit LES method in Section 4. Readers can refer to Table 4.2 for domain sizes and grid resolutions. We compare the non-Galilean invariant original model (equation (5.10), denoted as ’NG’) with the interscale energy transfer model (equation (5.47), denoted as ’INT’) and the new Galilean invariant modified model (equation (5.26), denoted as ’Model’) inaposteriori testing. The DNS benchmark are described in Section 3. Note that the original model violates Galilean invariance thus cannot be used generally as a SGS model, here we just use it as a benchmark to show how similar our modification is compared with the original idea. After some testings, a very weak filter is found to give the best results for the original model: b u n = 1 32 (u n−1 +u n+1 ) + 15 16 u n , (5.62) which corresponds to c = 1 8 in equation (4.4). On the other hand, the optimal filter found in Anderson & Domaradzki (2012) uses c = 1 3 . Since the new model removes energy more effectively, we are able to get good results using a weaker 107 filter. If the model is written in series formulation, the coefficients of nonlinear terms will be proportional to different powers of grid size, unlike the α models in which everything is absorbed intoα. Therefore, it is reasonable to expect the same filter to work consistently under different conditions, so we apply the same filter throughout this section. In Anderson & Domaradzki (2012), a local quadratic interpolation (Loh & Domaradzki, 1999) was applied in the vertical direction of the channel flow, so that thesamefilteringcanbeused. However,thiswillalsogeneratecommutationerrors, which means the Galilean invariance cannot be perfectly satisfied. The issue would be more severe when relatively high Reynolds number flows are simulated with coarse grid in the vertical. In the current work, the mesh in the vertical direction is based on the Chebyshev polynomials, with more points clustered near the wall thanthosefromtheLegendrepolynomialsusedinAnderson&Domaradzki(2012). Therefore, we choose not to apply vertical filtering in the model. In practice, due to the reduction of production in the horizontal planes, the smallest resolved scales are already less energetic, thus the application of the vertical filtering is found to have minimal impact on the results in a posteriori analysis. Note that when the interscale model is compared, we still used the same filtering in 3D as described in Anderson & Domaradzki (2012) to fully recover the their results. 5.3.1 Different Reynolds number Inthissection,themeanvelocityprofilesandtheReynoldsstressesforReynolds numbers Re τ = 180, 590, 1000, 2000 are compared. The benchmark DNS results are summarized in Table. 4.3. We show that our new model is capable of produc- ing very similar results as the original idea of energy removal (equation (5.10)). The model is also in general more accurate than the similarity-type model from 108 Anderson & Domaradzki (2012) when implemented in the same code, especially the predictions of the second order statistics. Figure 5.4: Mean velocity profiles: Re τ = 180, Mesh: 64× 65× 64. Results for Re τ = 180 with 64× 65× 64 mesh are shown in figures 5.4 and 5.5. The modified model in general produces comparable results as the original one. It is noteworthy that the optimal filter in the interscale energy transfer model (Anderson & Domaradzki (2012), equation (5.47)) is different. When the same filter is used (c = 1 8 ), the results from our new modification are much closer to the ones for the original model, indicating that the model can serve as a Galilean invariant surrogate for the original one. For the model of Anderson & Domaradzki (2012) the results are almost the same as for UDNS when this very weak filter is used. When it comes to lower resolution (case Re180) in Figs 5.6 and 5.7, the results from the original model are slightly less satisfactory due to the coarse resolution. Recall that the original form is derived entirely based on the removal of energy production through the interactions from the resolved scales in region 2. When much coarser grid is used than DNS, indicating less spectral support, this may 109 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.5: RMS velocities and Reynolds stresses: Re τ = 180, Mesh: 64× 65× 64. not be sufficient to fully represent the effect of subgrid scales. On the other hand, due to the additional modeling of SGS dissipation in our refined model, good agreement with DNS is still preserved. The interscale model with the optimal filter also accurately predicted the mean velocity profile, but the Reynolds stresses are still less satisfactory than the new model. For Re τ = 590, the previous observations still hold. When moderate grid resolution is used (case Re590, Figs 5.8 and 5.9), our new model produces very similar results as the original model. The agreement with DNS is good and much better than the no model simulation (UDNS). The interscale medel with optimal 110 Figure 5.6: Mean velocity profiles: Re τ = 180, Mesh: 48× 49× 48. filter also leads to accurate mean velocity profile, but the Reynolds stresses are consistently under-predicted. If the same weak filter is used, the interscale model shows no difference compared with UDNS. When it comes to a lower resolution (case 590l, Figs 5.10 and 5.11), our new model is still less sensitive to the reduced resolution than the original model. Nevertheless, the Reynolds stresses of our new modelisalsounder-predictedcomparedwithfinerresolutioncase. Thisisexpected since we are always comparing with the DNS results, while a more fair comparison should be made with filtered/sampled DNS data on LES grid. As a result, when coarser grid is used, a slight decrease of Reynolds stresses is observed, but it is still much better than the predictions from the interscale model. Similar to what we saw before, the results for the interscale model with the same filter are still almost the same as corresponding no model simulations, indicating that the model fails to effectively remove the nonlinear energy transfer. Since this phenomenon consistently occurs for other cases as well, it will not be shown for the rest of the cases for clarity. 111 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.7: RMS velocities and Reynolds stresses: Re τ = 180, Mesh: 48× 49× 48. For Re1000 case (Figs 5.12 and 5.13), all three models still generate satisfactory results despite the use of a much lower resolution in comparison with DNS. Overall the new model still performs the best and is slightly less dissipative than the inter- scale model, due to the application of a much weaker filter. It is interesting to see that the mean velocity profiles from the two approaches almost collapse on each other. This is most likely due to the fact that the leading order term for almost all the nonlinear models has a form close to a gradient model, if rewritten in series notation (see Section 5.1.3). The leading order term can become comparable by tuning the parameters in filtering. Since the first order statistics is relatively less 112 Figure 5.8: Mean velocity profiles: Re τ = 590, Mesh: 96× 97× 96. sensitive to the higher order terms in the SGS models, it makes sense for different nonlinear models with optimized filter to get rather similar mean velocity profiles. However, in order to effectively reduce the energy production for high wavenum- bers, the combination of higher order terms is really what makes the difference. Based on the discussions in Section 5.2.3, our model has much less impact on the well-resolved large scales in region 1. On the other hand, in order to compensate for the inability to reduce energy production, the interscale model needs to be used with a much stronger filter, which is more active at lower wavenumbers. As a consequence, the higher order statistics such as the Reynolds stresses are con- sistently under-predicated, while the results from the present model agree much better with the DNS benchmark. On the other hand, despite these second order statistics after other modifications being almost the same as the original model, the mean velocity profile from the present model is visibly more accurate. This observation in general holds for all other Reynolds numbers as well. When the grid resolution is fine or moderate, the modification hardly makes any difference compared with the original idea. When the resolution is reduced (all cases without 113 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.9: RMS velocities and Reynolds stresses: Re τ = 590, Mesh: 96× 97× 96. ‘*h’ in Table 4.2), our modification improves the results thanks to the additional compensation of subfiltered energy transfer (equation (5.33)). When Reynolds number increases to 2000, the interscale model failed to stabi- lize the flow due to the relatively high Reynolds number and low grid resolution. Therefore, the present model is compared only with the original model, the results are shown in figures 5.14 and 5.15. Both methods gives accurate perditions, but the present model is still slightly model accurate, which is observed consistently for all the cases we tested. 114 Figure 5.10: Mean velocity profiles: Re τ = 590, Mesh: 64× 65× 64. AmongdifferentReynoldsnumbers, thepresentmodelisabletooffersignificant improvement to corresponding no model simulations. The model is also very easy toimplement becausethethree-point filterinequation (5.62)offersnegligible extra computational cost. We believe the results in this section show that it is possible to provide sufficient SGS dissipation without extra support from linear dissipative models such as the eddy viscosity models. Unlike the α models, the automatic combination of mesh sizes through the filtering make the model less affected by the change of grid resolution. More detailed discussion about the grid convergence of the model is presented in the next section. 115 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.11: RMS velocities and Reynolds stresses: Re τ = 590, Mesh: 64×65×64. Figure 5.12: Mean velocity profiles: Re τ = 1000, Mesh: 128× 97× 128. 116 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.13: RMS velocities and Reynolds stresses: Re τ = 1000, Mesh: 128× 97× 128. 117 Figure 5.14: Mean velocity profiles: Re τ = 2000, Mesh: 256× 131× 128. 118 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.15: RMS velocities and Reynolds stresses: Re τ = 1000, Mesh: 256× 131× 128. 119 5.3.2 Different grid resolution Inthissection, weshowthatforrelativelyhighReynoldsnumbers, thestatistics gradually show better converge to DNS solution when the mesh is refined, which further justifies the consistent use of the weak 3-point filter in equation (5.62). Figure 5.16: Mean velocity profiles: Re τ = 1000, different resolution. ForRe τ = 1000 cases, the details about the grid can still be found in Table 4.2. Since the accuracy in the previous section for case Re1000 is already satisfactory, here we used one coarser and one finer grid to show convergence. For lower resolu- tion case Re1000l, the mean velocity profile in Fig. 5.16 already converged. Never- theless, theReynoldsstatisticsinFig.5.17areslightlylesssatisfactory. Theresults for case Re1000 we showed before are already almost the same as for case Re1000h, including the Reynolds stress. Therefore, no matter if resolution is increased or decreased, it is still not necessary to change to the filter; the grid sizes already implicitly made sufficient justifications for the present cases. When the friction Reynolds number doubles to 2000, the results for the second order statistics showed in Fig. 5.15 have some room of improvement. Therefore we choose to only refine the mesh. Corresponding grids are summarized in Table 5.1, 120 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.17: RMS velocities and Reynolds stresses: Re τ = 1000, different resolu- tion. in which Re1000l stands for the case we showed in the previous section. Recall that we did not apply filtering in the vertical direction. In test case Re2000 only the vertical resolution is increased, which improved on the under-predictions of the Reynolds stresses, as expected. In case Re2000h we doubled the spanwise resolution and slightly refined the vertical grid as well to be safe, and even better agreement with DNS is observed. These results confirm that our model is insensitive to the grid resolution, so that further improvements like a dynamic procedure are not necessary. For all the 121 Table 5.1: LES parameters for Re τ = 2000. Case Domain Resolution N/N DNS (%) Re2000h 2.5π× 2× 0.5π 256× 151× 256 0.102 Re2000 2.5π× 2× 0.5π 256× 181× 128 0.061 Re2000l 2.5π× 2× 0.5π 256× 131× 128 0.044 Figure 5.18: Mean velocity profiles: Re τ = 2000, different resolution. cases we tested, including the results not shown, we observe in general a gradually convergence to the DNS benchmarks when the mesh is refined. 122 (a) u + rms (b) v + rms (c) w + rms (d)−R + uv Figure 5.19: RMS velocities and Reynolds stresses: Re τ = 2000, different resolu- tion. 123 5.3.3 Energy spectra The ability to remove energy at the smallest resolved wavenumbers can be evaluated by checking the energy spectra directly. Similar to what we did for the implicit model, the spectrum is still sampled at the center of the channel and is normalized by u 2 τ . (a) E uu (b) E vv (c) E ww Figure 5.20: Mean velocity profile and streamwise energy spectra: Re τ = 590, Mesh: 64× 65× 64. For Re τ = 590 case in Fig. 5.20 (case Re590l in Table 4.2), the UDNS hardly follows the DNS profile, same for the interscale model with weak filter (c = 1 8 ). All other 3 models are capable to closely match the DNS spectra except relatively close 124 to the LES mesh cutoff (or the Nyquist wavenumber k c =π/Δ). This is why the low order statistics are significantly improved. Despite using a different filter, the spectra of the interscale model (c = 1 3 ) are close to those from the present model, which also explains the very similar mean velocity profiles we showed earlier for almost all test cases. For higher Reynolds number Re τ = 2000 case in Fig. 5.21 (case Re2000l in Table 5.1), we compare the energy spectra for different grid resolutions in Section 5.3.2. Since all of them have the same resolution in the streamwise direction, the 1D streamwise spectra all have exactly the same spectral support. Therefore, the spectra for no model simulations all look quite similar. However, when the grid resolution in the vertical and spanwise directions are inadequate (case Re2000l), the present model may not be able to completely remove the energy accumulation near the LES cutoff, indicating grid convergence has not yet been achieved and finer mesh is necessary. When the mesh is refined (cases Re2000 and Re2000h), the energy pileup for low and moderate wavenumbers has been removed, so the agreement of low order statistics with DNS also improved. For this relatively high Reynolds number, the energy spectra from LES is not very smooth, which is also observed for other nonlinear dispersive models. Nevertheless, for channel flow the present model is capable of reducing sufficiently energy at moderate and high wavenumbers to stabilize the simulations and improve the results. 125 (a) E uu (b) E uu (UDNS) (c) E vv (d) E vv (UDNS) (e) E ww (f) E ww (UDNS) Figure 5.21: One dimensional energy spectra the the center of the channel: Re τ = 2000, compare with corresponding no model simulations. 126 Chapter 6 Conclusions In LES due to the lack of small subgrid scales and insufficient spectral support from the mesh, it is necessary to provide adequate SGS dissipation in simulations. Failure to properly account for this would lead to mis-predicted flow quantities or even cause a blowup of kinetic energy in finite time. Purely dissipative models such as the eddy viscosity models are widely used in LES and give accurate results despite poor correlation with the real SGS stress. Nonlinear dispersive models such as the similarity model are normally reported to be not sufficiently dissipative and need to be combined with an eddy viscosity model, resulting in a mixed model. The present research emphasized the importance of reducing energy pileup in LES through simulations with a non-dissipative pseudo-spectral solver. We show that it is possible to overcome the problem in both explicit and implicit LES without any additional help from the eddy viscosity concept. Anewimplicitlargeeddysimulations(ILES)modelbasedonanadaptivefilter- ing criterion has been proposed. The approach was tested for a turbulent channel flow using a Fourier-Chebyshev pseudo-spectral code. The filtering criterion is based on the ratio of the kinetic energy removed by two different low-pass filters and is designed to remove the energy pileup near the cutoff wavenumber before it can contaminate larger scales. The design and testing of the method has been made possible by the recent availability of the procedure that allows to quantify the numerical dissipation in an arbitrary numerical code. The numerical dissipa- tion imparted by the filtering is found to linearly depend on the filtering period. 127 For a low-pass filter an optimal filtering period can be found, which remains con- stant in the statically steady state, making the ILES with a fixed filtering period almost as efficient as a no-model, under-resolved DNS. Low order statistics for Re τ between 180 and 2000 obtained using the method show a good agreement with DNS benchmarks, offering significant improvements over the corresponding UDNS and outperforming the dynamic Smagorinsky model. In general, different amounts of numerical dissipation are required for different cases, depending on the level of under-resolution of the mesh and an overall order of accuracy of a Navier-Stokes solver. The proposed filtering criterion can properly adjust to different Reynolds numbers and grid resolutions, without manually tun- ing any parameters, to provide an appropriate amount of numerical dissipation as a surrogate for explicit sub-grid (SGS) models. The filtering criterion in principle does not depend on the shape of the filter and a set of filters constructed by an approximate deconvolution method was tested in this work. When a filter removes more energy in one operation, fewer filtering operations are necessary over a given simulation time to have a comparable amount of time-averaged numerical dissipa- tion compared with an application of a weaker filter. However, to avoid negative effects on larger scales, a filter that weakly attenuates the small scales with mini- mal impact on the larger ones is safer, thus more advisable. The same procedure has also been applied to lower order schemes using the modified wavenumber con- cept. An increase of accuracy is also achieved in such cases but it is relatively less substantial due to the presence of the numerical dissipation generated by the solver itself. In summary, the proposed ILES method appears general, providing a simple and a computationally efficient tool to significantly improve the simulation accuracy of turbulence in a framework of arbitrary Navier-Stokes solvers. 128 The filtering in the ILES method modifies the kinetic energy spectrum, so it can still be regarded as a purely dissipative approach, similar to the eddy viscos- ity models. As a natural extension, an explicit SGS model is proposed, which reduces the energy production through the nonlinear term in the Navier-Stokes equations. The model is derived from the energy transfer in multi-level scales thus can be regarded as a direct removal of energy production in targeted regions. The resulting nonlinear model in general has a dissipative effect, but still allows for backscatter in a manner consistent with the action of real SGS stress. Despite the dispersive nature of nonlinear models, we think that for any SGS model to provide sufficient dissipation and stabilize flows, a partial removal and/or faster decay of energy at the smallest resolved scales is necessary. Nevertheless, a model that removes energy only at small scales cannot satisfy Galilean invariance, so it is combined with a gradient-type modification, which is found to have negligible impact on the original idea. No stabilization techniques such as adding dissipa- tion or controlling backscatter is used. Our new model is first compared with the regularization methods and similarity-type models. We show that all these models fail to properly redistribute energy transfer in desired regions, which motivates our model to remove energy production mostly at small scales, and decreasing gradually at larger scales through the application of a weak smooth filter. A posteriori analysis under the same condition as the implicit method is also conducted for the explicit model. Since the filtering operation already took grid size into account, the model was able to consistently produce accurate results. Therefore, a dynamic procedure may not be necessary, which makes the model easy to implement and computationally efficient. Due to the direct modeling of theSGSstressandtheallowanceforbackscatter, theexplicitmodelingeneralgives more accurate predictions of low order statistics than our implicit approach under 129 the same condition. The current research is focused on wall-bounded flows, which provides higher dissipation near the wall due to the no-slip boundary condition. Currentlynofilteringisimplementedintheverticaldirectionfortheexplicitmodel, so a further work could be extending the model to properly handle nonuniform mesh. Moreover, to further investigate the ability of energy removal by the present model, we plan to test it for homogeneous isotropic turbulence with high Reynolds numbers and low resolutions in the near future. 130 Reference List AluieH,EyinkGL(2009) Localnessofenergycascadeinhydrodynamicturbulence. II. Sharp spectral filter. Physics of Fluids 21:115108. Anderer M (2015) Explicit subgrid scale modeling for the simulation of turbu- lent flows with high order Discontinuous Galerkin methods. Bachelor’s thesis, Universität Stuttgart. Anderson BW, Domaradzki JA (2012) A subgrid-scale model for large-eddy sim- ulation based on the physics of interscale energy transfer in turbulence. Physics of Fluids 24:065104. Bardina J, Ferziger J, Reynolds W (1980) Improved subgrid-scale models for large-eddy simulation In 13th Fluid and PlasmaDynamics Conference. American Institute of Aeronautics and Astronautics. Bardina J, Ferziger J, Reynolds W (1983) Improved turbulence models based on large eddy simulation of homogeneous, incompressible turbulent flows. Stanford Univ. Report . Berland J, Lafon P, Daude F, Crouzet F, Bogey C, Bailly C (2011) Filter shape dependence and effective scale separation in large-eddy simulations based on relaxation filtering. Computers & Fluids 47:65–74. Berselli L, Iliescu T, Layton W (2006) Mathematics of Large Eddy Simulation of Turbulent Flows Springer-Verlag. Bogey C, Bailly C (2006) Computation of a high Reynolds number jet and its radiated noise using large eddy simulation based on explicit filtering. Computers & Fluids 35:1344–1358. Boris J, Grinstein F, Oran E, Kolbe R (1992) New insights into large eddy simu- lation. Fluid dynamics research 10:199–228. Borue V, Orszag SA (1998) Local energy flux and subgrid-scale statistics in three- dimensional turbulence. Journal of Fluid Mechanics 366:1–31. 131 Bose ST, Moin P, You D (2010) Grid-independent large-eddy simulation using explicit filtering. Physics of Fluids 22:105103. Bull JR, Jameson A (2016) Explicit filtering and exact reconstruction of the sub-filter stresses in large eddy simulation. Journal of Computational Physics 306:117–136. Buzzicotti M, Linkmann M, Aluie H, Biferale L, Brasseur J, Meneveau C (2018) Effect of filter type on the statistics of energy transfer between resolved and subfilter scales from a-priori analysis of direct numerical simulations of isotropic turbulence. Journal of Turbulence 19:167–197. Cadieux F, Domaradzki JA (2016) Periodic filtering as a subgrid-scale model for LES of laminar separation bubble flows. Journal of Turbulence 17:954–965. Cadieux F, Domaradzki JA (2015) Performance of subgrid-scale models in coarse large eddy simulations of a laminar separation bubble. Physics of Flu- ids 27:045112. Cadieux F, Sun G, Domaradzki JA (2017) Effects of numerical dissipation on the interpretation of simulation results in computational fluid dynamics. Computers & Fluids 154:256–272. Canuto C, Hussaini MY, Quarteroni A, Zang TA (2007) Spectral Methods : Fun- damentals in Single Domains Springer, Dordrecht. CaratiD,WinckelmansGS,JeanmartH(1999) Exactexpansionsforfiltered-scales modelling with a wide class of LES filters In Direct and Large-Eddy Simulation III, pp. 213–224. Springer Netherlands. Castiglioni G, Domaradzki J (2015) A numerical dissipation rate and viscos- ity in flow simulations with realistic geometry using low-order compressible Navier–Stokes solvers. Computers & Fluids 119:37–46. Chen S, Foias C, Holm DD, Olson E, Titi ES, Wynne S (1999a) A connection betweentheCamassa–Holmequationsandturbulentflowsinchannelsandpipes. Physics of Fluids 11:2343–2353. Chen S, Foias C, Holm D, Olson E, Titi E, Wynne S (1999b) The Camassa–Holm equations and turbulence. Physica D: Nonlinear Phenomena 133:49–65. Chen S, Foias C, Holm DD, Olson E, Titi ES, Wynne S (1998) Camassa-Holm equationsasaclosuremodelforturbulentchannelandpipeflow. Physical Review Letters 81:5338–5341. 132 Chen S, Holm DD, Margolin LG, Zhang R (1999) Direct numerical simulations of the Navier–Stokes alpha model. Physica D: Nonlinear Phenomena 133:66–83. Cheskidov A, Holm DD, Olson E, Titi ES (2005) On a Leray–α model of turbu- lence. Proceedings of the Royal Society A: Mathematical, Physical and Engineer- ing Sciences 461:629–649. Chumakov SG (2006) Statistics of subgrid-scale stress states in homogeneous isotropic turbulence. Journal of Fluid Mechanics 562:405. Cimarelli A, Angelis ED (2014) The physics of energy transfer toward improved subgrid-scale models. Physics of Fluids 26:055103. Clark RA, Ferziger JH, Reynolds WC (1979) Evaluation of subgrid-scale models using an accurately simulated turbulent flow. Journal of Fluid Mechanics 91:1. Cook AW (1997) Determination of the constant coefficient in scale similarity models of turbulence. Physics of Fluids 9:1485–1487. Dairay T, Lamballais E, Laizet S, Vassilicos JC (2017) Numerical dissipation vs. subgrid-scale modelling for large eddy simulation. Journal of Computational Physics 337:252–274. Domaradzki JA (2010) Large eddy simulations without explicit eddy viscosity models. International Journal of Computational Fluid Dynamics 24:435–447. Domaradzki JA, Adams NA (2002a) Direct modelling of subgrid scales of turbu- lence in large eddy simulations. Journal of Turbulence 3:N24. Domaradzki JA, Adams NA (2002b) Direct modelling of subgrid scales of turbu- lence in large eddy simulations. Journal of Turbulence 3:N24. Domaradzki JA, Yee PP (2000) The subgrid-scale estimation model for high Reynolds number turbulence. Physics of Fluids 12:193–196. Domaradzki JA, Carati D (2007a) An analysis of the energy transfer and the locality of nonlinear interactions in turbulence. Physics of Fluids 19:085112. Domaradzki JA, Carati D (2007b) A comparison of spectral sharp and smooth filters in the analysis of nonlinear interactions and energy transfer in turbulence. Physics of Fluids 19:085111. Domaradzki JA, Holm DD (2001) Navier-Stokes-alpha model: LES equations with nonlinear dispersion. arXiv preprint nlin/0103036 . 133 Domaradzki JA, Liu W (1995) Approximation of subgrid-scale energy trans- fer based on the dynamics of resolved scales of turbulence. Physics of Flu- ids 7:2025–2035. Domaradzki JA, Liu W, Brachet ME (1993) An analysis of subgrid-scale interac- tions in numerically simulated isotropic turbulence. Physics of Fluids A: Fluid Dynamics 5:1747–1759. Domaradzki JA, Liu W, Härtel C, Kleiser L (1994) Energy transfer in numerically simulated wall-bounded turbulent flows. Physics of Fluids 6:1583–1599. DomaradzkiJA,LohKC(1999) Thesubgrid-scaleestimationmodelinthephysical space representation. Physics of Fluids 11:2330–2342. DomaradzkiJA,LohKC,YeePP(2002) Largeeddysimulationsusingthesubgrid- scale estimation model and truncated Navier-Stokes dynamics. Theoretical and Computational Fluid Dynamics 15:421–450. Domaradzki JA, Metcalfe RW, Rogallo RS, Riley JJ (1987) Analysis of subgrid- scale eddy viscosity with use of results from direct numerical simulations. Phys- ical Review Letters 58:547–550. Domaradzki JA, Radhakrishnan S (2005) Effective eddy viscosities in implicit modeling of decaying high Reynolds number turbulence with and without rota- tion. Fluid Dynamics Research 36:385–406. Domaradzki JA, Rogallo RS (1990) Local energy transfer and nonlocal interac- tions in homogeneous, isotropic turbulence. Physics of Fluids A: Fluid Dynam- ics 2:413–426. Domaradzki JA, Saiki EM (1997) A subgrid-scale model based on the estimation of unresolved scales of turbulence. Physics of Fluids 9:2148–2164. Domaradzki JA, Xiao Z, Smolarkiewicz PK (2003) Effective eddy viscosi- ties in implicit large eddy simulations of turbulent flows. Physics of Flu- ids 15:3890–3893. Dunca A, Epshteyn Y (2006) On the Stolz–Adams deconvolution model for the large-eddy simulation of turbulent flows. SIAM Journal on Mathematical Anal- ysis 37:1890–1902. Eyink GL (1994) Energy dissipation without viscosity in ideal hydrodynamics i. Fourier analysis and local energy transfer. Physica D: Nonlinear Phenom- ena 78:222–240. 134 Eyink GL (2005) Locality of turbulent cascades. Physica D: Nonlinear Phenom- ena 207:91–116. EyinkGL,AluieH(2009) Localnessofenergycascadeinhydrodynamicturbulence. i. Smooth coarse graining. Physics of Fluids 21:115107. Flad D, Beck A, Munz CD (2016) Simulation of underresolved turbulent flows by adaptive filtering using the high order discontinuous Galerkin spectral element method. Journal of Computational Physics 313:1–12. Foias C, Holm DD, Titi ES (2001) The Navier–Stokes-alpha model of fluid turbu- lence. Physica D: Nonlinear Phenomena 152-153:505–519. Germano M, Piomelli U, Moin P, Cabot WH (1991) A dynamic subgrid-scale eddy viscosity model. Physics of Fluids A: Fluid Dynamics 3:1760–1765. Geurts BJ (1997) Inverse modeling for large-eddy simulation. Physics of Flu- ids 9:3585–3587. Geurts BJ, Fröhlich J (2002) A framework for predicting accuracy limitations in large-eddy simulation. Physics of Fluids 14:L41–L44. Geurts BJ, Holm DD (2003) Regularization modeling for large-eddy simulation. Physics of Fluids 15:L13–L16. Geurts BJ, Holm DD (2006) Leray and LANS-α modelling of turbulent mixing. Journal of Turbulence 7:N10. Gibbon J, Holm D (2006) Length-scale estimates for the LANS- equations in terms of the reynolds number. Physica D: Nonlinear Phenomena 220:69–78. Gottlieb D, Hesthaven JS (2001) Spectral methods for hyperbolic problems. Jour- nal of Computational and Applied Mathematics 128:83–131. Graham JP, Holm DD, Mininni PD, Pouquet A (2007) Highly turbulent solu- tions of the Lagrangian-averaged Navier-Stokes α model and their large-eddy- simulation potential. Physical Review E 76. GrahamJP,HolmDD,MininniPD,PouquetA(2008) Threeregularizationmodels of the Navier–Stokes equations. Physics of Fluids 20:035107. Grinstein FF, Margolin LG, Rider WJ (2007) Implicit large eddy simulation: com- puting turbulent fluid dynamics Cambridge university press. Guermond J, Oden J, Prudhomme S (2003) An interpretation of the Navier–Stokes-alpha model as a frame-indifferent Leray regularization. Phys- ica D: Nonlinear Phenomena 177:23–30. 135 Hickel S, Adams NA, Domaradzki JA (2006) An adaptive local deconvolution method for implicit LES. Journal of Computational Physics 213:413–436. Holm DD, Marsden JE, Ratiu TS (1998a) Euler-poincaré models of ideal fluids with nonlinear dispersion. Physical Review Letters 80:4173–4176. Holm DD, Marsden JE, Ratiu TS (1998b) The euler–poincaré equations and semidirect products with applications to continuum theories. Advances in Math- ematics 137:1–81. Horiuti K (1989) The role of the Bardina model in large eddy simulation of tur- bulent channel flow. Physics of Fluids A: Fluid Dynamics 1:426–428. Hughes TJ, Feijóo GR, Mazzei L, Quincy JB (1998) The variational multi- scale method—a paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering 166:3–24. Hughes TJ, Mazzei L, Jansen KE (2000) Large eddy simulation and the variational multiscale method. Computing and Visualization in Science 3:47–59. Hughes TJ, Oberai AA, Mazzei L (2001) Large eddy simulation of turbulent chan- nel flows by the variational multiscale method. Physics of Fluids 13:1784–1799. Iliescu T, John V, Layton W, Matthies G, Tobiska L (2003) A numerical study of a class of LES models. International Journal of Computational Fluid Dynam- ics 17:75–85. Ilyin AA, Lunasin EM, Titi ES (2006) A modified-Leray-α subgrid scale model of turbulence. Nonlinearity 19:879–897. Jeanmart H, Winckelmans G (2007) Investigation of eddy-viscosity models modi- fied using discrete filters: A simplified “regularized variational multiscale model” and an “enhanced field model”. Physics of Fluids 19:055110. Jiang GS, Shu CW (1996) Efficient implementation of weighted ENO schemes. Journal of Computational Physics 126:202–228. Karamanos GS, Karniadakis G (2000) A spectral vanishing viscosity method for large-eddy simulations. Journal of Computational Physics 163:22–50. Karniadakis GE, Israeli M, Orszag SA (1991) High-order splitting methods for the incompressible Navier-Stokes equations. Journal of Computational Physics 97:414–443. Kerr RM, Domaradzki JA, Barbier G (1996) Small-scale properties of nonlinear interactions and subgrid-scale energy transfer in isotropic turbulence. Physics of Fluids 8:197–208. 136 Kokkinakis I, Drikakis D (2015) Implicit large eddy simulation of weakly- compressible turbulent channel flow. Computer Methods in Applied Mechanics and Engineering 287:229–261. Kopriva DA (2009) Implementing Spectral Methods for Partial Differential Equa- tions Springer Netherlands. Kraichnan RH (1967) Inertial ranges in two-dimensional turbulence. Physics of Fluids 10:1417. Kravchenko A, Moin P (1997) On the effect of numerical errors in large eddy simulations of turbulent flows. Journal of Computational Physics 131:310–322. Lamballais E, Fortuné V, Laizet S (2011) Straightforward high-order numerical dissipation via the viscous term for direct and large eddy simulation. Journal of Computational Physics 230:3270–3275. Lamorgese AG, Caughey DA, Pope SB (2005) Direct numerical simulation of homogeneous turbulence with hyperviscosity. Physics of Fluids 17:015106. Layton W, Lewandowski R (2008) A high accuracy Leray-deconvolution model of turbulence and its limiting behavior. Analysis and Applications 06:23–49. Lee M, Moser RD (2015) Direct numerical simulation of turbulent channel flow up to Re τ ≈ 5200. Journal of Fluid Mechanics 774:395–415. Lele SK (1992) Compact finite difference schemes with spectral-like resolution. Journal of Computational Physics 103:16–42. Leonard A (1975) Energy cascade in large-eddy simulations of turbulent fluid flows In Advances in geophysics, pp. 237–248. Elsevier. Leray J (1934) Sur le mouvement d’un liquide visqueux emplissant l’espace. Acta mathematica 63:193–248. Leslie DC, Quarini GL (1979) The application of turbulence theory to the formu- lation of subgrid modelling procedures. Journal of Fluid Mechanics 91:65. Li CG, Tsubokura M (2017) An implicit turbulence model for low-mach Roe scheme using truncated Navier–Stokes equations. Journal of Computational Physics 345:462–474. Lilly DK (1992) A proposed modification of the Germano subgrid-scale closure method. Physics of Fluids A: Fluid Dynamics 4:633–635. 137 Liu S, Meneveau C, Katz J (1994a) On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet. Journal of Fluid Mechanics 275:83. Liu XD, Osher S, Chan T (1994b) Weighted Essentially Non-oscillatory schemes. Journal of Computational Physics 115:200–212. Liu Y, Tucker P, Kerr R (2008) Linear and nonlinear model large-eddy simulations of a plane jet. Computers & Fluids 37:439–449. LohKC,DomaradzkiJA(1999) Thesubgrid-scaleestimationmodelonnonuniform grids. Physics of Fluids 11:3786–3792. Marsden JE, Shkoller S (2003) The anisotropic Lagrangian averaged Euler and Navier-Stokes equations. Archive for Rational Mechanics and Analy- sis 166:27–46. Montgomery DC, Pouquet A (2002) An alternative interpretation for the Holm “alpha model”. Physics of Fluids 14:3365–3366. Moser RD, Kim J, Mansour NN (1999) Direct numerical simulation of turbulent channel flow up to Re τ =590. Physics of Fluids 11:943–945. Nicoud F, Ducros F (1999) Subgrid-scale stress modeling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion 62:183–200. Nicoud F, Toda HB, Cabrit O, Bose S, Lee J (2011) Using singular values to build a subgrid-scale model for large eddy simulations. Physics of Fluids 23:085106. Pasquetti R, Xu CJ (2002) High-order algorithms for large-eddy simulation of incompressible flows. Journal of Scientific Computing 17:273–284. Peyret R (2002) Spectral Methods for Incompressible Viscous Flow Springer New York. Piomelli U, Cabot WH, Moin P, Lee S (1991) Subgrid-scale backscatter in turbu- lent and transitional flows. Physics of Fluids A: Fluid Dynamics 3:1766–1771. Piomelli U, Moin P, Ferziger JH (1988) Model consistency in large eddy simulation of turbulent channel flows. Physics of Fluids 31:1884. Piomelli U, Zang TA, Speziale CG, Hussaini MY (1990) On the large-eddy sim- ulation of transitional wall-bounded flows. Physics of Fluids A: Fluid Dynam- ics 2:257–265. Rebholz LG (2007) Conservation laws of turbulence models. Journal of Mathe- matical Analysis and Applications 326:33–45. 138 Rebholz LG, Kim TY, Byon YJ (2017) On an accurate α model for coarse mesh turbulent channel flow simulation. Applied Mathematical Modelling 43:139–154. Sagaut P (2006) Large Eddy Simulation for Incompressible Flows Springer, Berlin, Germany. Schranner FS, Domaradzki JA, Hickel S, Adams NA (2015) Assessing the numer- ical dissipation rate and viscosity in numerical simulations of fluid flows. Com- puters & Fluids 114:84–97. Schultz MP, Flack KA (2013) Reynolds-number scaling of turbulent channel flow. Physics of Fluids 25:025104. Smagorinsky J (1963) General circulation experiments with the primitive equa- tions: I. The basic experiment. Monthly weather review 91:99–164. Stolz S, Adams NA (1999) An approximate deconvolution procedure for large-eddy simulation. Physics of Fluids 11:1699–1701. Stolz S, Adams NA, Kleiser L (2001) An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded flows. Physics of Fluids 13:997–1015. Sun G, Domaradzki JA, Yang X, Chen K (2017) Assessing accuracy of CFD simulations through quantification of a numerical dissipation rate. paper 167, TSFP-10 . Sun G, Domaradzki JA (2018) Implicit LES using adaptive filtering. Journal of Computational Physics 359:380–408. Tantikul T, Domaradzki JA (2010) Large eddy simulations using truncated Navier–Stokes equations with the automatic filtering criterion. Journal of Tur- bulence 11:N21. Thiry O, Winckelmans G (2016) A mixed multiscale model better accounting for the cross term of the subgrid-scale stress and for backscatter. Physics of Fluids 28:025111. Thornber B, Mosedale A, Drikakis D (2007) On the implicit large eddy simulations of homogeneous decaying turbulence. Journal of Computational Physics 226:1902–1929. Toosi S, Larsson J (2017) Anisotropic grid-adaptation in large eddy simulations. Computers & Fluids 156:146–161. 139 Van Driest ER (1956) On turbulent flow near a wall. Journal of the Aeronautical Sciences 23:1007–1011. van Reeuwijk M, Jonker HJJ, Hanjalić K (2006) Incompressibility of the Leray-α model for wall-bounded flows. Physics of Fluids 18:018103. van Reeuwijk M, Jonker HJ, Hanjalić K (2009) Leray-α simulations of wall-bounded turbulent flows. International Journal of Heat and Fluid Flow 30:1044–1053. Verstappen R (2008) On restraining the production of small scales of motion in a turbulent channel flow. Computers & Fluids 37:887–897. Vreman AW (2004) An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications. Physics of Fluids 16:3670–3681. Winckelmans GS, Wray AA, Vasilyev OV, Jeanmart H (2001) Explicit-filtering large-eddy simulation using the tensor-diffusivity model supplemented by a dynamic smagorinsky term. Physics of Fluids 13:1385–1403. Zang Y, Street RL, Koseff JR (1993) A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows. Physics of Fluids A: Fluid Dynam- ics 5:3186–3196. 140

Abstract (if available)

Abstract In large eddy simulations (LES) with coarse mesh, various methods can be used to model the unknown effects of small unresolved scales of turbulent flows. Due to the dominant average effect of forward energy cascade from large to small scales, unphysical energy pileup generally occurs for under-resolved LES meshes. The subgrid scale dissipation can be accounted for by explicit subgrid scale (SGS) stresses. Many widely used SGS models adopt the eddy viscosity concept, which models energy transfer in analogy to the molecular diffusion based on the Boussinesq hypothesis. It is well-known, however, that some eddy viscosity closures exhibit relatively weak correlation with real SGS stress, so it is desirable to seek some alternative ways that prevent the buildup of small scale energy. ❧ Implicit large eddy simulations (ILES) are increasingly popular in recent years due to ease of implementation and the relaxation of some physical assumptions, where the numerical dissipation works in a manner similar to the explicit SGS models. If spectral methods are used the numerical dissipation is negligible but it can be introduced by applying a low-pass filter in the physical space, resulting in an effective ILES. In the present work we provide a comprehensive analysis of the numerical dissipation produced by different filtering operations in a turbulent channel flow simulated using a non-dissipative, pseudo-spectral Navier-Stokes solver. The amount of numerical dissipation imparted by filtering can be easily adjusted by changing how often the filter is applied. We show that when the additional numerical dissipation is close to the subgrid-scale (SGS) dissipation of an explicit LES the overall accuracy of ILES is also comparable, indicating that periodic filtering can replace explicit SGS models. A new method is proposed, which does not require any prior knowledge of a flow, to determine the filtering period adaptively. Once an optimal filtering period is found, the accuracy of ILES is significantly improved at low implementation complexity and computational cost. The method is general, performing well for different Reynolds numbers, grid resolutions, and filter shapes. ❧ Another option to remove the unphysical energy buildup is through explicit SGS models based on the multi-scale method. By separating different scales with filtering, we can directly remove some of the energy production near the LES cutoff, whose collective effect is similar as adding subgrid dissipation. The scale separation is facilitated by a smooth low-pass filter, which becomes increasingly more active for smaller resolved wavenumbers to reduce more production. A gradient-type modification is included to model extra dissipative effect and make the model Galilean invariant, which is found to have negligible impact to the original model. The model is compared with other nonlinear models and regularization techniques both theoretically and numerically. In a posteriori analysis of wall-bounded flow, our new model were able to provide sufficient SGS dissipation without a help of extra dissipative terms such as an eddy viscosity model. This explicit structural model makes it possible to evaluate the nonlinear energy transfer in LES more rigorously. Due to the removal of energy production in desired scales, the model consistently produces more accurate results than the similarity-type model being compared.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

On the simulation of stratified turbulent flows

PDF

Large eddy simulations of laminar separation bubble flows

PDF

A subgrid-scale model for large-eddy simulation based on the physics of interscale energy transfer in turbulence

PDF

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

Supersonic flow study of a capsule geometry using large-eddy simulation

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Particle-in-cell simulations of kinetic-scale turbulence in the solar wind

PDF

Development and evaluation of high‐order methods in computational fluid dynamics

PDF

Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers

PDF

Numerical study of focusing effects generated by the coalescence of multiple shock waves

PDF

Dynamic modeling and simulation of flapping-wing micro air vehicles

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Understanding the transport potential of nearshore tsunami currents

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders

PDF

Numerical simulations of linearly stratified flow around a sphere

Asset Metadata

Creator Sun, Guangrui (author)

Core Title Large eddy simulations of turbulent flows without use of the eddy viscosity concept

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publication Date 12/06/2019

Defense Date 08/05/2019

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag CFD,DNS,LES,numerical dissipation,OAI-PMH Harvest,turbulence

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Domaradzki, Julian (committee chair), Bermejo-Moreno, Ivan (committee member), Nakano, Aiichiro (committee member), Oberai, Assad (committee member), Sadhal, Satwindar (committee member)

Creator Email guangrus@usc.edu,mpesung@nus.edu.sg

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-246852

Unique identifier UC11675313

Identifier etd-SunGuangru-8013.pdf (filename),usctheses-c89-246852 (legacy record id)

Legacy Identifier etd-SunGuangru-8013.pdf

Dmrecord 246852

Document Type Dissertation

Rights Sun, Guangrui

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA