# Triple Hill’s Vortex Synthetic Eddy Method

John S. Haywood email<EMAIL_ADDRESS>Mississippi State University,
MS 39762, USA  Adrian Sescu email<EMAIL_ADDRESS>Mississippi State
University, MS 39762, USA  Shanti Bhushan  Mississippi State University, MS
39762, USA  Chris Kees  US Army Engineer Research and Development Center
(ERDC), Vicksburg, MS

###### Abstract

The generation of initial or inflow synthetic turbulent velocity or scalar
fields reproducing statistical characteristics of realistic turbulence is
still a challenge. The synthetic eddy method, previously introduced in the
context of inflow conditions for large eddy simulations, is based on the
assumption that turbulence can be regarded as a superposition of coherent
structures. In this paper, a new type of synthetic eddy method is proposed,
where the fundamental eddy is constructed by superposing three Hill’s
vortices, with their axes orthogonal to each other. A distribution of Hill’s
vortices is used to synthesize an anisotropic turbulent velocity field that
satisfies the incompressibility condition and match a given Reynolds stress
tensor. The amplitudes of the three vortices that form the fundamental eddy
are calculated from known Reynolds stress profiles through a transformation
from the physical reference frame to the principal-axis reference frame. In
this way, divergence-free anisotropic turbulent velocity fields are obtained
that can reproduce a given Reynolds stress tensor. The model was tested on
both isotropic and anisotropic turbulent velocity fields, in the framework of
grid turbulence decay and turbulent channel flow, respectively. The transition
from artificial to realistic turbulence in the proximity to the inflow
boundary was found to be small in all test cases that were considered.

##  1 Introduction

There is a wide range of models for generating turbulent velocity fields at
the LES inflow boundary, roughly divided into several classes [41]: precursor
DNS or LES; turbulence recycling; synthetic turbulence; artificial forcing or
volume source terms; and vortex generating techniques. The most accurate
technique to impose turbulent inflow conditions is based on a precursor
simulation performed additionally in a separate flow domain with periodic
boundary conditions in the streamwise direction; data is extracted from a
cross-sectional plane and imposed at the inflow boundary of the main
simulation. This method is known to produce realistic turbulent structures at
the inflow boundary and is able to capture the correct flow physics in the
flow region of interest. However, this technique has two major drawbacks:
firstly, it is restricted to simple applications where the flow at the inflow
boundary can be regarded as a fully developed turbulent flow [20] or a
spatially developed turbulent boundary layer [30]. Secondly, this technique is
computationally expensive because two simulations have to be carried out.

Another choice is to employ recycling techniques by which the streamwise
periodic conditions of developing flows are modified to account for the
streamwise variations of the flow [30]. Spalart [45] developed such a method
to account for spatial growth in boundary layers, by using periodic boundary
conditions in the streamwise direction and source terms in the Navier-Stokes
equations to account for the growth. Spalart’s method is ideal to generate
inflow boundary conditions because it is capable of producing a spatially
evolving boundary layer. Methods based on proper orthogonal decomposition [12,
19] are less expensive than the other methods, but they require appropriate
DNS [19] or high-fidelity experimental data [12] to obtain the energetic
modes.

Synthetic turbulence modeling is an attractive and inexpensive approach that
is utilized to generate initial or inflow conditions for Large Eddy
Simulations (LES) or Direct Numerical Simulations (DNS), or to analyze
fundamental structures of turbulent flows. To be feasible, a synthetic
turbulence approach must ensure that specific statistics from the target
turbulent data are reproduced. A number of proposed synthetic turbulence
models [29, 2, 13] are based on spectral methods, which are mostly utilized to
generate isotropic turbulent flow-fields (the pioneering work of Kraichnan
[28] is worth mentioning here). The resulting turbulent kinetic energy in this
case is uniformly spread over all wavenumbers, and as a result the pseudo-
turbulence is rapidly dissipated, and a long distance is required to recover
the realistic turbulent flow. The model of Maxey [32] is a successful
application of Kraichnan’s method to applications involving anisotropic
turbulent flows. His technique is based on filtering and scaling operations
applied to the flow field. Smirnov et al. [44] extended the method proposed by
Maxey [32] by using scaling and simple coordinate transformation operations,
which were found to be efficient and fast. Based on the work of Smirnov et al.
[44], Batten et al. [1] reconstructed an inflow turbulent fluctuation flow
field by summing a set of Fourier modes with prescribed random phases and
amplitudes, and preserving the space/time correlations and the second order
moments. Keating et al. [21] used Batten’s model to generate inflow boundary
conditions for a plane channel flow and found that the development of the
turbulent flow was slow: after a distance of ten channel heights the turbulent
flow in the core region was clearly not yet fully developed. Billson et al.
[5] used the idea of scaling and coordinate transformation as proposed by
Smirnov et al. [44] to generate anisotropic synthesized velocity fields from
isotropic Reynolds stress tensors. Davidson and Billson [9] then applied this
method to fully-developed channel flow. Rosales and Meneveau [38, 39] proposed
a simple method, called multiscale minimal Lagrangian map (MMLM), to generate
synthetic vector fields, starting with an initial Gaussian field and using
deformed random-phase Fourier modes. Klein et al. [25] developed a synthetic
turbulence model that reproduced first and second order one-point statistics
and local correlations, based on three-dimensional digital filters. The
temporal correlations were well-predicted by shifting the random data in time.
Kempf et al. [22] showed that a method solving for an appropriate diffusion
equation using three-dimensional random data was equivalent to applying a
filter based on Gaussian functions. Using this idea, a new synthetic
turbulence model was devised which can be easily implemented in both
structured and unstructured grid codes. These previous methods approach the
modeling of realistic turbulence from a statistical point of view, modify
fields of random fluctuations as a means to match as many turbulent statistics
as possible. There is another class of method that is based in a physical
description of turbulence and tries to model the coherent structures present
in turbulence through the combination of discrete synthetic structures.

A synthetic eddy method (SEM) based on the assumption that turbulence can be
considered as a superposition of coherent structures was proposed by Jarrin et
al. [16]. The eddies were convected downstream with the mean flow, through the
inlet plane. Jarrin’s synthetic eddy method showed similar behavior as the
other existent models based on Fourier modes or digital filters, especially
for spatial decay of isotropic turbulence and turbulent channel flows. Pamiès
et al. [34] finely tuned the shape functions to mimic the elongated vortices
observed in the near-wall region. Roidl et al. [37] then tuned the shape
functions in the logarithmic and wake regions for a turbulent boundary layer.
Marked improvements in the reduction of the development length were seen as
more realistic modeling of the coherent structures was included in the imposed
synthetic eddies. In order to further reduce the development length, Poletto
et al. [35] applied the synthetic eddy framework to the vorticity field and
solved for the fluctuating velocity components using the Biot-Savart law. From
this solution process, restrictions on the shape functions were found that
were then used to select a function that satisfied the divergence-free
condition. This formulation only includes one shape function, based on a sine
function, and one characteristic eddy length scale, and as such, could only
match a limited range of anisotropic Reynolds stress tensors. To be able to
reproduce a wide range of anisotropic turbulence, Poletto et al. [36]
reformulated the method to allow for separate shape functions and lengths
scales for each dimension. This destroyed the satisfaction of the divergence-
free condition guaranteed by the previous formulation, but by taking the
divergence of the prospective eddy velocity field sufficient conditions were
found to be able to select divergence-free shape functions, based on a
quadratic polynomial. Poletto et al. [36] showed that pressure fluctuations at
the inlet, caused by the divergence-free condition not being satisfied, were
greatly reduced when compared with the original synthetic eddy method for
plane channel flow. Skillen et al. [43] introduced shape functions that are
scaled with the local synthetic eddy population density and also modified the
homogeneous eddy location distribution of the original SEM to allow for
clustering of smaller synthetic eddies near the wall. In the context of
computational aeroacoustics, a divergence-free synthetic eddy model was
constructed by Sescu and Hixon [40] based on the vector potential instead of
the velocity field.

In this work, the proposed Triple Hill’s Vortex (THV) synthetic eddy method
builds upon the foundation of the current SEM’s, but instead of using abstract
mathematical functions to determine the shape of eddies, the THV SEM uses a
superposition of Hill’s spherical vortices. Since the Hill’s spherical vortex
is an actual solution to the Euler equation, the synthetic inflow velocity
field will reproduce more realistic turbulent flow conditions. Specifically,
within the proposed THV model, three orthogonal Hill’s vortices are
superposed, with the amplitudes of the individual vortices determined by using
a transformation from the physical reference frame to the principal-axis
frame, as outlined in Billson et al. [5].

##  2 Governing Equations

The governing equations employed consist of the Favre-filtered full
compressible Navier-Stokes equations written in curvilinear coordinates and
conservative form. A generalized curvilinear coordinate transformation in the
three-dimensional form
$\xi=\xi\left(x,y,z\right),\eta=\eta\left(x,y,z\right),\zeta=\zeta\left(x,y,z\right)$,
is considered, where $\xi$, $\eta$, and $\zeta$ are the spatial coordinates in
the computational space, and $x$, $y$, and $z$ are the spatial coordinates in
physical space. As a note, a tilde ($\widetilde{\quad}$) represents Favre-
filtering at the grid level $\Delta$. In conservative form, the Navier-Stokes
equations are written as

$\displaystyle\mathbf{Q}_{t}+\mathbf{F}_{\xi}+\mathbf{G}_{\eta}+\mathbf{H}_{\zeta}=\mathbf{S}$
(1)

where the vector of conservative variables is given by

$\mathbf{Q}=\frac{1}{J}\\{\begin{array}[]{rrrrrr}\overline{\rho},\hskip
11.38109pt\overline{\rho}\widetilde{u}_{i},\hskip
11.38109pt\overline{\rho}\widetilde{E}\end{array}\\}^{T},\quad i=1,2,3$ (2)

$\overline{\rho}$ is the mean density of the fluid,
$\widetilde{u}_{i}=(\widetilde{u},\widetilde{v},\widetilde{w})$ is the
filtered velocity vector in physical space, and $\overline{\rho}\widetilde{E}$
is the total energy. The flux vectors, $\mathbf{F}$, $\mathbf{G}$, and
$\mathbf{H}$, are given by

$\displaystyle\mathbf{F}=\frac{1}{J}\left\\{\begin{array}[]{c}\overline{\rho}U,\hskip
11.38109pt\overline{\rho}\widetilde{u}_{i}U+\xi_{x_{i}}(\overline{p}+\tau_{i1}),\hskip
11.38109pt\overline{\rho}\widetilde{E}U+\overline{p}U+\xi_{x_{i}}\Theta_{i}\end{array}\right\\}^{T},$
(4)
$\displaystyle\mathbf{G}=\frac{1}{J}\left\\{\begin{array}[]{c}\overline{\rho}V,\hskip
11.38109pt\overline{\rho}\widetilde{u}_{i}V+\eta_{x_{i}}(\overline{p}+\tau_{i2}),\hskip
11.38109pt\overline{\rho}\widetilde{E}V+\overline{p}V+\eta_{x_{i}}\Theta_{i}\end{array}\right\\}^{T},$
(6)
$\displaystyle\mathbf{H}=\frac{1}{J}\left\\{\begin{array}[]{c}\overline{\rho}W,\hskip
11.38109pt\overline{\rho}\widetilde{u}_{i}W+\zeta_{x_{i}}(\overline{p}+\tau_{i3}),\hskip
11.38109pt\overline{\rho}\widetilde{E}W+\overline{p}W+\zeta_{x_{i}}\Theta_{i}\end{array}\right\\}^{T}$
(8)

where the contravariant velocity components are given by

$\displaystyle U=\xi_{x_{i}}\widetilde{u}_{i},\hskip
11.38109ptV=\eta_{x_{i}}\widetilde{u}_{i},\hskip
11.38109ptW=\zeta_{x_{i}}\widetilde{u}_{i}$ (9)

with the Einstein summation convention applied over $i=1,2,3$, the shear
stress tensor and the heat flux are given as

$\tau_{ij}=\frac{\widetilde{\mu}}{Re}\left[\left(\frac{\partial\xi_{k}}{\partial
x_{j}}\frac{\partial\widetilde{u}_{i}}{\partial\xi_{k}}+\frac{\partial\xi_{k}}{\partial
x_{i}}\frac{\partial\widetilde{u}_{j}}{\partial\xi_{k}}\right)-\frac{2}{3}\delta_{ij}\frac{\partial\xi_{l}}{\partial
x_{k}}\frac{\partial\widetilde{u}_{k}}{\partial\xi_{l}}\right]+\tau_{ij}^{sgs}$
(10)

$\Theta_{i}=\widetilde{u}_{j}\tau_{ij}+\frac{\widetilde{\mu}}{(\gamma-1)M_{\infty}^{2}RePr}\frac{\partial\xi_{l}}{\partial
x_{i}}\frac{\partial\widetilde{T}}{\partial\xi_{l}}+q_{i}^{sgs}$ (11)

respectively. $\tau_{ij}^{sgs}$ and $q_{i}^{sgs}$ are the subgrid scale stress
and subgrid scale heat flux terms, which are defined in Section 2.1. The
pressure $\overline{p}$, the temperature $\widetilde{T}$, and the density of
the fluid are combined in the equation of state,
$\overline{p}=\overline{\rho}\widetilde{T}/\gamma M_{\infty}^{2}$. Other
notations include the dynamic viscosity $\mu$, Reynold’s number
$Re=\rho_{\infty}V_{\infty}L/\mu$ based on a characteristic velocity
$V_{\infty}$, and a characteristic length $L$, the free-stream Mach number
$M_{\infty}=V_{\infty}/a$ (with $a$ being the speed of sound), Prandtl’s
number $Pr=C_{p}\mu/k$ (where $k$ is thermal conductivity), the specific heat
at constant pressure $C_{p}$, and the ratio between the specific heats
$\gamma$. The Jacobian of the curvilinear transformation from the physical
space to computational space is denoted by $J$. The derivatives $\xi_{x}$,
$\xi_{y}$, $\xi_{z}$, $\eta_{x}$, $\eta_{y}$, $\eta_{z}$, $\zeta_{x}$,
$\zeta_{y}$, and $\zeta_{z}$ represent grid metrics. The variables are non-
dimensionalized by their respective free-stream variables, except from
pressure which is non-dimensionalized by $\rho_{\infty}V_{\infty}$. The
dynamic viscosity and thermal conductivity are linked to the temperature using
the Sutherland’s equations in dimensionless form,

$\displaystyle\widetilde{\mu}=\widetilde{T}^{3/2}\frac{1+C_{1}/T_{\infty}}{\widetilde{T}+C_{1}/T_{\infty}};\hskip
11.38109pt\widetilde{k}=\widetilde{T}^{3/2}\frac{1+C_{2}/T_{\infty}}{\widetilde{T}+C_{2}/T_{\infty}},$
(12)

where for air at sea level, $C_{1}=110.4K$, $C_{2}=194K$, and $T_{\infty}$ is
a reference temperature.

###  2.1 Subgrid Scale Model

The SGS stress is modeled using the Coherent Structure Model (CSM) developed
by Kobayashi [26] for incompressible flow and later applied to compressible
flows by Hadjadj et al. [14] and Ben-Nasr et al. [3]. The CSM is based on the
assumption that the SGS dissipation is small at the center of a coherent eddy
and that the energy transfer between the resolved scales and the SGS occurs
around the edge of this coherent eddy [26, 27].

The SGS stress tensor,
$\tau_{ij}^{sgs}=\overline{\rho}\left(\widetilde{u_{i}u_{j}}-\widetilde{u_{i}}\widetilde{u_{j}}\right)$,
is defined using an eddy-viscosity model as as follows

$\tau_{ij}^{sgs}=-2\mu_{sgs}\left(\widetilde{S}_{ij}-\frac{1}{3}\widetilde{S}_{kk}\delta_{ij}\right)-\frac{1}{3}\tau_{kk}^{sgs}\delta_{ij}$
(13)

where $\widetilde{S}_{ij}$ is the resolved strain-rate tensor. The SGS
viscosity, $\mu_{sgs}$, is defined as

$\mu_{sgs}=\overline{\rho}C_{s}\Delta^{2}|\widetilde{S}|$ (14)

where $\Delta$ is the grid scale, $C_{s}$ is the dynamically calculated
Smagorinsky coefficient, and $|\widetilde{S}|$ is the strain-rate magnitude.
The isotropic part of the SGS tensor, $\tau_{kk}^{sgs}$, is modeled using the
following relationship proposed by Yoshisawa [48]

$\tau_{kk}^{sgs}=2\overline{\rho}C_{I}\Delta^{2}|\widetilde{S}|^{2}$ (15)

where $C_{I}$ is another dynamically calculated model coefficient.

For the CSM, the Smagorinsky coefficient is defined as

$C_{s}=C_{csm}|F_{cs}|^{3/2}\left(1-F_{cs}\right)$ (16)

where $C_{csm}$ is the CSM model coefficient, equal to $1/22$ [26, 27], and
$F_{cs}$ is the coherent structure function.

$F_{cs}=\frac{\widetilde{Q}}{\widetilde{E}}$ (17)

$\widetilde{Q}$ is the second invariant of the resolved velocity gradient and
$\widetilde{E}$ is the magnitude of the resolved velocity gradient tensor.

$\displaystyle\widetilde{Q}$
$\displaystyle=\frac{1}{2}\left(\widetilde{W}_{ij}\widetilde{W}_{ij}-\widetilde{S}_{ij}\widetilde{S}_{ij}\right)=$
$\displaystyle-\frac{1}{2}\frac{\partial\widetilde{u}_{j}}{\partial
x_{i}}\frac{\partial\widetilde{u}_{i}}{\partial x_{j}}$ (18)
$\displaystyle\widetilde{E}$
$\displaystyle=\frac{1}{2}\left(\widetilde{W}_{ij}\widetilde{W}_{ij}+\widetilde{S}_{ij}\widetilde{S}_{ij}\right)=$
$\displaystyle\ \ \ \frac{1}{2}\frac{\partial\widetilde{u}_{j}}{\partial
x_{i}}\frac{\partial\widetilde{u}_{j}}{\partial x_{i}}$ (19)

$\widetilde{S}_{ij}$ is the resolved strain-rate tensor and
$\widetilde{W}_{ij}$ is the resolved vorticity tensor.

The coefficient for the isotropic part of the SGS tensor, $C_{I}$, is
dynamically modeled in a similar fashion to the Smagorinsky coefficient.

$C_{I}=C_{csm_{I}}|F_{cs}|^{3/2}\left(1-F_{cs}\right)$ (20)

$C_{csm_{I}}$ is the CSM isotropic model coefficient and is chosen here such
that the maximum value of $C_{I}$ is $0.0066$ [49]

The SGS heat flux term, $q_{i}^{sgs}$, is also modeled using an eddy-viscosity
model.

$q_{i}^{sgs}=-\frac{\mu_{sgs}}{\left(\gamma-1\right)M_{\infty}^{2}RePr_{sgs}}\frac{\partial\widetilde{T}}{\partial
x_{i}}$ (21)

The SGS Prandtl number, $Pr_{sgs}$, is given the constant value of $0.7$ [49,
31].

##  3 Numerical Algorithm

The compressible Navier-Stokes equations are solved in the framework of Large
Eddy Simulations, where the Coherent Structure Model is applied to account for
the missing sub-grid scale energy. The numerical algorithm uses high-order
finite difference approximations for the spatial derivatives and explicit time
marching. The time integration is performed using a fully-explicit fourth-
order Adams-Bashforth scheme [6] of the form

$\displaystyle\mathbf{Q}^{n+1}=\mathbf{Q}^{n}+\Delta
t\left[\frac{55}{24}L\left(\mathbf{Q}^{n}\right)-\frac{59}{24}L\left(\mathbf{Q}^{n-1}\right)+\frac{37}{24}L\left(\mathbf{Q}^{n-2}\right)-\frac{9}{24}L\left(\mathbf{Q}^{n-3}\right)\right]$
(22)

where $L(\mathbf{Q})$ is the residual and $n$ is the current time level.

The spatial derivatives are discretized using the dispersion-relation-
preserving finite difference scheme of Tam and Webb [46]. The first derivative
at the $l$th node is approximated using $M$ values of $f$ to the right and $N$
values of $f$ to left of the node.

$\left(\frac{\partial f}{\partial x}\right)_{l}\simeq\frac{1}{\Delta
x}\sum_{j=-N}^{M}a_{j}f_{l+j}$ (23)

By taking the Fourier transform of the above equation, the coefficients
$a_{j}$ are found by minimizing the integrated error of the difference between
the wavenumber of the finite difference scheme and the wavenumber of the
Fourier transform of the finite difference scheme. The coefficients $a_{j}$
are given in Table 1.

Stencil | $a_{1}=-a_{-1}$ | $a_{2}=-a_{-2}$ | $a_{3}=-a_{-3}$ | $a_{4}=-a_{-4}$
---|---|---|---|---
$DRP$ | 0.77088238 | -0.16670590 | 0.02084314 | 0
Table 1: Weights of the centered stencil.

To damp out the unwanted high wavenumber waves from the solution, high-order
spatial filters, as developed by Kennedy and Carpenter [24], are used.

The Triple Hill’s Vortex Synthetic Eddy Method (THV SEM) was applied at the
inlet for all cases. No slip boundary condition for velocity and adiabatic
condition for temperature are imposed wherever a wall is present and periodic
boundary conditions are utilized in homogeneous directions. Standard outflow
boundary conditions are applied at the outlet along with a region of
artificially increased viscosity immediately upstream of the outlet to damp
spurious waves. The dynamic viscosity, $\mu_{m}$, in the exit region is
smoothly increased in the streamwise direction by a weighting function,
$w_{\mu}$, an artificial viscosity, $\mu_{BC}$, defined using a Smagorinsky
type eddy-viscosity model [47].

$\mu\left(x\right)=\mu_{m}+w_{\mu}\left(x\right)\mu_{BC}=w_{\mu}\left(x\right)\left[\overline{\rho}C_{s}\Delta^{2}|\widetilde{S}|\right]$
(24)

In practice, $\mu_{BC}$ is equal to the local SGS viscosity calculated by the
SGS model in Section 2.1.

##  4 Triple Hill’s Vortex Synthetic Eddy Method

Hill’s spherical vortex (Hill [15]) represents one of the best-known examples
of a steady rotational solution to the classical Euler equations, modeling an
inviscid incompressible flow. It is characterized by one amplitude, and by one
degree of freedom associated with the translation along its axis, which is a
fundamental property of the Hill’s spherical vortex as a consequence of its
definition in relation to a uniform flow. Thus, a synthetic turbulence model
that is based on a superposition of eddies representing single Hill’s vortices
would only have the freedom to reproduce isotropic turbulence. By considering
a new vortex structure, one that is composed of three superposed Hill’s
spherical vortices with the axes perpendicular to each other, two more
associated amplitudes are introduced and thus, two more degrees of freedom
that can be utilized to match a given Reynolds stress tensor are made
available. A synthetic turbulence model consisting of a superposition of such
Triple Hill’s Vortices (THV) will have enough freedom to reproduce anisotropic
turbulence.

###  4.1 Hill’s Spherical Vortex

The base component of the proposed synthetic eddy method is the Hill’s
spherical vortex. It is a steady, axisymmetric solution to the Euler equations
for an incompressible flow. Derived from the incompressible Euler equations,
the Helmhotz equation for vorticity is

$\frac{\partial\boldsymbol{\omega}}{\partial
t}+(\mathbf{u}\cdot\nabla)\boldsymbol{\omega}=(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}$
(25)

where $\boldsymbol{\omega}=(\omega_{r},\omega_{\theta},\omega_{z})$. It can be
combined with the continuity equation for an incompressible flow,
$\nabla\cdot\mathbf{u}=0$, and through using the definition of a
streamfunction in cylindrical coordinates for an axisymmetric flow,

$u_{r}(r,z)=\frac{1}{r}\frac{\partial\psi}{\partial r},\hskip
14.22636ptu_{z}(r,z)=-\frac{1}{r}\frac{\partial\psi}{\partial z}$ (26)

and the fact that $\omega_{\theta}/r$ is constant along a streamline and only
depends on the value of the stream function (i.e.
$\partial\left(\omega_{\theta}/r\right)/\partial t=0$),

$\frac{\omega_{\theta}}{r}=f(\psi),$ (27)

the following equation governing the streamfunction is obtained.

$\frac{\partial^{2}\psi}{\partial z^{2}}+\frac{\partial^{2}\psi}{\partial
r^{2}}-\frac{1}{r}\frac{\partial\psi}{\partial r}=-r^{2}f(\psi)$ (28)

Considering the assumption that $\omega_{\theta}/r=f(\psi)=A=const$ inside a
sphere of radius $a$, an exact solution to equation (28) can be found. This
solution is known as Hill’s spherical vortex [15]. Using the boundary
condition at the surface of the sphere, $\psi=0$, the streamfunction inside
the sphere ($r^{2}+z^{2}<a^{2}$) is

$\psi(r,z)=\frac{Ar^{2}}{10}\left(a^{2}-z^{2}-r^{2}\right)$ (29)

The streamfunction outside the sphere ($r^{2}+z^{2}>a^{2}$),

$\psi(r,z)=-\frac{u_{0}r^{2}}{2}\left[1-a^{3}\left(r^{2}+z^{2}\right)^{-3/2}\right],$
(30)

corresponds to the potential flow around a solid sphere of radius $a$ in a
uniform flow of speed $u_{0}$ in the negative z-axis direction. By matching
the two streamfunction solutions at the surface of the sphere, the constant
$A$ can be found as

$A=\frac{15u_{0}}{2a}$ (31)

From equation (26), the velocities in cylindrical coordinates ( $u_{r}$,
$u_{z}$ ) of the Hill’s spherical vortex inside the sphere

$u_{r}(r,z)=\frac{3}{2}u_{0}\frac{zr}{a^{2}},\hskip
17.07164ptu_{z}(r,z)=\frac{3}{2}u_{0}\left(1-\frac{z^{2}+2r^{2}}{a^{2}}\right)$
(32)

and outside the sphere

$u_{r}(r,z)=\frac{3}{2}u_{0}\frac{zr}{a^{2}}\left(\frac{a^{2}}{z^{2}+r^{2}}\right)^{5/2},\hskip
17.07164ptu_{z}(r,z)=u_{0}\left[\left(\frac{a^{2}}{z^{2}+r^{2}}\right)^{5/2}\frac{2z^{2}-r^{2}}{2a^{2}}-1\right]$
(33)

can be found. As a note, the axial velocity outside of the sphere, $u_{z}$, is
presented after the subtraction of the surrounding uniform flow in order to
isolate the velocity contribution of the Hill’s spherical vortex.

###  4.2 Triple Hill’s Vortex

Instead of viewing the fundamental synthetic eddy as a distinct Hill’s vortex
structure, this method considers a more complex flow structure: the Triple
Hill’s Vortex (THV). The THV is a superposition of three independent Hill’s
vortices that all share the same center ($x_{0}$, $y_{0}$, $z_{0}$) and the
same outer radius $a$. The rotation axis of each Hill’s vortex is orthogonal
to the rotation axis of the other two Hill’s vortices. In Cartesian
coordinates, the rotation axis of one Hill’s vortex with amplitude
$u_{0}=u_{0}^{x}$ is oriented in the $x$-direction, and has the associated
velocity components given as

$\displaystyle u^{x}(x,y,z)$ $\displaystyle=$ $\displaystyle u_{z}(r,z)$
$\displaystyle v^{x}(x,y,z)$ $\displaystyle=$ $\displaystyle
u_{r}(r,z)sin(\theta)$ (34) $\displaystyle w^{x}(x,y,z)$ $\displaystyle=$
$\displaystyle u_{r}(r,z)cos(\theta)$ $r=\sqrt{y^{2}+z^{2}},\hskip
17.07164ptz=x,\hskip 17.07164pt\theta=arctan\left(\frac{y}{z}\right)$

The rotation axis of another Hill’s vortex with amplitude $u_{0}=u_{0}^{y}$ is
oriented in the $y$-direction, and has the form

$\displaystyle u^{y}(x,y,z)$ $\displaystyle=$ $\displaystyle
u_{r}(r,z)cos(\theta)$ $\displaystyle v^{y}(x,y,z)$ $\displaystyle=$
$\displaystyle u_{z}(r,z)$ (35) $\displaystyle w^{y}(x,y,z)$ $\displaystyle=$
$\displaystyle u_{r}(r,z)sin(\theta)$ $r=\sqrt{x^{2}+z^{2}},\hskip
17.07164ptz=y,\hskip 17.07164pt\theta=arctan\left(\frac{z}{x}\right)$

The rotation axis of the third Hill’s vortex with amplitude $u_{0}=u_{0}^{z}$
is oriented in the $z$-direction, and has the form

$\displaystyle u^{z}(x,y,z)$ $\displaystyle=$ $\displaystyle
u_{r}(r,z)sin(\theta)$ $\displaystyle v^{z}(x,y,z)$ $\displaystyle=$
$\displaystyle u_{r}(r,z)cos(\theta)$ (36) $\displaystyle w^{z}(x,y,z)$
$\displaystyle=$ $\displaystyle u_{z}(r,z)$ $r=\sqrt{x^{2}+y^{2}},\hskip
17.07164ptz=z,\hskip 17.07164pt\theta=arctan\left(\frac{x}{y}\right)$

In equations (4.2)-(4.2), $r$ and $z$ are the local polar coordinates
associated with the vortex (the $z$ coordinate is in the direction of the
rotation axis), and $u_{r}$ and $u_{z}$ are the velocity components in the
local polar coordinate system. The velocity components $u_{r}$ and $u_{z}$ of
each of the vortices can be found from equations (32) and 33). The superscript
signifies the orientation direction of the Hill’s vortex rotation axis (for
example, $u^{y}$ is the $x$-component of velocity for the Hill’s vortex with
the rotation axis oriented in the $y$-direction). Thus, the velocity
components for the THV can be obtained by superposition according to

$\displaystyle\widehat{\widetilde{u}}$ $\displaystyle=$ $\displaystyle
u^{x}+u^{y}+u^{z}$ $\displaystyle\widehat{\widetilde{v}}$ $\displaystyle=$
$\displaystyle v^{x}+v^{y}+v^{z}$ (37) $\displaystyle\widehat{\widetilde{w}}$
$\displaystyle=$ $\displaystyle w^{x}+w^{y}+w^{z}$

where the tilde ($\widehat{\widetilde{\quad}}$) signifies the entire THV.

###  4.3 Divergence of a Triple Hill’s Vortex

Consider the divergence of the velocity field for a single THV,
$\mathbf{\widehat{\widetilde{u}}}$.

$\displaystyle\nabla\cdot\mathbf{\widehat{\widetilde{u}}}$ $\displaystyle=$
$\displaystyle\frac{\partial\widehat{\widetilde{u}}}{\partial
x}+\frac{\partial\widehat{\widetilde{v}}}{\partial
y}+\frac{\partial\widehat{\widetilde{w}}}{\partial z}$ $\displaystyle=$
$\displaystyle\frac{\partial\left(u^{x}+u^{y}+u^{z}\right)}{\partial
x}+\frac{\partial\left(v^{x}+v^{y}+v^{z}\right)}{\partial
y}+\frac{\partial\left(w^{x}+w^{y}+w^{z}\right)}{\partial z}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial u^{x}}{\partial x}+\frac{\partial
v^{x}}{\partial y}+\frac{\partial w^{x}}{\partial
z}\right)+\left(\frac{\partial u^{y}}{\partial x}+\frac{\partial
v^{y}}{\partial y}+\frac{\partial w^{y}}{\partial
z}\right)+\left(\frac{\partial u^{z}}{\partial x}+\frac{\partial
v^{z}}{\partial y}+\frac{\partial w^{z}}{\partial z}\right)$ $\displaystyle=$
$\displaystyle\left(\nabla\cdot\mathbf{u}^{x}\right)+\left(\nabla\cdot\mathbf{u}^{y}\right)+\left(\nabla\cdot\mathbf{u}^{z}\right)$

where, for example, $u^{x}$, $v^{x}$ and $w^{x}$ are the velocity components
in the local coordinate system that is associated with the vortex with its
axis aligned with the $x$ axis (of the global coordinate system). Therefore,
since each Hill’s vortex satisfies the divergence-free condition in its local
coordinate system,

$\nabla\cdot\mathbf{\widehat{\widetilde{u}}}=0$ (39)

the THV is divergence-free.

###  4.4 Convection of a Triple Hill’s Vortex

The THV’s are convected through the inlet using Taylor’s frozen turbulence
hypothesis

$\displaystyle x^{\prime}=(x-x_{0})-U(t-t0)$ $\displaystyle
y^{\prime}=(y-y_{0})-V(t-t0)$ (40) $\displaystyle
z^{\prime}=(z-z_{0})-W(t-t0)$

where ($U$,$V$,$W$) are the mean velocity components and
($x_{0}$,$y_{0}$,$z_{0}$,$t_{0}$) are the space and time coordinates of the
center of the THV. The divergence of velocity for a single THV is only
identically zero when ($U$,$V$,$W$) are constant over the entire THV. For
spatially varying inflow velocity profiles, such as in channel flow, the mean
inflow velocity vary continuously across the THV. So, the following assumption
was made: the ($U$,$V$,$W$) for a single THV are calculated at the center of
the THV and applied over the entire THV.

$\displaystyle x^{\prime}=(x-x_{0})-U(x_{0},y_{0},z_{0})(t-t0)$ $\displaystyle
y^{\prime}=(y-y_{0})-V(x_{0},y_{0},z_{0})(t-t0)$ (41) $\displaystyle
z^{\prime}=(z-z_{0})-W(x_{0},y_{0},z_{0})(t-t0)$

###  4.5 Inflow Velocity

The calculation of the THV Synthetic Eddy Method is similar to the method
proposed by Jarrin et al. [16], except the matching between the imposed and
given Reynolds stress tensors is performed differently, in order to maintain
the divergence-free condition. Also, the previous methods used Gaussian and
tent functions to define the shape of these eddies, which may not reproduce
realistic flow fields. Since the synthetic eddies (THV) proposed in this model
are composed of Hill’s vortices, which are actual solutions to the Euler
equation, it is expected that they will resemble more realistic turbulent flow
conditions. The cases presented in Section 5 will investigate whether this
expectation is met. The turbulent inflow is viewed as a collection of eddies
added to a mean flow. In the following, the velocities in equations (4.2),
(4.2), and (4.2) are assumed to be products of a constant amplitude
$\left(u_{0}^{x},u_{0}^{y},u_{0}^{z}\right)$ multiplied by a shape function
$\left(f_{i}^{x},f_{i}^{y},f_{i}^{z}\right)$,

$\left[\begin{array}[]{c}u^{x}\\\ v^{x}\\\
w^{x}\end{array}\right]=\left[\begin{array}[]{c}u_{0}^{x}f_{u}^{x}\\\
u_{0}^{x}f_{v}^{x}\\\ u_{0}^{x}f_{w}^{x}\end{array}\right],\hskip
17.07164pt\left[\begin{array}[]{c}u^{y}\\\ v^{y}\\\
w^{y}\end{array}\right]=\left[\begin{array}[]{c}u_{0}^{y}f_{u}^{y}\\\
u_{0}^{y}f_{v}^{y}\\\ u_{0}^{y}f_{w}^{y}\end{array}\right],\hskip
17.07164pt\left[\begin{array}[]{c}u^{z}\\\ v^{z}\\\
w^{z}\end{array}\right]=\left[\begin{array}[]{c}u_{0}^{z}f_{u}^{z}\\\
u_{0}^{z}f_{v}^{z}\\\ u_{0}^{z}f_{w}^{z}\end{array}\right]$ (42)

where, instead of using Gaussian shape functions or other simple functions (as
in Jarrin et al. [16] and Poletto [35]), THV shape functions are utilized
(extracted from equations (32) and (33)). The three amplitudes associated with
the THV are the means by which the Reynolds stress tensor components of a
prescribed inflow velocity are matched. The inflow velocity components,
($u_{in}$, $v_{in}$, $w_{in}$), are composed of the mean base flow and the
fluctuating components

$\displaystyle u_{in}$ $\displaystyle=$ $\displaystyle
U+\sum_{j=1}^{N}\widehat{\widetilde{u}}_{j}$ $\displaystyle v_{in}$
$\displaystyle=$ $\displaystyle V+\sum_{j=1}^{N}\widehat{\widetilde{v}}_{j}$
(43) $\displaystyle w_{in}$ $\displaystyle=$ $\displaystyle
W+\sum_{j=1}^{N}\widehat{\widetilde{w}}_{j}$

where ($U$,$V$,$W$) is the mean velocity vector, and
($\widehat{\widetilde{u}}_{j}$,$\widehat{\widetilde{v}}_{j}$,$\widehat{\widetilde{w}}_{j}$)
is the velocity vector of the $j$th THV. Each THV has an independent random
center, radius, and amplitude.

###  4.6 Determination of the Amplitudes

In the original SEM of Jarrin et al. [16], the fluctuating velocity field is
rescaled using the Cholesky decomposition of the Reynolds stress tensor in
order to ensure reproduction of the desired Reynolds stresses. Since those
synthetic eddies already violate the divergence-free condition, the fact that
rescaling using the Cholesky decomposition generally violates the divergence-
free condition is of little consequence. Because the THV is divergence-free,
the amplitudes of the THV’s need to be calculated using a different method in
order to match the desired Reynolds stresses and also preserve the divergence-
free nature of the synthetic field.

Smirnov et al. [44] proposed a method, which was later expanded upon by
Davidson and Bilson [9], to match the Reynolds stresses for anisotropic
turbulence. The idea behind the method is to find a local reference system
where the normal Reynolds stress terms are non-zero, while the off-diagonal
Reynolds stress terms are zero, calculate the associated amplitudes, and then
transform the signal back to the global reference system. To be able to apply
this method in the current THV framework, the creation of the THV’s has to
take place in the local principal-axis reference system. This begins by
calculating the principal-axis Reynolds stresses and eigenvectors of the local
Reynolds stress tensor. Examining the principal-axis Reynolds stresses for a
single stream of THV’s can give insight into how to set each of the principal-
axis amplitudes, ($u_{0}^{x,p}$, $u_{0}^{y,p}$, $u_{0}^{z,p}$), where
superscript p denotes variables in the principal-axis coordinate system. Since
the velocity components are at a maximum, consider the principal-axis Reynolds
stresses at the center of a single stream of THV’s. A single stream of THV’s
is defined as an infinite number of identical THV’s whose centers pass through
the same point in space on the inflow plane. The frequency with which the
THV’s pass through the inflow plane is such that the preceeding and subsequent
THV’s do not interact with the THV currently at the inlet.

$\displaystyle\left<\widehat{\widetilde{u}}\widehat{\widetilde{u}}\right>^{p}$
$\displaystyle=$
$\displaystyle\left<\left(u_{0}^{x,p}\right)^{2}\left(f_{u}^{x,p}\right)^{2}\right>+\left<\left(u_{0}^{y,p}\right)^{2}\left(f_{u}^{y,p}\right)^{2}\right>+\left<\left(u_{0}^{z,p}\right)^{2}\left(f_{u}^{z,p}\right)^{2}\right>$
$\displaystyle\left<\widehat{\widetilde{v}}\widehat{\widetilde{v}}\right>^{p}$
$\displaystyle=$
$\displaystyle\left<\left(u_{0}^{x,p}\right)^{2}\left(f_{v}^{x,p}\right)^{2}\right>+\left<\left(u_{0}^{y,p}\right)^{2}\left(f_{v}^{y,p}\right)^{2}\right>+\left<\left(u_{0}^{z,p}\right)^{2}\left(f_{v}^{z,p}\right)^{2}\right>$
(44)
$\displaystyle\left<\widehat{\widetilde{w}}\widehat{\widetilde{w}}\right>^{p}$
$\displaystyle=$
$\displaystyle\left<\left(u_{0}^{x,p}\right)^{2}\left(f_{w}^{x,p}\right)^{2}\right>+\left<\left(u_{0}^{y,p}\right)^{2}\left(f_{w}^{y,p}\right)^{2}\right>+\left<\left(u_{0}^{z,p}\right)^{2}\left(f_{w}^{z,p}\right)^{2}\right>$
$\left<\widehat{\widetilde{u}}\widehat{\widetilde{v}}\right>^{p}=\left<\widehat{\widetilde{u}}\widehat{\widetilde{w}}\right>^{p}=\left<\widehat{\widetilde{v}}\widehat{\widetilde{w}}\right>^{p}=0$

Let the amplitudes in the principal-axis coordinate system associated with a
THV be defined as a product of a constant amplitude and a random number,

$\displaystyle\left[\begin{array}[]{c}u_{0}^{x,p}\\\ u_{0}^{y,p}\\\
u_{0}^{z,p}\end{array}\right]=\left[\begin{array}[]{c}\epsilon^{x}\
\hat{u}_{0}^{x,p}\\\ \epsilon^{y}\ \hat{u}_{0}^{y,p}\\\ \epsilon^{z}\
\hat{u}_{0}^{z,p}\end{array}\right]$ (51)

where $\epsilon^{x}$, $\epsilon^{y}$, and $\epsilon^{z}$ are independent
random numbers such that $\epsilon=\pm 1$, $\left<\epsilon\right>=0$,
$\left<\epsilon^{2}\right>=1$, and
$\left<\epsilon^{x}\epsilon^{y}\right>=\left<\epsilon^{x}\epsilon^{z}\right>=\left<\epsilon^{y}\epsilon^{z}\right>=0$
. Since the amplitudes and the constant shape functions are independent, after
inserting equation (51) and using the properties of $\epsilon$, equation (4.6)
becomes

$\displaystyle\left[\begin{array}[]{c}\left<\widehat{\widetilde{u}}\widehat{\widetilde{u}}\right>^{p}\\\
\left<\widehat{\widetilde{v}}\widehat{\widetilde{v}}\right>^{p}\\\
\left<\widehat{\widetilde{w}}\widehat{\widetilde{w}}\right>^{p}\end{array}\right]=\left[\begin{array}[]{ccc}\left(f_{u}^{x,p}\right)^{2}\left(f_{u}^{y,p}\right)^{2}\left(f_{u}^{z,p}\right)^{2}\\\
\left(f_{v}^{x,p}\right)^{2}\left(f_{v}^{y,p}\right)^{2}\left(f_{v}^{z,p}\right)^{2}\\\
\left(f_{w}^{x,p}\right)^{2}\left(f_{w}^{y,p}\right)^{2}\left(f_{w}^{z,p}\right)^{2}\end{array}\right]\left[\begin{array}[]{c}\left(\hat{u}_{0}^{x,p}\right)^{2}\\\
\left(\hat{u}_{0}^{y,p}\right)^{2}\\\
\left(\hat{u}_{0}^{z,p}\right)^{2}\end{array}\right]=F_{uvw}^{xyz}\left[\begin{array}[]{c}\left(\hat{u}_{0}^{x,p}\right)^{2}\\\
\left(\hat{u}_{0}^{y,p}\right)^{2}\\\
\left(\hat{u}_{0}^{z,p}\right)^{2}\end{array}\right]$ (64)

$F_{uvw}^{xyz}$ is constant and invertible. Thus, the amplitudes of a THV can
be calculated from the given principal-axis Reynolds stresses. By combining
the amplitudes calculated in equation (64) with equations (51) and (42),

$\left[\begin{array}[]{c}u^{x,p}\\\ v^{x,p}\\\
w^{x,p}\end{array}\right]=\left[\begin{array}[]{c}\epsilon^{x}\
\hat{u}_{0}^{x,p}f_{u}^{x,p}\\\ \epsilon^{x}\ \hat{u}_{0}^{x,p}f_{v}^{x,p}\\\
\epsilon^{x}\
\hat{u}_{0}^{x,p}f_{w}^{x,p}\end{array}\right],\left[\begin{array}[]{c}u^{y,p}\\\
v^{y,p}\\\ w^{y,p}\end{array}\right]=\left[\begin{array}[]{c}\epsilon^{y}\
\hat{u}_{0}^{y,p}f_{u}^{y,p}\\\ \epsilon^{y}\ \hat{u}_{0}^{y,p}f_{v}^{y,p}\\\
\epsilon^{y}\
\hat{u}_{0}^{y,p}f_{w}^{y,p}\end{array}\right],\left[\begin{array}[]{c}u^{z,p}\\\
v^{z,p}\\\ w^{z,p}\end{array}\right]=\left[\begin{array}[]{c}\epsilon^{z}\
\hat{u}_{0}^{z,p}f_{u}^{z,p}\\\ \epsilon^{z}\ \hat{u}_{0}^{z,p}f_{v}^{z,p}\\\
\epsilon^{z}\ \hat{u}_{0}^{z,p}f_{w}^{z,p}\end{array}\right]$ (65)

and then inserting the result into equation (4.2), the velocity contribution
of a single THV in the local principal-axis reference system can be found as

$\displaystyle\widehat{\widetilde{u}}^{p}$ $\displaystyle=$ $\displaystyle
u^{x,p}+u^{y,p}+u^{z,p}$ $\displaystyle\widehat{\widetilde{v}}^{p}$
$\displaystyle=$ $\displaystyle v^{x,p}+v^{y,p}+v^{z,p}$ (66)
$\displaystyle\widehat{\widetilde{w}}^{p}$ $\displaystyle=$ $\displaystyle
w^{x,p}+w^{y,p}+w^{z,p}$

The last step is to transform the velocities from the local principal-axis
reference system to the global reference system using the transformation
matrix created from the eigenvectors of the local Reynolds stress tensor,
$T_{p}^{G}$.

$\displaystyle\left[\begin{array}[]{c}\widehat{\widetilde{u}}\\\
\widehat{\widetilde{v}}\\\ \widehat{\widetilde{w}}\end{array}\right]$
$\displaystyle=$ $\displaystyle
T_{p}^{G}\left[\begin{array}[]{c}\widehat{\widetilde{u}}^{p}\\\
\widehat{\widetilde{v}}^{p}\\\ \widehat{\widetilde{w}}^{p}\end{array}\right]$
(73)

This is the velocity contribution of a single THV in the global reference
system that is combined with all the other THV’s in equation (4.5) to create
the imposed inflow.

###  4.7 Modification of the Target Reynolds Stresses

As is defined in Section 4.6, the imposed Reynolds stress tensor at any point
in space is reproduced by a single stream of THV’s moving through that point
in time. To account for that fact that the imposed turbulent inflow is a
superposition of many THV’s over a range of different sizes, the method for
determining the amplitudes of each THV needs to include the influence of the
THV’s currently at the inlet when a new THV is created in order to ensure that
Reynolds stress matching is recovered. The target Reynolds stress modification
procedure is proposed to address these issue.

Consider the inlet plane shown in Figure 1. The black circles represent THV’s
that are currently passing through the inlet and the red circle represents a
new THV that is about to be created. If the surrounding THV’s did not exist,
the calculated amplitudes of the new THV ensure Reynolds stress matching at
the center of the THV. Since the surrounding THV’s already exist and are
influencing the flow field at the center of the new THV, the amplitudes of the
new THV need to be modified to account for the contributions of its neighbors.

new$1$$2$$3$$4$$5$$6$$8$ Figure 1:  An example of inlet plane when a new THV
is being created: black) The current THV’s being convected through the inflow
plane; red) The random location on the inflow plane where a new THV is being
created.

Looking at equation (33), the velocity components for a Hill’s vortex approach
zero quickly when moving away from the surface of the vortex. Thus, only the
contribution of THV’s that encircle the center of the new THV will be
considered (the center of THV new falls within the boundaries of THV 1,2, and
3). In the global coordinate system, the Reynolds stresses at the point on the
inlet where the center of the new THV will pass through are defined as

$\displaystyle\langle u^{\prime}_{i}u^{\prime}_{j}\rangle_{existing}$
$\displaystyle=$
$\displaystyle\left\langle\left(\sum_{m=1}^{M_{s}}\left(\widetilde{u}_{i}\right)_{m}\right)\left(\sum_{n=1}^{N_{s}}\left(\widehat{\widetilde{u}}_{j}\right)_{n}\right)\right\rangle$
(74) $\displaystyle=$
$\displaystyle\sum_{m=1}^{M_{s}}\sum_{n=1}^{N_{s}}\langle\left(\widehat{\widetilde{u}}_{i}\right)_{m}\left(\widehat{\widetilde{u}}_{j}\right)_{n}\rangle$
(75)

where $M_{s}$ and $N_{s}$ are both equal to the number of encircling THV’s.
Since the random numbers associated with each THV are independent,

$\displaystyle\langle\left(\widehat{\widetilde{u}}_{i}\right)_{m}\left(\widehat{\widetilde{u}}_{j}\right)_{n}\rangle=0\quad\text{for}\quad
m\neq n$ (76)

Thus, the Reynolds stress contribution from the surrounding THV’s is

$\displaystyle\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{existing}=\sum_{n=1}^{N_{s}}\langle\widehat{\widetilde{u}}_{i}\widehat{\widetilde{u}}_{j}\rangle_{n}$
(77)

where $N_{s}$ is again the number of encircling THV’s. This contribution
Reynolds stress tensor is then used to modify the given Reynolds stress tensor
that is being matched.

$\displaystyle\langle u^{\prime}_{i}u^{\prime}_{j}\rangle_{target}=\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{given}-\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{existing}$ (78)

The target Reynolds stress tensor can then be used in the THV creation
framework to determine the amplitudes of the new THV. As a note, if any of the
eigenvalues of the target Reynolds stress tensor are negative, that means that
the given Reynolds stresses are already reproduced by the present THV’s and
the new THV does not need to be created at that point on the inlet.

####  4.7.1 Reynolds Stresses for a Single Stream of THV

The Reynolds stress tensors of the individual THV’s needed by equation (77)
can be found by first considering the velocity components of a single THV in
the global coordinate system.

$\displaystyle\widehat{\widetilde{u}}_{i}=u_{0}^{x}f_{i}^{x}+u_{0}^{2}f_{i}^{y}+u_{0}^{3}f_{i}^{z}$
(79)

$\left(f_{i}^{x},f_{i}^{y},f_{i}^{z}\right)$ are the shape functions for the
three Hill’s vortices associated with the $\widehat{\widetilde{u}}_{i}$
velocity component. $\left(u_{0}^{x},u_{0}^{y},u_{0}^{z}\right)$ are the
randomized amplitudes in the global coordinate system and are defined as

$\displaystyle\left[\begin{array}[]{c}u_{0}^{x}\\\ u_{0}^{y}\\\
u_{0}^{z}\end{array}\right]$ $\displaystyle=$ $\displaystyle
T_{p}^{G}\left[\begin{array}[]{c}\epsilon^{x}\thinspace\hat{u}_{0}^{x,p}\\\
\epsilon^{y}\thinspace\hat{u}_{0}^{y,p}\\\
\epsilon^{z}\thinspace\hat{u}_{0}^{z,p}\end{array}\right]$ (86)

where $T_{p}^{G}$ is the eigenvector transformation matrix,
$\left(\hat{u}_{0}^{x,p},\hat{u}_{0}^{y,p},\hat{u}_{0}^{z,p}\right)$ are the
constant amplitudes in the principal coordinate system, and
$\left(\epsilon^{x},\epsilon^{y},\epsilon^{z}\right)$ are independent random
numbers that have the following properties

$\displaystyle\epsilon^{x}=\pm 1\quad;\quad\epsilon^{y}=\pm
1\quad;\quad\epsilon^{z}=\pm 1$ (87)
$\displaystyle\langle\epsilon^{x}\rangle=\langle\epsilon^{y}\rangle=\langle\epsilon^{z}\rangle=0$
$\displaystyle\langle\left(\epsilon^{x}\right)^{2}\rangle=\langle\left(\epsilon^{y}\right)^{2}\rangle=\langle\left(\epsilon^{z}\right)^{2}\rangle=1$
$\displaystyle\langle\epsilon^{x}\epsilon^{y}\rangle=\langle\epsilon^{x}\epsilon^{z}\rangle=\langle\epsilon^{y}\epsilon^{z}\rangle=0$

By taking the time-average of the products of the velocity components in
equation (79) and recognizing that the amplitudes and shape functions are
independent from each other, expressions can be found for the Reynolds
stresses of an infinite stream of THV’s being convected through the same
location at the inlet.

$\displaystyle\langle\widehat{\widetilde{u}}_{i}\widehat{\widetilde{u}}_{j}\rangle=\langle\left(u_{0}^{x}\right)^{2}\rangle\langle
f_{i}^{x}f_{j}^{x}\rangle$ $\displaystyle+$
$\displaystyle\langle\left(u_{0}^{y}\right)^{2}\rangle\langle
f_{i}^{y}f_{j}^{y}\rangle+\langle\left(u_{0}^{z}\right)^{2}\rangle\langle
f_{i}^{z}f_{j}^{z}\rangle$ $\displaystyle+$ $\displaystyle\left(\langle
f_{i}^{x}f_{j}^{y}\rangle+\langle f_{j}^{x}f_{i}^{y}\rangle\right)\langle
u_{0}^{x}u_{0}^{y}\rangle+\left(\langle f_{i}^{x}f_{j}^{z}\rangle+\langle
f_{j}^{x}f_{i}^{z}\rangle\right)\langle u_{0}^{x}u_{0}^{z}\rangle$
$\displaystyle+$ $\displaystyle\left(\langle f_{i}^{y}f_{j}^{z}\rangle+\langle
f_{j}^{y}f_{i}^{z}\rangle\right)\langle u_{0}^{y}u_{0}^{z}\rangle$

The time-average of the products of the shape functions can be solved for
analytically by assuming the time period for a single THV to flow through the
inlet is

$\displaystyle T_{in}=\frac{4\thinspace a}{U}$ (89)

where $a$ is the radius of the THV and $U$ is the mean velocity at the center
of the THV.

Through the multiplying of the global amplitudes in equation (86), taking the
time-average of the resulting products, recognizing that the eigenvector
transformation matrix and the principal coordinate system amplitudes are all
constant, and finally inserting the properties of the random numbers seen in
equation (87); the time-averaged global coordinate system amplitude
correlations can be determined.

$\displaystyle\langle
u_{0}^{i}u_{0}^{j}\rangle=\left(t_{p}^{G}\right)_{i1}\left(t_{p}^{G}\right)_{j1}\left(\hat{u}_{0}^{x,p}\right)^{2}+\left(t_{p}^{G}\right)_{i2}\left(t_{p}^{G}\right)_{j2}\left(\hat{u}_{0}^{y,p}\right)^{2}+\left(t_{p}^{G}\right)_{i3}\left(t_{p}^{G}\right)_{j3}\left(\hat{u}_{0}^{z,p}\right)^{2}$
(90)

$t_{p}^{G}$ are the individual elements of the eigenvector transformation
matrix, $T_{p}^{G}$.

###  4.8 Generations of Triple Hill’s Vortices

To be able to recreate turbulent flow, the Reynolds stress tensor is not the
only quantity that needs to be matched. Turbulent structures occur over a wide
range of length scales, and the population varies over a length scale
distribution which can be taken from the turbulent kinetic energy spectrum.
Thus, the number of larger structures is less than the number of medium sized
structures, which is less than the number of smaller structures. This type of
distribution also needs to be replicated to make sure the space is fully
covered by structures in the whole range of wavenumbers. The minimum,
$a_{min}$, and maximum, $a_{max}$, radii of the THV allowed are dictated by
the grid resolution and the integral length scale or the domain size,
respectively. As shown in Jarrin et al. [17], the overall number of imposed
THV’s are calculated as follows,

$\displaystyle N$ $\displaystyle=$ $\displaystyle C\frac{A_{in}}{A_{THV}}$
(91)

where $N$ is the total number of THV’s, $A_{in}$ is the area of the inlet
plane, $A_{THV}$ is the projected area of a THV onto the inlet plane, and $C$
is a proportionality constant. This number of THV’s will ensure that the
inflow plane is fully populated with THV’s at any moment in time. Jarrin et
al. [17] recommends that $C=1$ to minimize the computational overhead required
with any unnecessary THV’s. Given the limits for the smallest and the largest
eddy, the idea of defining multiple generations of THV’s is introduced here.
The range of radii is divided discretely into several subranges, or
generations. Based on a given distribution, the number of THV’s in each
generation is set. Table 2 shows an example of five generations of THV’s,
where $a_{min}<a_{4}<a_{3}<a_{2}<a_{1}<a_{max}$ and
$N=N_{1}+N_{2}+N_{3}+N_{4}+N_{5}$

Generation | Range of Radii | Number of THV’s
---|---|---
$1$ | $a_{1}<a<a_{max}$ | $N_{1}$
$2$ | $a_{2}<a<a_{1}$ | $N_{2}$
$3$ | $a_{3}<a<a_{2}$ | $N_{3}$
$4$ | $a_{4}<a<a_{3}$ | $N_{4}$
$5$ | $a_{min}<a<a_{4}$ | $N_{5}$
Table 2: Example of five generations of THV’s

The multiple generations of THV’s ensure that the inlet is covered with enough
multiscale structures. The range of locations for these generations can also
be restricted, for example, to cluster smaller THV’s near a wall for wall-
bounded flows.

Each generation of THV can either match all or only a portion of the overall
Reynolds stress tensor. The given Reynolds stress tensor in equation (78) is
multiplied by a factor, $b_{m}$, based on the ratio of turbulent kinetic
energy (TKE) contained in that specific generation,$TKE_{m}$, to the total
turbulent kinetic energy, $TKE$.

$\displaystyle
b_{m}=\frac{TKE_{m}}{TKE}\quad\quad;\quad\quad\sum_{m=1}^{M}b_{m}=1$ (92)

$M$ is the total number of THV generations and $m$ represents the current
generation. Thus, the equation for the target Reynolds stresses for a
particular THV, equation (78), is modified as follows.

$\displaystyle\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{target}=b_{m}\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{given}-\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{existing}$ (93)

The Reynolds stress contribution from the surrounding THV’s, $\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{existing}$, is only calculated over the
THV’s in the same generation as the THV being created.

$\displaystyle\langle
u^{\prime}_{i}u^{\prime}_{j}\rangle_{existing}=\sum_{n=1}^{N_{M}}\langle\widehat{\widetilde{u}}_{i}\widehat{\widetilde{u}}_{j}\rangle_{n}$
(94)

$N_{m}$ is the number of encircling THV’s in the same generation. By only
matching a portion of the TKE, THV’s from multiple generations can exist on
top of each other. This enables energy to be contributed over a wider range of
scales. The number of generations, the radii bounds of the generations, and
the factor, $b_{m}$ all provide an element of control over how the TKE is
distributed over the resolvable lengthscales. More generations and tighter
radii bounds for the generations allow for a finer tuning of the spectra.

###  4.9 Near-wall THV Stretching

A spherical Triple Hill’s Vortex created near a wall does not agree with the
observed elongated flow structures observed in the near-wall region. A
stretching operation is introduced to produce angled, elongated THV’s in the
near-wall region consistent with the observations of Jeong et al. [18] and
Sibilla and Beretta [42]. The physical coordinates system associated with the
center of a near-wall THV is stretched according to a streamwise length scale,
rotated about the spanwise coordinate axis according to a inclination angle,
and rotated about the vertical coordinate axis according to a tilting angle.
Figure 2 shows a diagram illustrating the stretched and rotated THV.

$\left[\begin{array}[]{c}x^{\prime}_{s}\\\ y^{\prime}_{s}\\\
z^{\prime}_{s}\end{array}\right]=\left[\begin{array}[]{ccc}\frac{1}{2}\frac{l_{x}}{2a}&0&0\\\
0&1&0\\\
0&0&1\end{array}\right]\left[\begin{array}[]{ccc}cos(\alpha)&-sin(\alpha)&0\\\
sin(\alpha)&cos(\alpha)&0\\\
0&0&1\end{array}\right]\left[\begin{array}[]{ccc}cos(\beta)&0&-sin(\beta)\\\
0&1&0\\\
sin(\beta)&0&cos(\beta)\end{array}\right]\left[\begin{array}[]{c}x^{\prime}\\\
y^{\prime}\\\ z^{\prime}\end{array}\right]$ (95)

$(x^{\prime},y^{\prime},z^{\prime})$ are the physical coordinates associated
with the center of a THV and $(x_{s}^{\prime},y_{s}^{\prime},z_{s}^{\prime})$
are the stretched and rotated coordinates that are then used to calculate the
velocity components of the THV. $l_{x}/2a$ is the length to diameter ratio and
it is multiplied by $1/2$ in equation (95) because the spherical THV is
stretched equally in the positive and negative streamwise direction. The
inclination angle is $\alpha$ and $\beta$ is the tilting angle. Jeong et al.
[18] found that the inclination angle was constant and that the tilting angle
varied between $-\beta_{max}$ and $\beta_{max}$ depending on the local
vertical component of vorticity. Here,

$\beta=\epsilon^{\beta}\beta_{max}$ (96)

where $\epsilon^{\beta}$ is a random number defined such that
$\langle\epsilon^{\beta}\rangle=0$ and
$\langle\epsilon^{\beta}\thinspace{}^{2}\rangle=1$. The constant stretching
and rotation procedure allows for the continued satisfaction of the
divergence-free condition for a THV.

a)xy$\alpha$b)xz$\beta$ Figure 2: A diagram of the stretched THV imposed in
the near-wall region: a) inclination angle; b) tilting angle.

##  5 Results

Three test cases were considered: convection of an isolated Triple Hill’s
Vortex, homogeneous isotropic turbulence, and turbulent channel flow.

###  5.1 Single Triple Hill’s Vortex

The convection of a single Triple Hill’s Vortex was investigated. The THV was
generated at the inlet plane and convected downstream by a uniform mean flow.
The dimensions of the domain were $12a\times 6a\times 6a$ in the streamwise,
vertical, and spanwise directions and it was discretized using a uniform
Cartesian grid with $120\times 60\times 60$ grid points. Far field boundary
conditions were used on the vertical and spanwise boundaries.

Contours of the velocity magnitude on planes through the center of a single
convected THV are shown in Figure 3. The contours clearly shows the spherical,
and thus symmetric, nature of the THV. The xy- and xz-planes also show the
stretching of the THV by the convecting flow.

xy-plane) xz-plane) yz-plane)

Figure 3: Contours of the velocity magnitude on xy (left), xz (middle), and yz
(right) planes through the center of a THV.

Figure 4 depicts the generation in time of a single THV from the inlet. Each
of the frames is a constant time step apart. In the first frame, the THV can
be seen just beginning to emerge from the inlet. Notice that there are no
spurious waves originating from the front of the THV. Then, moving forward in
time, the next two frames show the THV as half of it is generated and then as
it just leaves the inlet. In the last frame, the THV has fully released from
the inlet. Notice that the THV passed through the inflow boundary cleanly,
with no spurious waves originating from the back of the THV, which is an
indication that the divergence-free condition is satisfied (otherwise,
spurious waves may be generated in the downstream of the eddy). In the
synthetic turbulence model, analyzed next, when a THV is released from the
inlet, it is no longer acted upon by the THV SEM and the generation process
begins again with a new THV at a new random location.

a) b) c) d)

Figure 4: A single THV being generated at the inlet. Time is increasing with a
constant time step from a) to d) and the inlet is on the left side of each
frame.

###  5.2 Homogeneous Isotropic Turbulent Flow

Large Eddy Simulation of isotropic turbulence was performed and the results
were compared to experimental data collected by Comte-Bellot and Corrsin [7,
8]. Comte-Bellot and Corrsin measured the temporal decay of the turbulent
kinetic energy of mesh-generated isotropic turbulence in a wind tunnel. For
this case, a study of the spatial decay of the turbulent kinetic energy was
used to validate the THV synthetic eddy method.

The numerical domain corresponds to the contracted test section of the wind
tunnel downstream of the mesh. The dimensions of the domain were
$120L_{mesh}\times 10L_{mesh}\times 10L_{mesh}$ in the streamwise, vertical,
and spanwise directions, where $L_{mesh}$ is the size of the experimental
mesh, which equals $0.0508$m [7]. In both the vertical and spanwise
directions, the numerical domain size was chosen to be ten experimental mesh
sizes with periodic boundary conditions, instead of using the entire cross-
section of the wind tunnel test section. A uniform Cartesian grid consisting
of $1200\times 100\times 100$ grid points was used. The increased artificial
viscosity region was added for the last $20L_{mesh}$ in the streamwise
direction at the outflow. A $M_{\infty}=0.05$ uniform mean flow was imposed at
the inflow, which is slightly larger than the experimental mean velocity
($12.7$ m/s); this is because the numerical algorithm is not able to handle
very low Mach number flows. The non-dimensional streamwise turbulence
intensity of the synthetic eddies imposed at the inlet corresponds to level
reported at $42L_{mesh}$ in Comte-Bellot and Corrsin [8]. The THV’s were
broken into ten generations, as given in Table 3, where the number of THV’s in
each generation was determined using equation (91). The maximum radius of a
THV was in the order of the experimental mesh size, $1.5L_{mesh}$, and the
minimum radius was based on the grid resolution (five to six grid points
across the smallest THV). Each THV generation is given an equal portion of the
TKE, $b_{m}=1/M$.

Generation | Range of Radii | Number of THV’s
---|---|---
$1$ | $1.375\ L_{mesh}<a<1.5\ L_{mesh}$ | $8$
$2$ | $1.25\ L_{mesh}<a<1.375\ L_{mesh}$ | $12$
$3$ | $1.125\ L_{mesh}<a<1.25\ L_{mesh}$ | $22$
$4$ | $1.0\ L_{mesh}<a<1.125\ L_{mesh}$ | $28$
$5$ | $0.875\ L_{mesh}<a<1.0\ L_{mesh}$ | $36$
$6$ | $0.75\ L_{mesh}<a<0.875\ L_{mesh}$ | $48$
$7$ | $0.625\ L_{mesh}<a<0.75\ L_{mesh}$ | $68$
$8$ | $0.5\ L_{mesh}<a<0.625\ L_{mesh}$ | $104$
$9$ | $0.375\ L_{mesh}<a<0.5\ L_{mesh}$ | $180$
$10$ | $0.25\ L_{mesh}<a<0.375\ L_{mesh}$ | $400$
Table 3: THV generations for the homogeneous cases

Isosurfaces of Q-criterion are shown in Figure 5 . The THV’s are generated at
the inlet on the left and are convected downstream to the right. It can be
seen from the isosurfaces that the shape of the generated THV’s change from
larger and compact structures near the inflow to smaller and stretched eddies
farther downstream. In the proximity to the inflow boundary, one can hardly
notice any transition from artificial to realistic turbulence. Comparing the
regions near the inlet and outlet, there has been a noticeable dissipation of
the eddies.

Figure 5: Isosurfaces of Q-criterion for homogeneous isotropic turbulent flow
(inlet is on the left).

Two-point spatial velocity correlations and non-dimensionalized turbulent
kinetic energy spectrum of the synthetic fluctuations at the inlet plane are
shown in Figure 6 along with experimental data collected at $42L_{mesh}$ by
Comte-Bellot and Corrsin [8]. Looking at the spatial correlations in Figure
6(a), very good agreement is found for both the longitudinal and transverse
correlations. Other than providing insight into what the maximum radius of the
largest sized THV generation should be, the spatial correlations are not used
as a control on the creation of the THV’s. The energy spectrum in Figure 6(b)
also shows the same good agreement with the experimental data. Like with the
spatial correlations, the only influence exerted over the energy spectrum of
the created THV’s is through the definition of the THV generations. While each
THV generation is given an equal fraction of the total TKE, the wavenumber
range inherently present in each THV generation is different. Because each THV
is a coherent structure, each one contributes to the energy spectrum over the
range of scales from the overall size of the THV to the smallest scale
resolved by the grid.

a) b)

Figure 6: Two-point spatial correlations (a) and turbulent kinetic energy
spectrum as a function of wavenumber (b), captured at the inlet, compared with
the experimental data collected by Comte-Bellot and Corrsin [8]. a) —–)
numerical transverse; - - -) numerical longitudinal; O) experimental
transverse; X) experimental longitudinal. b) O) numerical; —–) experimental.

Figure 7(a) shows the spatial decay of the turbulent kinetic energy for the
numerical results compared with the $k^{-1.25}$ power law fit to experimental
data collected by Comte-Bellot and Corrsin [7]. The THV SEM was able to
reproduce the rate of turbulent kinetic energy decay seen in the experimental
data. The lack of a recovery region just after the inlet further reinforces
that the THV SEM is divergence-free. The slower decay of TKE present near the
inlet is caused by the larger THV’s at the inlet containing more of the TKE.
This is seen in the energy spectrum at the inlet shown in Figure 6(b), where
the energy of the lower wavenumbers is slightly greater than the experimental
results and the energy of the higher wavenumbers is slightly less than the
experimental results. Even though the filtering required by the numerical
method for stability was keep to a minimum, this added dissipation contributed
to the slightly quicker TKE decay.

The downstream development of the velocity derivative skewness is shown in
Figure 7(b) along with LES data from Jarrin et al. [16] of spatially decaying
homogeneous isotropic turbulence at a similar integral length scale Reynolds
number. At the inlet, the skewness of the synthetic fluctuations generated
using THV’s is zero, which is expected because the THV’s are symmetric. Moving
shortly downstream, the skewness rapidly increases in magnitude for the flow
generated by the THV SEM, while the skewness of the original SEM flow requires
a much longer development length. Even though the Gaussian and tent function
synthetic eddies of the original SEM and the THV’s are all symmetric, the more
physically based THV allows for the synthetic field to transition to realistic
turbulence much quicker than the more abstract Gaussian shape.

a) b)

Figure 7: Turbulent kinetic energy (a) and velocity derivative skewness (b) in
the streamwise direction compared with the experimental data collected by
Comte-Bellot and Corrsin [7] and LES data from Jarrin et al. [16] : —–)
numerical; O) experimental; - - -) LES utilizing the original SEM.

A comparison of the turbulent kinetic energy decay and the longitudinal two-
point spatial correlation at the inlet (already presented in Figure 7(a) and
Figure 6(a), respectively) with the LES data from Dietzel et al. [11] is shown
in Figure 8. Dietzel et al. [11] simulated the temporal decay of homogeneous
isotropic turbulence using the following three synthetic turbulence methods to
generate synthetic initial conditions: the digital filtering method proposed
by Klein et al. [25] and improved by Kempf et al. [23], the diffusion method
of Kempf et al. [22], and the Fourier method of Billson et al. [4] and
Davidson [10]. The decay of the TKE for the THV SEM exhibits similar behavior
to the divergence-free (for isotropic turbulence) Fourier method, whereas the
digital filtering and diffusion methods both show the significant dissipation
immediately after synthetic fluctuations are imposed that is characteristic of
methods that do not satisfy the divergence-free condition.

a) b)

Figure 8: Turbulent kinetic energy in the streamwise direction (a) and
longitudinal two-point spatial correlation at the inlet compared with the
experimental data collected by Comte-Bellot and Corrsin [7] and LES data from
Dietzel et al. [11]: black) THV SEM; O) experimental; blue) digital filtering
method [25, 23]; red) diffusion method [22]; yellow) Fourier method [4, 10].

Two-point spatial velocity correlations and non-dimensionalized turbulent
kinetic energy spectrum $56L_{mesh}$ downstream of the inlet plane are shown
in Figure 9 along with experimental data collected at the corresponding
downstream location by Comte-Bellot and Corrsin [8]. Again, very good
agreement with the experimental results is seen for both the spatial
correlations in Figure 9(a) and the energy spectrum in Figure 9(b). The well-
modeled synthetic fluctuations imposed by the THV SEM at the inlet quickly
developed into realistic turbulence and correctly reproduced downstream
statistics.

a) b)

Figure 9: Two-point spatial correlations (a) and turbulent kinetic energy
spectrum as a function of wavenumber (b), captured $56L_{mesh}$ downstream of
the inlet, compared with the experimental data collected by Comte-Bellot and
Corrsin [8]. a) —–) transverse; - - -) longitudinal; O) experimental
transverse; X) experimental longitudinal. b) O) numerical; —–) experimental.

###  5.3 Turbulent Channel Flow

To further test the capability of the THV SEM to model anisotropic non-
homogeneous turbulence, a turbulent channel flow case was considered. Reynolds
stress tensor mean velocity profiles from the DNS of a fully turbulent channel
at $Re_{\tau}=395$ collected by Moser et al. [33] were imposed at the inlet
plane. The dimensions of the domain were $15\delta\times 2\delta\times
3\delta$ in streamwise, vertical, and spanwise directions, where $\delta$ is
the channel half height. A Cartesian grid composed of $200\times 100\times
100$ grid points was used. The grid was uniformly spaced in the streamwise and
spanwise directions, and stretched in the vertical direction in order to
increase resolution at the walls. This allowed for grid spacings of $\Delta
x^{+}=30$, $\Delta z^{+}=12$, $\Delta y_{min}^{+}=1$, and $\Delta
y_{max}^{+}=17.5$. No-slip wall boundary conditions were used on the bottom
and top boundaries, and periodic boundary conditions were applied in the
spanwise direction. The increased artificial viscosity region was added for
the last $3\delta$ in the streamwise direction at the outflow. $6218$ THV’s
were imposed at any one instant in time on to the mean flow at the inflow
boundary. The maximum radius allowed was $0.6\delta$, while the minimum radius
was $0.006\delta$. The THV’s were split into five generations, as given in
Table 4, with the first generation being the elongated eddies described in
Section 4.9 and the fifth generation being clustered very near the walls. This
clustering of THV’s replicates the smallest turbulent structures very near the
wall. Because the energy spectra are not global quantities and a relationship
between the factor $b_{m}$ and distance from the wall was not readily
apparent, $b_{m}=1$.

Generation | Range of Radii | Range of y | Number of THV’s
---|---|---|---
$1$ | $0.24\ \delta<a<0.32\ \delta$ | $0.64\ \delta<|y|<0.68\ \delta$ | $30$
$2$ | $0.4\ \delta<a<0.6\ \delta$ | $0<|y|<(\delta-a)$ | $8$
$3$ | $0.1\ \delta<a<0.2\ \delta$ | $0<|y|<(\delta-a)$ | $180$
$4$ | $0.02\ \delta<a<0.1\ \delta$ | $0<|y|<(\delta-a)$ | $4000$
$5$ | $0.006\ \delta<a<0.01\ \delta$ | $0.9\ \delta<|y|<0.99\ \delta$ | $2000$
Table 4: THV generations for the channel case, where $-\delta<y<\delta$.

Contours of streamwise velocity and vorticity magnitude at the inflow plane
are shown in Figure 10. These contours show the synthetic inflow, the
superposed distorted THV’s, that is imposed by the THV SEM. Notice how the
concept of the generations of THV’s manifests itself in the synthetic inflow
with a few large eddies surrounded by ever smaller structures. Moving from the
center of the channel to the walls, the size the imposed eddies decreases. The
vorticity magnitude contour especially illuminates the effect of the clustered
eddies near the wall, where the vorticity magnitude is larger as compared to
the center of the channel. The development of the elongated turbulent
structures stretching downstream from the walls can be seen from the
Q-criterion isosurfaces in Figure 11. While the synthetic eddies are released
into the the domain and quickly evolve into realistic turbulent structures,
there is still a small region near the inlet where the non-stretched THV’s of
generations two through five need to develop.

a) b)

Figure 10: Contours for turbulent channel flow at the inlet plane(from the
synthetic turbulence model): a) streamwise velocity; b) vorticity magnitude.

Figure 11: Isosurfaces of Q-criterion for turbulent channel flow.

Vertical profiles of the mean streamwise velocity and Reynolds stresses are
plotted in Figure 12, where averaging in both time and the homogeneous
spanwise directions has been taken. As seen in Figure 12(a), the DNS mean
streamwise velocity profile is well reproduced. When looking at the Reynolds
stresses in Figure 12(b), there is excellent agreement between the heights of
$-0.9\delta$ and $0.9\delta$. This height range corresponds to the effective
region where most of the THV’s are imposed. A portion of the smaller THV’s of
the fourth generation and the fifth clustered generation are imposed at
heights closer to the wall, where their effect is clearly seen in the bulge of
the $\langle u^{\prime}u^{\prime}\rangle$ profile very close to the walls, but
those THV’s are not enough to reproduce the DNS profiles. With more smaller
and smaller THV’s in those nearest wall regions, those Reynolds stress
profiles will be better matched.

a) b)

Figure 12: Vertical profiles at the inlet compared with the DNS data collected
by Moser et al. [33] of: a) mean streamwise velocity; b) Reynolds stresses;
—–) numerical; - - -) DNS data.

Two-point spanwise spatial velocity correlations of the synthetic fluctuations
at the inlet plane and at three different heights from the wall are shown in
Figure 13 along with the corresponding DNS data from Moser et al. [33]. There
is generally good agreement at spanwise distances less than $0.4\delta$.
Closer to the wall, Figure 13(a), the synthetic fluctuations are composed of
the smaller THV’s. They are only locally correlated, so the spatial
correlations quickly go to zero. Moving away from the wall, Figure 13(b), and
towards the center of the channel, Figure 13(c), the largest THV’s become more
present. This is seen in the higher correlations at larger spanwise distances,
which might be a hint that the largest THV’s are too large in size. As a
general trend across all three heights from the wall, the THV’s do not
reproduce negative spatial correlations, especially at larger spanwise
distances. The transverse correlation $0.2\delta$ from the wall does dip
slightly below zero, which was also seen for the transverse correlation at the
inlet of the homogeneous isotropic turbulence (Figure 6(a)), but it does not
reproduce the negative correlation.

a) b) c)

Figure 13: Two-point spanwise spatial correlations, captured at the inlet and
at three different heights from the wall , compared with the DNS data
collected by Moser et al. [33]. blue) transverse; red) longitudinal; - - -)
DNS transverse; -.-.-) DNS longitudinal.

Figure 14 shows the non-dimensionalized spanwise turbulent kinetic energy
spectra of the synthetic fluctuations at the inlet plane and at three
different heights from the wall along with the corresponding DNS data from
Moser et al. [33]. There is good agreement at all three heights in the higher
wavenumbers, with better agreement when moving towards the wall. Unlike for
homogeneous turbulence where each THV generation matches an equal fraction of
the TKE at any point, because of the non-homogeneous nature of the wall normal
direction, each individual THV matches the total TKE at any point. So, as seen
in the darker large circular areas in the inlet vorticity contour of Figure
10(b), the generations of smaller THV’s are not created inside as many of the
largest THV’s which does not allow for the same sort of fine tuning of the
energy spectrum as for the homogeneous turbulence. This causes the general
under-prediction of the peaks of the energy spectra. Looking back again to
Figure 10(b) and considering the larger THV present towards the middle top of
the contour, these larger THV’s match the Reynolds stresses away from the
walls, but can contribute disproportionately to energy at the lowest
wavenumbers. This is why the spectra of the synthetic fluctuations levels off
at the lowest wavenumbers.

a) b) c)

Figure 14: Turbulent kinetic energy spectra, captured at the inlet and at
three different heights from the wall , compared with the DNS data collected
by Moser et al. [33]. blue) transverse; red) longitudinal; - - -) DNS
transverse; -.-.-) DNS longitudinal.

The development of the friction coefficient at the walls is presented in
Figure 15. Results from the THV SEM are compared with data reported by Jarrin
et al. [16] using the original SEM (OSEM), Poletto et al. [36] using a
divergence-free SEM (DFSEM), and Skillen et al. [43] using a SEM with improved
eddy positioning and amplitude calculation (SSEM). The THV SEM shows the same
spike in friction coefficient as the DFSEM and the SSEM immediately downstream
of the inlet as the synthetic fluctuations begin to interact with the wall.
$C_{F}$ then quickly settles to within $\pm 2\%$ of the fully-developed value
by $2\delta$ downstream. The THV SEM shows a significant reduction in the
development distance over the DFSEM and a similar development distance with
the SSEM.

Figure 15: Downstream development of the friction coefficient compared with
the DNS data collected by Moser et al. [33] and three different synthetic eddy
methods: black) THV SEM; blue) Original SEM (OSEM) of Jarrin et al. [16]; red)
Divergence-Free SEM (DFSEM) of Poletto et al. [36]; yellow) SEM (SSEM) of
Skillen et al. [43]; - - -) DNS of Moser et al. [33].

Vertical profiles of the turbulent kinetic energy and Reynolds shear stress at
four streamwise locations are plotted in Figure 16. Both the TKE and Reynolds
shear stress have recovered their given values away from the wall
($-0.5\delta<y<0.5\delta$) by one channel height downstream of the inlet
($x=2\delta$). The decrease in both quantities near the wall is because the
grid in the streamwise direction is too coarse to properly resolve the smaller
THV’s in the fourth and fifth THV generations. Returning to Figure 10, the
larger THV’s away from the wall do influence this nearer-wall region; which is
why the TKE and Reynolds shear stress both decrease to values consistent with
the Reynolds stresses at the center of the larger THV’s. By $10\delta$
downstream of the inlet, the desired values are recovered in the near-wall
region.

a) b)

Figure 16: Vertical profiles at four different streamwise locations compared
with the DNS data collected by Moser et al. [33] of: a) turbulent kinetic
energy; b) Reynolds shear stress; blue) $x=0$ (inlet); red) $x=2\delta$ ;
yellow) $x=4\delta$ ; purple) $x=10\delta$ ; - - -) DNS of Moser et al. [33].

Two-point streamwise spatial velocity correlations at two different heights
from the wall, calculated from the inlet, are presented in Figure 17. The
influence of stretched THV’s defined in Section 4.9 can clearly be seen in the
differences between the longitudinal correlation near the wall in Figure 17(a)
and near the height where the stretched THV’s are clustered in Figure 17(b).
Close to the wall, only the smaller spherical THV’s are imposed. Although
those small THV’s give good agreement at distances less than $0.5\delta$ from
the inlet, their spherical nature is unable to model the larger downstream
coherence. This is in stark contrast to the longitudinal correlation
$0.5\delta$ away from the wall in Figure 17(b). The spatial correlation is in
agreement much farther downstream. This is because the stretched THV’s are
being created at the inlet at this general height. The stretched THV’s
introduce larger correlated scales which provides a better model of the
physical structures.

a) b)

Figure 17: Two-point streamwise spatial correlations, captured at two
different heights from the wall , compared with the DNS data collected by
Moser et al. [33]. blue) transverse; red) longitudinal; - - -) DNS transverse;
-.-.-) DNS longitudinal.

The downstream evolution of the two-point spanwise transverse spatial velocity
correlation is shown in Figure 18 and compared with DNS data from Moser et al.
[33] at two representative heights. Again, at the inlet, the THV’s are well
correlated at small distances, but they are unable to reproduce the larger
negative correlations. As the synthetic fluctuations are transformed by the
Navier-Stokes equations, the correct negative correlations quickly develop by
$2\delta$ downstream of the inlet. This rapid development of the physical
coherence of the synthetic fluctuations agrees with the rapid settling of the
friction coefficient, the turbulent kinetic energy, and the Reynolds shear
stress.

a) b)

Figure 18: Downstream evolution of the two-point spanwise spatial
correlations, captured at two different heights from the wall, compared with
the DNS data collected by Moser et al. [33]. blue) $x=0$ (inlet); red)
$x=2\delta$ ; yellow) $x=4\delta$ ; - - -) DNS transverse.

##  6 Conclusion

The Triple Hill’s Vortex was proposed as a superposition of three Hill’s
vortices with their axes perpendicular to each other. The proposed new
synthetic eddy was applied in the framework of the synthetic eddy method in
order to generate synthetic turbulent inflow velocity fields that satisfied
the divergence-free condition and matched given Reynolds stress profiles.
Simulation of homogeneous isotropic turbulence and turbulent channel flow were
both able to reproduce given Reynolds stress tensors. The transition from
artificial to realistic turbulence in the proximity to the inflow boundary was
found to be small in all test cases that were considered. Excellent agreement
between the turbulent kinetic energy spectrum and the two-point spatial
correlations were found for the homogeneous case both at the inlet and far
downstream. The spatial decay of the turbulent kinetic energy for isotropic
homogeneous turbulence was shown to be in agreement with experimental data of
isotropic turbulence in a wind tunnel and the synthetic fluctuations quickly
developed skewness downstream of the inlet.

For the channel flow, Reynolds stress profiles taken from DNS data were able
to be reproduced away from the wall, but the distribution of Triple Hill’s
Vortices was not sufficient to provide good agreement very close to the walls.
The recovery of the friction coefficient was shown to be as just quick as for
one of the latest synthetic eddy methods. Although good agreement with the DNS
data was found at the inlet for the spanwise two-point correlations, recovery
of the correct correlations occurred quickly downstream.

For the future, the THV SEM needs to be expanded to allow for the use of a
wider range of arbitrary turbulent statistics as inputs. The near-wall
stretching and the uniform division of the total TKE across the THV
generations both showed that rudimentary indirect control over the spatial
correlations and energy spectra can be exerted. Through a systematic approach,
general relationships might be able to be formulated to connect arbitrary
energy spectra to how the target TKE for each THV generation is allocated and
to connect arbitrary spatial correlations to the stretching of the THV’s. The
intrinsic coherence of a synthetic eddy, where all of the resolved scales
smaller than the size of the eddy are contained within the eddy, introduces an
additional layer of complexity that could hinder the formulation of the
general relationships described above. To introduce non-zero higher order
statistics (for example, skewness or flatness) into the final synthetic
turbulence field, the streamfunction for the THV could be distorted such that
the THV is no longer symmetric while preserving the divergence-free condition.

## Funding Sources

Effort sponsored by the Engineering Research & Development Center under
Cooperative Agreement number W912HZ-17-C-0021. The views and conclusions
contained herein are those of the authors and should not be interpreted as
necessarily representing the official policies or endorsements, either
expressed or implied, of the Engineering Research & Development Center.

## Disclaimer

Reference herein to any specific commercial company, product, process, or
service by trade name, trademark, manufacturer, or otherwise, does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or the Department of the Army (DoA). The
opinions of the authors expressed herein do not necessarily state or reflect
those of the United States Government or the DoA, and shall not be used for
advertising or product endorsement purposes.

##  Appendix A Triple Hill’s Vortex Synthetic Eddy Method Outline

This appendix serves as an outline for a general implementation of the THV SEM
at an inflow plane with static Reynolds stresses.

Before time marching,

  1. 1.

Set the radius ranges and populations of the THV generations. Also, set which
generations are stretched or clustered is applicable.

  2. 2.

Set the target Reynolds stress profiles

  3. 3.

Create all of the THV’s. For each THV:

    1. (a)

Assign a random radius within the bounds of the generation size

    2. (b)

Assign a random center location within the bounds of the generation clustering
location, if applicable

    3. (c)

Calculate the three amplitudes

      1. i.

Calculate the Reynolds stress contribution of the THV’s already present on the
inflow plane at the location of the new THV

      2. ii.

Modify the target Reynolds stresses at the center of the new THV with the
calculated contribution Reynolds stresses

      3. iii.

Calculate the principal Reynolds stresses and eigenvectors

      4. iv.

If any of the principal Reynolds stresses are negative, assign a new random
center location and recalculate the amplitudes at the new location. Negative
principal Reynolds stresses indicate that the target Reynolds stresses are
already satisfied.

    4. (d)

Multiply each of the principal Reynolds stresses by an independent random
number

During time marching within the calculation of boundary conditions,

  1. 1.

Check if current THV’s need to be released (if they are over two radii
downstream of the inflow plane, for example)

  2. 2.

Create new THV’s to replace released THV’s using the above creation steps

  3. 3.

Calculate the velocity contribution from each eddy at each grid point on the
inflow plane (contributions of THV’s a certain distance away, for example two
radii, from a grid point are assumed to be zero)

  4. 4.

Add the fluctuating velocity to a given mean flow and impose the combination
at the inflow

## References

  * [1] Batten, P., Goldberg, U., and Chakravarthy, S.: Interfacing Statistical Turbulence Closures with Large-Eddy Simulation. AIAA Journal 42, 3, 485-492 (2004)
  * [2] Bechara, M.F., Christophe, B., Lafon, P. and Candel, S.: A numerical technique for trailing edge acoustic scattering problems. AIAA Journal, 32, 455-463 (1994)
  * [3] Ben-Nasr, O., Hadjadj, A., Chaudhuri, A., Shadloo, M.S.: Assessment of subgrid-scale modeling for large-eddy simulation of a spatially-evolving compressible turbulent boundary layer. Comput. Fluids 151, 144-158 (2017)
  * [4] Billson, M., Eriksson, L,-E., and Davidson, L.: Jet Noise Prediction Using Stochastic Turbulence Modeling. AIAA Paper 2003-3282, (2003)
  * [5] Billson, M., Eriksson, L,-E., and Davidson, L.: Modeling of Synthetic Anisotropic Turbulence and Its Sound Emission. AIAA Journal 2004-2857, (2004)
  * [6] Butcher, J.C.: Numerical methods for ordinary differential equations, John Wiley and Sons Inc, (2003)
  * [7] Comte-Bellot, G. and Corrsin, S.: The use of a contraction to improve the isotropy of grid-generated turbulence. J. Fluid Mech. 25, 657-682 (1966)
  * [8] Comte-Bellot, G. and Corrsin, S.: Simple Eulerian time correlation of full and narrow-band velocity signals in grid-generated, isotropic turbulence. J. Fluid Mech. 48, 273-337 (1971)
  * [9] Davidson, L. and Billson, M.: Hybrid LES-RANS using synthesized turbulent fluctuations for forcing in the interface region. Int. J. Heat Fluid Fl. 27, 1028-1042 (2006)
  * [10] Davidson L.: Using isotropic synthetic fluctuations as inlet boundary conditions for unsteady simulations. Advances and Applications in Fluid Mechanics 1, 1, 1-35 (2007)
  * [11] Dietzel, D., Messig, D., Piscaglia, F., Montorfano, A., Olenik, G., Stein, O.T., Kronenburg, A., Onorati, A., and Hasse, C.: Evaluation of scale resolving turbulence generation methods for Large Eddy Simulation of turbulent flows. Comput. Fluids 93, 116-128 (2014)
  * [12] Druault, P., Lardeau, S., Bonnet, J.-P., Coiffet, F., Delville, J., Lamballais, E., Largeau, J.F. and Perret, L.: Generation of three-dimensional turbulent inlet conditions for large-eddy simulation. AIAA Journal 42, 447-456 (2004)
  * [13] Fung, J., Hunt, J., Malik, N., and Perkins, R.: Kinematic simulation of homogeneous turbulence by unsteady random Fourier modes. J. Fluid Mech. 236, 281-318 (1992)
  * [14] Hadjadj, A., Ben-Nasr, O., Shadloo, M.S., and Chaudhuri, A., 2015: Effect of wall temperature in supersonic turbulent boundary layers: A numerical study. Int. J. Heat Mass Transf. 81, 426-438 (2015)
  * [15] Hill, M.J.M.: On a spherical vortex. Philosophical Transactions of the Royal Society London 185, 213 (1894)
  * [16] Jarrin,N., Benhamadouche, S., Laurence, D. and Prosser, R.: A synthetic-eddy method for generating inflow conditions for large-eddy simulations. Int. J. Heat Fluid Fl. 27, 585-593 (2006)
  * [17] Jarrin, N., Prosser, R., Uribe, J.-C., Benhamadouche, S., and Laurence, D.: Reconstruction of turbulent fluctuations for hybrid RANS/LES simulations using a Synthetic-Eddy method. Int. J. Heat Fluid Fl. 30, 435-442 (2009)
  * [18] Jeong, J., Hussain, F., Schoppa, W., and Kim, J.: Coherent structures near the wall in a turbulent channel flow. J. Fluid Mech. 332, 185-214 (1997)
  * [19] Johansson, P.S. and Andersson, H.I.: Generation of inflow data for inhomogeneous turbulence. Theor. Comput. Fluid Dyn. 18, 371-389 (2004)
  * [20] Kaltenbach, H.-J., Fatica, M., Mittal, R., Lund, T.S. and Moin, P.: Study of flow in a planar asymmetric diffuser using large-eddy simulation. J. Fluid Mech. 390, 151-185 (1999)
  * [21] Keating, A., Piomelli, U., Balaras, E., Kaltenbach, H.-J.: A priori and a posteriori tests of inflow conditions for large-eddy simulation. Phys. Fluids 16, 4696-4712 (2004)
  * [22] Kempf, A.M., Klein, M., and Janicka, J.: Efficient Generation of Initial- and Inflow-Conditions for Transient Turbulent Flows in Arbitrary Geometries. Flow Turbul. Comb. 81 449-470 (2005)
  * [23] Kempf, A.M., Wysock, S., and Pettit, M.: An efficient, parallel low-storage implementation of Klein;s turbulence generator for LES and DNS. Comput. Fluids 60, 58-60 (2012)
  * [24] Kennedy, C.A. and Carpenter, M.H.: Comparison of several numerical methods for simulation of compressible shear layers. NASA Technical Report NASA-97-TP3484 (1997)
  * [25] Klein, M., Sadiki, A. and Janicka, J.: A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations. J. Comput. Phys. 86, 652-665 (2003)
  * [26] Kobayashi, H.: The subgrid-scale models based on coherent structures for rotating homogeneous turbulence and turbulent channel flow. Phys. Fluids 17, 045104 (2005)
  * [27] Kobayashi, H., Ham, F., and Wu, X.: Application of a local SGS model based on coherent structures to complex geometries. Int. J. Heat Fluid Fl. 29, 640-653 (2008)
  * [28] Kraichnan, R.H.: Diffusion by a random velocity field. Phys. Fluids 13, No. 1, 22-31 (1970)
  * [29] Li, A., Ahmadi, G., Bayer, R. and Gaynes, M.: Aerosol particle deposition in an obstructed turbulent duct flow. J. Aerosol Sci. 25, 91-112 (1994)
  * [30] Lund, T.S., Wu, X. and Squires, D.: Generation of turbulent inflow data for spatially developing boundary layer simulations. J. Comput. Phys. 140, 233-258 (1998)
  * [31] Martin, M.P., Piomelli, U., and Chandler, G.V.: Subgrid-Scale Models for Compressible Large-Eddy Simulations. Theor. Comput. Fluid Dyn. 13, 361-376 (2000)
  * [32] Maxey, M.: The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields. J. Fluid Mech. 174, 441 (1987)
  * [33] Moser, R.D., Kim, J., and Mansour, N.N.: Direct numerical simulation of turbulent channel flow up to $Re_{\tau}=590$. Phys. Fluids 11, 4, 943-945 (1999)
  * [34] Pamiés, M., Weiss, P.-É., Deck, S., and Sagaut, P.: Generation of synthetic turbulent inflow data for large eddy simulation of spatially evolving wall-bounded flows. Phys. Fluids 21, 045103 (2009)
  * [35] Poletto, R., Revell, A., Craft, T., and Jarrin, N.: divergence-free synthetic eddy method for embedded les inflow boundary conditions. Seventh International Symposium On Turbulence and Shear Flow Phenomena (TSFP-7), Ottawa, (2011)
  * [36] Poletto, R., Craft, T., and Revell, A.: A New divergence-free Synthetic Eddy Method for the Reproduction of Inlet Flow Conditions for LES. Flow Turbul. Comb. 91, 519-539 (2013)
  * [37] Roidl, B., Meinke, M., and Schröder, W.: A reformulated synthetic turbulence generation method for a zonal RANS-LES method and its application to zero-pressure gradient boundary layers. Int. J. Heat Fluid Fl. 44, 28-40 (2013)
  * [38] Rosales, C. and Meneveau, C.: A minimal multiscale Lagrangian map approach to synthesize non-Gaussian turbulent vector fields. Phys. Fluids 18, 075104 (2006)
  * [39] Rosales, C. and Meneveau, C.: Anomalous scaling and intermittency in three-dimensional synthetic turbulence. Phys. Rev. E 78, 016313 (2008)
  * [40] Sescu, A. and Hixon, R.: Towards Low-Noise Synthetic Turbulent Inflow Conditions for Aeroacoustic Calculations. Int. J. Numer. Methods Fluids 73, 1001-1010 (2013)
  * [41] Shur, L.M., Spalart, P.M., Strelets, M.K. and Travin, A.K.: Synthetic Turbulence Generators for RANS-LES Interfaces in Zonal Simulations of Aerodynamic and Aeroacoustic Problems. Flow Turbul. Comb. 93, 63-92 (2014).
  * [42] Sibilla, S. and Beretta, C.P.: Near-wall coherent structures in the turbulent channel flow of a dilute polymer solution. Fluid Dyn. Res. 37, 183-202 (2005)
  * [43] Skillen, A., Revell, A., and Craft, T.: Accuracy and efficiency improvements in synthetic eddy methods. Int. J. Heat Fluid Fl. 62, 386-394 (2016)
  * [44] Smirnov, A., Shi, S. and Celik, I.: Random flow generation technique for large eddy simulations and particle-dynamics modeling. J. Fluids Eng. 123, 359-371 (2001)
  * [45] Spalart, P.R.: Direct simulation of a turbulent boundary layer up to $Re_{\theta}=1410$. J. Fluid Mech. 187, 61-98 (1988)
  * [46] Tam, C.K.W. and Webb, J.C.: Dispersion-relation-preserving finite difference schemes for Computational Aeroacoustics. J. Comput. Phys. 107, 262-281 (1993)
  * [47] Walchshofer, C., Steiner, H., and Günter, B.: Robust Outflow Boundary Conditions for Strongly Buoyant Turbulent Jet Flames. Flow Turbul. Comb. 86, 713-734 (2011)
  * [48] Yoshisawa, A: Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling. Phys. Fluids 29, No. 7, 2152-2164 (1986)
  * [49] Zang, T.A., Dahlburg, R.B., and Dahlburg, J.P.: Direct and large-eddy simulations of three-dimensional compressible Navier-Stokes turbulence. Phys. Fluids A 4, 1, 127-140 (1992)