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.
INTRODUCTION
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
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
According to Kraichnan, a sufficient condition for the divergence to be zero is
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
and this operation returns
where
It is clear that, unless (Δx, Δy, Δz) → 0, . This means that vn ⋅ kn≠0 and wn ⋅ kn≠0 and therefore the discrete velocity divergence, ∇d ⋅ u, 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
rather than Eq. (3). Therefore, one can set and 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, 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 is
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.
IMPACT OF BOUNDARY CONDITIONS
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.
MODIFICATION OF KINETIC ENERGY CONTENT VIA PRESSURE PROJECTION
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
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
A straightforward way to evaluate ΔE is to use Fourier series. Let
where , , and are complex Fourier coefficients. Also note that
This gives us the following solution for the Poisson equation:
Finally, we have
which implies that
Now by virtue of Parseval’s identity, we have
where a starred quantity denotes a complex conjugate. Upon careful substitutions and evaluations, we have the following:
where , , and n are unit vectors such that , , and . The symbol “•” denotes scalar multiplication.
The modal energy loss, at wavenumber k, can be easily calculated as
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 . Recognizing that leads to
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.
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 vortex . | Isotropic turbulence . | ||||
---|---|---|---|---|---|---|
Resolution . | LHS of (17) . | RHS of (17) . | Relative error . | LHS 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 vortex . | Isotropic turbulence . | ||||
---|---|---|---|---|---|---|
Resolution . | LHS of (17) . | RHS of (17) . | Relative error . | LHS 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 |
Acknowledgments
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.”