appendix - added chapter on the synthetic turbulence generator

    framesep=2\fboxsep,

    framerule=0.01pt,

]{python}{data/lbm/lattice_weights.py}

\chapter{Synthetic Turbulence Generator}

\label{appendix:turbgen}

The procedure for generating isotropic synthetic turbulence is based on \cite{davidson2006,davidson2007,davidson2008,davidson:lecture}.

It can be used to generate velocity fluctuations for an initial condition $\vec v(\vec x, 0)$ for $\vec x \in \Omega$, as well as for an inflow boundary condition $\vec v_{\mathrm{in}}(\vec x, t)$ for $\vec x \in \Gamma_{\mathrm{in}} \subset \partial \Omega$, $t \ge 0$.

The velocity field $\vec v(\vec x, t)$ is first decomposed as

\begin{equation}

    \vec v(\vec x, t) = \overline{\vec v}(\vec x) + \vec v'(\vec x, t),

\end{equation}

where $\overline{\vec v}(\vec x)$ is the mean (time-averaged) field and $\vec v'(\vec x, t)$ is the velocity fluctuation.

The mean field $\overline{\vec v}$ can be set as desired for the initial or boundary condition and the fluctuation $\vec v'$ must be generated.

The algorithm described below is implemented in the \C++ language with CUDA acceleration and available as part of the TNL library \cite{oberhuber:2021tnl,TNL:documentation}.

\section{Computing the fluctuations}

Turbulent fluctuations consist of a wide range of scales and their analysis is based on using Fourier series.

A general form of a fluctuating velocity field $\vec v'$ at a fixed time level with zero average ($\overline{\vec v'} = 0$) can be expressed as \cite[Eq. (27.3)]{davidson:lecture}

\begin{equation} \label{eq:turbgen:Fourier series}

    \vec v'(\vec x) = 2 \sum\limits_{n = 1}^{N_{\mathrm{modes}}} \hat{v}^n \cos(\vec \kappa^n \cdot \vec x + \psi^n) \vec \sigma^n,

\end{equation}

where $\hat{v}^n$ and $\psi^n$ are the amplitude and phase of the $n$-th mode, $\vec \sigma^n$ is a unit vector that determines the direction of the $n$-th mode, and $\vec \kappa^n$ is the wave-number vector of the $n$-th mode.

Note that only the first $N_{\mathrm{modes}}$ terms of the full series are considered in the generator.

The wave-number vector $\vec \kappa^n$ is written as

\begin{align} \label{eq:turbgen:kappa vector}

    \kappa_x^n &= \rho^n \sin(\theta^n) \cos(\phi^n), \\

    \kappa_y^n &= \rho^n \sin(\theta^n) \sin(\phi^n), \\

    \kappa_z^n &= \rho^n \cos(\theta^n),

\end{align}

where $\theta^n$ and $\phi^n$ are randomly selected angles and $\rho^n$ is the magnitude of the $n$-th wave number.

In order to have $\nabla \cdot \vec v' = 0$, the direction vector $\vec \sigma^n$ must be orthogonal to the wave-number vector $\vec \kappa^n$ \cite{davidson:lecture}.

Hence, $\vec \sigma^n$ can be written as

\begin{align} \label{eq:turbgen:sigma vector}

    \sigma_x^n &= \cos(\phi^n) \cos(\theta^n) \cos(\alpha^n) - \sin(\phi^n) \sin(\alpha^n), \\

    \sigma_y^n &= \sin(\phi^n) \cos(\theta^n) \cos(\alpha^n) + \cos(\phi^n) \sin(\alpha^n), \\

    \sigma_z^n &= - \sin(\theta^n) \cos(\alpha^n),

\end{align}

where $\alpha^n$ is a randomly selected angle that gives the direction of $\vec \sigma^n$ in the plane orthogonal to $\vec \kappa^n$.

The phase $\psi^n$ in \cref{eq:turbgen:Fourier series} is selected randomly as well.

The magnitudes $\rho^n$ of the discrete wave numbers in \cref{eq:turbgen:kappa vector} are selected by sampling an interval $[\kappa_{\min}, \kappa_{\max}]$ with $N_{\mathrm{modes}}$ points, i.e.

\begin{equation} \label{eq:turbgen:discrete wave number}

    \rho^n = \kappa_{\min} + \left( n - \frac{1}{2} \right) \Delta \kappa

\end{equation}

for $n \in \{1, \ldots, N_{\mathrm{modes}}\}$, where $\Delta \kappa = (\kappa_{\max} - \kappa_{\min}) / N_{\mathrm{modes}}$.

The highest wave number $\kappa_{\max} = 2 \pi / \delta_x$ is set based on the grid spacing $\delta_x$ and the smallest wave number is set as $\kappa_{\min} = \kappa_e / p$, where the factor $p=5$ is chosen as suggested in \cite{davidson:lecture} and $\kappa_e$ is specified below in \cref{eq:turbgen:kappa_e} based on the turbulence model.

The modified von Kármán spectrum \cite{davidson2007,davidson:lecture} is used to model the turbulent energy spectrum.

The amplitude $\hat{v}^n$ of the $n$-th mode in \cref{eq:turbgen:Fourier series} is obtained from

\begin{subequations}

    \begin{align}

        \hat{v}^n &= v_{\mathrm{RMS}} \sqrt{E(\rho^n) \Delta k}, \\

        E(\kappa) &= \frac{c_E}{\kappa_e} \frac{ \left( \frac{\kappa}{\kappa_e} \right)^4 }{ \left[ 1 + \left( \frac{\kappa}{\kappa_e} \right)^2 \right]^{17/6} }

        \exp\left( -2 \left( \frac{\kappa}{\kappa_{\eta}} \right)^2 \right), \\

        c_E &= \frac{4}{\sqrt{\pi}} \frac{\Gamma\left(\frac{17}{6}\right)}{\Gamma\left(\frac{1}{3}\right)} \approx 1.452762113, \\

        \label{eq:turbgen:kappa_e}

        \kappa_e &= \frac{9 \pi c_E}{55 \mathcal L_{\mathrm{int}}}, \quad\quad

        \kappa_{\eta} = \left( \frac{\epsilon}{\nu^3} \right)^{1/4},

    \end{align}

\end{subequations}

where $v_{\mathrm{RMS}} = \sqrt{\frac{2}{3}k}$~[\si{\metre\per\second}] is the turbulent velocity scale (root mean square), $k$~[\si{\metre\squared\per\second\squared}] is the turbulent kinetic energy, $\mathcal L_{\mathrm{int}}$~[\si{\metre}] is the turbulent integral length scale, $\epsilon = \frac{k^{3/2}}{\mathcal L_{\mathrm{int}}}$~[\si{\metre\squared\per\second\cubed}] is the dissipation rate of the turbulent kinetic energy, and $\nu$~[\si{\metre\squared\per\second}] is the kinematic viscosity of the fluid.

The input parameters of the turbulence generator are the quantities $N_{\mathrm{modes}}$, $k$, and $\mathcal L_{\mathrm{int}}$.

For channel flow applications, $\mathcal L_{\mathrm{int}}$ may be estimated as a fraction (e.g., $1/10$ to $1/2$) of the expected boundary layer height \cite{davidson:lecture}.

\section{Introducing time correlation}

To generate velocity fluctuations at multiple discrete time levels $t_n = n \delta_t$, where $n$ is an integer denoting the time level and $\delta_t$ is the time step used in a numerical simulation, additional computations are needed.

Firstly, independent realizations of random fluctuations $\hat{\vec v}'$ are generated for each time level using the procedure described above.

Then, time correlation between the realizations is introduced using an asymmetric time filter

\begin{equation}

    (\vec v')^n = a (\vec v')^{n-1} + b (\hat{\vec v}')^n,

\end{equation}

where $\vec v'$ denotes the time-correlated field, $\hat{\vec v}'$ denotes the time-independent field, superscripts denote the time levels and the coefficients are chosen as $a = \exp(-\Delta t / \mathcal T_{\mathrm{int}})$ and $b = \sqrt{1 - a^2}$, where $\mathcal T_{\mathrm{int}} = \mathcal L_{\mathrm{int}} / V_b$ is set using the Taylor's hypothesis based on the desired bulk velocity $V_b \approx \abs{\overline{\vec v}}$.

The time filter ensures that $\mathcal T_{\mathrm{int}}$ corresponds to the turbulent integral time scale and that the variance of the generated fluctuations is preserved~\cite{davidson:lecture}.

\section{Example}

\Cref{fig:turbgen} shows an example of the synthetic fluctuating velocity field $v'_x$ computed on a unit square discretized by $128 \times 256$ points.

The parameters used in the generator are: grid spacing $\delta_x = 1/\SI{255}{\metre}$, integral length scale $\mathcal L_{\mathrm{int}} = \SI{0.01}{\metre}$, number of discrete modes $N_{\mathrm{modes}} = 3000$, turbulent kinetic energy $k = \SI{e-2}{\metre\squared\per\second\squared}$, and kinematic viscosity $\nu = \SI{1.5e-5}{\metre\squared\per\second}$.

\begin{figure}[h]

    \centering

    \includegraphics[width=\textwidth]{figures/turbgen/turbgen_vx_2.png}

    \caption{First component of synthetic fluctuating velocity field ($v'_x$) generated on a unit square discretized by $128 \times 256$ points.}

    \label{fig:turbgen}

\end{figure}