Synthetic turbulence generator according to Shur et al \(\text{.}\) [49] was implemented by Morsbach and Franke [37] and tested by Matha et al \(\text{.}\) [31] .

Input from panels

For panels (faces), the mesh size \( \vec{h} \) is calculated as:

\[ h_x = 2 \cdot \left| {\mathbf x}_{inner} - {\mathbf x}_{face} \right| \]

with \( {\mathbf x}_{inner} \) as cell center and \( {\mathbf x}_{face} \) as face center coordinates.

The mesh size in y- and z-direction are calculated as:

\[ h_y = \underset{j}{\max} \left| y_0 - y_j \right| \]

with \( y_j\) as y-coordinate of cell node \( j\). Same for z-direction.

The wall distance \(y_n\) is copied from the cell center value and the turbulent length scale \( L_t\) has to be given via input file.

Random vectors

The vector \(\mathbf d\) is constructed as:

\[ \mathbf d = \begin{pmatrix}r(-0.5, 0.5) \\ r(-0.5, 0.5) \\ r(-0.5, 0.5)\end{pmatrix} \]

as long \( \left|{\mathbf d}\right|>1 \) . If \( \left|{\mathbf d}\right|<=1 \) the vector is scaled to 1 and stored. This approach provides the best results, tested by Christian M. The first step to create the vector \({\mathbf \sigma}\) in plane orthogonal to \(\mathbf d\) rotated by angle \(\phi = r(0, 2\pi)\) is to derive arbitrary normal vector

\[ \begin{pmatrix}d_1 \\ d_2 \\ d_3\end{pmatrix} \times \begin{pmatrix}1 \\ 0 \\ 0\end{pmatrix} = \begin{pmatrix}0 \\ d_3 \\ -d_2\end{pmatrix} \]

and rotating it by \(\phi\). This can fail if \(d_2 = d_3 = 0\). The resulting \({\mathbf \sigma}\) is then given as:

\[ {\mathbf \sigma} = \frac{1}{\sqrt{d_2^2 + d_3^2}} \begin{pmatrix} - (d_2^2 + d_3^2) \sin \phi \\ d_1 d_2 \sin \phi + d_3 \cos \phi \\ d_1 d_3 \sin \phi - d_2 \cos \phi \end{pmatrix} \]

For periodic boundaries, the extension formulated by Morsbach and Franke [37] can be applied using the control file command STG_PERIODIC_FIX.

Time-constant variables

The following variables are constant in time, hence they are calculated at the simulation start.

Local variables:

\[k_e = \frac{2 \pi}{l_e} \qquad \text{with} \qquad l_e = \min(2 y_n, C_l L_T) \quad , \quad C_l=3.0\]

\[ k_{\text{cut}} = \frac{2\pi}{l_{cut}} \qquad \text{with} \qquad l_{\text{cut}} = 2 \min\left( \left[\max\left(h_y,h_z,0.3h_\max\right)+0.1y_n\right], h_\max\right) \]

with \( h_{\max} = \max \left( h_x, h_y, h_z\right) \) .

Global variables:

\[ l_{e,\max} = \underset{\mathbf x}\max\ l_e\left(\mathbf x\right) \]

\[ k_{\text{cut,max}} = \underset{\mathbf x}{\max}\ k_{cut}\left(\mathbf x\right) \]

Number of modes is given by the minimum integer \( N\) satisfying:

\[ k^N \geq k_\max = 1.5 \max\left(k_{\text{cut}}\left(\mathbf x\right)\right) \]

with

\[ k^n = k_{\min}\left(1+\alpha\right)^{n-1} \quad , \quad k_{\min}=\beta\frac{2\pi}{l_{e,\max}} \quad , \quad \alpha = 0.01 \quad , \quad \beta=0.5 , \quad n = 1, 2, \dots, N; \]

and differences

\[ \Delta k^n = k^{n+1} - k^n \]

leading to following relation for \( N\):

\[ N = 1 + \frac{\log \frac{1.5 k_\text{cut,max}}{k_\text{min}}}{\log (1 + \alpha)} \]

With known number of modes the uniform random value for angle \(\phi^n\in\left[0,2\pi\right)\), the unit vector \(\mathbf d\) uniformly randomUnitVectoristributed over a sphere and the random unit vector \(\mathbf \sigma\) with \(\mathbf{d}^n \cdot \mathbf{\sigma}^n = 0\) can be calculated.

Computation

During simulation the following calculations have to be made:

To calculate the normalized wave amlitudes:

\[ q^n = \frac{E(k^n) \Delta k^n}{\sum_{n=1}^N E(k^n) \Delta k^n} \]

the prescribed spatial spectrum of the kinetic energy of turbulence is defined by a modified von Kármán [60] spectrum:

\[ E(k) = \frac{(k / k_e)^4}{\left[ 1 + 2.4 (k / k_e)^2 \right]^{17/6}} f_\eta f_\text{cut} \]

with:

\[ f_\eta = \exp \left[ - \left( \frac{12 k}{k_\eta} \right)^2 \right] \]

\[ f_\text{cut} = \exp \left[ - \left( \frac{4 \max \left( k - 0.9 k_\text{cut}, 0 \right)}{k_\text{cut}} \right)^3 \right] \]

and \(k_\eta = \frac{2 \pi}{l_\eta}\) with \(l_\eta = \left( \frac{\nu^3}{\epsilon} \right)^{\frac 1 4}\) with viscosity \(\mu=\nu\rho\) calculated according to Sutherland's law [56] (see Sutherland's law for calorically perfect gas ). The kinetic energy of turbulence \( k\) is calculated as half the trace of the Reynolds stress tensor. According to Wilcox [64] \(\epsilon=\beta k\omega\) (eq.A4b) and \(L_T=\sqrt{k}/\omega\) (eq.24):

\begin{align} \epsilon&=\beta k\omega\\ &=\beta k \frac{\sqrt{k}}{L_T}\\ &=\beta \frac{\sqrt{k^2}\sqrt{k}}{L_T}\\ &=\beta \frac{\sqrt{k^3}}{L_T}\\ \frac{1}{\epsilon}&= \frac{L_T}{\beta k^\frac{3}{2}} \end{align}

Thus, for non-dimensionalization can be written:

\begin{align} l_{\eta}&= \left( \frac{\nu^3}{\epsilon} \right)^{\frac 1 4} = \left( \frac{\mu^3L_T}{\rho^3\beta k^{3/2}} \right)^{\frac 1 4}\\ \Rightarrow l_0l^*_{\eta}&= \left( \frac{\mu_0^3\left(\mu^*\right)^3l_0L_T^*}{\left(\rho_0 \rho^*\right)^3\beta\left(a_0^2k^*\right)^{3/2}}\right)^{\frac 1 4}\\ l^*_{\eta}&= \left(\left(\frac{\mu^*}{\rho^*}\right)^3 \frac{L_T^*\mu_0^3l_0}{\beta\left(k^*\right)^{3/2}\rho^3_0a_0^3} \right)^{\frac 1 4}\frac{1}{l_0} &= \left( \left(\frac{\mu^*}{\rho^*}\right)^3\frac{L_T^*\nu_0^3l_0}{\beta\left(k^*\right)^{3/2}a_0^3}\frac{1}{l_0^4} \right)^{\frac 1 4}\\ &= \left( \left(\frac{\mu^*}{\rho^*}\right)^3\frac{L_T^*}{\beta\left(k^*\right)^{\frac 3 2}Re_0^3} \right)^{\frac{1}{4}} \end{align}

The ''pseudo-position vector'' [49] is defined as:

\[\mathbf{r}'=\left(\begin{matrix}\frac{2\pi}{k^nl_{e,\max}}\left(x-U_0t\right)\\y\\z\end{matrix}\right)\]

If the command STG_REMOVE_PSEUDOPOSITION_FACTOR is used, the factor \(\frac{2\pi}{k^nl_{e,\max}}\) is replaced by \(1\).

The auxiliary fluctuation vector:

\[ {\mathbf v}^{\prime}\left(\mathbf x, t\right) = 2\sqrt{3/2} \sum\limits_{n=1}^{N} \sqrt{q^n}\left[{\mathbf \sigma}^n \cos\left(k^n{\mathbf d}^n\cdot {\mathbf r}^{\prime} + \phi^n\right)\right]\]

and the Cholesky decomposition \(\hat{\mathbf A}\) of the Reynolds stress tensor \( \mathbf R = \hat{\mathbf A}^T\hat{\mathbf A}\) with \(a_{ij}=\left(\hat{\mathbf A}\right)a_{ij}\):

\[ \hat{\mathbf A} = \left[ \begin{matrix} \sqrt{r_{11}} &0&0\\ r_{21}/a_{11} &\sqrt{r_{22}-a_{21}^2}&0\\ r_{31}/a_{11} &\left(r_{32}-a_{21}a_{31}\right)/a_{22}&\sqrt{r_{33}-a_{31}^2-a_{32}^2} \end{matrix}\right]\]

are multiplied to get the velocity fluctuations:

\[ \mathbf{u}^{\prime}\left(\mathbf x, t\right) = \hat{\mathbf{A}} \mathbf{v}^{\prime}\left(\mathbf x, t\right) \]

which are then added to the mean velocity:

\[ {\mathbf u} = \overline{{\mathbf u}} + {\mathbf u}^{\prime}\left(\mathbf x, t\right) \]