This comment aims at addressing a mass conservation issue in a paper published in the physics of fluids. The paper [R. H. Kraichnan, “Diffusion by a random velocity field,” Phys. Fluids 13(1), 22 (1970)] introduces a novel method to generate synthetic isotropic turbulence for computational purposes. The method has been used in the literature to generate inlet boundary conditions and to model aeroacoustic noise as well as for validation and verification purposes. However, the technique uses a continuous formulation to derive the mass conservation constraint. In this comment, we argue that the continuous constraint is invalid on a discrete grid and provide an alternative derivation using the discrete divergence. In addition, we present an analysis to quantify the impact of a pressure projection on the kinetic energy of a non-solenoidal velocity field.

In 1970, Kraichnan published a seminal paper entitled “Diffusion by a Random Velocity Field.”1 Among the many foundational contributions found in that paper is the definition of a method to generate a synthetic, statistically isotropic turbulent velocity field. Kraichnan’s turbulence generator was used in many subsequent studies to model aeroacoustic noise,2–4 inlet boundary conditions,5–10 verification and validation of turbulence solvers,11,12 and a priori and a posteriori testing of turbulence models.13,14

To enforce a divergence-free velocity field, Kraichnan’s method requires that two vectors be perpendicular. In this comment we show that this condition is only valid in the infinitesimal limit and should be replaced with an otherwise simple modification based on the discrete divergence of the velocity field.

The essence of Kraichnan’s technique is to generate a velocity field in the form

u(x)=n=1Nvncos(knx+ωnt)+n=1Nwnsin(knx+ωnt),
(1)

where t is time, ωn is a frequency, x ≡ (x, y, z) is a spatial coordinate, kn ≡ (kx,n, ky,n, kz,n) is a wave vector, and vn and wn are amplitudes to be determined using the velocity divergence constraint.

The (spatially continuous) divergence of this velocity field is

u(x)=n=1Nvnknsin(knx+ωnt)+n=1Nwnkncos(knx+ωnt).
(2)

According to Kraichnan, a sufficient condition for the divergence to be zero is

vnkn=0,wnkn=0
(3)

and hence, one sets vn = ζn × kn and wn = ξn × kn, where ζn and ξn are randomly chosen from a Gaussian distribution. However, when used to initialize a discrete velocity field, (3) does not necessarily guarantee a divergence-free velocity field. Indeed, for example, on a staggered structured uniform grid with grid spacing Δx ≡ (Δx, Δy, Δz), the discrete velocity divergence is

du(x)=u(x+12Δx)u(x12Δx)Δx,
(4)

and this operation returns

du(x)=n=1Nvnk̃nsin(knx+ωnt)+n=1Nwnk̃ncos(knx+ωnt),
(5)

where

k̃n=2Δxsin(12kx,nΔx),2Δysin(12ky,nΔy),2Δzsin(12kz,nΔz).
(6)

It is clear that, unless (Δx, Δy, Δz) → 0, k̃nkn. This means that vnkn≠0 and wnkn≠0 and therefore the discrete velocity divergence, ∇du, is not zero. Bailly and Juvé allude to this defect in Ref. 4 and attribute it to inherent inhomogeneity in the generated turbulence. The fact is that the solenoidal condition derived by Kraichnan (and others) is only valid in the infinitesimal limit.

A simple remedy can be used to fix this. Rather than using the spatially continuous divergence condition to compute vn and wn, one should instead use the discrete divergence constraint (7). Again, using a staggered grid arrangement on a uniform structured grid, a sufficient condition for a (discrete) divergence-free velocity field is

vnk̃n=0;wnk̃n=0
(7)

rather than Eq. (3). Therefore, one can set vn=ζn×k̃n and wn=ξn×k̃n where, again, ζn and ξn are randomly sampled from a Gaussian (or other) distribution. It is also noteworthy that, in the limit of a vanishingly small grid size, k̃n approaches the value of kn thus recovering the original continuous divergence-free requirement.

It is possible to derive this solenoidal condition for other grid arrangements as well as higher-order discretization. For example, for a collocated grid, where the velocity field is stored at cell centers, the solenoidal wave vector k̃m is

k̃n=1Δxsin(kx,nΔx),1Δysin(ky,nΔy),1Δzsin(kz,nΔz).
(8)

An open-source, isotropic turbulence generator based on Eqs. (1) and (7) has been developed in Ref. 15 and can be found online at http://turbulence.utah.edu.

When periodic or specified boundary conditions are used in conjunction with Eqs. (1) and (7), then one will lose the divergence-free constraint in the cells near the boundaries. However, using the discrete constraint, the number of cells with divergence is only a fraction of those resulting from using the continuous constraint.

For example, under periodic conditions on a staggered grid with N cells in each direction, the continuous constraint will produce N3 cells with divergence while the discrete constraint results in only 3N2 cells—a 3/N fraction of cells with divergence.

For an incompressible flow, the impact of a divergence-free initial condition can be quantified by evaluating the impact of the mass correction scheme on the kinetic energy. Given a diverging velocity vector field q such that ∇ ⋅ q≠0, define the mass conserving velocity field via the projection

u=qp
(9)

such that ∇ ⋅ u = 0. This implies that ∇2p = ∇ ⋅ q. What is the impact of the projection on the kinetic energy of the system? Mathematically, we are interested in calculating

ΔE=Vu2dVVq2dV.
(10)

A straightforward way to evaluate ΔE is to use Fourier series. Let

q=qˆeikx,u=uˆeikx,p=pˆeikx
(11)

where qˆ, uˆ, and pˆ are complex Fourier coefficients. Also note that

p=ipˆkeikx,2p=pˆk2eikx,q=iqˆkeikx.
(12)

This gives us the following solution for the Poisson equation:

pˆ=iqˆkk2.
(13)

Finally, we have

uˆeikx=qˆeikxipˆkeikx
(14)

which implies that

uˆ=qˆipˆk.
(15)

Now by virtue of Parseval’s identity, we have

Vu2dV=uˆ2=uˆuˆ,
(16)

where a starred quantity denotes a complex conjugate. Upon careful substitutions and evaluations, we have the following:

ΔE=Vu2dVVq2dV=qˆkqˆkk2=qˆ2eˆneˆn,
(17)

where eˆ, eˆ, and n are unit vectors such that qˆ=qˆeˆ, qˆ=qˆeˆ, and k=kn. The symbol “•” denotes scalar multiplication.

The modal energy loss, at wavenumber k, can be easily calculated as

uˆ2qˆ2=qˆkqˆkk2=qˆ2eˆneˆn.
(18)

This tells us that the energy deficit, at the kth-wavenumber, between the uncorrected and correct velocity fields is independent of the wavenumber but is proportional to the magnitude of the uncorrected velocity.

Another way to look at this is to divide (18) by qˆ2. Recognizing that 1eˆneˆn1 leads to

0uˆ2qˆ2=1eˆneˆn2.
(19)

This gives us a bound on the modal energy at wave number k.

By way of validation, (17) has been discretely computed for a number of cases with a diverging initial condition. The first set of cases consisted of a single projection for a perturbed three dimensional Taylor-Green vortex on a domain of size L = 2π. The perturbation was intended to create a diverging initial condition while maintaining the periodic properties of the Taylor-Green vortex. The second set of cases consisted of using Kraichnan’s turbulence generator with divergence in the initial condition fitted to produce the von Kármán-Pao spectrum4 on a domain size of length L = 18π/100. Three resolutions were chosen: 323, 643, and 1283. The results for this validation are shown in Table I. We found that the relative error for the Taylor-Green vortex is several orders of magnitude smaller than that produced by the isotropic turbulence data. We believe that this is due to the periodicity of the signals and the role that it plays in making sure derivatives of Fourier series are consistent with spatial derivatives in real space. Overall, the results are quite satisfactory. The validation could also be further improved by paying close attention to correcting terms in the Fourier series that correspond to spatial derivatives.

TABLE I.

Validation of (17) for velocity fields corresponding to the 3D Taylor-Green vortex and isotropic turbulence. The initial conditions for both cases had non-zero divergence.

Taylor-Green vortexIsotropic turbulence
ResolutionLHS of (17)RHS of (17)Relative errorLHS of (17)RHS of (17)Relative error
32 −1.974 × 10−4 −1.974 × 10−4 2.122 × 10−13 −4.358 × 10−3 −4.371 × 10−3 3.033 × 10−3 
64 −1.974 × 10−4 −1.974 × 10−4 4.913 × 10−13 −4.852 × 10−3 −4.916 × 10−3 1.314 × 10−2 
128 −1.974 × 10−4 −1.974 × 10−4 2.137 × 10−13 −4.449 × 10−3 −4.443 × 10−3 1.458 × 10−3 
Taylor-Green vortexIsotropic turbulence
ResolutionLHS of (17)RHS of (17)Relative errorLHS of (17)RHS of (17)Relative error
32 −1.974 × 10−4 −1.974 × 10−4 2.122 × 10−13 −4.358 × 10−3 −4.371 × 10−3 3.033 × 10−3 
64 −1.974 × 10−4 −1.974 × 10−4 4.913 × 10−13 −4.852 × 10−3 −4.916 × 10−3 1.314 × 10−2 
128 −1.974 × 10−4 −1.974 × 10−4 2.137 × 10−13 −4.449 × 10−3 −4.443 × 10−3 1.458 × 10−3 

The authors gratefully acknowledge support from the National Science Foundation PetaApps Award No. 0904631 and XPS Award No. 1337145 as well as the Department of Energy Award No. DE-SC0008998. The authors also wish to thank the anonymous reviewers for their insightful remarks and for their suggestions on the analysis presented in the section titled “Modification of kinetic energy content via pressure projection.”

1.
R. H.
Kraichnan
, “
Diffusion by a random velocity field
,”
Phys. Fluids
13
(
1
),
22
(
1970
).
2.
M.
Karweit
,
P.
Blanc-Benon
,
D.
Juvé
, and
G.
Comte-Bellot
, “
Simulation of the propagation of an acoustic wave through a turbulent velocity field: A study of phase variance
,”
J. Acoust. Soc. Am.
89
(
1
),
52
62
(
1991
).
3.
W.
Bechara
,
C.
Bailly
,
P.
Lafon
, and
S. M.
Candel
, “
Stochastic approach to noise modeling for free turbulent flows
,”
AIAA J.
32
(
3
),
455
463
(
1994
).
4.
C.
Bailly
and
D.
Juve
, “
A stochastic approach to compute subsonic noise using linearized Euler’s equations
,” presented at the
5th AIAA/CEAS Aeroacoustics Conference and Exhibit
,
Bellevue, WA
,
1999
.
5.
M.
Billson
,
Computational Techniques for Turbulence Generated Noise
(
Chalmers University of Technology
,
2004
).
6.
M.
Billson
,
L.-E.
Eriksson
, and
L.
Davidson
, “Jet noise modeling using synthetic anisotropic turbulence,” AIAA Paper 2004-3028, 2004.
7.
L.
Davidson
, “
Hybrid les-rans: Inlet boundary conditions for flows including recirculation
,” in
Proceedings of the 5th International Symposium on Turbulence and Shear Flow Phenomena
(
Citeseer
,
Munich, Germany
,
2007
).
8.
L.
Davidson
, “
Using isotropic synthetic fluctuations as inlet boundary conditions for unsteady simulations
,”
Adv. Appl. Fluid Mech.
1
(
1
),
1
35
(
2007
), available online at http://www.pphmj.com/abstract/2288.htm.
9.
L.
Davidson
, “
Hybrid LES-RANS: Inlet boundary conditions for flows with recirculation
,” in
Advances in Hybrid RANS-LES Modelling
(
Springer
,
2008
), pp.
55
66
.
10.
L.
Davidson
and
S.-H.
Peng
, “
Embedded large-eddy simulation using the partially averaged Navier–Stokes model
,”
AIAA J.
51
(
5
),
1066
1079
(
2013
).
11.
S. M. D. B.
Kops
, “
Numerical simulation of non-premixed turbulent combustion
,” Ph.D. thesis,
University of Washington
,
1999
.
12.
R. J.
McDermott
, “
Toward one-dimensional turbulence subgrid closure for large-eddy simulation
,” Ph.D. thesis,
The University of Utah
,
2005
.
13.
A.
Keating
,
U.
Piomelli
,
E.
Balaras
, and
H.-J.
Kaltenbach
, “
A priori and a posteriori tests of inflow conditions for large-eddy simulation
,”
Phys. Fluids
16
(
12
),
4696
4712
(
2004
).
14.
Y.
Li
and
C.
Rosales
, “
Constrained multi-scale turnover Lagrangian map for anisotropic synthetic turbulence: A priori tests
,”
Phys. Fluids
26
(
7
),
075102
(
2014
).
15.
T.
Saad
,
D.
Cline
,
R.
Stoll
, and
J. C.
Sutherland
, “
Scalable tools for generating synthetic isotropic turbulence with arbitrary spectra
,”
AIAA J.
(published online).