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,

**k**

_{n}≡ (

*k*

_{x,n},

*k*

_{y,n},

*k*

_{z,n}) is a wave vector, and

**v**

_{n}and

**w**

_{n}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 **v**_{n} = *ζ*_{n} × **k**_{n} and **w**_{n} = *ξ*_{n} × **k**_{n}, 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, $k\u0303n\u2260kn$. This means that **v**_{n} ⋅ **k**_{n}≠0 and **w**_{n} ⋅ **k**_{n}≠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 **v**_{n} and **w**_{n}, 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 $vn=\zeta n\xd7k\u0303n$ and $wn=\xi n\xd7k\u0303n$ 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\u0303n$ approaches the value of **k**_{n} 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\u0303m$ 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 *N*^{3} cells with divergence while the discrete constraint results in only 3*N*^{2} 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 ∇^{2}*p* = ∇ ⋅ **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 $q\u02c6$, $u\u02c6$, and $p\u02c6$ 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 $e\u02c6$, $e\u02c6\u2217$, and **n** are unit vectors such that $q\u02c6=q\u02c6e\u02c6$, $q\u02c6\u2217=q\u02c6e\u02c6\u2217$, and $k=kn$. 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 *k*th-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\u02c62$. Recognizing that $\u22121\u2264e\u02c6\u22c5n\u2022e\u02c6\u2217\u22c5n\u22641$ 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 spectrum^{4} on a domain size of length *L* = 18*π*/100. Three resolutions were chosen: 32^{3}, 64^{3}, and 128^{3}. 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.

. | 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.”