- Email: [email protected]
Large-eddy simulation of turbulent channel flows with conservative IDO scheme
Journal of Computational Physics 230 (2011) 5787–5805
Contents lists available at ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Large-eddy simulation of turbulent channel flows with conservative IDO scheme Naoyuki Onodera a,⇑, Takayuki Aoki b, Hiromichi Kobayashi c,d a
Department of Energy Sciences, Tokyo Institute of Technology, 2-12-1 O-okayama, Meguro-ku, Tokyo 152-8550, Japan Global Scientific Information and Computing Center, Tokyo Institute of Technology, 2-12-1 O-okayama, Meguro-ku, Tokyo 152-8550, Japan c Department of Physics, Keio University, 4-1-1 Hiyoshi, Kouhoku-ku, Yokohama 223-8521, Japan d Research and Education Center for Natural Sciences, Keio University, 4-1-1 Hiyoshi, Kouhoku-ku, Yokohama 223-8521, Japan b
a r t i c l e
i n f o
Article history: Received 3 August 2010 Received in revised form 2 April 2011 Accepted 3 April 2011 Available online 12 April 2011 Keywords: Multi-moment scheme IDO scheme Large-eddy simulation Coherent structure model
a b s t r a c t The resolution of a numerical scheme in both physical and Fourier spaces is one of the most important requirements to calculate turbulent flows. A conservative form of the interpolated differential operator (IDO-CF) scheme is a multi-moment Eulerian scheme in which point values and integrated average values are separately defined in one cell. Since the IDO-CF scheme using high-order interpolation functions is constructed with compact stencils, the boundary conditions are able to be treated as easy as the 2nd-order finite difference method (FDM). It is unique that the first-order spatial derivative of the point value is derived from the interpolation function with 4th-order accuracy and the volume averaged value is based on the exact finite volume formulation, so that the IDO-CF scheme has higher spectral resolution than conventional FDMs with 4th-order accuracy. The computational cost to calculate the first-order spatial derivative with non-uniform grid spacing is one-third of the 4th-order FDM. For a large-eddy simulation (LES), we use the coherent structure model (CSM) in which the model coefficient is locally obtained from a turbulent structure extracted from a second invariant of the velocity gradient tensor, and the model coefficient correctly satisfies asymptotic behaviors to walls. The results of the IDO-CF scheme with the CSM for turbulent channel flows are compared to the FDM with the CSM and dynamic Smagorinsky model as well as the direct numerical simulation (DNS) by Moser et al. Adding the sub-grid scale stress tensor of LES to the IDOCF scheme improves the profile of the mean velocity in comparison with an implicit eddy viscosity of the IDO-CF upwind scheme. The IDO-CF scheme with the CSM gives better turbulent intensities than conventional FDMs with the same number of grid points. The turbulent statistics calculated by IDO-CF scheme are in good agreement with the DNS at the various values of Reynolds number Res = 180,395, and 590. It is found that the IDO-CF scheme is suitable for the turbulent flow computation and improves the turbulent statistics with compact stencils. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Large-eddy simulation (LES) is one of successful approaches for unsteady turbulent flows [1,2]. LES resolves the dynamics of large-scale structures on the grid scale (GS), and the effect of smaller-scale turbulent structures are taken into account by using the subgrid-scale (SGS) model. A number of SGS models have been proposed as described below. In order to apply the ⇑ Corresponding author. Tel.: +81 3 5734 3575. E-mail address: [email protected] (N. Onodera). 0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.04.004
5788
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
SGS model to complex flows, however, it is preferable that the SGS model can be applied to wall bounded flows with local determination of its model coefficient. The Smagorinsky model (SM) [3] is a prominent SGS model based on a concept of eddy viscosity. The SGS stress tensor Þ2 S and the resolved strain rate tensor Sij , where CS is sij ¼ 2mSGS Sij is explicitly modeled by the eddy viscosity mSGS ¼ ðC S M ffiffiffiffiffiffiffiffiffiffiffiffi q is the filter width, and S ¼ 2Sij Sij . In the SM, the undetermined parameter CS is kept constant the Smagorinsky constant, M in a computational domain, although the model parameter should be changed among various turbulent flows. The SM gives a non-zero eddy viscosity in laminar-flow regions, and requires a wall function to damp its SGS viscosity in wall-bounded flows [4,5]. The dynamic Smagorinsky model (DSM) [6,7] is proposed to overcome these deficiencies by locally calculating the model coefficient. In addition to the grid filter, the test filter is applied to the incompressible Navier–Stokes equations. The model parameter is determined only by the ratio of the grid filter width and the test filter width. This model can be applied to various kinds of flows without empirical parameters. The DSM is the most notable breakthrough in LES, and many SGS models based on the DSM are proposed [8,9]. As one of the weak points, the model parameter introduces negative SGS eddy viscosity, which causes numerical instability. Since the clipping and averaging of the model parameter are required to avoid breakdown, the DSM is difficult to apply to complex geometries. The Lagrangian dynamic models [10,11] calculate the model parameter without averaging in homogeneous directions. The model parameter is time-averaged along fluid path lines. These models are applicable to complex geometries, however it takes more computational time in comparison with the DSM. The coherent structure model (CSM) [12] is one of remedies for these problems. A turbulent structure, which is extracted from a second invariant of the velocity gradient tensor, determines a model parameter locally. As a result, the SGS eddy viscosity becomes zero near the wall without a wall-damping function. In addition, the model parameter is guaranteed to be positive, and its variance is small, so that averaging is not necessary. The CSM shows good performance in complex geometries [13]. In this paper, we apply the CSM to computations of turbulent phenomena. 4th- or higher-order numerical schemes are necessary for turbulent simulations because low-order discretization error exceeds the amount of the SGS stress [14]. The Morinishi scheme [15], one of high-order finite difference methods (FDMs), is often applied to LES. The 4th-order one uses 7-point stencil in staggered grid, which makes it difficult to handle the wall boundary condition. On the other hand, multi-moment schemes achieve high-order accuracy with compact stencils. The multi-moment schemes define different moments in addition to a point value for a dependent variable like spectral difference methods [16]. A high-order interpolation function is constructed locally by using several neighbor moments. The cubic-interpolated pseudo-particle (CIP) scheme [17,18] is the originally proposed multi-moment scheme. First-order spatial derivative values and point values are defined at the same grid point, and the advection equation is solved with high accuracy using a semiLagrangian formulation. The interpolated differential operator (IDO) scheme [19,20] is a multi-moment Eulerian scheme. High-order accuracy is obtained not only in the advection equation but also in more general equations, such as a diffusion, pressure and Poisson equation [21]. In addition, it is confirmed that those multi-moment schemes have better spectral resolution than the conventional FDM with the same order accuracy [22,23]. The conservation property is necessary to carry out numerical experiments for stable long-time integration. A conservative semi-Lagrangian CIP (CIP-CSL) scheme [24,25] and a conservative formulation of IDO (IDO-CF) scheme [26] guarantee exactly mass conservation. It is confirmed that the truncation error of the IDO-CF scheme is smaller than that of conventional FDMs [26]. In this paper, the IDO-CF scheme with the CSM is proposed and studied for turbulent channel flows which has been studied for several decades as one of the most basic benchmarks [27,28]. The IDO-CF scheme can handle the wall boundary condition with high-order accuracy. The IDO-CF scheme has twice the number of degrees of freedom compared to conventional FDMs in every direction. Therefore, numerical tests are examined with the same total number of degrees of freedom. The turbulent statistics of the IDO-CF scheme with the CSM are assessed and compared to those of the FDM with the DSM and the CSM as well as the direct numerical simulation (DNS). This paper is organized as follows. The IDO-CF scheme is reviewed, and easiness to handle the boundary condition and the computational cost with non-uniform grid spacing are discussed in Section 2. The accuracy and stability are studied in comparison with the FDM in Section 3. The governing equation and SGS models are presented in Section 4. Numerical results for the turbulent channel flows are given in Section 5. The paper ends with some concluding remarks in Section 6.
2. Conservative form of IDO scheme The turbulent flow calculation requires high-order discretization, and the given partial differential equations (PDEs) have been generally discretized by the FDM, FVM, or finite element method (FEM). Multi-moment schemes formulated as the different approach realize high-order accuracy with compact stencils. Unlike conventional FDMs, they define several kinds of moments in one cell for each dependent variable. In this section, the IDO-CF [26], one of multi-moment schemes, is compared to the FDM with the same order of accuracy. The three major differences are the configuration of variable moments, the computational cost with non-uniform grid spacing, and the spectral resolution in Fourier space.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5789
2.1. Formulation and spatial discretization In the section, we give discussions of the IDO-CF scheme in one-dimensional case for simplicity and they are available for two- and three-dimensional cases. The IDO-CF scheme defines the value fi = f(xi) and its integrated average value R x fiþ1=2 ¼ D1x f ðxÞdx. Fig. 1 illustrates the locations of the variables in a one-dimensional case. Both the variables are time-integrated simultaneously by solving different discretized equations. For solving PDEs, spatial derivatives fx(x), fxx(x), . . . , fx(n)(x) are obtained by differentiating of an interpolation function. One of the interpolation functions is the following fourth-order polynomial for the local domain xi1 6 x 6 xi+1.
F c ðxÞ ¼ ac4 X 4 þ ac3 X 3 þ ac2 X 2 þ ac1 X þ ac0 ;
fxc
ðkÞ
¼ k!ack ;
ð1Þ
where X = x xi and the superscript ‘‘c’’ represents a central scheme. The interpolation function is called the fourth-order central interpolation function. For simplicity, we demonstrate the case of the uniform grid spacing Dx = xx+1 xi. The coefficients of the interpolation function are determined with the five matching conditions of Fc(xi) = fi, Fc(xi±1) = fi±1, and R Dx c 1 F ðxÞdx ¼ x fi1=2 . The spatial derivatives at x = xi are represented as follows Dx 0
ðfiþ1 fi1 Þ 4 x fiþ1=2 x fi1=2 ¼ ; 2 Dx 3ðfiþ1 þ 8f i þ fi1 Þ 15 x fiþ1=2 þ x fi1=2 c fxx;i ¼ ; 2 Dx2 6ðfiþ1 fi1 Þ 12 x fiþ1=2 x fi1=2 c ð3Þ fx;i ¼ ; Dx3 30ðfiþ1 þ 4f i þ fi1 Þ 90 x fiþ1=2 þ x fi1=2 c ð4Þ fx;i ¼ : Dx4
c fx;i
ð2Þ
ð3Þ
ð4Þ
ð5Þ
The first- and second-order spatial derivatives have 4th-order accuracy (O(Dx4)). c fx;i ¼ fx;i
1 1 ð5Þ 4 f Dx ; 3 5! x;i
c ¼ fxx;i fxx;i
ð6Þ
6 1 ð6Þ 4 f Dx : 7 6! x;i
ð7Þ
The major error of the first-order derivative comes from dispersion errors (fx(5)Dx4), and does not contain numerical dissipation. In solving the incompressible Navier–Stokes equations, a convection term, stress tensor, and the Poisson equation are discretized by those expressions (Eqs. (2) and (3)). Another interpolation function is the upwind interpolation function (IDO-CF upwind) constructed for the local domain xi1 6 x 6 xi (ui > 0) up up 2 F up ðxÞ ¼ aup 2 X þ a1 X þ a0
ð8Þ up
up
1 Dx
R0
up
x
with the three matching conditions of F (xi) = fi, F (xi1) = fi1, and F ðxÞdx ¼ fi1=2 . The superscript ‘‘up’’ represents Dx a upwind scheme. The spatial derivatives are also represented by the following approximations, up fx;i ¼
ð4f i þ 2f i1 6x fi1=2 Þ ; Dx
up ¼ fxx;i
ð6f i þ 6f i1 12x fi1=2 Þ : Dx2
ð9Þ
Fig. 1. Fourth-order interpolation function of the IDO-CF scheme.
ð10Þ
5790
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
The discretization errors are found to be up fx;i ¼ fx;i
1 1 ð3Þ 2 4 1 ð4Þ 3 f Dx þ f Dx ; 2 3! x;i 5 4! x;i
up fxx;i ¼ fxx;i 3
ð11Þ
1 ð3Þ 18 1 ð4Þ 2 f Dx þ f Dx : 3! x;i 5 4! x;i
ð12Þ
The first-order spatial derivative, derived from the upwind interpolation function, has 2nd-order accuracy (O(Dx2)) and is rewritten by the combination of the spatial derivative of the central interpolations (Eqs. (2), (4), and (5)). up c fx;i ¼ fx;i
1 1 c ð3Þ 2 4 1 c ð4Þ 3 Dx þ Dx f f 2 3! x;i 5 4! x;i
ðu > 0Þ:
ð13Þ
The last term of the fourth-order spatial derivative works as a numerical viscosity. The first-order spatial derivative of the cell-integrated average value is derived as follows. It consists of the cell surface values. x
fx;iþ1=2
@ 1 ¼ @x Dx
Z
!
xiþ1
FðxÞdx
¼
xi
fiþ1 fi : Dx
ð14Þ
This is exactly the same as the FVM formulation. The second-order spatial derivative is also the similar formulation, x
@2 1 ¼ 2 @x Dx
fxx;iþ1=2
Z
xiþ1
! FðxÞdx
xi
¼
fx;iþ1 fx;i ; Dx
ð15Þ
where we apply Eq. (2) to the first-order spatial derivatives. The conservative property of the cell-integrated average is completely guaranteed. For the FDM, the variables are defined at the grid points illustrated in Fig. 3. The first-order spatial derivative with 4thorder accuracy is represented as c fx;i ¼
fiþ2 8f iþ1 þ 8f i1 fi2 12h
ðFDMÞ;
ð16Þ
where h = xi+1 xi. The grid space h represents the distance between neighbor values. In the two-dimensional case, four kinds of variable moments are defined such as fi;j ¼ f ðxi ; yj Þ; x fiþ1=2;j ¼ R R RR 1 f ðx; yÞdx; y fi;jþ1=2 ¼ D1y f ðx; yÞdy, and xy fiþ1=2;jþ1=2 ¼ Dx1Dy f ðx; yÞdx dy in Fig. 5. Each variable is time-integrated simultaDx neously. The eight kinds of variable moments are defined in the three-dimensional case. The formulation for system equations in multi-dimensional cases was presented in the paper [26]. 2.2. Boundary treatment For the IDO-CF scheme, the boundary is imposed on xwall = x1 as shown in Fig. 2. The values f0, xf1/2, and f1 are given by the boundary condition. The first-order spatial derivatives xfx,3/2 and fx,2 are represented as follows, x
fx;3 ¼ 2
f2 f1 ; Dx
c fx;2 ¼
ð17Þ
ðf3 f1 Þ 4 x f5=2 x f3=2 : 2 Dx
ð18Þ
Fig. 2. Boundary treatment of the IDO-CF scheme.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5791
Eqs. (17) and (18) do not refer the outside values f0 and xf1/2. The boundary condition of the FDM is imposed on xwall = x2 as shown in Fig. 4. The corresponding first-order spatial derivative with 4th-order accuracy at x3 is represented as c fx;3 ¼
f5 8f 4 þ 8f 2 f1 12h
ðFDMÞ:
ð19Þ
The value f1 is located outside the domain. Eq. (19) is subjected to the boundary conditions. For that reason, low-order (2ndorder) or less stable one-side discretization is often used at the wall boundaries in FDMs. 2.3. Spatial discretization with non-uniform grid spacing Turbulent shear stress augments near the wall, and the spatial- and time-scales become much smaller than those away from the wall. In order to resolve the near-wall phenomena, the non-uniform grid spacing is commonly used in association with complicated derivation and computational cost. The procedure to derive the IDO-CF scheme with non-uniform grid spacing is the same as that with uniform grid spacing. Rx The point value is defined as fi = F(xi), and the volume-integrated average value is also defined as x fiþ1=2 ¼ x 1x xiiþ1 FðxÞdx. iþ1 i The first-order product of coefficients vector and variables vector as spatial derivative atTx =xi is expressed by theinner T fx;i ¼ ~ ci ~ f i ¼ ci1 ; ci1=2 ; ci ; ciþ1=2 ; ciþ1 fi1 ; x fi1=2 ; fi ; x fiþ1=2 ; fiþ1 . The coefficients of the central interpolation function are written as
ci1 ¼
2Dx2i
Dxi1 ðDxi þ Dxi1 Þ2
ci1=2 ¼
ci ¼
ð20Þ
;
Dx2i ð6Dxi þ 10Dxi1 Þ Dxi1 ðDxi þ Dxi1 Þ3
ð21Þ
;
4ðDxi Dxi1 Þ ; Dxi1 Dxi
ciþ1=2 ¼
ð22Þ
Dx2i1 ð10Dxi þ 6Dxi1 Þ
ciþ1 ¼
Dxi ðDxi þ Dxi1 Þ3 2Dx2i1
Dxi ðDxi þ Dxi1 Þ2
ð23Þ
;
ð24Þ
;
where the grid spacing is denoted as Dxi = xi+1 xi. The first-order spatial derivative is found to be composed of two types of grid spacing (Dxi1 and Dxi). The computational costs are evaluated as the number of floating-point operations (FLOP) for the first-order derivative at the point as
FLOPIDOPV ci1 ; ci1=2 ; ci ; ciþ1=2 ; ciþ1 ¼ ð7; 13; 4; 12; 8Þ ¼ 51: fx
ð25Þ
For the integrated average value of the IDO-CF scheme, Eq. (14) does not change. The FLOP of the first-order derivative on the cell is shown as
FLOPIDOIA ¼ 2: fx
ð26Þ
Each variable is time-integrated simultaneously. The total FLOP with the number of cell N is counted as
Total FLOPIDO ¼ FLOPIDOPV N þ FLOPIDOIA N: fx fx fx
Fig. 3. Fourth-order interpolation function of the FDM.
ð27Þ
5792
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
Fig. 4. Boundary treatment of the FDM.
Fig. 5. Definition of moments in two dimension.
Since Eq. (27) contains 2N independent variables, the FLOP per variable moment of the IDO-CF scheme is simply estimated as
FLOPIDO ¼ fx
Total FLOPIDO FLOPIDOPV þ FLOPIDOIA fx fx fx ¼ ¼ 26:5: 2N 2
ð28Þ
In the case of the spatial discretization of the FDM, all the variables are defined at grid points as fi = F(xi). The variables T vector and the coefficient vector are similarly expressed as ~ ci ¼ ðci2 ; ci1 ; ci ; ciþ1 ; ciþ2 ÞT and ~ f i ¼ ðfi2 ; fi1 ; fi ; fiþ1 ; fiþ2 Þ . The coefficients of the central interpolation function are written as
ci2 ¼
hi1 hi ðhiþ1 þ hi Þ ; hi2 ðhi1 þ hi2 Þðhi þ hi1 þ hi2 Þðhiþ1 þ hi þ hi1 þ hi2 Þ
ci1 ¼
ci ¼
ð29Þ
ðhi1 þ hi2 Þhi ðhiþ1 þ hi Þ ; hi2 hi1 ðhi þ hi1 Þðhiþ1 þ hi þ hi1 Þ
ð30Þ
2 2 2 2 hi2 hi hiþ1 hi1 hiþ1 hi2 hi1 hiþ1 þ hi2 hi þ 2 hi1 hi hiþ1 þ hi1 hi hi1 hi hi2 hi1 hi
ciþ1 ¼
hi1 ðhi1 þ hi2 Þhi ðhiþ1 þ hi Þ hi1 ðhi1 þ hi2 Þðhiþ1 þ hi Þ ; hi ðhi þ hi1 Þðhi þ hi1 þ hi2 Þhiþ1
ciþ2 ¼
hi1 ðhi1 þ hi2 Þhi hiþ1 ðhiþ1 þ hi Þðhiþ1 þ hi þ hi1 Þðhiþ1 þ hi þ hi1 þ hi2 Þ
;
ð31Þ
ð32Þ
ðFDMÞ;
ð33Þ
where the grid spacing is hi = xi+1 xi. The FLOP of the first-order derivative is counted as
FLOPFDM ðci2 ; ci1 ; ci ; ciþ1 ; ciþ2 Þ ¼ ð13; 12; 30; 11; 14Þ ¼ 87: fx
ð34Þ
The spatial derivative of the FDM consists of four types of the grid spacing (hi2, hi1, hi, and hi+1), and requires more computational cost than the IDO-CF scheme. The Flop advantage of IDO-CF scheme leads to shorter computational time: 775 ms for IDO-CF scheme and 1660 ms for FDM per 10,000 steps on Intel Core i7 2.8 GHz machine.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5793
The second-order spatial derivatives of the IDO-CF scheme and the 4th-order FDM are evaluated in a similar way.
FLOPIDOPV ð73Þ þ FLOPIDOIA ð53Þ fxx fxx ¼ 63; 2
FLOPIDO fxx ¼
¼ 107: FLOPFDM fxx
ð35Þ ð36Þ
3. Fourier analysis Turbulent phenomena contain a lot of small-scale vortices, because the scale of the turbulent vortex spreads proportional to Re3/4. Since the error of the numerical scheme is mainly caused by such small vortices, it is necessary to use high-order discretization for turbulent simulations in Fourier space. Although the spectral method [29] has the highest accuracy, it is only applicable to periodic boundary conditions. In this section, the accuracy and stability in the one-dimensional IDO-CF scheme are analyzed in comparison with the FDM in Fourier space, and the overall characteristics do not change in twoand three-dimensional stability analyses. 3.1. Basic formulation The spatial profile of the variable f(x) defined over the domain [0, L] with an uniform grid spacing Dx = L/N is decomposed into Fourier series
f ðxÞ ¼
X
^f ðkÞeiWx=Dx ;
ð37Þ
k
where i ¼
fj ¼
pffiffiffiffiffiffiffi 1, and W = 2pkDx/L is a scaled wavenumber. The values at the points xj and xj±1 are also decomposed as
X
^f ðkÞeiWxj =Dx ;
f j1 ¼
k
X
^f ðkÞeiWxj =Dx eiW :
ð38Þ
k
Using Eq. (38), the value at the point xj+m is simply decomposed as
fjþm ¼ fj eiWm :
ð39Þ x
The integrated average value fj of the IDO-CF scheme is also decomposed as x
fj ¼
1 Dx
Z
iW
þDx=2
f ðxj þ xÞdx;
¼ fj
Dx=2
iW
e 2 e 2 : iW
ð40Þ
Since this equation expresses the relation between the point value and its integrated average value, the stability and accuracy of the IDO-CF scheme can be examined by using Eqs. (39) and (40). 3.2. First-order spatial derivative In the IDO-CF scheme having the point value and integrated average value, we define a grid spacing h = Dx/2 and a scaled wavenumber w = 2pkh/L, and we compare the IDO-CF scheme with the FDM with the same total number of independent variables. The first-order derivatives of the IDO-CF scheme and FDM reduce to the following approximations, respectively, IDO: 4th-order central (point value and corresponding integrated average value)
fx ðxj Þ ¼ x
ðfjþ1 fj1 Þ 4 x fjþ1=2 x fj1=2 ; 2 Dx
fx ðxjþ1=2 Þ ¼
fjþ1 fj : Dx
ð41Þ
ð42Þ
FDM: 4th-order central (point value and corresponding integrated average value)
fx ðxj Þ ¼ x
fx ðxj Þ ¼
fjþ2 8f jþ1 þ 8f j1 fj2 ; 12h
fjþ1=2 fj1=2 ; h
x
¼
fjþ2 8x fjþ1 þ 8x fj1 x fj2 : 12h
ð43Þ
ð44Þ
The first-order derivatives of the point value and integrated average value are expressed as the different formulation in the IDO-CF scheme and are described in the following matrix form,
F nx ¼ A F n ;
ð45Þ
5794
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
Fn ¼
fn ; x n f
F nx ¼
fxn ; x n fx
A¼
a11
a12
a21
a22
:
ð46Þ
The notation F = (f, xf)T is a variable vector which consists of the point value and its integrated average value. The vector Fx is its first-order spatial derivative. The superscript n indicates the time step. The elements a11 and a12 for the point value (PV) of the matrix A vary with the accuracy of the first-order spatial derivatives and are obtained from Eqs. (41), (45) and (46)
a11 ¼
i sinð2wÞ ; h 2
a12 ¼
i 2 sinðwÞ ðIDO-CF PVÞ: h
ð47Þ
The elements a21 and a22 for the integrated average value (IA) do not change their values because of their definitions,
! Z þDx=2 @ 1 f ðxÞdx ; @x Dx Dx=2 f Dx f D2x ¼ 2 ; Dx i sinðwÞ: ¼f h
ðx f Þx ¼
ð48Þ
We obtain the following a21 and a22 from Eq. (48),
a21 ¼
i sinðwÞ; h
a22 ¼ 0 ðIDO-CF IAÞ:
ð49Þ
This is the exactly same as the FVM formulation. Fig. 6 shows the resolution of the first-order spatial derivative in Fourier space. The accuracy of the IDO-CF scheme at the point (IDO-CF PV) is almost the same as that of the FDM with 4th-order accuracy. The first-order spatial derivative on the cell (IDO-CF IA), defined as Eq. (49), is the analytically exact solution. 3.3. Fourier analysis for linear advection equation In this subsection, we examine a linear advection equation. In order to discuss the time-integration method, we consider the following Taylor expansion around the time equal to tn,
f nþ1 ¼ f n þ
! 2 ! k n @f @ 2 f n dt @ k f n dt dt þ þ þ : @t 2! k! @t2 @t k
ð50Þ
In the advection equation, a time derivative is written by the following equation,
@f @f ¼ u ; @t @x
ð51Þ
where the velocity u is assumed to be constant in time and space. The time derivatives in Eq. (50) are replaced by using the relation Eq. (52)
! ! ! @ @ 2 ðu dtÞ2 @ k ðu dtÞk ðu dtÞ þ 1þ þ þ f n: @x @x2 2! @xk k!
¼
3 2.5
Imaginary
10
IDO-CF PV FDM (PV & IA) Exact (IDO-CF IA)
Error of spatial derivative
f
nþ1
2 1.5 1 0.5
0.5
1
1.5
2
Wavenumber
2.5
3
ð52Þ
0
10
-2
10
-4
10
-6
10
-8
10
-10
10
-12
10
IDO-CF PV FDM (PV & IA)
-3
10
-2
10
-1
10
0
Wavenumber
Fig. 6. Fourier analysis for the first-order spatial derivative: (left) imaginary part and (right) error; IDO-CF IA is not drawn in the right figure as there does not exist any error.
5795
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
The Runge–Kutta method is based on the Taylor expansion to obtain the accuracy in time. The resolutions for the advection equation are shown in Fig. 7. The accuracy of the point value (IDO-CF PV) is almost the same as that of the 4th-order FDM. The integrated average values (IDO-CF IA) are quite better than those of the FDM in highwavenumber area. 3.4. The stability analysis for the advection equation The eigenvalues of the numerical scheme show the stability of the time-integration. The IDO-CF scheme has various kinds of moments. The matrix eigenvalue problem is solved as
F nþ1 ¼ DT;k F n ; DT;k ¼ PKk P1 ;
Kk ¼
kþ
0
0
k
;
ð53Þ
where kk is an eigenvalue matrix and the components k+ and k are eigenvalues of the point value and integrated average value. Fig. 8 shows the contour plot of the eigenvalue k = max(jk+j, jkj) as functions of Courant–Friedrichs–Lewy (CFL) number and wavenumber, and the region where k 6 1 is stable. For the 2nd-order Runge–Kutta method, the stable area is too small (CFL < 0.05) to solve the linear advection equation. When we use the 3rd-order Runge–Kutta time integration, the IDO-CF scheme is stable for CFL < 1.15. The stable region of the CFL number for the 4th-order Runge–Kutta method is enlarged to CFL < 1.87. Therefore the IDO-CF scheme requires at least 3rd-order accuracy in time. In the following section of this paper, numerical simulations are carried out by using the 3rd-order Runge–Kutta method.
4. Governing equations and SGS models 4.1. Governing equations The governing equations for incompressible flows consist of the Navier–Stokes equations and the continuity equation,
@ui @ui uj @p 1 @ 2 ui ¼ þ ; þ @xi Res @xj @xj @t @xj
@ui ¼ 0: @xi
ð54Þ
Here ui, p and Res denote the velocity, pressure and the Reynolds number which is based on the friction velocity and the channel half-width. In LES, governing equations are respectively filtered in physical space, and filtered variables f denoted by an overbar are defined as
f x; D; t ¼
Z
0 f ðx0 ; tÞG x; x0 ; D dx ;
ð55Þ
where G is the filter function and D is the filter width. We adopt that the filter width is the same as the grid spacing. The grid filter function G is normalized to satisfy the following condition,
0 G x; x0 ; D dx ¼ 1:
ð56Þ
3
10
IDO-CF PV IDO-CF IA FDM (PV & IA) Exact
Error of phase
2.5 2
Phase
Z
1.5 1 0.5
0.5
1
1.5
2
Wavenumber
2.5
3
0
10
-2
10
-4
10
-6
10
-8
10
-10
10
-12
10
IDO-CF PV IDO-CF IA FDM (PV & IA)
-3
10
-2
10
-1
10
0
Wavenumber
Fig. 7. Fourier analysis for the advection equation with 4th-order Runge–Kutta method: (left) phase and (right) phase error.
5796
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
2
2
1.5
1.5
Courantnumber
1 1.
Courantnumber
Unstable
1
1.0
5
0.95
1
CFL=1.15
Stable
0.5
0.5 1 0
1 1.0 5 1.1
Unstable
Stable 0
CFL=0.05 0.5
1
0
1.5
0
0.5
1.5
Unstable
2
0.9 1 1.1
CFL=1.87
0. 0.9 95
1.5
Courantnumber
1
Wavenumber
Wavenumber
Stable 1
0.5
0
0
0.5
1
1.5
Wavenumber : (top-left) 2nd-, (top-right) 3rd-, (bottom) 4th-order Fig. 8. Contour plots of eigenvalues for time integration of the advection equation. CFL ¼ udt h Runge–Kutta method.
Applying the grid filter to Eq. (54), we obtain the filtered Navier–Stokes equations and the filtered continuous equation,
i u j i i @ u @ sij @u @p 1 @2u ¼ þ ; þ @xi @xj Res @xj @xj @t @xj
i @u ¼ 0: @xi
ð57Þ
j which is derived from the non-linear coni u The difference between Eqs. (54) and (57) is the SGS stress tensors sij ¼ ui uj u vection term. Since the filtered Navier–Stokes equations are not closed by themselves, the SGS stress tensor should be modeled. 4.2. Subgrid-scale models 4.2.1. Dynamic Smagorinsky model In the Smagorinsky eddy viscosity model, the SGS stress tensor sij is modeled with the filter width D as
sij ¼ 2C D2 jSjSij ¼ 2mSGS Sij ; Sij ¼
j @ u i 1 @u ; þ 2 @xi @xj
1=2 jSj ¼ 2Sij Sij ;
ð58Þ
where Sij and S are the velocity-strain tensor and its magnitude in the GS, respectively. In the Smagorinsky model (SM) [3], the model coefficient C is constant in an entire computational domain, and the SM does not satisfy a correct asymptotic behavior to a wall. The dynamic Smagorinsky model (DSM) [6,7] improves those defects by calculating the model parameter dynamically. The model parameter of the DSM is determined by using a least square procedure h,i with an average in homogeneous directions as
C DSM ¼
hLij M ij i ; hMij M ij i
ð59Þ
b and the grid scale D. They are given by where Lij and Mij are the resolved stress tensors between the test scale D
^ i u ^ j ; d j u Lij ¼ u iu
ð60Þ
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
c d b 2 b Mij ¼ 2D2 j SjSij 2 D S Sij :
5797
ð61Þ
b ¼ 2D, and the test-filtered value ^ f in homogeneous directions is calculated by using SimpThe test filter width is chosen as D son’s rule
1 ^ f i ¼ ðf i1 þ 4f i þ f iþ1 Þ: 6
ð62Þ
4.2.2. Coherent-structure model The model parameter of the DSM is determined dynamically by the width ratio of grid filter and test filter. This procedure is the most notable breakthrough in LES. However, the averaging in homogeneous directions is required for numerical stability. That makes it difficult to treat a complex geometry with the DSM. In the CSM [12], the model parameter C is dynamically calculated from a flow based on turbulence structures. No averaging in homogeneous directions is needed, so that this eddy viscosity model can treat a complex geometry. The model parameter of the CSM is defined as
C ¼ C CSM jF CS j3=2 ;
ð63Þ
where CCSM = 1/30 is a fixed model parameter and FCS is the coherent structure function defined as
F CS ¼
Q E
ð1 6 F CS 6 1Þ;
ð64Þ
Q¼
j i 1 1 @u @u ; W ij W ij Sij Sij ¼ 2 2 @xi @xj
ð65Þ
E¼
1 @ u j j 1 @u : W ij W ij þ Sij Sij ¼ 2 2 @xi @xi
ð66Þ
The coherent structure function FCS, composed of the second invariant Q and the magnitude of a velocity gradient tensor E, plays a role of a wall function, and its magnitude is less than unity. Therefore the model parameter can be locally determined without averaging in homogeneous directions. The CSM shows almost the same performance as the DSM and is applicable to complex geometries [13]. In the following section, the SGS stress tensor is calculated by using the CSM. 5. Turbulent channel flow simulations The turbulent channel flows have been studied for several decades as one of the most basic benchmarks of wall turbulence and are often used for validation of numerical schemes. The flows between two parallel walls are driven by a mean pressure gradient. The periodic boundary condition is considered in the streamwise and spanwise directions (x1 and x3). The computational grid is uniform in the periodic directions. The non-slip condition is imposed on the wall. The computational grids in this section are concentrated near the wall by using a hyperbolic-tangent function. Reynolds number is based on the channel half width d and friction velocity us (Res = usd/m). The turbulent statistics are obtained with time–space averaging after the turbulence is well developed to a statistically stable state. The IDO-CF scheme has better accuracy than conventional FDMs with 4th-order accuracy for the linear advection equation. The spectral like resolution has been confirmed by solving two-dimensional homogeneous isotropic turbulence [26]. We solve the incompressible Navier–Stokes equations to satisfy the velocity divergence free condition by using the SMAC method. The BiCGSTAB method is employed to solve the sparse matrix of the Poisson equation. The 3rd-order Runge–Kutta method is used for the time-integration with the time step of CFL = 0.1. In this section, the results of the IDO-CF scheme with the CSM are compared to those of the 4th-order FDM with the DSM and CSM as well as DNS data in Moser et al. [27] at Res = 180. We also carried out the cases of different Reynolds numbers Res = 395 and 590 with larger grid numbers and found almost the same results. 5.1. Effect of grid resolution for the IDO-CF scheme LES of the turbulent channel flows is conducted at Res = 180 in various cases of the grid resolution. In the IDO-CF scheme, the number of degrees of freedom is twice compared to conventional FDMs in every direction. Numerical tests are diligently accumulated with the same total number of degrees of freedom. The governing equations are discretized by the central formulation (Eqs. (2) and (3)). The influence of wall-normal grid resolution is examined with changing the grid resolution in the periodic directions. Figs. 9 and 10 show the maximum values of the mean velocity and the turbulent intensity in the streamwise direction. þ The horizontal axis indicates a grid spacing in periodic directions defined as xþ 1 ¼ 3x3 ¼ Res 4p=N 1 . The number of grid þ min point in the wall direction is defined as 16, 24, 32, and 64 which correspond to x2 ¼ 14:9, 5.8, 2.7, and 1.3.
5798
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
þ min The mean velocity is overestimated about 4% at xþ ¼ 2:7, i.e. (N1, N2, N3) = (8, 32, 8) and comes closer to 1 ¼ 282 and x2 the DNS from the higher value with increasing the number of grid point. However, the maximum value of the mean velocity þ min is still underestimated at xþ ¼ 1:3, i.e. (N1, N2, N3) = (32, 64, 32) in Fig. 9. It is well known that the mean veloc1 ¼ 70 and x2 ity of the 4th-order scheme without SGS models is below that of DNS because of their truncation errors and the lack of grid resolution. For that reason, the SGS stress should be applied to compensate for the lack of grid resolution. The maximum value of the turbulent intensity approaches the result of DNS in each case and is mainly related to the resolution of grid spacing in both two periodic directions. It is confirmed that the maximum values of the mean velocity and the turbulent intensity are getting closer to the DNS in much finer grid spacing xþ 1 ¼ 47. The effect of grid resolution on the turbulent statistics is similar to the FDM [30,31]. This is why a proper SGS model should be applied to LES using the IDO-CF scheme.
5.2. The effectiveness of SGS model In LES, the flow field in GS is directly calculated, and the effect of the SGS is evaluated with the SGS stress tensor
sij ¼ 2mSGS Sij with the CSM. The SGS eddy viscosity becomes zero near the wall without a wall-damping function. The averaging in homogeneous direction is not conducted. In this subsection, the IDO-CF scheme with the CSM is compared to the IDO-CF scheme without the SGS model and the IDO-CF upwind scheme (Eq. (9)). The grid resolution is shown in Table 1. Fig. 11 shows the mean velocity profiles at (N1, N2, N3) = (16, 32, 16) (case 1). In the previous subsection, the IDO-CF scheme without the SGS model underestimates the maximum value of the mean velocity. The SGS stress tensor modifies the underestimation, and the IDO-CF scheme with the CSM is in good agreement with DNS. The IDO-CF upwind scheme implicitly introduces numerical viscosity (Eq. (13)) too much, so that it overestimates the profile of the mean velocity. Fig. 12 shows the turbulent intensity in the streamwise direction. The IDO-CF schemes with and without the CSM obtain almost the same profiles and have good convergences to the result of DNS. The IDO-CF upwind scheme overestimates the turbulent intensity as well because of the numerical viscosity.
20
N1=8
N1=12
N1=48 N =32 1
N1=16
U 1+
15
dx +2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS + dx1 = 3 dx 3+
10
5
(N2=16) (N2=24) (N2=32) (N2=64)
0 0
50
100
150
200
250
300
dx1+ Fig. 9. The maximum values of the mean velocity at various grid resolution without the SGS model.
N1=8 4
N1=12 N1=16
u’ 1
3
N1=48 N =32 1 +
2
dx 2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS + dx1 = 3 dx 3+
1
(N2=16) (N2=24) (N2=32) (N2=64)
0 0
50
100
150
200
250
300
dx1+ Fig. 10. The maximum values of the turbulent intensity at various grid resolution without the SGS model.
5799
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805 Table 1 The numerical condition of turbulent channel flows.
Case Case Case Case
1 2 3 4
þ
þ
þ
Res
L1 L2 L3
N1 N2 N3
dx1
dx2 min
dx3
180 180 395 590
4p 2 4/3p 4p 2 4/3p 2p 2 p 2p 2 p
16 32 16 32 64 32 32 64 32 32 64 32
141 70 77 116
2.7 1.3 2.9 4.3
47 23 39 58
25
IDO with CSM w/o CSM upwind DNS
20
U+
15
10
5
0 0 10
10
1
10
x2 +
2
Fig. 11. The profiles of the mean velocity using the CSM, upwind scheme and without CSM (case 1).
IDO with CSM w/o CSM upwind DNS
4
u’ 1
3
2
1
0 0
50
100
x 2+
150
Fig. 12. The profiles of the turbulent intensity using the CSM, upwind scheme and without CSM (case 1).
Figs. 13 and 14 show the SGS eddy viscosity and the shear stress profile of the IDO-CF scheme with the CSM. In the turbulent channel flows, the pressure gradient is balanced against the sum of total shear stress and viscosity stress. The total shear stress is composed of Reynolds stress u01 u02 and SGS stress s12 in LES. In Fig. 14, the SGS model gives the additional shear stress to the flow field and total shear stress is in good agreement with DNS. Figs. 15 and 16 show the maximum value of the mean velocity and the turbulent intensity at Res = 180. The spatial resolution is the same as Figs. 9 and 10. The profile of the mean velocity is improved well to compare with the profile of the IDOCF scheme without the SGS model. The profile of the turbulent intensity converges to the value of DNS. 5.3. The comparison with finite difference scheme Numerical results of the IDO-CF scheme are compared to the 4th-order FDM [15]. Note that the 2nd-order FDM is used in the wall direction because of facilitating the treatment of the boundary condition on the wall. The 4th-order FDM is often used to calculate turbulent flows because of their conservation of kinetic energy. The DSM and CSM are respectively applied to the FDM. The IDO-CF scheme has 4th-order accuracy in every direction and simulations are carried out at the same number of grid point (case 1). For example, the number of grid (N1, N2, N3) = (16, 32, 16) with the IDO-CF scheme corresponds to that of grid (N1, N2, N3) = (32, 64, 32) with the FDM.
5800
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
0.005
CSM
0.003
ν
SGS
0.004
0.002
0.001
0
0
50
100
x2
150
+
Fig. 13. The profile of the SGS eddy viscosity with the CSM (case 1).
1
u’ 1 u’ 2 τ 12 u’ 1 u’ 2 + τ 12 DNS u’ 1 u’ 2
2
0.6
u’ 1 u’
0.8
0.4
0.2
0 0
50
100
150
x 2+
Fig. 14. The profiles of the Reynolds shear stress with the CSM (case 1).
N1=48
20
N1=32
N1=8
N1=12
N1=16
U 1+
15
dx +2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS + dx1 = 3 dx 3+
10
5
0
0
50
100
150
dx1
200
(N2=16) (N2=24) (N2=32) (N2=64)
250
300
+
Fig. 15. The maximum values of the mean velocity at various grid resolution with the SGS model.
Fig. 17 shows the mean velocity profiles of the IDO-CF scheme with the CSM and the FDM with the DSM and CSM. All profiles are in good agreement with the profile of DNS. Figs. 18–20 show the profiles of the turbulent intensity in the streamwise, the wall, and the spanwise directions. The turbulent intensity in the streamwise direction strongly depends on the accuracy of a numerical scheme and the grid resolution. The streamwise turbulent intensity of the IDO-CF scheme with
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5801
N1=8
4
N1=12 N1=16
u’ 1
3
N1=48 N1=32 +
2
dx 2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS dx1+ = 3 dx 3+
1
0
0
50
100
150
dx1
200
(N2=16) (N2=24) (N2=32) (N2=64)
250
300
+
Fig. 16. The maximum values of the turbulent intensity at various grid resolution with the SGS model.
20
U
+
15
IDO with CSM FDM with CSM FDM with DSM DNS
10
5
0 0 10
10
1
10
x2 +
2
Fig. 17. The profiles of the mean velocity using the IDO-CF and FDM.
IDO with CSM FDM with CSM FDM with DSM DNS
3
u’ 1
2
1
0 0
50
100
x 2+
150
Fig. 18. The profiles of the turbulent intensity in the streamwise direction using the IDO-CF and FDM.
the CSM is lower than that of FDM with the CSM and DSM. The present scheme achieves better agreement with the DNS. In Figs. 19 and 20, it is found that the wall and spanwise turbulent intensities of the IDO-CF scheme are higher than those of FDM. The turbulence field of the IDO-CF scheme is more isotropic than that of FDM. The IDO-CF scheme gives good results for the incompressible Navier–Stokes equations with less number of grid points.
5802
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
IDO with CSM FDM with CSM FDM with DSM DNS
u’ 2
1
0.5
0 0
50
100
150
x 2+
Fig. 19. The profiles of the turbulent intensity in the wall direction using the IDO-CF and FDM.
1.5
IDO with CSM FDM with CSM FDM with DSM DNS
u’ 3
1
0.5
0 0
50
100
150
x 2+
Fig. 20. The profiles of the turbulent intensity in the spanwise direction using the IDO-CF and FDM.
25
U+
20
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
15
10
5
0 0 10
10
1
x2 +
10
2
Fig. 21. The profiles of the mean velocity at Res = 180, 395, and 590. (case 2–4).
5.4. The flow at various Reynolds number (Res = 180, 395, and 590) The IDO-CF scheme obtains better results in comparison with conventional FDMs of the same order accuracy. This is because the present scheme has higher spectral resolution in Fourier space. Further LES is conducted at higher Reynolds number as shown in Table 1. Fig. 21 shows the mean velocity profiles at (N1, N2, N3) = (32, 64, 32) (case 2–4). The profiles at each Reynolds number are in good agreement with the DNS. Figs. 22–24 show the turbulent intensity profiles in streamwise, wall, and spanwise directions at (N1, N2, N3) = (32, 64, 32) (case 2–4). The profiles of the mean velocity and the turbulent intensities at Res = 180 are in
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
3
5803
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
u’ 1
2
1
0 0
100
200
300
x 2+
400
500
Fig. 22. The profiles of the turbulent intensity in the streamwise direction at Res = 180, 395, and 590. (case 2–4).
1.5
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
u’ 2
1
0.5
0 0
100
200
300
x2
+
400
500
Fig. 23. The profiles of the turbulent intensity in the wall direction at Res = 180, 395, and 590. (case 2–4).
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
u’ 3
1.5
1
0.5
0 0
100
200
300
x 2+
400
500
Fig. 24. The profiles of the turbulent intensity in the spanwise direction at Res = 180, 395, and 590. (case 2–4).
good agreement with the DNS. At Res = 395 and 590, the profiles of the mean velocity and the turbulent intensities in streamwise, wall, and spanwise directions show good convergence to the DNS in spite of less number of grid points, and the maximum values of the turbulent intensities in the streamwise direction are in good agreement with the DNS in Fig. 22. The IDO-CF scheme with the CSM can be applied to various Reynolds number with a small grid number, and is one of good numerical schemes for turbulent phenomena.
5804
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
6. Summary and conclusions The IDO-CF scheme with the subgrid-scale model is applied to turbulent phenomena, possessing the following three advantages; (a) the implementation of the wall boundary condition with high-order accuracy is as easy as that with the 2nd-order FDM, (b) the computational cost for non-uniform grid spacing is much smaller than that of the 4th-order FDM, and (c) the IDO-CF scheme has better spectral resolution than the same order FDM. The first-order spatial derivative of the integrated average value is exactly the same as the FVM. Fourier analysis shows that the spatial derivative of the IDO-CF scheme has better spectral resolution in comparison with the 4th-order FDM. The IDO-CF scheme is suitable for the turbulent simulations because LES is sensitive to the accuracy and resolution in Fourier space. In order to evaluate the IDO-CF scheme with the CSM, the fully developed turbulent channel flows are examined at various grid resolutions. Numerical simulations at three Reynolds numbers Res = 180, 395, and 590 are compared to the reference of DNS data in Moser et al. [27]. We have three major conclusions. First, it is found that the IDO-CF scheme without the CSM underestimates the profile of the mean velocity in relatively coarse grid resolution, so that the SGS model is needed. The IDO-CF upwind scheme overestimates the profiles of the mean velocity and the turbulent intensity with the excessive numerical viscosity introduced implicitly. The CSM evaluates the SGS eddy viscosity correctly, and numerical results of the IDO-CF scheme with the CSM are in good agreement with the DNS. Second, the IDO-CF scheme is compared to the FDM with the DSM and the CSM. The IDO-CF scheme with the CSM achieves a better profile of the turbulent intensity in the streamwise direction. In LES of turbulent channel flows, the statistical value of the turbulent intensity mainly depends on the grid resolution and discretized accuracy in both two periodic directions. The IDO-CF scheme has higher spectral resolution for the incompressible Navier–Stokes equations than conventional FDMs with the same order of accuracy. Finally, we validate the dependence of the IDO-CF scheme with the CSM on various Reynolds numbers. The statistical values at Res = 180 are in good agreement with the DNS at (N1, N2, N3) = (32, 64, 32). At Res = 395 and 590, the mean velocity and the turbulent intensity show good convergence to the DNS in spite of less number of grid points. It is concluded that the present scheme is one of most promising candidates to simulate turbulence. In the future work, we apply the present scheme to turbulent phenomena with the complex geometry. Acknowledgements This research was supported in part by KAKENHI, Grant-in-Aid for Scientific Research (B) 19360043 from The Ministry of Education, Culture, Sports, Science and Technology (MEXT), JST-CREST project ‘‘ULP-HPC: Ultra Low-Power, High Performance Computing via Modeling and Optimization of Next Generation HPC Technologies’’ from Japan Science and Technology Agency (JST) and JSPS Global COE program ‘‘Computationism as a Foundation for the Sciences’’ from Japan Society for the Promotion of Science (JSPS). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
M. Lesieur, O. Metais, New trends in large-eddy simulations of turbulence, Annu. Rev. Fluid Mech. 28 (1) (1996) 45–82. C. Meneveau, J. Katz, Scale-invariance and turbulence models for large-eddy simulation, Annu. Rev. Fluid Mech. 32 (1) (2000) 1–32. J. Smagorinsky, General circulation model of the atmosphere, Mon. Weather Rev. 164 (1963) 91–99. E.R. Van Driest, On turbulent flow near a wall, J. Aeronaut. Sci 23 (1956) 1007. C. Haertel, L. Kleiser, Analysis and modelling of subgrid-scale motions in near-wall turbulence, J. Fluid Mech. 356 (1998) 327–352. M. Germano, U. Piomelli, P. Moin, W. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A: Fluid Dyn. 3 (1991) 1760. D. Lilly, A proposed modification of the germano subgrid-scale closure method, Phys. Fluids A: Fluid Dyn. 4 (1992) 633. Y. Zang, R. Street, J. Koseff, A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows, Phys. Fluids A: Fluid Dyn. 5 (1993) 3186. B. Vreman, B. Geurts, H. Kuerten, On the formulation of the dynamic mixed subgrid-scale model, Phys. Fluids 6 (12) (1994). C. Meneveau, T.S. Lund, W.H. Cabot, A Lagrangian dynamic subgrid-scale model of turbulence, J. Fluid Mech. 319 (353) (1996). E. Bou-Zeid, C. Meneveau, M. Parlange, A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows, Phys. Fluids 17 (2005) 025105. H. Kobayashi, The subgrid-scale models based on coherent structures for rotating homogeneous turbulence and turbulent channel flow, Phys. Fluids 17 (2005) 045104. H. Kobayashi, F. Ham, X. Wu, Application of a local sgs model based on coherent structures to complex geometries, Inter. J. Heat Fluid Flow 29 (3) (2008). S. Ghosal, An analysis of numerical errors in large-eddy simulations of turbulence, J. Comput. Phys. 125 (1) (1996) 187–206. Y. Morinishi, T. Lund, O. Vasilyev, P. Moin, Fully conservative higher order finite difference schemes for incompressible flow, J. Comput. Phys. 143 (1) (1998) 90–124. Y. Liu, M. Vinokur, Z. Wang, Spectral difference method for unstructured grids. I: basic formulation, J. Comput. Phys. 216 (2) (2006) 780–801. T. Yabe, T. Aoki, A universal solver for hyperbolic equations by cubic-polynomial interpolation. I: One-dimensional solver, Comput. Phys. Commun. 66 (2–3) (1991) 219–232. T. Yabe, T. Ishikawa, P. Wang, T. Aoki, Y. Kadota, F. Ikeda, A universal solver for hyperbolic equations by cubic-polynomial interpolation. II: Two-and threedimensional solvers, Comput. Phys. Commun. 66 (1991) 233–242. T. Aoki, Interpolated differential operator (IDO) scheme for solving partial differential equations, Comput. Phys. Commun. 102 (1–3) (1997) 132–146. K. Sakurai, T. Aoki, W. Lee, K. Kato, Poisson equation solver with fourth-order accuracy by using interpolated differential operator scheme, Comput. Math. Appl. 43 (6–7) (2002) 621–630. H. Yoshida, T. Aoki, T. Utsumi, Improvement of accuracy and stability in numerically solving hyperbolic equations by IDO (interpolated differential operator) scheme with Runge–Kutta time integration, Electron. Commun. Jpn. (Part III: Fundam. Electron. Sci.) 87 (2) (2004). T. Utsumi, T. Kunugi, T. Aoki, Stability and accuracy of the cubic interpolated propagation scheme, Comput. Phys. Commun. 101 (1–2) (1997) 9–20.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5805
[23] Y. Imai, T. Aoki, Accuracy study of the IDO scheme by Fourier analysis, J. Comput. Phys. 217 (2) (2006) 453–472. [24] T. Yabe, R. Tanaka, T. Nakamura, F. Xiao, An exactly conservative semi-Lagrangian scheme (CIP–CSL) in one dimension, Mon. Weather Rev. 129 (2) (2001) 332–344. [25] F. Xiao, T. Yabe, Completely conservative and oscillationless semi-Lagrangian schemes for advection transportation, J. Comput. Phys. 170 (2) (2001) 498–522. [26] Y. Imai, T. Aoki, K. Takizawa, Conservative form of interpolated differential operator scheme for compressible and incompressible fluid dynamics, J. Comput. Phys. 227 (4) (2008) 2263–2285. [27] R. Moser, J. Kim, N. Mansour, Direct numerical simulation of turbulent channel flow up to re = 590, Phys. Fluids 11 (1999) 943. [28] J. Kim, P. Moin, R. Moser, Turbulence statistics in fully developed channel flow at low reynolds number, J. Fluid Mech. 177 (1) (1987) 133–166. [29] D. Gottlieb, A. Orszag, Numerical Analysis of Spectral Methods, SIAM, Philadelphia, 1977. [30] Y. Morinishi, O. Vasilyev, A recommended modification to the dynamic two-parameter mixed subgrid scale model for large eddy simulation of wall bounded turbulent flow, Phys. Fluids 13 (2001) 3400. [31] Y. Morinishi, O. Vasilyev, Vector level identity for dynamic subgrid scale modeling in large eddy simulation, Phys. Fluids 14 (2002) 3616.
Contents lists available at ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Large-eddy simulation of turbulent channel flows with conservative IDO scheme Naoyuki Onodera a,⇑, Takayuki Aoki b, Hiromichi Kobayashi c,d a
Department of Energy Sciences, Tokyo Institute of Technology, 2-12-1 O-okayama, Meguro-ku, Tokyo 152-8550, Japan Global Scientific Information and Computing Center, Tokyo Institute of Technology, 2-12-1 O-okayama, Meguro-ku, Tokyo 152-8550, Japan c Department of Physics, Keio University, 4-1-1 Hiyoshi, Kouhoku-ku, Yokohama 223-8521, Japan d Research and Education Center for Natural Sciences, Keio University, 4-1-1 Hiyoshi, Kouhoku-ku, Yokohama 223-8521, Japan b
a r t i c l e
i n f o
Article history: Received 3 August 2010 Received in revised form 2 April 2011 Accepted 3 April 2011 Available online 12 April 2011 Keywords: Multi-moment scheme IDO scheme Large-eddy simulation Coherent structure model
a b s t r a c t The resolution of a numerical scheme in both physical and Fourier spaces is one of the most important requirements to calculate turbulent flows. A conservative form of the interpolated differential operator (IDO-CF) scheme is a multi-moment Eulerian scheme in which point values and integrated average values are separately defined in one cell. Since the IDO-CF scheme using high-order interpolation functions is constructed with compact stencils, the boundary conditions are able to be treated as easy as the 2nd-order finite difference method (FDM). It is unique that the first-order spatial derivative of the point value is derived from the interpolation function with 4th-order accuracy and the volume averaged value is based on the exact finite volume formulation, so that the IDO-CF scheme has higher spectral resolution than conventional FDMs with 4th-order accuracy. The computational cost to calculate the first-order spatial derivative with non-uniform grid spacing is one-third of the 4th-order FDM. For a large-eddy simulation (LES), we use the coherent structure model (CSM) in which the model coefficient is locally obtained from a turbulent structure extracted from a second invariant of the velocity gradient tensor, and the model coefficient correctly satisfies asymptotic behaviors to walls. The results of the IDO-CF scheme with the CSM for turbulent channel flows are compared to the FDM with the CSM and dynamic Smagorinsky model as well as the direct numerical simulation (DNS) by Moser et al. Adding the sub-grid scale stress tensor of LES to the IDOCF scheme improves the profile of the mean velocity in comparison with an implicit eddy viscosity of the IDO-CF upwind scheme. The IDO-CF scheme with the CSM gives better turbulent intensities than conventional FDMs with the same number of grid points. The turbulent statistics calculated by IDO-CF scheme are in good agreement with the DNS at the various values of Reynolds number Res = 180,395, and 590. It is found that the IDO-CF scheme is suitable for the turbulent flow computation and improves the turbulent statistics with compact stencils. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Large-eddy simulation (LES) is one of successful approaches for unsteady turbulent flows [1,2]. LES resolves the dynamics of large-scale structures on the grid scale (GS), and the effect of smaller-scale turbulent structures are taken into account by using the subgrid-scale (SGS) model. A number of SGS models have been proposed as described below. In order to apply the ⇑ Corresponding author. Tel.: +81 3 5734 3575. E-mail address: [email protected] (N. Onodera). 0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.04.004
5788
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
SGS model to complex flows, however, it is preferable that the SGS model can be applied to wall bounded flows with local determination of its model coefficient. The Smagorinsky model (SM) [3] is a prominent SGS model based on a concept of eddy viscosity. The SGS stress tensor Þ2 S and the resolved strain rate tensor Sij , where CS is sij ¼ 2mSGS Sij is explicitly modeled by the eddy viscosity mSGS ¼ ðC S M ffiffiffiffiffiffiffiffiffiffiffiffi q is the filter width, and S ¼ 2Sij Sij . In the SM, the undetermined parameter CS is kept constant the Smagorinsky constant, M in a computational domain, although the model parameter should be changed among various turbulent flows. The SM gives a non-zero eddy viscosity in laminar-flow regions, and requires a wall function to damp its SGS viscosity in wall-bounded flows [4,5]. The dynamic Smagorinsky model (DSM) [6,7] is proposed to overcome these deficiencies by locally calculating the model coefficient. In addition to the grid filter, the test filter is applied to the incompressible Navier–Stokes equations. The model parameter is determined only by the ratio of the grid filter width and the test filter width. This model can be applied to various kinds of flows without empirical parameters. The DSM is the most notable breakthrough in LES, and many SGS models based on the DSM are proposed [8,9]. As one of the weak points, the model parameter introduces negative SGS eddy viscosity, which causes numerical instability. Since the clipping and averaging of the model parameter are required to avoid breakdown, the DSM is difficult to apply to complex geometries. The Lagrangian dynamic models [10,11] calculate the model parameter without averaging in homogeneous directions. The model parameter is time-averaged along fluid path lines. These models are applicable to complex geometries, however it takes more computational time in comparison with the DSM. The coherent structure model (CSM) [12] is one of remedies for these problems. A turbulent structure, which is extracted from a second invariant of the velocity gradient tensor, determines a model parameter locally. As a result, the SGS eddy viscosity becomes zero near the wall without a wall-damping function. In addition, the model parameter is guaranteed to be positive, and its variance is small, so that averaging is not necessary. The CSM shows good performance in complex geometries [13]. In this paper, we apply the CSM to computations of turbulent phenomena. 4th- or higher-order numerical schemes are necessary for turbulent simulations because low-order discretization error exceeds the amount of the SGS stress [14]. The Morinishi scheme [15], one of high-order finite difference methods (FDMs), is often applied to LES. The 4th-order one uses 7-point stencil in staggered grid, which makes it difficult to handle the wall boundary condition. On the other hand, multi-moment schemes achieve high-order accuracy with compact stencils. The multi-moment schemes define different moments in addition to a point value for a dependent variable like spectral difference methods [16]. A high-order interpolation function is constructed locally by using several neighbor moments. The cubic-interpolated pseudo-particle (CIP) scheme [17,18] is the originally proposed multi-moment scheme. First-order spatial derivative values and point values are defined at the same grid point, and the advection equation is solved with high accuracy using a semiLagrangian formulation. The interpolated differential operator (IDO) scheme [19,20] is a multi-moment Eulerian scheme. High-order accuracy is obtained not only in the advection equation but also in more general equations, such as a diffusion, pressure and Poisson equation [21]. In addition, it is confirmed that those multi-moment schemes have better spectral resolution than the conventional FDM with the same order accuracy [22,23]. The conservation property is necessary to carry out numerical experiments for stable long-time integration. A conservative semi-Lagrangian CIP (CIP-CSL) scheme [24,25] and a conservative formulation of IDO (IDO-CF) scheme [26] guarantee exactly mass conservation. It is confirmed that the truncation error of the IDO-CF scheme is smaller than that of conventional FDMs [26]. In this paper, the IDO-CF scheme with the CSM is proposed and studied for turbulent channel flows which has been studied for several decades as one of the most basic benchmarks [27,28]. The IDO-CF scheme can handle the wall boundary condition with high-order accuracy. The IDO-CF scheme has twice the number of degrees of freedom compared to conventional FDMs in every direction. Therefore, numerical tests are examined with the same total number of degrees of freedom. The turbulent statistics of the IDO-CF scheme with the CSM are assessed and compared to those of the FDM with the DSM and the CSM as well as the direct numerical simulation (DNS). This paper is organized as follows. The IDO-CF scheme is reviewed, and easiness to handle the boundary condition and the computational cost with non-uniform grid spacing are discussed in Section 2. The accuracy and stability are studied in comparison with the FDM in Section 3. The governing equation and SGS models are presented in Section 4. Numerical results for the turbulent channel flows are given in Section 5. The paper ends with some concluding remarks in Section 6.
2. Conservative form of IDO scheme The turbulent flow calculation requires high-order discretization, and the given partial differential equations (PDEs) have been generally discretized by the FDM, FVM, or finite element method (FEM). Multi-moment schemes formulated as the different approach realize high-order accuracy with compact stencils. Unlike conventional FDMs, they define several kinds of moments in one cell for each dependent variable. In this section, the IDO-CF [26], one of multi-moment schemes, is compared to the FDM with the same order of accuracy. The three major differences are the configuration of variable moments, the computational cost with non-uniform grid spacing, and the spectral resolution in Fourier space.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5789
2.1. Formulation and spatial discretization In the section, we give discussions of the IDO-CF scheme in one-dimensional case for simplicity and they are available for two- and three-dimensional cases. The IDO-CF scheme defines the value fi = f(xi) and its integrated average value R x fiþ1=2 ¼ D1x f ðxÞdx. Fig. 1 illustrates the locations of the variables in a one-dimensional case. Both the variables are time-integrated simultaneously by solving different discretized equations. For solving PDEs, spatial derivatives fx(x), fxx(x), . . . , fx(n)(x) are obtained by differentiating of an interpolation function. One of the interpolation functions is the following fourth-order polynomial for the local domain xi1 6 x 6 xi+1.
F c ðxÞ ¼ ac4 X 4 þ ac3 X 3 þ ac2 X 2 þ ac1 X þ ac0 ;
fxc
ðkÞ
¼ k!ack ;
ð1Þ
where X = x xi and the superscript ‘‘c’’ represents a central scheme. The interpolation function is called the fourth-order central interpolation function. For simplicity, we demonstrate the case of the uniform grid spacing Dx = xx+1 xi. The coefficients of the interpolation function are determined with the five matching conditions of Fc(xi) = fi, Fc(xi±1) = fi±1, and R Dx c 1 F ðxÞdx ¼ x fi1=2 . The spatial derivatives at x = xi are represented as follows Dx 0
ðfiþ1 fi1 Þ 4 x fiþ1=2 x fi1=2 ¼ ; 2 Dx 3ðfiþ1 þ 8f i þ fi1 Þ 15 x fiþ1=2 þ x fi1=2 c fxx;i ¼ ; 2 Dx2 6ðfiþ1 fi1 Þ 12 x fiþ1=2 x fi1=2 c ð3Þ fx;i ¼ ; Dx3 30ðfiþ1 þ 4f i þ fi1 Þ 90 x fiþ1=2 þ x fi1=2 c ð4Þ fx;i ¼ : Dx4
c fx;i
ð2Þ
ð3Þ
ð4Þ
ð5Þ
The first- and second-order spatial derivatives have 4th-order accuracy (O(Dx4)). c fx;i ¼ fx;i
1 1 ð5Þ 4 f Dx ; 3 5! x;i
c ¼ fxx;i fxx;i
ð6Þ
6 1 ð6Þ 4 f Dx : 7 6! x;i
ð7Þ
The major error of the first-order derivative comes from dispersion errors (fx(5)Dx4), and does not contain numerical dissipation. In solving the incompressible Navier–Stokes equations, a convection term, stress tensor, and the Poisson equation are discretized by those expressions (Eqs. (2) and (3)). Another interpolation function is the upwind interpolation function (IDO-CF upwind) constructed for the local domain xi1 6 x 6 xi (ui > 0) up up 2 F up ðxÞ ¼ aup 2 X þ a1 X þ a0
ð8Þ up
up
1 Dx
R0
up
x
with the three matching conditions of F (xi) = fi, F (xi1) = fi1, and F ðxÞdx ¼ fi1=2 . The superscript ‘‘up’’ represents Dx a upwind scheme. The spatial derivatives are also represented by the following approximations, up fx;i ¼
ð4f i þ 2f i1 6x fi1=2 Þ ; Dx
up ¼ fxx;i
ð6f i þ 6f i1 12x fi1=2 Þ : Dx2
ð9Þ
Fig. 1. Fourth-order interpolation function of the IDO-CF scheme.
ð10Þ
5790
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
The discretization errors are found to be up fx;i ¼ fx;i
1 1 ð3Þ 2 4 1 ð4Þ 3 f Dx þ f Dx ; 2 3! x;i 5 4! x;i
up fxx;i ¼ fxx;i 3
ð11Þ
1 ð3Þ 18 1 ð4Þ 2 f Dx þ f Dx : 3! x;i 5 4! x;i
ð12Þ
The first-order spatial derivative, derived from the upwind interpolation function, has 2nd-order accuracy (O(Dx2)) and is rewritten by the combination of the spatial derivative of the central interpolations (Eqs. (2), (4), and (5)). up c fx;i ¼ fx;i
1 1 c ð3Þ 2 4 1 c ð4Þ 3 Dx þ Dx f f 2 3! x;i 5 4! x;i
ðu > 0Þ:
ð13Þ
The last term of the fourth-order spatial derivative works as a numerical viscosity. The first-order spatial derivative of the cell-integrated average value is derived as follows. It consists of the cell surface values. x
fx;iþ1=2
@ 1 ¼ @x Dx
Z
!
xiþ1
FðxÞdx
¼
xi
fiþ1 fi : Dx
ð14Þ
This is exactly the same as the FVM formulation. The second-order spatial derivative is also the similar formulation, x
@2 1 ¼ 2 @x Dx
fxx;iþ1=2
Z
xiþ1
! FðxÞdx
xi
¼
fx;iþ1 fx;i ; Dx
ð15Þ
where we apply Eq. (2) to the first-order spatial derivatives. The conservative property of the cell-integrated average is completely guaranteed. For the FDM, the variables are defined at the grid points illustrated in Fig. 3. The first-order spatial derivative with 4thorder accuracy is represented as c fx;i ¼
fiþ2 8f iþ1 þ 8f i1 fi2 12h
ðFDMÞ;
ð16Þ
where h = xi+1 xi. The grid space h represents the distance between neighbor values. In the two-dimensional case, four kinds of variable moments are defined such as fi;j ¼ f ðxi ; yj Þ; x fiþ1=2;j ¼ R R RR 1 f ðx; yÞdx; y fi;jþ1=2 ¼ D1y f ðx; yÞdy, and xy fiþ1=2;jþ1=2 ¼ Dx1Dy f ðx; yÞdx dy in Fig. 5. Each variable is time-integrated simultaDx neously. The eight kinds of variable moments are defined in the three-dimensional case. The formulation for system equations in multi-dimensional cases was presented in the paper [26]. 2.2. Boundary treatment For the IDO-CF scheme, the boundary is imposed on xwall = x1 as shown in Fig. 2. The values f0, xf1/2, and f1 are given by the boundary condition. The first-order spatial derivatives xfx,3/2 and fx,2 are represented as follows, x
fx;3 ¼ 2
f2 f1 ; Dx
c fx;2 ¼
ð17Þ
ðf3 f1 Þ 4 x f5=2 x f3=2 : 2 Dx
ð18Þ
Fig. 2. Boundary treatment of the IDO-CF scheme.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5791
Eqs. (17) and (18) do not refer the outside values f0 and xf1/2. The boundary condition of the FDM is imposed on xwall = x2 as shown in Fig. 4. The corresponding first-order spatial derivative with 4th-order accuracy at x3 is represented as c fx;3 ¼
f5 8f 4 þ 8f 2 f1 12h
ðFDMÞ:
ð19Þ
The value f1 is located outside the domain. Eq. (19) is subjected to the boundary conditions. For that reason, low-order (2ndorder) or less stable one-side discretization is often used at the wall boundaries in FDMs. 2.3. Spatial discretization with non-uniform grid spacing Turbulent shear stress augments near the wall, and the spatial- and time-scales become much smaller than those away from the wall. In order to resolve the near-wall phenomena, the non-uniform grid spacing is commonly used in association with complicated derivation and computational cost. The procedure to derive the IDO-CF scheme with non-uniform grid spacing is the same as that with uniform grid spacing. Rx The point value is defined as fi = F(xi), and the volume-integrated average value is also defined as x fiþ1=2 ¼ x 1x xiiþ1 FðxÞdx. iþ1 i The first-order product of coefficients vector and variables vector as spatial derivative atTx =xi is expressed by theinner T fx;i ¼ ~ ci ~ f i ¼ ci1 ; ci1=2 ; ci ; ciþ1=2 ; ciþ1 fi1 ; x fi1=2 ; fi ; x fiþ1=2 ; fiþ1 . The coefficients of the central interpolation function are written as
ci1 ¼
2Dx2i
Dxi1 ðDxi þ Dxi1 Þ2
ci1=2 ¼
ci ¼
ð20Þ
;
Dx2i ð6Dxi þ 10Dxi1 Þ Dxi1 ðDxi þ Dxi1 Þ3
ð21Þ
;
4ðDxi Dxi1 Þ ; Dxi1 Dxi
ciþ1=2 ¼
ð22Þ
Dx2i1 ð10Dxi þ 6Dxi1 Þ
ciþ1 ¼
Dxi ðDxi þ Dxi1 Þ3 2Dx2i1
Dxi ðDxi þ Dxi1 Þ2
ð23Þ
;
ð24Þ
;
where the grid spacing is denoted as Dxi = xi+1 xi. The first-order spatial derivative is found to be composed of two types of grid spacing (Dxi1 and Dxi). The computational costs are evaluated as the number of floating-point operations (FLOP) for the first-order derivative at the point as
FLOPIDOPV ci1 ; ci1=2 ; ci ; ciþ1=2 ; ciþ1 ¼ ð7; 13; 4; 12; 8Þ ¼ 51: fx
ð25Þ
For the integrated average value of the IDO-CF scheme, Eq. (14) does not change. The FLOP of the first-order derivative on the cell is shown as
FLOPIDOIA ¼ 2: fx
ð26Þ
Each variable is time-integrated simultaneously. The total FLOP with the number of cell N is counted as
Total FLOPIDO ¼ FLOPIDOPV N þ FLOPIDOIA N: fx fx fx
Fig. 3. Fourth-order interpolation function of the FDM.
ð27Þ
5792
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
Fig. 4. Boundary treatment of the FDM.
Fig. 5. Definition of moments in two dimension.
Since Eq. (27) contains 2N independent variables, the FLOP per variable moment of the IDO-CF scheme is simply estimated as
FLOPIDO ¼ fx
Total FLOPIDO FLOPIDOPV þ FLOPIDOIA fx fx fx ¼ ¼ 26:5: 2N 2
ð28Þ
In the case of the spatial discretization of the FDM, all the variables are defined at grid points as fi = F(xi). The variables T vector and the coefficient vector are similarly expressed as ~ ci ¼ ðci2 ; ci1 ; ci ; ciþ1 ; ciþ2 ÞT and ~ f i ¼ ðfi2 ; fi1 ; fi ; fiþ1 ; fiþ2 Þ . The coefficients of the central interpolation function are written as
ci2 ¼
hi1 hi ðhiþ1 þ hi Þ ; hi2 ðhi1 þ hi2 Þðhi þ hi1 þ hi2 Þðhiþ1 þ hi þ hi1 þ hi2 Þ
ci1 ¼
ci ¼
ð29Þ
ðhi1 þ hi2 Þhi ðhiþ1 þ hi Þ ; hi2 hi1 ðhi þ hi1 Þðhiþ1 þ hi þ hi1 Þ
ð30Þ
2 2 2 2 hi2 hi hiþ1 hi1 hiþ1 hi2 hi1 hiþ1 þ hi2 hi þ 2 hi1 hi hiþ1 þ hi1 hi hi1 hi hi2 hi1 hi
ciþ1 ¼
hi1 ðhi1 þ hi2 Þhi ðhiþ1 þ hi Þ hi1 ðhi1 þ hi2 Þðhiþ1 þ hi Þ ; hi ðhi þ hi1 Þðhi þ hi1 þ hi2 Þhiþ1
ciþ2 ¼
hi1 ðhi1 þ hi2 Þhi hiþ1 ðhiþ1 þ hi Þðhiþ1 þ hi þ hi1 Þðhiþ1 þ hi þ hi1 þ hi2 Þ
;
ð31Þ
ð32Þ
ðFDMÞ;
ð33Þ
where the grid spacing is hi = xi+1 xi. The FLOP of the first-order derivative is counted as
FLOPFDM ðci2 ; ci1 ; ci ; ciþ1 ; ciþ2 Þ ¼ ð13; 12; 30; 11; 14Þ ¼ 87: fx
ð34Þ
The spatial derivative of the FDM consists of four types of the grid spacing (hi2, hi1, hi, and hi+1), and requires more computational cost than the IDO-CF scheme. The Flop advantage of IDO-CF scheme leads to shorter computational time: 775 ms for IDO-CF scheme and 1660 ms for FDM per 10,000 steps on Intel Core i7 2.8 GHz machine.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5793
The second-order spatial derivatives of the IDO-CF scheme and the 4th-order FDM are evaluated in a similar way.
FLOPIDOPV ð73Þ þ FLOPIDOIA ð53Þ fxx fxx ¼ 63; 2
FLOPIDO fxx ¼
¼ 107: FLOPFDM fxx
ð35Þ ð36Þ
3. Fourier analysis Turbulent phenomena contain a lot of small-scale vortices, because the scale of the turbulent vortex spreads proportional to Re3/4. Since the error of the numerical scheme is mainly caused by such small vortices, it is necessary to use high-order discretization for turbulent simulations in Fourier space. Although the spectral method [29] has the highest accuracy, it is only applicable to periodic boundary conditions. In this section, the accuracy and stability in the one-dimensional IDO-CF scheme are analyzed in comparison with the FDM in Fourier space, and the overall characteristics do not change in twoand three-dimensional stability analyses. 3.1. Basic formulation The spatial profile of the variable f(x) defined over the domain [0, L] with an uniform grid spacing Dx = L/N is decomposed into Fourier series
f ðxÞ ¼
X
^f ðkÞeiWx=Dx ;
ð37Þ
k
where i ¼
fj ¼
pffiffiffiffiffiffiffi 1, and W = 2pkDx/L is a scaled wavenumber. The values at the points xj and xj±1 are also decomposed as
X
^f ðkÞeiWxj =Dx ;
f j1 ¼
k
X
^f ðkÞeiWxj =Dx eiW :
ð38Þ
k
Using Eq. (38), the value at the point xj+m is simply decomposed as
fjþm ¼ fj eiWm :
ð39Þ x
The integrated average value fj of the IDO-CF scheme is also decomposed as x
fj ¼
1 Dx
Z
iW
þDx=2
f ðxj þ xÞdx;
¼ fj
Dx=2
iW
e 2 e 2 : iW
ð40Þ
Since this equation expresses the relation between the point value and its integrated average value, the stability and accuracy of the IDO-CF scheme can be examined by using Eqs. (39) and (40). 3.2. First-order spatial derivative In the IDO-CF scheme having the point value and integrated average value, we define a grid spacing h = Dx/2 and a scaled wavenumber w = 2pkh/L, and we compare the IDO-CF scheme with the FDM with the same total number of independent variables. The first-order derivatives of the IDO-CF scheme and FDM reduce to the following approximations, respectively, IDO: 4th-order central (point value and corresponding integrated average value)
fx ðxj Þ ¼ x
ðfjþ1 fj1 Þ 4 x fjþ1=2 x fj1=2 ; 2 Dx
fx ðxjþ1=2 Þ ¼
fjþ1 fj : Dx
ð41Þ
ð42Þ
FDM: 4th-order central (point value and corresponding integrated average value)
fx ðxj Þ ¼ x
fx ðxj Þ ¼
fjþ2 8f jþ1 þ 8f j1 fj2 ; 12h
fjþ1=2 fj1=2 ; h
x
¼
fjþ2 8x fjþ1 þ 8x fj1 x fj2 : 12h
ð43Þ
ð44Þ
The first-order derivatives of the point value and integrated average value are expressed as the different formulation in the IDO-CF scheme and are described in the following matrix form,
F nx ¼ A F n ;
ð45Þ
5794
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
Fn ¼
fn ; x n f
F nx ¼
fxn ; x n fx
A¼
a11
a12
a21
a22
:
ð46Þ
The notation F = (f, xf)T is a variable vector which consists of the point value and its integrated average value. The vector Fx is its first-order spatial derivative. The superscript n indicates the time step. The elements a11 and a12 for the point value (PV) of the matrix A vary with the accuracy of the first-order spatial derivatives and are obtained from Eqs. (41), (45) and (46)
a11 ¼
i sinð2wÞ ; h 2
a12 ¼
i 2 sinðwÞ ðIDO-CF PVÞ: h
ð47Þ
The elements a21 and a22 for the integrated average value (IA) do not change their values because of their definitions,
! Z þDx=2 @ 1 f ðxÞdx ; @x Dx Dx=2 f Dx f D2x ¼ 2 ; Dx i sinðwÞ: ¼f h
ðx f Þx ¼
ð48Þ
We obtain the following a21 and a22 from Eq. (48),
a21 ¼
i sinðwÞ; h
a22 ¼ 0 ðIDO-CF IAÞ:
ð49Þ
This is the exactly same as the FVM formulation. Fig. 6 shows the resolution of the first-order spatial derivative in Fourier space. The accuracy of the IDO-CF scheme at the point (IDO-CF PV) is almost the same as that of the FDM with 4th-order accuracy. The first-order spatial derivative on the cell (IDO-CF IA), defined as Eq. (49), is the analytically exact solution. 3.3. Fourier analysis for linear advection equation In this subsection, we examine a linear advection equation. In order to discuss the time-integration method, we consider the following Taylor expansion around the time equal to tn,
f nþ1 ¼ f n þ
! 2 ! k n @f @ 2 f n dt @ k f n dt dt þ þ þ : @t 2! k! @t2 @t k
ð50Þ
In the advection equation, a time derivative is written by the following equation,
@f @f ¼ u ; @t @x
ð51Þ
where the velocity u is assumed to be constant in time and space. The time derivatives in Eq. (50) are replaced by using the relation Eq. (52)
! ! ! @ @ 2 ðu dtÞ2 @ k ðu dtÞk ðu dtÞ þ 1þ þ þ f n: @x @x2 2! @xk k!
¼
3 2.5
Imaginary
10
IDO-CF PV FDM (PV & IA) Exact (IDO-CF IA)
Error of spatial derivative
f
nþ1
2 1.5 1 0.5
0.5
1
1.5
2
Wavenumber
2.5
3
ð52Þ
0
10
-2
10
-4
10
-6
10
-8
10
-10
10
-12
10
IDO-CF PV FDM (PV & IA)
-3
10
-2
10
-1
10
0
Wavenumber
Fig. 6. Fourier analysis for the first-order spatial derivative: (left) imaginary part and (right) error; IDO-CF IA is not drawn in the right figure as there does not exist any error.
5795
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
The Runge–Kutta method is based on the Taylor expansion to obtain the accuracy in time. The resolutions for the advection equation are shown in Fig. 7. The accuracy of the point value (IDO-CF PV) is almost the same as that of the 4th-order FDM. The integrated average values (IDO-CF IA) are quite better than those of the FDM in highwavenumber area. 3.4. The stability analysis for the advection equation The eigenvalues of the numerical scheme show the stability of the time-integration. The IDO-CF scheme has various kinds of moments. The matrix eigenvalue problem is solved as
F nþ1 ¼ DT;k F n ; DT;k ¼ PKk P1 ;
Kk ¼
kþ
0
0
k
;
ð53Þ
where kk is an eigenvalue matrix and the components k+ and k are eigenvalues of the point value and integrated average value. Fig. 8 shows the contour plot of the eigenvalue k = max(jk+j, jkj) as functions of Courant–Friedrichs–Lewy (CFL) number and wavenumber, and the region where k 6 1 is stable. For the 2nd-order Runge–Kutta method, the stable area is too small (CFL < 0.05) to solve the linear advection equation. When we use the 3rd-order Runge–Kutta time integration, the IDO-CF scheme is stable for CFL < 1.15. The stable region of the CFL number for the 4th-order Runge–Kutta method is enlarged to CFL < 1.87. Therefore the IDO-CF scheme requires at least 3rd-order accuracy in time. In the following section of this paper, numerical simulations are carried out by using the 3rd-order Runge–Kutta method.
4. Governing equations and SGS models 4.1. Governing equations The governing equations for incompressible flows consist of the Navier–Stokes equations and the continuity equation,
@ui @ui uj @p 1 @ 2 ui ¼ þ ; þ @xi Res @xj @xj @t @xj
@ui ¼ 0: @xi
ð54Þ
Here ui, p and Res denote the velocity, pressure and the Reynolds number which is based on the friction velocity and the channel half-width. In LES, governing equations are respectively filtered in physical space, and filtered variables f denoted by an overbar are defined as
f x; D; t ¼
Z
0 f ðx0 ; tÞG x; x0 ; D dx ;
ð55Þ
where G is the filter function and D is the filter width. We adopt that the filter width is the same as the grid spacing. The grid filter function G is normalized to satisfy the following condition,
0 G x; x0 ; D dx ¼ 1:
ð56Þ
3
10
IDO-CF PV IDO-CF IA FDM (PV & IA) Exact
Error of phase
2.5 2
Phase
Z
1.5 1 0.5
0.5
1
1.5
2
Wavenumber
2.5
3
0
10
-2
10
-4
10
-6
10
-8
10
-10
10
-12
10
IDO-CF PV IDO-CF IA FDM (PV & IA)
-3
10
-2
10
-1
10
0
Wavenumber
Fig. 7. Fourier analysis for the advection equation with 4th-order Runge–Kutta method: (left) phase and (right) phase error.
5796
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
2
2
1.5
1.5
Courantnumber
1 1.
Courantnumber
Unstable
1
1.0
5
0.95
1
CFL=1.15
Stable
0.5
0.5 1 0
1 1.0 5 1.1
Unstable
Stable 0
CFL=0.05 0.5
1
0
1.5
0
0.5
1.5
Unstable
2
0.9 1 1.1
CFL=1.87
0. 0.9 95
1.5
Courantnumber
1
Wavenumber
Wavenumber
Stable 1
0.5
0
0
0.5
1
1.5
Wavenumber : (top-left) 2nd-, (top-right) 3rd-, (bottom) 4th-order Fig. 8. Contour plots of eigenvalues for time integration of the advection equation. CFL ¼ udt h Runge–Kutta method.
Applying the grid filter to Eq. (54), we obtain the filtered Navier–Stokes equations and the filtered continuous equation,
i u j i i @ u @ sij @u @p 1 @2u ¼ þ ; þ @xi @xj Res @xj @xj @t @xj
i @u ¼ 0: @xi
ð57Þ
j which is derived from the non-linear coni u The difference between Eqs. (54) and (57) is the SGS stress tensors sij ¼ ui uj u vection term. Since the filtered Navier–Stokes equations are not closed by themselves, the SGS stress tensor should be modeled. 4.2. Subgrid-scale models 4.2.1. Dynamic Smagorinsky model In the Smagorinsky eddy viscosity model, the SGS stress tensor sij is modeled with the filter width D as
sij ¼ 2C D2 jSjSij ¼ 2mSGS Sij ; Sij ¼
j @ u i 1 @u ; þ 2 @xi @xj
1=2 jSj ¼ 2Sij Sij ;
ð58Þ
where Sij and S are the velocity-strain tensor and its magnitude in the GS, respectively. In the Smagorinsky model (SM) [3], the model coefficient C is constant in an entire computational domain, and the SM does not satisfy a correct asymptotic behavior to a wall. The dynamic Smagorinsky model (DSM) [6,7] improves those defects by calculating the model parameter dynamically. The model parameter of the DSM is determined by using a least square procedure h,i with an average in homogeneous directions as
C DSM ¼
hLij M ij i ; hMij M ij i
ð59Þ
b and the grid scale D. They are given by where Lij and Mij are the resolved stress tensors between the test scale D
^ i u ^ j ; d j u Lij ¼ u iu
ð60Þ
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
c d b 2 b Mij ¼ 2D2 j SjSij 2 D S Sij :
5797
ð61Þ
b ¼ 2D, and the test-filtered value ^ f in homogeneous directions is calculated by using SimpThe test filter width is chosen as D son’s rule
1 ^ f i ¼ ðf i1 þ 4f i þ f iþ1 Þ: 6
ð62Þ
4.2.2. Coherent-structure model The model parameter of the DSM is determined dynamically by the width ratio of grid filter and test filter. This procedure is the most notable breakthrough in LES. However, the averaging in homogeneous directions is required for numerical stability. That makes it difficult to treat a complex geometry with the DSM. In the CSM [12], the model parameter C is dynamically calculated from a flow based on turbulence structures. No averaging in homogeneous directions is needed, so that this eddy viscosity model can treat a complex geometry. The model parameter of the CSM is defined as
C ¼ C CSM jF CS j3=2 ;
ð63Þ
where CCSM = 1/30 is a fixed model parameter and FCS is the coherent structure function defined as
F CS ¼
Q E
ð1 6 F CS 6 1Þ;
ð64Þ
Q¼
j i 1 1 @u @u ; W ij W ij Sij Sij ¼ 2 2 @xi @xj
ð65Þ
E¼
1 @ u j j 1 @u : W ij W ij þ Sij Sij ¼ 2 2 @xi @xi
ð66Þ
The coherent structure function FCS, composed of the second invariant Q and the magnitude of a velocity gradient tensor E, plays a role of a wall function, and its magnitude is less than unity. Therefore the model parameter can be locally determined without averaging in homogeneous directions. The CSM shows almost the same performance as the DSM and is applicable to complex geometries [13]. In the following section, the SGS stress tensor is calculated by using the CSM. 5. Turbulent channel flow simulations The turbulent channel flows have been studied for several decades as one of the most basic benchmarks of wall turbulence and are often used for validation of numerical schemes. The flows between two parallel walls are driven by a mean pressure gradient. The periodic boundary condition is considered in the streamwise and spanwise directions (x1 and x3). The computational grid is uniform in the periodic directions. The non-slip condition is imposed on the wall. The computational grids in this section are concentrated near the wall by using a hyperbolic-tangent function. Reynolds number is based on the channel half width d and friction velocity us (Res = usd/m). The turbulent statistics are obtained with time–space averaging after the turbulence is well developed to a statistically stable state. The IDO-CF scheme has better accuracy than conventional FDMs with 4th-order accuracy for the linear advection equation. The spectral like resolution has been confirmed by solving two-dimensional homogeneous isotropic turbulence [26]. We solve the incompressible Navier–Stokes equations to satisfy the velocity divergence free condition by using the SMAC method. The BiCGSTAB method is employed to solve the sparse matrix of the Poisson equation. The 3rd-order Runge–Kutta method is used for the time-integration with the time step of CFL = 0.1. In this section, the results of the IDO-CF scheme with the CSM are compared to those of the 4th-order FDM with the DSM and CSM as well as DNS data in Moser et al. [27] at Res = 180. We also carried out the cases of different Reynolds numbers Res = 395 and 590 with larger grid numbers and found almost the same results. 5.1. Effect of grid resolution for the IDO-CF scheme LES of the turbulent channel flows is conducted at Res = 180 in various cases of the grid resolution. In the IDO-CF scheme, the number of degrees of freedom is twice compared to conventional FDMs in every direction. Numerical tests are diligently accumulated with the same total number of degrees of freedom. The governing equations are discretized by the central formulation (Eqs. (2) and (3)). The influence of wall-normal grid resolution is examined with changing the grid resolution in the periodic directions. Figs. 9 and 10 show the maximum values of the mean velocity and the turbulent intensity in the streamwise direction. þ The horizontal axis indicates a grid spacing in periodic directions defined as xþ 1 ¼ 3x3 ¼ Res 4p=N 1 . The number of grid þ min point in the wall direction is defined as 16, 24, 32, and 64 which correspond to x2 ¼ 14:9, 5.8, 2.7, and 1.3.
5798
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
þ min The mean velocity is overestimated about 4% at xþ ¼ 2:7, i.e. (N1, N2, N3) = (8, 32, 8) and comes closer to 1 ¼ 282 and x2 the DNS from the higher value with increasing the number of grid point. However, the maximum value of the mean velocity þ min is still underestimated at xþ ¼ 1:3, i.e. (N1, N2, N3) = (32, 64, 32) in Fig. 9. It is well known that the mean veloc1 ¼ 70 and x2 ity of the 4th-order scheme without SGS models is below that of DNS because of their truncation errors and the lack of grid resolution. For that reason, the SGS stress should be applied to compensate for the lack of grid resolution. The maximum value of the turbulent intensity approaches the result of DNS in each case and is mainly related to the resolution of grid spacing in both two periodic directions. It is confirmed that the maximum values of the mean velocity and the turbulent intensity are getting closer to the DNS in much finer grid spacing xþ 1 ¼ 47. The effect of grid resolution on the turbulent statistics is similar to the FDM [30,31]. This is why a proper SGS model should be applied to LES using the IDO-CF scheme.
5.2. The effectiveness of SGS model In LES, the flow field in GS is directly calculated, and the effect of the SGS is evaluated with the SGS stress tensor
sij ¼ 2mSGS Sij with the CSM. The SGS eddy viscosity becomes zero near the wall without a wall-damping function. The averaging in homogeneous direction is not conducted. In this subsection, the IDO-CF scheme with the CSM is compared to the IDO-CF scheme without the SGS model and the IDO-CF upwind scheme (Eq. (9)). The grid resolution is shown in Table 1. Fig. 11 shows the mean velocity profiles at (N1, N2, N3) = (16, 32, 16) (case 1). In the previous subsection, the IDO-CF scheme without the SGS model underestimates the maximum value of the mean velocity. The SGS stress tensor modifies the underestimation, and the IDO-CF scheme with the CSM is in good agreement with DNS. The IDO-CF upwind scheme implicitly introduces numerical viscosity (Eq. (13)) too much, so that it overestimates the profile of the mean velocity. Fig. 12 shows the turbulent intensity in the streamwise direction. The IDO-CF schemes with and without the CSM obtain almost the same profiles and have good convergences to the result of DNS. The IDO-CF upwind scheme overestimates the turbulent intensity as well because of the numerical viscosity.
20
N1=8
N1=12
N1=48 N =32 1
N1=16
U 1+
15
dx +2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS + dx1 = 3 dx 3+
10
5
(N2=16) (N2=24) (N2=32) (N2=64)
0 0
50
100
150
200
250
300
dx1+ Fig. 9. The maximum values of the mean velocity at various grid resolution without the SGS model.
N1=8 4
N1=12 N1=16
u’ 1
3
N1=48 N =32 1 +
2
dx 2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS + dx1 = 3 dx 3+
1
(N2=16) (N2=24) (N2=32) (N2=64)
0 0
50
100
150
200
250
300
dx1+ Fig. 10. The maximum values of the turbulent intensity at various grid resolution without the SGS model.
5799
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805 Table 1 The numerical condition of turbulent channel flows.
Case Case Case Case
1 2 3 4
þ
þ
þ
Res
L1 L2 L3
N1 N2 N3
dx1
dx2 min
dx3
180 180 395 590
4p 2 4/3p 4p 2 4/3p 2p 2 p 2p 2 p
16 32 16 32 64 32 32 64 32 32 64 32
141 70 77 116
2.7 1.3 2.9 4.3
47 23 39 58
25
IDO with CSM w/o CSM upwind DNS
20
U+
15
10
5
0 0 10
10
1
10
x2 +
2
Fig. 11. The profiles of the mean velocity using the CSM, upwind scheme and without CSM (case 1).
IDO with CSM w/o CSM upwind DNS
4
u’ 1
3
2
1
0 0
50
100
x 2+
150
Fig. 12. The profiles of the turbulent intensity using the CSM, upwind scheme and without CSM (case 1).
Figs. 13 and 14 show the SGS eddy viscosity and the shear stress profile of the IDO-CF scheme with the CSM. In the turbulent channel flows, the pressure gradient is balanced against the sum of total shear stress and viscosity stress. The total shear stress is composed of Reynolds stress u01 u02 and SGS stress s12 in LES. In Fig. 14, the SGS model gives the additional shear stress to the flow field and total shear stress is in good agreement with DNS. Figs. 15 and 16 show the maximum value of the mean velocity and the turbulent intensity at Res = 180. The spatial resolution is the same as Figs. 9 and 10. The profile of the mean velocity is improved well to compare with the profile of the IDOCF scheme without the SGS model. The profile of the turbulent intensity converges to the value of DNS. 5.3. The comparison with finite difference scheme Numerical results of the IDO-CF scheme are compared to the 4th-order FDM [15]. Note that the 2nd-order FDM is used in the wall direction because of facilitating the treatment of the boundary condition on the wall. The 4th-order FDM is often used to calculate turbulent flows because of their conservation of kinetic energy. The DSM and CSM are respectively applied to the FDM. The IDO-CF scheme has 4th-order accuracy in every direction and simulations are carried out at the same number of grid point (case 1). For example, the number of grid (N1, N2, N3) = (16, 32, 16) with the IDO-CF scheme corresponds to that of grid (N1, N2, N3) = (32, 64, 32) with the FDM.
5800
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
0.005
CSM
0.003
ν
SGS
0.004
0.002
0.001
0
0
50
100
x2
150
+
Fig. 13. The profile of the SGS eddy viscosity with the CSM (case 1).
1
u’ 1 u’ 2 τ 12 u’ 1 u’ 2 + τ 12 DNS u’ 1 u’ 2
2
0.6
u’ 1 u’
0.8
0.4
0.2
0 0
50
100
150
x 2+
Fig. 14. The profiles of the Reynolds shear stress with the CSM (case 1).
N1=48
20
N1=32
N1=8
N1=12
N1=16
U 1+
15
dx +2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS + dx1 = 3 dx 3+
10
5
0
0
50
100
150
dx1
200
(N2=16) (N2=24) (N2=32) (N2=64)
250
300
+
Fig. 15. The maximum values of the mean velocity at various grid resolution with the SGS model.
Fig. 17 shows the mean velocity profiles of the IDO-CF scheme with the CSM and the FDM with the DSM and CSM. All profiles are in good agreement with the profile of DNS. Figs. 18–20 show the profiles of the turbulent intensity in the streamwise, the wall, and the spanwise directions. The turbulent intensity in the streamwise direction strongly depends on the accuracy of a numerical scheme and the grid resolution. The streamwise turbulent intensity of the IDO-CF scheme with
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5801
N1=8
4
N1=12 N1=16
u’ 1
3
N1=48 N1=32 +
2
dx 2 min= 14.9 = 5.8 = 2.7 = 1.3 DNS dx1+ = 3 dx 3+
1
0
0
50
100
150
dx1
200
(N2=16) (N2=24) (N2=32) (N2=64)
250
300
+
Fig. 16. The maximum values of the turbulent intensity at various grid resolution with the SGS model.
20
U
+
15
IDO with CSM FDM with CSM FDM with DSM DNS
10
5
0 0 10
10
1
10
x2 +
2
Fig. 17. The profiles of the mean velocity using the IDO-CF and FDM.
IDO with CSM FDM with CSM FDM with DSM DNS
3
u’ 1
2
1
0 0
50
100
x 2+
150
Fig. 18. The profiles of the turbulent intensity in the streamwise direction using the IDO-CF and FDM.
the CSM is lower than that of FDM with the CSM and DSM. The present scheme achieves better agreement with the DNS. In Figs. 19 and 20, it is found that the wall and spanwise turbulent intensities of the IDO-CF scheme are higher than those of FDM. The turbulence field of the IDO-CF scheme is more isotropic than that of FDM. The IDO-CF scheme gives good results for the incompressible Navier–Stokes equations with less number of grid points.
5802
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
IDO with CSM FDM with CSM FDM with DSM DNS
u’ 2
1
0.5
0 0
50
100
150
x 2+
Fig. 19. The profiles of the turbulent intensity in the wall direction using the IDO-CF and FDM.
1.5
IDO with CSM FDM with CSM FDM with DSM DNS
u’ 3
1
0.5
0 0
50
100
150
x 2+
Fig. 20. The profiles of the turbulent intensity in the spanwise direction using the IDO-CF and FDM.
25
U+
20
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
15
10
5
0 0 10
10
1
x2 +
10
2
Fig. 21. The profiles of the mean velocity at Res = 180, 395, and 590. (case 2–4).
5.4. The flow at various Reynolds number (Res = 180, 395, and 590) The IDO-CF scheme obtains better results in comparison with conventional FDMs of the same order accuracy. This is because the present scheme has higher spectral resolution in Fourier space. Further LES is conducted at higher Reynolds number as shown in Table 1. Fig. 21 shows the mean velocity profiles at (N1, N2, N3) = (32, 64, 32) (case 2–4). The profiles at each Reynolds number are in good agreement with the DNS. Figs. 22–24 show the turbulent intensity profiles in streamwise, wall, and spanwise directions at (N1, N2, N3) = (32, 64, 32) (case 2–4). The profiles of the mean velocity and the turbulent intensities at Res = 180 are in
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
3
5803
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
u’ 1
2
1
0 0
100
200
300
x 2+
400
500
Fig. 22. The profiles of the turbulent intensity in the streamwise direction at Res = 180, 395, and 590. (case 2–4).
1.5
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
u’ 2
1
0.5
0 0
100
200
300
x2
+
400
500
Fig. 23. The profiles of the turbulent intensity in the wall direction at Res = 180, 395, and 590. (case 2–4).
IDO with CSM (Re τ = 180) (Re τ = 395) (Re τ = 590) DNS (Re τ = 180, 395, 590)
u’ 3
1.5
1
0.5
0 0
100
200
300
x 2+
400
500
Fig. 24. The profiles of the turbulent intensity in the spanwise direction at Res = 180, 395, and 590. (case 2–4).
good agreement with the DNS. At Res = 395 and 590, the profiles of the mean velocity and the turbulent intensities in streamwise, wall, and spanwise directions show good convergence to the DNS in spite of less number of grid points, and the maximum values of the turbulent intensities in the streamwise direction are in good agreement with the DNS in Fig. 22. The IDO-CF scheme with the CSM can be applied to various Reynolds number with a small grid number, and is one of good numerical schemes for turbulent phenomena.
5804
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
6. Summary and conclusions The IDO-CF scheme with the subgrid-scale model is applied to turbulent phenomena, possessing the following three advantages; (a) the implementation of the wall boundary condition with high-order accuracy is as easy as that with the 2nd-order FDM, (b) the computational cost for non-uniform grid spacing is much smaller than that of the 4th-order FDM, and (c) the IDO-CF scheme has better spectral resolution than the same order FDM. The first-order spatial derivative of the integrated average value is exactly the same as the FVM. Fourier analysis shows that the spatial derivative of the IDO-CF scheme has better spectral resolution in comparison with the 4th-order FDM. The IDO-CF scheme is suitable for the turbulent simulations because LES is sensitive to the accuracy and resolution in Fourier space. In order to evaluate the IDO-CF scheme with the CSM, the fully developed turbulent channel flows are examined at various grid resolutions. Numerical simulations at three Reynolds numbers Res = 180, 395, and 590 are compared to the reference of DNS data in Moser et al. [27]. We have three major conclusions. First, it is found that the IDO-CF scheme without the CSM underestimates the profile of the mean velocity in relatively coarse grid resolution, so that the SGS model is needed. The IDO-CF upwind scheme overestimates the profiles of the mean velocity and the turbulent intensity with the excessive numerical viscosity introduced implicitly. The CSM evaluates the SGS eddy viscosity correctly, and numerical results of the IDO-CF scheme with the CSM are in good agreement with the DNS. Second, the IDO-CF scheme is compared to the FDM with the DSM and the CSM. The IDO-CF scheme with the CSM achieves a better profile of the turbulent intensity in the streamwise direction. In LES of turbulent channel flows, the statistical value of the turbulent intensity mainly depends on the grid resolution and discretized accuracy in both two periodic directions. The IDO-CF scheme has higher spectral resolution for the incompressible Navier–Stokes equations than conventional FDMs with the same order of accuracy. Finally, we validate the dependence of the IDO-CF scheme with the CSM on various Reynolds numbers. The statistical values at Res = 180 are in good agreement with the DNS at (N1, N2, N3) = (32, 64, 32). At Res = 395 and 590, the mean velocity and the turbulent intensity show good convergence to the DNS in spite of less number of grid points. It is concluded that the present scheme is one of most promising candidates to simulate turbulence. In the future work, we apply the present scheme to turbulent phenomena with the complex geometry. Acknowledgements This research was supported in part by KAKENHI, Grant-in-Aid for Scientific Research (B) 19360043 from The Ministry of Education, Culture, Sports, Science and Technology (MEXT), JST-CREST project ‘‘ULP-HPC: Ultra Low-Power, High Performance Computing via Modeling and Optimization of Next Generation HPC Technologies’’ from Japan Science and Technology Agency (JST) and JSPS Global COE program ‘‘Computationism as a Foundation for the Sciences’’ from Japan Society for the Promotion of Science (JSPS). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
M. Lesieur, O. Metais, New trends in large-eddy simulations of turbulence, Annu. Rev. Fluid Mech. 28 (1) (1996) 45–82. C. Meneveau, J. Katz, Scale-invariance and turbulence models for large-eddy simulation, Annu. Rev. Fluid Mech. 32 (1) (2000) 1–32. J. Smagorinsky, General circulation model of the atmosphere, Mon. Weather Rev. 164 (1963) 91–99. E.R. Van Driest, On turbulent flow near a wall, J. Aeronaut. Sci 23 (1956) 1007. C. Haertel, L. Kleiser, Analysis and modelling of subgrid-scale motions in near-wall turbulence, J. Fluid Mech. 356 (1998) 327–352. M. Germano, U. Piomelli, P. Moin, W. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A: Fluid Dyn. 3 (1991) 1760. D. Lilly, A proposed modification of the germano subgrid-scale closure method, Phys. Fluids A: Fluid Dyn. 4 (1992) 633. Y. Zang, R. Street, J. Koseff, A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows, Phys. Fluids A: Fluid Dyn. 5 (1993) 3186. B. Vreman, B. Geurts, H. Kuerten, On the formulation of the dynamic mixed subgrid-scale model, Phys. Fluids 6 (12) (1994). C. Meneveau, T.S. Lund, W.H. Cabot, A Lagrangian dynamic subgrid-scale model of turbulence, J. Fluid Mech. 319 (353) (1996). E. Bou-Zeid, C. Meneveau, M. Parlange, A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows, Phys. Fluids 17 (2005) 025105. H. Kobayashi, The subgrid-scale models based on coherent structures for rotating homogeneous turbulence and turbulent channel flow, Phys. Fluids 17 (2005) 045104. H. Kobayashi, F. Ham, X. Wu, Application of a local sgs model based on coherent structures to complex geometries, Inter. J. Heat Fluid Flow 29 (3) (2008). S. Ghosal, An analysis of numerical errors in large-eddy simulations of turbulence, J. Comput. Phys. 125 (1) (1996) 187–206. Y. Morinishi, T. Lund, O. Vasilyev, P. Moin, Fully conservative higher order finite difference schemes for incompressible flow, J. Comput. Phys. 143 (1) (1998) 90–124. Y. Liu, M. Vinokur, Z. Wang, Spectral difference method for unstructured grids. I: basic formulation, J. Comput. Phys. 216 (2) (2006) 780–801. T. Yabe, T. Aoki, A universal solver for hyperbolic equations by cubic-polynomial interpolation. I: One-dimensional solver, Comput. Phys. Commun. 66 (2–3) (1991) 219–232. T. Yabe, T. Ishikawa, P. Wang, T. Aoki, Y. Kadota, F. Ikeda, A universal solver for hyperbolic equations by cubic-polynomial interpolation. II: Two-and threedimensional solvers, Comput. Phys. Commun. 66 (1991) 233–242. T. Aoki, Interpolated differential operator (IDO) scheme for solving partial differential equations, Comput. Phys. Commun. 102 (1–3) (1997) 132–146. K. Sakurai, T. Aoki, W. Lee, K. Kato, Poisson equation solver with fourth-order accuracy by using interpolated differential operator scheme, Comput. Math. Appl. 43 (6–7) (2002) 621–630. H. Yoshida, T. Aoki, T. Utsumi, Improvement of accuracy and stability in numerically solving hyperbolic equations by IDO (interpolated differential operator) scheme with Runge–Kutta time integration, Electron. Commun. Jpn. (Part III: Fundam. Electron. Sci.) 87 (2) (2004). T. Utsumi, T. Kunugi, T. Aoki, Stability and accuracy of the cubic interpolated propagation scheme, Comput. Phys. Commun. 101 (1–2) (1997) 9–20.
N. Onodera et al. / Journal of Computational Physics 230 (2011) 5787–5805
5805
[23] Y. Imai, T. Aoki, Accuracy study of the IDO scheme by Fourier analysis, J. Comput. Phys. 217 (2) (2006) 453–472. [24] T. Yabe, R. Tanaka, T. Nakamura, F. Xiao, An exactly conservative semi-Lagrangian scheme (CIP–CSL) in one dimension, Mon. Weather Rev. 129 (2) (2001) 332–344. [25] F. Xiao, T. Yabe, Completely conservative and oscillationless semi-Lagrangian schemes for advection transportation, J. Comput. Phys. 170 (2) (2001) 498–522. [26] Y. Imai, T. Aoki, K. Takizawa, Conservative form of interpolated differential operator scheme for compressible and incompressible fluid dynamics, J. Comput. Phys. 227 (4) (2008) 2263–2285. [27] R. Moser, J. Kim, N. Mansour, Direct numerical simulation of turbulent channel flow up to re = 590, Phys. Fluids 11 (1999) 943. [28] J. Kim, P. Moin, R. Moser, Turbulence statistics in fully developed channel flow at low reynolds number, J. Fluid Mech. 177 (1) (1987) 133–166. [29] D. Gottlieb, A. Orszag, Numerical Analysis of Spectral Methods, SIAM, Philadelphia, 1977. [30] Y. Morinishi, O. Vasilyev, A recommended modification to the dynamic two-parameter mixed subgrid scale model for large eddy simulation of wall bounded turbulent flow, Phys. Fluids 13 (2001) 3400. [31] Y. Morinishi, O. Vasilyev, Vector level identity for dynamic subgrid scale modeling in large eddy simulation, Phys. Fluids 14 (2002) 3616.