Characteristics of 1D spectra in finite-volume LES with ODT subgrid closure R. McDermott, A. Kerstein, R. Schmidt, and P. Smith November, 2005
1
Introduction
These are notes that accompany a talk to be given at the 2005 APS-DFD meeting in Chicago. The focus of the talk is to highlight the connection between the implied filter for energy-conserving LES numerics (e.g., the Harlow and Welch [3] staggered-grid scheme) and resolution of the spectral dip problem, which has been reported previously in LES/ODT simulations of isotropic turbulence. We will first present the finite-volume large-eddy simulation (FV-LES) formulation and highlight the equations to be solved. We next give a brief overview of the one-dimensional turbulence model (ODT) of Kerstein [5] and show how it is used as a subgrid-scale (SGS) model for FV-LES. We then review the “spectral dip problem” encountered in previous LES/ODT simulations [7, 9]. The reasons behind this problem are explained and a means of correction is proposed that turns out to work well for the transverse ODT spectrum. The correction leads to a statement of the “implied filter” transfer function for energy-conserving schemes, which is nontrivial and is claimed to be generally applicable regardless of the SGS model employed.
2
Finte-volume large-eddy simulation
In FV-LES we seek to model the evolution of “cell-average” fields. For the incompressible flow considered here the transport equation for the cell average is an integro-differential equation obtained by substitution of the following 1
definitions into an integral momentum balance around an Eulerian control volume (CV): 1 Z (k) U (x, t) ≡ U (x, t) dSk , (1) Sk Sk 1 Z 1 Z (k) Ü U (x, t) ≡ U (x, t) dx , (2) U (x, t) dV = V V h h where h is the LES grid spacing, which is herein taken to be uniform. Thus, Sk = h2 is the surface area of face k of a cubic CV of volume V = h3 . The volume-filtered momentum transport equation is Ü (k) dU Sk i = − nkj U i U j + P δij + τ ij + τijsgs . dt V
(3)
Here P is the modified pressure, τij is the viscous stress tensor, and the SGS stress tensor, obtained from decomposition of the nonlinear advection term, is defined to be τijsgs ≡ Ui U j − U i U j . (4) Note that this equation is valid for cell averages defined continuously in space. The “FV-LES equations” are obtained by sampling the continuous space at discrete intervals on the computational grid. Discrete continuity is enforced through the relation: X
(k)
Sk U j nj = 0 .
(5)
k
3
One-dimensional turbulence
A detailed discussion of ODT is beyond the scope of this paper. The interested reader is referred to [5]. Here we give a brief summary of the ODT concept. An “ODT field” represents a notional 1D line-of-sight through a fully-resolved (space and time) 3D flow field. Viscous transport along the ODT line direction is treated using molecular transport coefficients and the inviscid process of turbulent stirring along the the 1D line is modelled by measure-preserving stochastic mapping events that represent notional eddies in the turbulent flow. The location, size, and frequency of the eddy events are sampled from an eddy event-rate distribution that evolves in time with the flow. Hence, the distribution responds to the real-time conditions of the flow field and generates a cascade of length scales characteristic of turbulence. 2
4
LES/ODT coupling
When using ODT as the SGS closure model for LES we envision a lattice of ODT lines that intersect at LES CV centers. Note that for an N 3 LES there will be 3N 2 ODT lines (N 2 for each of the three coordinate directions), each resolved down to the Kolmogorov length scale. Due to this high resolution, the cost of LES/ODT scales as Re3/2 , which is a dramatic improvement over direct simulation (Re3 ). Eddy events along a line in a given direction contribute to the SGS momentum fluxes across the LES cell face normal to that direction. From a consistent initial state, all the ODT lines evolve independently for the period of one LES time step. At this point the flux divergence due to eddy events is compiled and used as the SGS force on the right hand side of a semi-implicit (i.e., fractional step) LES advancement. The LES field and ODT fields will no longer be consistent with each other because, with the exception of the local mean velocity gradient, the ODT field is oblivious to large-scale, 3D forces (e.g., mean advection) during the course of its time advancement cycle. The splitting error is corrected via a reconstruction procedure which adjusts the filtered ODT field to be consistent with the LES cell average. This is a correction to the local mean of the velocity field. Reconstruction of the variance is a work in progress. A key difference between stand-alone ODT and LES/ODT is that in the latter we truncate the eddy size distribution at roughly two times the LES grid spacing. This is because, in principle, LES can resolve eddies of this size. The exact value of the maximum eddy size, `max , is an adjustable parameter in the model and is analogous to the filter width, ∆, used in conventional LES modelling.
5
The spectral dip problem
The reconstruction step mentioned in the previous section is a key component of the algorithm and significantly effects the way the energy in the ODT field is distributed in wavespace near the LES Nyquist limit. To be clear, when we refer to the ODT field we mean a field, u(x), which contains low-wavenumber Ü (x), and high-wavenumber contributions contributions from the LES field, U from the SGS field, u0 (x). The reconstruction step forces the filtered value Ü (x ). Hence, of the ODT field to match the LES cell average: ue(xj ) = U j 3
u(x) = ue(x) + u0 (x). Note that the SGS field is defined such that uÜ0 (xj ) = 0 and hence the SGS field is undisturbed in the reconstruction. Figure 1 shows the 3D LES spectrum for an initial condition corresponding to a sharp spectral cutoff of the isotropic turbulence data of Comte-Bellot and Corrsin [1]. Figure 2 shows the corresponding longitudinal LES spectrum and the resulting longitudinal ODT spectrum after reconstruction. By forcing the ODT and LES cell means to match we imprint the LES 1D spectrum onto the ODT 1D spectrum up to the neighborhood of the LES Nyquist limit. For reasons that are discussed in the next section, if the 3D LES field is filtered with a spectral cutoff, tophat, or Gaussian filter, as is the common practice, then the LES 1D spectrum falls to zero too quickly and causes a dip in the spectrum obtained from the ODT field.
6
Relating 1D and 3D energy spectra
The “energy spectrum”, E(κ), so often discussed in turbulence, is obtained by integrating the spectral energy density over a spherical shell of radius √ κ ≡ κi κi in Cartesian wavenumber space. The longitudinal 1D energy spectrum, E11 (κ1 ), and transverse 1D energy spectrum, E22 (κ1 ), are obtained by integrating the appropriate energies over a plane perpendicular to the κ1 axis. For isotropic vector fields, the 1D and 3D spectra are related as follows (see, e.g., [8]): Z ∞ E(κ) κ21 E11 (κ1 ) = 1 − 2 dκ , (6) κ1 κ κ
κ2 1 Z ∞ E(κ) 1 + 12 E22 (κ1 ) = 2 κ1 κ κ
dκ .
(7)
For our purposes the important observations to be made about these relations are: 1. The 1D spectra must be monotonically decreasing functions of κ1 . 2. If the 3D spectrum is truncated (i.e., filtered) then it cannot contribute energy to the 1D spectrum above the wavenumber where E(κ) falls off. 3. Due to the (±) sign difference in the integrands of (6) and (7), we will find that we cannot force smooth transitions across the Nyquist limit for both longitudinal and transverse ODT spectra simultaneously. 4
4. Again, due to the sign difference, the dip in the longitudinal spectrum (E11 [κ1 ]) will be more pronounced than the dip in the transverse spectrum (E22 [κ2 ]). By examining the coupling behavior of the conventional filter kernels (i.e., cutoff, tophat, and Gaussian) it is apparent that none of these will suffice for generating a smooth transition across the LES Nyquist limit for the ODT spectra.
7
Behavior of energy-conserving schemes in the inviscid limit
Consider simulation of an unforced flow in a periodic cubic domain. The flow may be initialized arbitrarily subject to the periodic boundary condition. In practice we simply use an initial condition from one of our decaying isotropic turbulence simulations. We turn off the SGS model and set the viscosity to zero. What is the long-time behavior of such a simulation, assuming the numerical scheme identically conserves energy? The answer is: that the energy becomes uniformly distributed in the 3D Cartesian wavespace. Given the definition of the 3D energy spectrum (an integral over a spherical shell of radius κ) a uniform distribution implies that, up to the grid Nyquist limit (κc = π/h), the [implied] filtered energy density simply goes as the surface area of the spherical shell. That is, E(κ) ∼ κ2 . Beyond κc the energy of the shell is clipped by the Cartesian box that is inscribed by the sphere of radius κc (see Figure 4). Thus, the filtered √ energy decays, and, as is shown below, between the limits κc ≤ κ ≤ 2κc , the surface area of the clipped sphere, and hence the √filtered energy √ density, goes 2 as E(κ) ∼ κ (3κc /κ − 2). Between the limits 2κc < κ ≤ 3κc the kernel formula is much more complicated, but it is a rapidly decreasing function of κ,√and so for practical purposes we will at present assume E(κ) = 0 for κ > 2κc . Now (sidestepping the issue of infinite energy on an infinite domain), the unfiltered energy density should scale like E(κ) ∼ κ2 as κ becomes large. Noting that the relationship between the filtered energy spectrum and the Ò2 (κ)E(κ) where G Ò2 (κ) is the “attenuation unfiltered spectrum is E(κ) = G factor” for the filter and is the square of the filter “transfer function” (see, e.g., [8] for a detailed discussion), we have that the transfer function for the 5
energy-conserving numerical scheme is 8 > 1 > < r
Ò G(κ) => > :
for κ ≤ κc , √ κc 3 − 2 for κc < κ ≤ 2κc , κ √ 0 for 2κc < κ .
(8)
Again, the first two parts of the piecewise function (8) are exact and the last is an approximation, which we will see has no adverse effect. The derivation of the attenuation factor for the middle section of the implied filter is completely geometrical and proceeds as follows. As the sphere begins to grow beyond the bounding box in wavespace, there are six (one for each face of the box) equal solid angles that are generated and must be subtracted from the total surface area of the sphere. Consider one of the solid angles. Here we use as θ the polar angle and φ is the azimuthal angle. Note that for a given spectral radius, κ, the maximum value of the azimuthal angle is given by φmax = cos−1 (κc /κ). The surface area is then given by Z
Ssolid angle = =
S
dS ,
Z 2π Z φmax 0
= κ2
0
κ dφ × κ sin(φ) dθ ,
Z 2π Z φmax 0
0
sin(φ) dφ dθ ,
= κ2 2π [− cos(φmax ) + cos(0)] , = κ2 2π − cos(cos−1 {κc /κ}) + 1 , κc . = κ2 2π 1 − κ
(9)
The area of the truncated sphere is now given by SE(κ) = Ssphere − 6 × Ssolid angle , κc = 4πκ2 − 6 × 2πκ2 1 − , κ κc = 4πκ2 1 − 3 1 − , κ κc = 4πκ2 3 − 2 . κ
(10)
The Fourier transform of the transfer function gives the physical space filter kernel, G(r; ∆), that one is accustom to seeing used to define the LES 6
equations. The velocity field is typically defined by the convolution ZZ∞ Z
U (x, t) ≡
G(r; ∆)U (x − r, t) dr .
(11)
−∞
For the filter conceptualized here the kernel becomes ∞
3 1 ZZZ −iκ·x Y G(r; ∆) ≡ e H(κc − |κ(j) |) dκ . (2π)3 −∞ j=1 κc Zκc Zκc
1 Z = (2π)3 −κ
c
=
3 Y
1 j=1 2π 3 Y
1 = j=1 2π =
e−i(κ1 r1 +κ2 r2 +κ3 r3 ) dκ1 dκ2 dκ3 .
−κc −κc
Z κc −κc
e−iκ(j) r(j) dκ(j) ,
−1 ir(j)
!
e−iκ(j) r(j)
3 Y
sin(κc r(j) ) , πr(j) j=1
κc −κc
, (12)
where κc = π/hLES and the last step makes use of trigonometric function, sin(x) = (eix −e−ix )/(2i) . It is to be emphasized that this is not the standard “sharp spectral cutoff” filter typically employed, which clips energy at the spectral radius κ = κc . The above formula is anisotropic and clips energy at the edge of the box centered at the origin in wavespace with side 2κc . It should be further noted that this is not the exact transform pair of the proposed transfer function due to the approximation made in (8). A key question becomes: Is the implied filter a projector? That is, does the filtered function change upon repeated application of the attenuation factor? Well, from the very source of our analysis we can see that repeated application of the numerical procedure does not change the shape of the filtered energy spectrum. Therefore, we know that the true implied filter is a projector. On the other hand, even once the third section of the proposed attenuation factor is worked out, it seems that (8) is not a projector. It is not entirely clear what is going on here. The answer most likely has to do with aliasing, and a continuous flux of spectral energy, even at the steady state. 7
8
Behavior of the dynamic model in FV-LES
As we will see, the proposed filter is not model specific. In this section we present results from an FV-LES using the dynamic model [2] as the SGS closure. The spectra are compared against the decaying isotropic turbulence data of Kang et al. [4]. Two different initializations are considered. The first is using a spectral cutoff filter. The other is using the proposed implied filter of energy-conserving schemes. Recall that no explicit filtering operation is employed with FV-LES. The simulation initialized with the cutoff filter clearly tries to approach the shape of the implied filter as time progresses and the nonlinear term begins to fill the corners of the spectral space. The second case, where the field is initialized with the implied filter, clearly retains its shape.
9
Results for FV-LES with ODT SGS closure
In this section we show how the implied filter impacts the spectral transition for the longitudinal and transverse ODT spectra. We also show 1D and 3D spectral results for simulations of decaying isotropic turbulence using ODT as the SGS closure. The ODT field is initialized to match the spectral data of Kang et al. [4]. The LES field is initialized to match the filtered field implied by the transfer function (8). At this point the LES and ODT fields share the same low-wavenumber spectral fields but their physical space fields are completely uncorrelated. Therefore, we perform a mean reconstruction procedure as outlined in Section 9. The resulting initialized transverse and longitudinal ODT spectra can be seen as the green curves in Figures 7 and 8, respectively. Clearly, the transverse spectrum is accurate, while the longitudinal spectrum still dips near the LES Nyquist limit. At the present time there seems to be no remedy for this dilemma. Also shown in Figures 7 and 8 are the results for the simulations compared with the last time station of the experiment. There are several things to note. First, the high-wavenumber portion of the spectrum could not be resolved with ODT and so we used a scaled viscosity that causes the ODT spectra to roll off faster than in the actual experiment. In this simulation the number of ODT points per line was NODT = 4096. The LES resolution was 323 , putting the total number of ODT points in the simulation at 3(32)2 × 4096 = 8
12, 582, 912 (This simulation runs on a desktop PC in a matter of minutes). The next thing to point out is that the experimental data for the transverse spectra are not monotonically decreasing. This is a problem with the data, not the simulation. So, one should not be alarmed that the data and simulation do not match in this range. With regard to the quality of the ODT results, it should be noted that the smoothness of the transition across the LES Nyquist limit is consistent for the duration of the simulation. That is, the longitudinal spectrum, which starts out with a dip, retains its dip throughout the run. The transverse spectrum starts out very smooth and retains this smoothness. It may be seen that in both cases the energy of the simulated spectra is lower than that of the data in the inertial range. We suspect this to be due to an energy conservation problem in the LES/ODT algorithm. This is currently being corrected. The resulting LES spectra can be seen in Figure 9. The results are quite encouraging. Clearly, as was the case with the dynamic model, the 3D spectra retain the implied filter shape throughout the simulation. The model coefficient used for this case was determined through an equilibrium analysis as shown in [6]. The maximum eddy size for the ODT model was set to be `max = 2hLES .
10
Conclusions
Here we have shown that, due to the mean reconstruction procedure, the shape of the LES 3D spectrum, E(κ), significantly impacts the ODT spectrum across the region near the LES Nyquist limit. The implied filter shape for energy-conserving numerical schemes, introduced here for the first time, adds just the right amount of energy on the high-wavenumber side of the Nyquist limit to rebuild the transverse ODT spectrum and form a smooth transition. Due to the (±) difference in the formulas (6) and (7) the transition is still not smooth for the longitudinal spectrum. This observation leads one to conclude that for a given LES 3D spectrum, the shape of which is determined a priori for energy-conserving numerics, it will be extremely difficult to design a reconstruction procedure that simultaneously produces smooth transitions for both longitudinal and transverse ODT spectra. Perhaps the key will be to limit the impact of the reconstruction. This is currently a research issue for LES/ODT. 9
The preceding sections bring to bear a key observation about energy conserving finite-volume schemes, and that is that the implied filter is not a box filter. The box filter is poor at attenuating high-wavenumber modes and attenuates low-wavenumber modes too strongly. The implied filter is much closer in character to a spectral cutoff filter. As a final note, it is to be emphasized that the implied filter presented here was deduced purely from idea that the energy becomes uniform in spectral space at long times for energy-conserving numerical schemes. It follows that any scheme that conserves energy in the inviscid limit should obey a similar type of filter. The exact form of the kernel will, of course, depend on the geometry of the domain. The form presented here is based on a cubic domain. It has been demonstrated that the shape does not depend upon the SGS closure for the cases of the dynamic procedure and ODT, two vastly different models. It is suspected, therefore, that the concept applies generally.
Acknowledgements Funding for this work was provided through the Department of Energy Computational Science Graduate Fellowship (DE-FG02-97ER25308).
References [1] G. Comte-Bellot and S. Corrsin. Simple Eularian time correlation of fulland narrow-band velocity signals in grid-generated, ‘isotropic’ turbulence. J. Fluid Mech., 48:273–337, 1971. [2] M. Germano, U. Piomelli, P. Moin, and W. Cabot. A dynamic subgridscale eddy viscosity model. Phys. Fluids A, 3(7):1760–1765, 1991. [3] F.H. Harlow and J.E. Welch. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. Phys. Fluids, 8:2182, 1965. [4] H.S. Kang, S. Chester, and C. Meneveau. Decaying turbulence in an active-grid-generated flow and comparisons with large-eddy simulation. J. Fluid Mech., 480:129–160, 2003.
10
[5] A.R. Kerstein, W.T. Ashurst, S. Wunsch, and V. Nilsen. One-dimensional turbulence: vector formulation and application to free shear flows. J. Fluid Mech., 447:85–109, 2001. [6] R. McDermott, A. Kerstein, R. Schmidt, and P. Smith. The ensemble mean limit of the one-dimensional turbulence model and application to residual stress closure in finite-volume large-eddy simulation. Journal of Turbulence, in press. [7] Randall J. McDermott. Toward one-dimensional turbulence subgrid closure for large-eddy simulation. PhD thesis, The University of Utah, 2005. [8] Stephen B. Pope. Turbulent Flows. Cambridge, 2000. [9] R.C. Schmidt, R. McDermott, and A. Kerstein. ODTLES: A Model for 3D Turbulent Flow Based on One-dimensional Turbulence Modeling Concepts. Sandia National Laboratories Report 2005–0206, 2005.
11
−3
10
−4
E(k), m3/s2
10
−5
10
CBC data, t=0 Cutoff LES LES Nyquist limit
−6
10
1
10
2
10 k, 1/m
3
10
Figure 1: LES initialization with a sharp spectral cutoff at the grid Nyquist limit.
12
−3
10
−4
−5
10
11
1
E (k ), m3/s2
10
−6
10
−7
CBC data LESODT LES LES Nyquist limit
10
−8
10
1
10
2
10 k1, 1/m
3
10
Figure 2: Longitudinal spectrum of the LES and ODT fields after the initial reconstruction.
13
Figure 3: The green connected dots are the initial spectrum. The red connected dots are the 3D spectrum for energy-conserving scheme in the inviscid limit at long times. The black line is the curve ∼ κ2 and the blue curve is ∼ 3κc /κ − 2. The Nyquist limit is at wavenumber κ = 16 1/m. √ Notice that the blue curve departs from the red curve precisely at κ = 2(16) = 22.3 1/m.
14
Figure 4: Spherical shell of radius κ =
15
√
2κc .
Figure 5: DSM with sharp spectral cutoff initialization
16
Figure 6: DSM with implied filter initialization
17
0
10
−1
10
−2
−3
10
22
1
E (k ), m3/s2
10
−4
10
Kang data (t=0.42s) LESODT LES Nyquist limit min eddy size
−5
10
−6
10
0
10
1
2
10
10
3
10
k1, 1/m
Figure 7: Transverse 1D spectra for ODT field for FV-LES/ODT using `max = 2 and implied filter LES initialization.
18
Figure 8: Longitudinal 1D spectrum for initial ODT field with implied filter imposed on 3D LES field.
19
0
10
−1
E(k), m3/s2
10
−2
10
Kang data LES LES Nyquist limit max eddy size
−3
10
0
1
10
10 k, 1/m
Figure 9: 3D LES spectra for FV-LES/ODT using `max = 2.
20
