A novel numerical forcing for homogeneous stratified turbulence in full energy equilibrium

Computers & Fluids 39 (2010) 1789–1795

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A novel numerical forcing for homogeneous stratified turbulence in full energy equilibrium S. Hirabayashi, T. Sato * Department of Ocean Technology, Policy, and Environment, University of Tokyo, 5-1-5 Kashiwa-no-ha, Kashiwa 277-8583, Japan

a r t i c l e

i n f o

Article history: Received 2 December 2009 Received in revised form 2 June 2010 Accepted 2 June 2010 Available online 11 June 2010 Keywords: Linear forcing Homogeneous stratified turbulence Full equilibrium Stationary turbulence Direct numerical simulation

a b s t r a c t A new numerical method to simulate stationary, homogeneous, and stratified turbulence was developed with an arbitrary energy dissipation rate. Conventional numerical simulations for homogeneous stratified turbulence were accomplished by given mean shears and generated turbulence field either developed or decayed with time. In order to keep stationary turbulence, a linear forcing method for anisotropic turbulence was developed to apply stratified turbulence, where shear was calculated dynamically to satisfy the condition of energy equilibrium. The present method was implemented to direct numerical simulations to realize the generation of small-scale eddy fields, in which the Reynolds number based on the Taylor microscale was about 50. It was confirmed that simulated flow field successfully maintained energy equilibrium with keeping the given dissipation rate. The gradient and flux Richardson numbers determined by the balance of turbulence intensity and stratification were also found to be constant with time. As an application, diffusion of massless dye was numerically simulated in the reproduced turbulent field. The calculated vertical diffusivity was found to be well comparable with that of active heat estimated by a conventional model. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Stratification is one of the main characteristics of ocean turbulence and it is known to suppress the vertical diffusion of mass and heat in the ocean. Numerical simulation is useful to extract the statistics of stratified turbulence since it is capable to control turbulence intensity and stratification systematically. When thinking of numerical simulations of stably stratified turbulence, a driving force is necessary to keep turbulent motions. Gerz et al. [1] carried out direct numerical simulations (DNS) of homogeneous turbulence driven by mean shear under linear stratifications and reproduced turbulence fields that develop or decay depending on given Richardson numbers. They found the stationary Richardson number, RS, of about 0.1 at which the turbulence was nearly stationary without developing or decaying. Similar methods were used by Holt et al. [2] and Shih et al. [3] in their DNS, while Kaltenbach et al. [4] adopted this method in the large eddy simulations. However, in their studies, RS varies depending on the combination of mean shear and stratification and the determination of RS has to rely on trial and errors. Sato et al. [5] numerically generated a small-scale ocean turbulence by forcing low-wavenumber components of velocity, which were measured in the real ocean. However, because their forcing components change with time, it is still not suitable to apply the method to study turbulent behavior * Corresponding author. Tel.: +81 4 7136 4726; fax: +81 4 7136 4727. E-mail address: [email protected] (T. Sato). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.06.001

statistically. For the statistical evaluation of physical processes, turbulence properties should be kept constant, in other words, in full equilibrium. Lundgren [6] developed a linear forcing method under the condition of the energy equilibrium to generate isotropic turbulence for constant energy dissipation rates by DNS. The capability to control the energy dissipation rate of the reproduced turbulence fields in Lundgren [6] is effective to investigate the relationship between the elementary parameters of flow field and resultant turbulence statistics. Based on the assumption that the energy equilibrium is kept in the small-scale ocean turbulence, a similar forcing concept to that of Lundgren [6] may have a possibility to reproduce anisotropic turbulence in stratified fluids. In this study, a linear forcing method is developed to reproduce homogeneous shear-driven turbulence with arbitrary stratifcations and energy dissipation rates. The proposed method is implemented to DNS and the verification of this method is demonstrated. As an application of the present method, diffusion of massless dye was numerically simulated in the generated stationary turbulence field and its vertical diffusivity was estimated. 2. Materials and methods 2.1. Linear forcing of Lundgren [6] for isotropic turbulence At first, the original linear forcing of Lundgren [6] for isotropic turbulence is reviewed. The governing equations for isotropic tur-

1790

S. Hirabayashi, T. Sato / Computers & Fluids 39 (2010) 1789–1795

bulence are the equation of continuity and the Navier–Stokes (NS) equation:

@uj ¼ 0; @xj

ð1Þ

@ui @ui 1 @p @ 2 ui ¼ þm 2 ; þ uj @t @xj q0 @xi @xj

ð2Þ

where ui is the velocity, p is the pressure from which the hydrostatic pressure was subtracted, q0 is the standard density, m is the kinematic viscosity of water. Indices i, j = 1, 2, and 3 indicate the components in the x, y, and z directions, respectively, of the Cartesial coordinates. Repeated indices mean summation. Now we divide primitive variables into mean and fluctuating components,

u0i ; 0

ui ¼ hui i þ p ¼ hpi þ p ;

ð3Þ ð4Þ

where an angle bracket denotes volumetric average and primes mean fluctuations. By applying volumetric average to Eqs. (1) and (2), we obtain

@huj i ¼ 0; @xj   @hui i @ui 1 @hpi @ 2 hui i ¼ þ uj þm ; @t @xj q0 @xi @x2j

@xj

¼ 0;

@u0i @u0 @u0 1 @p0 þm þ huj i i þ u0j i ¼  @t @xj @xj q0 @xi

ð18Þ

ð6Þ

@T @T @2T þ uj ¼ jT 2 ; @t @xj @xj

ð19Þ

ð9Þ

where fi is the forcing term, which is given by

fi ¼

ð10Þ u0i

where Q is a constant. By multiplying to Eq. (9) and taking the volumetric average, one obtains the transport equation of kinetic energy:

@hEi ¼ hei þ hPi; @t

ð11Þ

where E, e, and P are the kinetic energy, the energy dissipation rate, and the energy production rate, respectively, expressed as

1 0 0 uu; 2 i i @ 2 u0 e ¼ mu0i 2i ; @xj E¼

P ¼ Qu0i u0i :

ð12Þ ð13Þ ð14Þ

For stationary turbulence, the left-hand side (LHS) of Eq. (11) should be zero, and we obtain

Q ¼

e  :

u0i u0i

Here we try to apply the idea of linear forcing of Lundgren [6] to stratified turbulence where turbulence structure is anisotropic. The governing equations for stratified turbulence are the equation of continuity, and the NS equations with the Boussinesq approximation, and the heat transfer equation:

@ui @ui 1 @p @ 2 ui ¼ þ m 2 þ gbTdi3 ; þ uj @t @xj q0 @xi @xj

The last term of the right-hand side (RHS) is corresponding to the production of turbulent kinetic energy, u0i u0j @hui i=@xj , and Eq. (8) indicates that the fluctuating velocity is driven by mean shear. Although the gradient of mean velocity is zero in homogeneous turbulence, it is inferred that the appropriate forcing term for isotropic and stationary turbulence is also isotropic and proportional to the fluctuating velocity, u0i [7]. Therefore, the NS equation for the fluctuations become

Qu0i ;

2.2. Modification for stratified turbulence

ð5Þ

  @ui @hui i  u0j þ uj : ð8Þ @xj @xj

@u0i @u0 @ 2 u0 1 @p0 þ m 2i þ fi ; þ u0j i ¼  @t @xj q0 @xi @xj

Eq. (16) is the momentum equation for stationary turbulence where the forcing term is dynamically calculated so that the energy dissipation rate, e0, is always constant. Since the forcing term is linear to u0i , this method was named the linear forcing by Lundgren [6]. Rosales and Meneveau [7] performed DNS of isotropic turbulence with the linear forcing and showed good agreement with the result of a conventional forcing in the wavenumber domain.

ð17Þ

ð7Þ @ 2 u0i @x2j

ð16Þ

j j

@uj ¼ 0; @xj

Subtracting Eqs. (5) and (6) from Eqs. (1) and (2) gives

@u0j

Then Eq. (9) is rewritten as

@u0i @u0 @ 2 u0 1 @p0 e E 0 ui : þ m 2i þ D þ u0j i ¼  @t @xj q0 @xi @xj u0 u0

where T is the temperature, g is the gravity acceleration, b is the coefficient of thermal expansion, jT is the molecular thermal diffusivity in water, and dij is the Kronecker delta. In this study, it is assumed that linear density stratification is made only by temperature and the effect of salinity is neglected. Then we divide temperature into linear stratification and fluctuating terms:

T ¼ hTih þ T 0 ;

where hih is the horizontal average and s = d hTih/dx3 becomes constant because of the assumption of linear stratification. Then we can rewrite Eqs. (18) and (19) as

  @ui @ui 1 @p @ 2 ui ¼  q0 gbhTih þ m 2 þ gbT 0 di3 ; þ uj @t @xj q0 @xi @xj

ð21Þ

@T 0 @T 0 @2T 0 ¼ jT 2  u3 s: þ uj @t @xj @xj

ð22Þ

For generating stationary turbulence, Lundgren [6] took the assumption of isotropy, in which the magnitude of shear components are equivalent in all directions. However, this assumption cannot be applied to stratified turbulence, because the main source of energy for turbulence is assumed to be vertical shear in anisotropic turbulence. Based on this assumption, velocity is decomposed into a vertical shear and a fluctuating components:

ui ¼ hui ih þ u0i ;

ð23Þ

where the vertical gradient of huiih corresponds to the vertical shear in the x1 direction, Shdi1. Then Eqs. (17), (21) and (22) are rewritten as

@u0j @xj

¼ 0;

ð24Þ 



@u0i @u0 @u0 @ 2 u0 1 @p  q0 gbhTih þ m 2i þ u0j i þ huj ih i ¼  @t @xj @xj q0 @xi @xj þ gbT 0 di3  u03 Sh di1 ; 0

ð15Þ

ð20Þ

0

0

ð25Þ 2 0

@T @T @T @ T þ huj ih ¼ jT 2  u03 s; þ u0j @t @xj @xj @xj

ð26Þ

1791

S. Hirabayashi, T. Sato / Computers & Fluids 39 (2010) 1789–1795

where the following relations are used:

xM;i ¼ xM;i L1 ;

ð37Þ

@hui ih ¼ huj ih Sh di1 dj3 ¼ hu3 ih di1 ¼ /i ; @xj @hTih ¼ huj ih sdj3 ¼ hu3 ih ¼ 0; huj ih @xj

u0i

¼ u0i ð 0 LÞ1=3 ; t ¼ tð 0 L2 Þ1=3 ;

ð38Þ

huj ih

ð27Þ ð28Þ

where /i is the zero vector. The time derivative of h uiih is removed in Eq. (25) because of the assumption of constant shear in stationary turbulence in a sufficiently long time period. Following Gerz et al. [1], we introduce the coordinate system, xM,i, which moves with the mean shear component of velocity, huiih, and the fourth term on the LHS of Eq. (25) and the third term on the LHS of Eq. (26) are eliminated. Then, the Eqs. (24)–(26) are rewritten in the new coordinate system, xM,i, as

@u0j @xM;j

e

e

p



p ¼

q0

ð29Þ

  @u0i @u0i @ 2 u0 1 @p þ u0j ¼  q0 gbhTih þ m 2 i @t @xM;j q0 @xM;i @xM;j þ gbT 0 di3  u03 Sh di1 ; @T 0 @T 0 @2T 0 ¼ jT 2  u03 s: þ u0j @t @xM;j @xM;j

ð30Þ ð31Þ

ð32Þ

ð41Þ

Sh ¼ Sh ðe0 L2 Þ1=3 ;

ð42Þ

e e ¼ ; e0

ð43Þ

where the upper asterisk denotes the nondimensionalization. Sh corresponds to the dimensionless shear number. Here Eqs. (29)–(31) are nondimensionalized:

Eqs. (29)–(31) are equivalent to the conventional governing equations of homogeneous shear turbulence under density stratification (e.g.[1–4,8]). While all those studies coped with developing or decaying turbulence under given shears and stratifications, the main focus of the present study is to control the shear dynamically in order to maintain stationary turbulence. Similar to Eq. (11), one can obtain the transport equation of energy by taking the inner product with u0i and the volumetric average to Eq. (30):

ð33Þ

where B is the buoyancy flux expressed by

B ¼ gbT 0 u03 ;

ð34Þ

and P is changed from that in Eq. (14) to

P ¼ u01 u03 Sh :

ð35Þ

When the flow field is under the equilibrium condition, P, e, and B balance. Therefore, Sh is determined by

Sh ¼ 

hB þ ei : hu01 u03 i

¼ 0;

@xM;j

ð44Þ

ð36Þ

By Eq. (36), the mean vertical shear, Sh, which is required to maintain turbulence, can be dynamically obtained at every time step in a numerical calculation with a given constant energy dissipation rate. 2.4. Nondimensionalization By introducing representative length, L, velocity, (e0L)1/3, and temperature, sL, primitive variables are nondimensionalized as follows:

! þ

1 @ 2 u0i ReL @x2 M;j

þ RiL T 0 di3  u03 Sh di1 ; 0

0

ð45Þ 2 0

@T @T 1 @ T þ u0j ¼  u03 ; @t @xM;j PrReL @x2 M;j

ð46Þ

where ReL, RiL, and Pr, are the Reynolds, Richardson, and Prandtl numbers, respectively, defined by

ReL ¼ RiL ¼

2.3. External force term

@hEi ¼ hei  hBi þ hPi; @t

ð40Þ

;

0 @u0i @p 0 @ui ¼  RiL hT  ih  þ uj  @t @xM;j @xM;i

The last term on the RHS of Eq. (30) represents the driving force by the mean shear:

fi ¼ u03 Sh di1 :

ðe0 LÞ

ð39Þ

T 0 ¼ T 0 ðsLÞ1 ;

@u0j ¼ 0;

2=3

Pr ¼

4=3 e1=3 0 L ; m

N2 L4=3

e2=3 0

;

m ; jT

ð47Þ ð48Þ ð49Þ

where N is the Brunt–Väisälä frequency:

N ¼ ðgbsÞ1=2 :

ð50Þ

2.5. Numerical methods, boundary and initial conditions Eqs. (44)–(46) were solved numerically in a time-marching way in a finite difference method. Pressure was solved by means of Marker and Cell method developed by Harlow and Welch [9], where the Poisson equation for pressure is solved iteratively to satisfy the continuity equation. The primitive variables were in staggered allocation, in which the velocity components are defined on the cell face and the pressure and the temperature are at the cell center. For spatial differentiation, a 4th-order energy- and momentum-conservative scheme [10] was applied to the convection term of the NS equation and a 2nd-order central difference scheme was adopted for the other terms. For time accuracy, the 3rd-order Adams–Bashforth method was used. The domain used in the simulation was a cubic, the edge length of which was 1 nondimensionalized by its own length, L. A schematic drawing of the simulation domain was shown in Fig. 1. As was mentioned before, Eqs. (44)–(46) do not include the effect of convection of the mean shear explicitly and they are expressed in a Lagrangian coordinate system, in which the simulation domain deforms in the direction of shear with time as was proposed by Rogallo [11]. Since it is easy to perform the simulation in a rectangular domain, an interpolation was applied to the intermediate velocity at each time step in order to remove the deformation effect as Gerz et al. [1] suggested: im uip i ðxM;i Þ ¼ ui ðxM;i  DtSh x3 di1 Þ;

ð51Þ

1792

S. Hirabayashi, T. Sato / Computers & Fluids 39 (2010) 1789–1795

LK ¼ Shear-periodic B. C.

*

0.5

x2 *

-0.5 1

0

x1

im where uip is the intermediate i is the interpolated velocity and ui velocity, which is described by

#

@ sij @u0i 1 @ 2 u0i   þ þ fi ;  @xM;j @xM;j ReL @x2 M;j

ð52Þ

where Dt is the interval of time steps. To obtain uip i accurately, the 2-dimensional Fourier interpolation was applied to each horizontal plane. Based on the assumption of homogeneity, fluctuating components were regarded as periodic in the simulation domain. On the horizontal boundaries, normal periodic boundary condition was adopted. For the top and bottom boundaries, the shear-periodic boundary condition of Schumann [12] was adopted in order to take the effect of the deformation into consideration. The initial condition for ui was given by the superposition of waves, the amplitude and the direction of which were determined randomly to satisfy the following equations: 5=3

EðkÞ ¼ 1:6e2=3 0 k

ð56Þ

:

3. Results and discussions

Fig. 1. Schematic diagram of simulation domain and boundary conditions.

0 uim i ¼ ui þ Dt uj

N3

ð55Þ

0.5 * -0.5 0

"

m ; e0  1=2 e0

The size of D is about twice as large as LK and it is expected that the grid size is adequate for DNS. As shown in Table 2, the fact that the grid size is sufficiently smaller than the Ozmidov scale also supports the capability to express the inertial subrange appropriately in the present DNS.

0

. B. C odic Peri

B. C.

LO ¼

x3

Pe rio dic

0.5

 3 1=4

ð53Þ

;

kj xj ¼ 0:

ð54Þ

3.1. Energy spectra Fig. 2a shows the instantaneous energy spectrum of simulated flow field at t* = 20.8, when the turbulence is fully developed, as shown later in the time series of turbulence kinetic energy in Fig. 3. It is found that the spectrum follows the slope of k5/3 in 4 6 k* 6 7 and it can be said that the inertial subrange was properly expressed in the present DNS. The largest scale in the inertial subrange was corresponding to the Ozmidov wavenumber,  kO ¼ L1 ¼ 4:2. This is reasonable considering that the Ozmidov O scale is a threshold between the buoyancy-dominating range by large-scale motions and the turbulence-dominating range by small-scale motions (e.g. [13]). The smallest scale resolved in the present simulation seems to be within the dissipation range and the spatial resolution is considered to be sufficient. Fig. 2b is the shear spectra calculated by using the results in Fig. 2a. It is noted that the contribution of high-wavenumber components to the total energy dissipation is small compared with that of inertial subrange.

(a)

*-1

*-1

LO

100

LK

10-1

The other components were set to be 0 at the initial state.

-2

Table 1 shows the constant parameters used in this study. The total grid numbers used in the simulation was 1283. Time step interval, Dt, was set to be 2.5  104. Table 2 shows the conditions where the Kolmogorov scale, LK, and the Ozmidov scale, LO, were calculated by

E* (k*)

10

2.6. Simulation case

10-3 -4

10

-5

10

-6

10

0

Unit

Value

m

m2 s1 m s2 kg m3 – m m3 s2 s1 – – –

1  106 9.817 103 1.0 1.5 1  109 2.0  103 1700 6.9 250

g

q0 Pr L

e0 N ReL RiL ReB

Table 2 Simulation conditions for stationary, homogeneous shear turbulence in stratification. Grids

D*

D/LK

LO/D

1283

7.8  103

2.1

30

2

10

10 *

k

(b) 10

2

*-1

*-1

LO

LK

101

k*2 E* (k*)

Variable

1

10

Table 1 Simulation constants for stationary, homogeneous shear turbulence in stratification.

100 10

-1

10

-2

10

-3

100

101

102 *

k

Fig. 2. (a) Wavenumber-based instantaneous spectrum of kinetic energy at t* = 20.8 and (b) shear spectra obtained from (a).

1793

S. Hirabayashi, T. Sato / Computers & Fluids 39 (2010) 1789–1795

(a)

Fig. 3a shows the time series of grid-scale kinetic energies, E ¼ hui ui i=2. It is found that the kinetic energy was maintained without further development or decay when t* > 2. Fig. 3b shows the time series of energy dissipation rates. The energy dissipation rate is calculated by

  @ui @ui @uj : þ @xM;j @xM;j @xM;i

0.1

ð57Þ

It is shown that the constant energy dissipation can be achieved under the given e0 in the present forcing method. Fig. 4 shows the time series of (a) the gradient Richardson number, Rg, and (b) the flux Richardson number, Rf, respectively, where Rg and Rf are defined by

Rg ¼

N2

0.2

0

0

5

10

t

(b)

20

15

20

0.3

ð58Þ

;

S2h hBi : Rf ¼ hPi

0.2

ð59Þ

Although the former is the ratio of forces and the latter is the ratio of energies, the both show a similar trend. Since the simulation was performed under the condition of energy equilibrium, the timeaverage of Rg in the range of stationary turbulence (t* > 8) is equivalent to the stationary Richardson number, RS.

0.1

0

0

5

10

3.3. Turbulence statistics

t

The temporal and volumetric averages of turbulence statistics are listed in Table 3 where the time average was taken over 8 < t* < 24 for all the parameters, such as the kinetic energy, energy dissipation rate, production rate, buoyancy flux, the Reynolds

(a)

15 *

Rf

e¼m

0.3

Rg

3.2. Temporal variation

1 0.8

*

Fig. 4. Time series of (a) the gradient Richardson number and (b) the flux Richardson number.

Table 3 Time- and volumetric averages of turbulence statistics of simulated flow field; kinetic energy, energy dissipation and production rates, buoyancy flux, Reynolds number based on the Taylor microscale, and gradient and flux Richardson numbers. Temporal average was taken over 8 < t* < 24. Rek

hE*i

h e*i

hP*i

h B*i

Rg

Rf

C

50

0.46

0.96

1.12

0.12

0.088

0.10

0.12

E*

0.6 0.4

number based on the Taylor microscale, and the gradient and the flux Richardson numbers. The Reynolds number based on the Taylor microscale was about 50, which is defined by

0.2 0

0

5

10

15

t

20

*

(b) 1.5

k ¼ ½15mq2 =e1=2 :

ð60Þ ð61Þ ð62Þ

where q and k are the representative velocity and the Taylor microscale, respectively. The temporal average of e* in the range of 8 < t* < 24 was 0.96 and this indicates that the error from the given dissipation rate was within 4% in the present DNS. Rg = 0.093 is comparable with RS = 0.088 in the results of Holt et al. [2], which is reasonable considering that Rek = 52 in their study was similar to that of the our present DNS.

ε*

1

0.5

0

qk Rek ¼ ; m  1=2 2E ; q¼ 3

4. Diffusion simulation of passive scalar 0

5

10

15

20

t* Fig. 3. Time series of volumetric-averaged (a) kinetic energy and (b) energy dissipation rate, he*i.

4.1. Governing equations Diffusion of a passive scalar in stationary, homogeneous, and stratified turbulence was numerically simulated to evidence the

1794

S. Hirabayashi, T. Sato / Computers & Fluids 39 (2010) 1789–1795

applicability of the present forcing method. A massless dye was adopted as a passive scalar in this study. The transport equation of the dye concentration is solved by

@C 1 @2C 0 @C ¼ ;  þ uj  @t @xM;j ReL Sc @x2 M;j

ð63Þ

where C is the nondimensional concentration of the dye and Sc is the Schmidt number defined by

Sc ¼

x3

m ; jC

original domain x3 = 0

ð64Þ

where jC is the molecular diffusivity for the dye. In this study, Sc was set to be 1. Fig. 6. Schematic diagram of initial dye distribution.

4.2. Boundary conditions and dummy domain Simulation conditions are the same as those adopted in Section 3. For the dye, periodic boundary conditions were imposed in the horizontal directions. In the vertical direction, however, the periodic condition cannot be applied since the concentration should diffuse infinitely. In order to pursue this requirement, dummy domains were introduced. Fig. 5 is the schematic diagram of the allocation of the dummy domains. In the dummy domains, flow fields are the same as that in the original domain locating at the vertical center. Because the Ozmidov scale is smaller than the size of the original domain, the eddies that are largest in the vertical direction can be accommodated in the domain. The dye, on the other hand, diffuses across the interfaces between the domains. The number of the total domains for the dye diffusion was 2n + 1 including the original domain. At the upper boundary of the top domain and the lower boundary of the bottom domain, a Neumann-type boundary condition of zero-gradient was adopted. From several test runs, it was confirmed that n = 8 was

adequate for having little effect on the distribution of the dye concentration in the time period of the present interest. 4.3. Initial conditions The initial condition for the dye was given by the Gaussian distribution in the vertical direction:

  x2 C ¼ exp  3 2 ; 2r 0

ð65Þ

where r0 is the initial standard deviation which was set to be 0.5, and x3 ¼ 0 corresponds to the vertical center of the original domain. The initial distribution of the dye concentration is illustrated in Fig. 6. Fully developed flow field as the result in Section 3 at t* = 8 was used at the start of the dye diffusion. 4.4. Distribution of dye patch

C

Fig. 7 shows the contours of dye concentration on a vertical cross-section of the vertically central domain at t* = 16 after the start of diffusion simulation. It is found that the dye was mixed up vertically by small eddies. Fig. 8 shows the temporal change of the vertical distributions of the dye concentration averaged in

dummy domain Pn

. . .

dummy domain P1 eddies original domain

distribution of dye

dummy domain M1

. . .

dummy domain Mn

Fig. 5. Schematic diagrams of dummy domains for dye diffusion simulation.

Fig. 7. Distributions of dye concentration in a vertical cross-section (y = 0) in the original domain at t* = 16 after start of diffusion calculation of dye.

S. Hirabayashi, T. Sato / Computers & Fluids 39 (2010) 1789–1795

1

Vertical diffusivity for an active scaler, such as heat, Kq, is modelled by Osborn [14] for a patch of stratified stationary turbulence:

t* = 0.5 * t =8 * t = 16

0.8

Kq ¼ C

0.6

C 0.2 0 -4

-2

0

2

A novel forcing method was developed to generate anisotropic turbulence in thermal stratification under the condition of energy equilibrium. This linear forcing method was implemented to DNS to simulate stationary homongeneous turbulence with arbitrary Reynolds and Richardson numbers. Turbulence statistics calculated from the generated flow field showed that the turbulence is properly maintained under the given conditions of energy dissipation rate and stratification. As an application of the present method, the vertical diffusivity of a passive scalar was estimated by simulating the diffusion of a massless dye in the DNS of stationary, homogeneous stratified turbulence with a given energy dissipation rate. The estimated diffusivity for the passive scalar was well comparable with that for heat obtained by a conventional model.

Fig. 8. Vertical distributions of horizontal-mean dye concentration.

1 0.8

σ∗2

0.6 0.4 0.2

0

2

4

6

8

t

10

12

14

16

the horizontal plane. It is seen that the dye diffusion keeps following the Gaussian distribution well while its variance increases with time. 4.5. Diffusivity of dye concentration The diffusivity for the dye concentration can be estimated from the change of its distribution with time. In order to quantify the size of dye patch, the variance of its distribution was calculated by the following formula:

x3 hCih

;

ð66Þ

where r is the standard deviation of the dye distribution. Fig. 9 is the time variations of the variance, r2. It is noted that the variance increases almost linearly with time. The diffusion coefficient, KC, in the vertical direction is then calculated by

KC ¼

1 dr2 : 2 dt

References

*

Fig. 9. Time variations of the variance of the dye distribution. The gradient of the dashed line was obtained by the least square method to fit the simulation results.

2 x x3 hCih

ð68Þ

5. Concluding remarks

x

P

;

4

* 3

r2 ¼ P3

e N2

where C = B/e is the mixing efficiency. From the result of DNS, the time average of C in 0 < t* < 16 was 0.12 and the resultant K q becomes 1.71  102, which is almost identical to K C obtained by the present DNS, 1.67  102. This does not contradict the result of Lindborg and Fedina [15], in which turbulent flux of passive scalar is the same as that of active scalar.

0.4

0

1795

ð67Þ

By means of the least square method, K C was calculated as the half of the slope of variance as 1.67  102, as shown by the dashed line in Fig. 9.

[1] Gerz T, Schumann U, Elghobashi SE. Direct numerical simulation of stratified homogeneous turbulent shear flows. J. Fluid Mech. 1989;200:563–94. [2] Holt SE, Koseff JR, Ferziger JH. A numerical study of the evolution and structure of homogeneous stably stratified sheared turbulence. J. Fluid Mech. 1992;237:499–539. [3] Shih LH, Koseff JR, Ferziger JH, Rehmann CR. Scaling and parameterization of stratified homogeneous turbulent shear flow. J. Fluid Mech. 2000;412:1–20. [4] Kaltenbach H-J, Gerz T, Schumann U. Large-eddy simulation of homogeneous turbulence and diffusion in stably stratified shear flow. J. Fluid Mech. 1994;280:1–40. [5] Sato T, Sato K, Jeong S-M. Numerical reproduction of turbulent flow field in a small patch of open ocean by low-wavenumber forcing. Comput. Fluids 2007;36:540–8. [6] T.S. Lundgren, Linearly forced isotropic turbulence, in: Annual Research Briefs, Center for Turbulence Research, Stanford, 2003, pp. 461–473. [7] Rosales C, Meneveau C. Linear forcing in numerical simulations of isotropic turbulence: physical space implementations and convergence properties. Phys. Fluids 2005;17. doi:10.1063/1.2047568. [8] Itsweire EC, Koseff JR, Briggs DA, Ferziger JH. Turbulence in stratified shear flows: implications for interpreting shear-induced mixing in the ocean. J. Phys. Oceanogr. 1993;23:1508–22. [9] Harlow FH, Welch JE. Numerical calculation of time dependent viscous incompressible flow of fluid with free surface. Phys. Fluids 1965;8:2182–9. [10] Morinishi Y, Luid TS, Vasilyev OV, Moin P. Fully conservative higher order finite difference schemes for incompressible flow. J. Comp. Phys. 1998;143:90–124. [11] Rogallo RS. Numerical experiments in homogeneous turbulence. NASA Tech. Memorandum 1981:81315. [12] U. Schumann, Algorithms for direct numerical simulation of shear-periodic turbulence, in: Proceedings of 9th International Conference on Numerical Methods in Fluid Dynamics. Lecture Notes in Physics, vol. 218, 1985, pp. 492– 496. [13] Carnevale GF, Briscolini M, Orland P. Buoyancy- to inertial- range transition in forced stratified turbulence. J. Fluid Mech. 2001;427:205–39. [14] Osborn TR. Estimates of the local rate of vertical diffusion from dissipation measurements. J. Phys. Oceanogr. 1980;10:83–9. [15] Lindborg E, Fedina E. Vertical turbulent diffusion in stably stratified flows. Geophys. Res. Lett. 2009;36:L01605. doi:10.1029/2008GL036437.