1.4. Turbulence Inflow¶

Introduction¶

In computational wind engineering (CWE), the generation of inflow turbulence satisfying prescribed mean-velocity profiles, turbulence spectra, and spatial and temporal correlations is of great importance for the accurate evaluation of wind effects on buildings and structures. More specifically, the task is to generate a turbulent velocity field \(\boldsymbol{u}(\boldsymbol{x},t)\) with the form

(1.4.1)¶\[\boldsymbol{u}(\boldsymbol{x},t) = \boldsymbol{U}(\boldsymbol{x},t)+\boldsymbol{u}'(\boldsymbol{x},t)\]

where \(\boldsymbol{U}(\boldsymbol{x},t)\) and \(\boldsymbol{u}'(\boldsymbol{x},t)\) are the mean and fluctuating velocities at the position \(\boldsymbol{x}\). The turbulent velocity field \(\boldsymbol{u}\) and its fluctuation \(\boldsymbol{u}'\) need to satisfy a number of properties which are list below:

\(\boldsymbol{u}'\) should be spatially and temporally correlated.

item \(\boldsymbol{u}'\) needs to have prescribed Reynolds stresses tensor \(R_{ij}(\boldsymbol{x}) = \overline{u_i'u_j'}(\boldsymbol{x})\) where \(u_i'\) (\(i=1,2,3\)) is the \(i\)-th component of \(\boldsymbol{u}'\) and the over line denotes the time average.
\(\boldsymbol{u}'\) needs to have prescribed integral length scales \(L_{ij}(\boldsymbol{x},\boldsymbol{e})\)

(1.4.2)¶\[L_{ij}(\boldsymbol{x},\boldsymbol{e}) = \int_{0}^{\infty} \rho_{ij}(\boldsymbol{x},r\boldsymbol{e})\ \mathrm{d}r,\]

where \(\rho_{ij}(\boldsymbol{x},\boldsymbol{e})\) is the correlation function given by

(1.4.3)¶\[\rho_{ij}(\boldsymbol{x},\boldsymbol{e}) = \frac{\overline{u_i'(\boldsymbol{x},t)u_j'(\boldsymbol{x}+\boldsymbol{e},t)}}{\overline{u_i'(\boldsymbol{x},t)u_j'(\boldsymbol{x},t)}}.\]
\(\boldsymbol{u}\) should fulfil the divergence free constraint \(\nabla \cdot \boldsymbol{u} = 0\).
\(\boldsymbol{u}'\) should have prescribed correlation functions \(\rho_{ij}(\boldsymbol{x},\boldsymbol{e})\) or spectra.

Several methodologies have been proposed for this purpose which can be classified into three general categories: precursor simulation methods, recycling methods and synthetic methods. Compared with precursor simulation and recycling methods, synthetic methods in general offer a more practical and relatively efficient approach to generate inflow turbulence, and is therefore chosen as the subject of this section. Research activities on synthetic turbulence generation have been vigorous over the past decades and have branched out into several categories of techniques [Wu17], including the synthetic random Fourier method [Kra70, Hos72, MB21], the synthetic digital filtering method [KSJ03] and the synthetic eddy methods [JBLP06]. A brief introduction regarding to these techniques is given below and emphasis is placed on their abilities to capture the statistical characteristics as well as the spatial and temporal coherence of turbulence. Also note that since real turbulence is very complex, in most cases, not all of the above listed features can be fulfilled. There is always some adaptation time required for artificial turbulence to evolve into real turbulence. Fulfilling the properties above with the synthetic turbulence is important to minimize the adaptation time or length.

Synthetic Random Fourier Method¶

The so-called synthetic random Fourier method (SRFM) attempts to model turbulent flow field indirectly by imposing constraints on uncorrelated random fields through an energy spectrum to account for the spatial and temporal correlations, which can be further classified into two groups. The first group of the SRFM was based on the pioneering work in [Hos72] and [SJ72] on the simulation of multi-correlated random processes using a weighted amplitude wave superposition (WAWS) method. This approach has an advantage that both the targeted power- and cross-spectra can be imposed in the generation process so that the prescribed target characteristics can be maintained. A major drawback of this method is that the generated turbulence does not satisfy the continuity equation of the flow, or in other words, the divergence-free condition is not guaranteed. As a consequence, it would take enormous effort for the solver to enforce the continuity by correcting the turbulence inflow inserted into the computational domain, and the statistical characteristics of the corrected flow field differ from the target values. More recently, [MB21] proposed a computationally efficient implementation of the method in [SJ72] to generate inflow turbulence by enforcing continuity requirements posteriorly.

The second group of the SRFM was initiated by the work in [Kra70] on divergence-free homogeneous isotropic turbulence synthesis through the superposition of random harmonic functions. [SSC01] took a step forward by combining Kraichnan’s technique with scaling and orthogonal transformation operations in a procedure known as the random flow generation (RFG) which allows to generate inhomogeneous and anisotropic turbulence. However, the scaling operation introduced in the RFG technique can result in a velocity field that is not divergence-free for inhomogeneous turbulence. Modifications to enforce the divergence-free the constraint for inhomogeneous turbulence was discussed in [YB14]. A major drawback of RFG technique is that the power spectra of the generated turbulence only follow Gaussian spectra model, so it is not suitable for simulating flows in the atmospheric boundary layer. [HLW10] revisited Kraichnan’s method and proposed a technique called DSRFG (for discretizing and synthesizing random flow generation) which allows to generate turbulent inflow from any prescribed spectrum. Instead of using the scaling and orthogonal transformation, the anisotropy of turbulence is realized by modifying the distribution strategy of the wave vector in Kraichnan’s original method. A drawback of the DSRFG technique is that it produces fluctuating velocities with high correlation due to the fact that in this method the spatial correlation is modeled by a parameter which is not a function of frequency but a constant value. Inspired by the DSRFG method, [CPMS17] proposed some modifications to this technique to obtain the velocity field that had a better match with the target turbulent statistics. This method, known as modified discretizing and synthesizing random flow generation (MDSRFG), is capable of removing the dependence of statistical quantities of synthetic turbulence on spectra discretization resolution. [AEBED15] also proposed a technique called consistent discrete RFG (CDRFG) to accurately model the target spectra and the coherence function. In both two methods mentioned above, the parameter that characterizes the spatial correlation is expressed as a function of frequency to account for the damping of coherence with the increase of frequency. An attractive feature of second group of SRFM is that the generation procedures are usually independent at each point and each time-instant so that it can be easily accelerated by conducting parallel computation, although the generated random flow may not satisfy the continuity equation.

Synthetic Eddy Method¶

The synthetic eddy method (SEM) initiated by [JBLP06] is based on the classical view of turbulence as a superposition of the representative coherent eddies. In the SEM, the flow is assumed to consist of randomly distributed turbulent spots, and each turbulent spot is modeled by a three-dimensional shape function with compact support and satisfies a proper normalization condition. The spots are then assumed to be convected through an inlet plane with a reference velocity using Taylor’s frozen turbulence hypothesis. The resulting inflow turbulence is then reconstructed using the method proposed to recover the desired statistical characteristics and to account for the conditions of inhomogeneity and anisotropy. The choice of the shape function plays an important role in the SEM since it is directly related to the two-point autocorrelation function, and consequently the power spectrum of the synthetic turbulence. Enforcement of the continuity condition in the SEM was discussed in [PCR13].

A brief introduction on the SEM presented by [JBLP06] is given as follows. To start with, the turbulent spot mentioned above can be represented as eddies defined by shape function \(f\) which has a compact support on \([-1,1]\) and has the normalization

(1.4.4)¶\[\int_{-1}^1 f^2(x) \mathrm{d}x = 1\]

The inflow plane on which we want to generate the synthetic turbulence with the SEM is basically a finite set of points \(S =\{\boldsymbol{x}_1,\boldsymbol{x}_2,\ldots,\boldsymbol{x}_s\}\). The first step is to create a box of eddies \(B\) surrounding \(S\) which is going to contain the synthetic eddies. It is defined by

(1.4.5)¶\[B = \big\{(x_1,x_2,x_3)\in \mathbb{R}^3: x_{i,\text{min}}<x_i<x_{i,\text{max}}\big\}\]

where

(1.4.6)¶\[x_{i,\text{min}} = \text{min}(x_i-\sigma_i(\boldsymbol{x})), \quad x_{i,\text{max}} = \text{max}(x_i+\sigma_i(\boldsymbol{x})), \quad \boldsymbol{x}\in S\]

The volume of the box of eddies is noted by \(V_B\). In the synthetic eddy method, the velocity signal generated by \(N\) eddies has the representation

(1.4.7)¶\[u_i(\boldsymbol{x}) = U_i(\boldsymbol{x}) + \frac{1}{\sqrt{N}}\sum_{k=1}^N a_{ij} \epsilon_j^k f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k)\]

where \(\boldsymbol{x}\) represent the coordinates of computational points and \(\boldsymbol{x}^k\) represent the coordinates of eddies. The coefficient \(a_{ij}\) results from the Cholesky decomposition of a prescribed Reynolds stress tensor \(R_{ij}\)

(1.4.8)¶\[\begin{split}\left(\begin{matrix} \sqrt{R_{11}} & 0 & 0 \\ R_{21}/a_{11} & \sqrt{R_{22}-a_{21}^2} & 0 \\ R_{31}/a_{11} & (R_{32}-a_{21}a_{31})/a_{22} & \sqrt{R_{33}-a_{31}^2-a_{32}^2} \end{matrix}\right)\end{split}\]

The coefficient \(\epsilon_j^k\) (\(j=1,2,3\)) is is the uniformly random intensity factor of values \(+1\) or \(-1\), and \(f_{\boldsymbol{\sigma}(\boldsymbol{x})} (\boldsymbol{x}-\boldsymbol{x}^k)\) is the velocity distribution at \(\boldsymbol{x}\) of the eddy located at \(\boldsymbol{x}^k\) defined as follows:

(1.4.9)¶\[f_{\boldsymbol{\sigma}(\boldsymbol{x})} (\boldsymbol{x}-\boldsymbol{x}^k) = \sqrt{\frac{V_B}{\sigma_1\sigma_2\sigma_3}}f\left(\frac{x_1-x_1^k}{\sigma_1}\right)f\left(\frac{x_2-x_2^k}{\sigma_2}\right)f\left(\frac{x_3-x_3^k}{\sigma_3}\right)\]

where \(\boldsymbol{\sigma}=(\sigma_1,\sigma_2,\sigma_3)^T\). The position of the eddies \(\boldsymbol{x}^k\) before the first time step are independent from each other and taken from a uniform distribution over the box of eddies \(B\). The eddies are convected through the box of eddies \(B\) with the mean velocity \(\boldsymbol{U}(\boldsymbol{x})\). At each time step, the new position of eddy \(k\) is given by

(1.4.10)¶\[\boldsymbol{x}^k(t+\varDelta t) = \boldsymbol{x}^k(t)+\boldsymbol{U}(\boldsymbol{x}^k)\varDelta t\]

where \(\varDelta t\) is the time step of the simulation. If an eddy \(k\) is convected out of the box \(B\), then it is immediately regenerated randomly with in the region

(1.4.11)¶\[B_{\varDelta t} = \left\{ \boldsymbol{x}\notin B, \ \boldsymbol{x}+\boldsymbol{U}(\boldsymbol{x})\varDelta t \in B \right\}\]

with a new random intensity vector \(\epsilon_j^k\). \(B_{\varDelta t}\) denotes the region in which regenerated eddy \(\boldsymbol{x}^k(t) \in B_{\varDelta t}\) dose not effect the synthetic velocity at the inflow plane until the next time-step.

Mean flow and Reynolds stresses¶

The mean value of the velocity signal (1.4.7) can be expressed as

(1.4.12)¶\[\left\langle u_i \right\rangle = U_i(\boldsymbol{x}) + \frac{1}{\sqrt{N}}\sum_{k=1}^N \left\langle a_{ij} \varepsilon_j^k f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k) \right\rangle\]

where the angles denote the mean operator. The independence between the random variables \(\boldsymbol{x}^k\) and \(\varepsilon_j^k\) in the mean operator implies that

(1.4.13)¶\[\left\langle a_{ij} \varepsilon_j^k f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k) \right\rangle = a_{ij} \left\langle\varepsilon_j^k\right\rangle \left\langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k) \right\rangle\]

The term \(\langle\varepsilon_j^k\rangle = 0\) since the intensities of the eddies is either \(1\) or \(-1\) with equal probability. Consequently, we obtain

(1.4.14)¶\[\left\langle u_i \right\rangle = U_i(\boldsymbol{x}).\]

The Reynolds stresses \(\langle u_i u_j \rangle\) of the synthesized write

(1.4.15)¶\[\langle u_i u_j \rangle = \frac{1}{N}\sum_{k=1}^N\sum_{k=1}^N a_{im}a_{jn} \langle \varepsilon_m^k \varepsilon_n^l \rangle \langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k) f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^l) \rangle\]

Using again the independence between the random variables \(\boldsymbol{x}^k\) and \(\varepsilon_j^k\), the above equation reduces to

(1.4.16)¶\[\langle u_i u_j \rangle = \frac{1}{N}\sum_{k=1}^N a_{im}a_{jm} \langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}^2(\boldsymbol{x}-\boldsymbol{x}^k)\]

The term

(1.4.17)¶\[\langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}^2(\boldsymbol{x}-\boldsymbol{x}^k) \rangle = \int_{\mathbb{R}^3} p(\boldsymbol{y}) f_{\boldsymbol{\sigma}(\boldsymbol{x})}^2(\boldsymbol{x}-\boldsymbol{x}^k) = 1\]

follows from the fact that \(\boldsymbol{x}^k\) follows a uniform distribution over \(B\), i.e.

(1.4.18)¶\[\begin{split}p(\boldsymbol{y}) = \begin{cases} \frac{1}{V_B} & \boldsymbol{y} \in B \\ 0 & \boldsymbol{y} \notin B \end{cases}.\end{split}\]

Finally, we arrive at

(1.4.19)¶\[\langle u_i u_j \rangle = \frac{1}{N}\sum_{k=1}^N a_{im}a_{jm} = R_{ij}\]

Hence the Reynolds stresses of the velocity fluctuations generated by the SEM reproduce exactly the input Reynolds stresses.

Two-point correlation¶

The two-point cross-correlation of the velocity fluctuations writes

(1.4.20)¶\[R_{ij}(\boldsymbol{x},\boldsymbol{r}) = \langle u_i(\boldsymbol{x},t) u_j(\boldsymbol{x}+\boldsymbol{r},t) \rangle\]

where \(\boldsymbol{r} = (r_1,r_2,r_3)\) is a vector defining the relative positions between the two points at which the velocity correlations are computed. By (1.4.7) and the linearity of the statistical mean, we obtain

(1.4.21)¶\[R_{ij}(\boldsymbol{x},\boldsymbol{r}) = \frac{1}{N}\sum_{k=1}^N\sum_{k=1}^N a_{im}a_{jn} \langle \varepsilon_m^k \varepsilon_n^l \rangle \langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k) f_{\boldsymbol{\sigma}(\boldsymbol{x}+\boldsymbol{r})}(\boldsymbol{x}+\boldsymbol{r}-\boldsymbol{x}^l) \rangle\]

Using again the independence between the positions \(\boldsymbol{x}^k\) and the intensities \(\varepsilon^k\) of the eddies, this yields

(1.4.22)¶\[R_{ij}(\boldsymbol{x},\boldsymbol{r}) = \frac{1}{N}\sum_{k=1}^N a_{im}a_{jm} \langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k) f_{\boldsymbol{\sigma}(\boldsymbol{x}+\boldsymbol{r})}(\boldsymbol{x}+\boldsymbol{r}-\boldsymbol{x}^k) \rangle\]

By (1.4.18), the term in the mean operator writes

(1.4.23)¶\[\langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k) f_{\boldsymbol{\sigma}(\boldsymbol{x}+\boldsymbol{r})}(\boldsymbol{x}+\boldsymbol{r}-\boldsymbol{x}^k) \rangle = \frac{1}{V_B} \int_B f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{y}) f_{\boldsymbol{\sigma}(\boldsymbol{x}+\boldsymbol{r})}(\boldsymbol{x}+\boldsymbol{r}-\boldsymbol{y}) \mathrm{d}\boldsymbol{y}\]

Inserting (1.4.23) back to (1.4.22) and using (1.4.9), this yields

(1.4.24)¶\[R_{ij}(\boldsymbol{x},\boldsymbol{r}) = R_{ij} \cdot \prod_{l=1}^3 \left[f_{\boldsymbol{\sigma}(\boldsymbol{x})} *f_{\boldsymbol{\sigma}(\boldsymbol{x}+\boldsymbol{r})} \right](r_l)\]

where \(∗\) denotes the convolution product. For homogeneous turbulence where integral length scales \(\boldsymbol{\sigma}(\boldsymbol{x}) = \boldsymbol{\sigma}(\boldsymbol{x}+\boldsymbol{r}) =(\sigma,\sigma,\sigma)^T\), the two-point cross-correlation tensor \(R_{ij}(\boldsymbol{x},\boldsymbol{r})\) only depends on \(\boldsymbol{r}\) and consequently (1.4.24) simplifies to

(1.4.25)¶\[R_{ij}(\boldsymbol{r}) = R_{ij} \cdot \prod_{l=1}^3 \left[f*f\right]\left(\frac{r_l}{\sigma}\right)\]

Recall the integral length scale \(L_{ij}\) is defined as the integral of the two-point correlation \(R_{ij}(\boldsymbol{x},\boldsymbol{r})\) in a particular direction and is thus proportional to \(\sigma\). By integrating (1.4.25), one easily verifies that (for homogeneous turbulence) \(L_{ij}=C_f\sigma\) in every direction where \(C_f\) only depends on the choice of \(f\).

Fourier analysis can also be used to obtain the spectra of the synthetic turbulence. Note that the velocity spectrum tensor \(\phi_{ij}(k)\) is the Fourier transform of the two-point correlation tensor

(1.4.26)¶\[\phi_{ij}(\boldsymbol{k}) = \mathcal{F}_{\boldsymbol{k}}\left\{R_{ij}(\boldsymbol{r})\right\}\]

Recall the convolution theorem for cross-correlation states that

(1.4.27)¶\[\mathcal{F}_{\boldsymbol{k}}\left\{f * f\right\} = |\mathcal{F}_{\boldsymbol{k}}\left\{f\right\}|^2\]

Hence the spatial velocity spectrum tensor can be expressed as

(1.4.28)¶\[\phi_{ij}(\boldsymbol{k}) = R_{ij}\sigma^3 \cdot \prod_{l=1}^3|\mathcal{F}_{k_l\sigma}\left\{f\right\}|^2\]

where \(\boldsymbol{k} = (k_1,k_2,k_3)\). More specifically for instance, the one-dimensional spectra in the \(x\) direction is

(1.4.29)¶\[E_{ij}(k) = R_{ij}\sigma^3 \cdot |\mathcal{F}_{k_l\sigma}\left\{f\right\}|^2\]

Two-time correlation¶

The two-time correlation tensor of the velocity, denoted by \(R_{ij}(\boldsymbol{x},\tau)\), is the correlation between \(u_i(\boldsymbol{x},t)\) and \(u_j(\boldsymbol{x},t+\tau)\) at times \(t\) and \(t + \tau\) respectively, i.e.,

(1.4.30)¶\[R_{ij}(\boldsymbol{x},\tau) = \langle u_i(\boldsymbol{x},t) u_j(\boldsymbol{x},t+\tau) \rangle.\]

By (1.4.7) and the linearity of the statistical mean, we have

(1.4.31)¶\[R_{ij}(\boldsymbol{x},\tau) = \frac{1}{N}\sum_{k=1}^N\sum_{k=1}^N a_{im}a_{jn} \langle \varepsilon_m^k(t) \varepsilon_n^l(t+\tau) f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k(t)) f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^l(t+\tau)) \rangle\]

The independence between the position \(\boldsymbol{x}^k\) and intensity \(\varepsilon_m^k\) of different eddies implies that, for \(k \neq l\), the statistical mean in (1.4.31) can be split as follows

(1.4.32)¶\[\langle \varepsilon_m^k(t) \rangle \langle \varepsilon_n^l(t+\tau) \rangle \langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k(t)) \rangle \langle f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^l(t+\tau)) \rangle = 0\]

Consequently (1.4.31) reduces to

(1.4.33)¶\[R_{ij}(\boldsymbol{x},\tau) = \frac{1}{N}\sum_{k=1}^N a_{im}a_{jn} \langle \varepsilon_m^k(t) \varepsilon_n^k(t+\tau) f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k(t)) f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{x}^k(t+\tau)) \rangle\]

Before computing the term in the angles, we define \(B_{\tau} \in B\) such that all eddies that present in \(B_{\tau}\) at time \(t\) will be convected far enough so that they will be recycled at least once before time \(t+\tau\)

(1.4.34)¶\[B_{\tau} = \left\{\boldsymbol{x}\in B, \ \boldsymbol{x}+\tau \boldsymbol{U}(\boldsymbol{x}) \in B \right\}\]

If \(\boldsymbol{x}^k(t)\in B_{\tau}\), then it is going to be recycled between time \(t\) and \(t+\tau\) and hence both \(\boldsymbol{x}^k(t+\tau)\) and \(\varepsilon_m^k(t+\tau)\) will be independent of their previous values. The contribution of an eddy \(k\) located within the region where \(\boldsymbol{x}^k(t) \in B_{\tau}\) to the term in the angles of (1.4.33) is thus zero. On the contrary if \(\boldsymbol{x}^k(t) \in B_{\tau}\), the eddy \(k\) will remain inside of the box \(B\) at time \(t + \tau\) and hence \(\varepsilon_m^k(t+\tau) = \varepsilon_m^k(t)\) and \(\boldsymbol{x}^k(t+\tau) =\boldsymbol{x}^k(t)+\tau\boldsymbol{U}(\boldsymbol{x}^k)\). Thus both \(\varepsilon_n^k(t+\tau) = \varepsilon_n^k(t)\) and \(\boldsymbol{x}^k(t+\tau)\) depend on the previous position \(\boldsymbol{x}^k(t)\) of eddy \(k\) relative to \(B_{\tau}\). By (1.4.19) and the definition of \(B_{\tau}\), (1.4.20) can then be replaced by

(1.4.35)¶\[R_{ij}(\boldsymbol{x},\tau) = R_{ij} \int_{B/B_{\tau}}f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{y}) f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-(\boldsymbol{y}+\tau\boldsymbol{U}_c)) \ \mathrm{d}\boldsymbol{y}\]

Since \(\boldsymbol{y}\in B_{\tau}\) leads to \(f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-(\boldsymbol{y}+\tau\boldsymbol{U}))=0\), the integral over \(B/B_{\tau}\) in the above expression can be extended to an integral over \(B\). Besides \(\boldsymbol{y}\in B\) suggests \(f_{\boldsymbol{\sigma}(\boldsymbol{x})}(\boldsymbol{x}-\boldsymbol{y})=0\) as previously demonstrated, therefore the integral in (1.4.35) can be further extended to an integral over \(\mathbb{R}^3\). Using (1.4.9), we finally arrive at

(1.4.36)¶\[R_{ij}(\boldsymbol{x},\tau) = R_{ij} \cdot \prod_{l=1}^3[f*f]\left(\frac{\tau U_{l}(\boldsymbol{x})}{\sigma_l(\boldsymbol{x})}\right)\]

In the case where the mean velocity is in the x-direction only \(\boldsymbol{U} = (U,0,0)\) and the target turbulence is homogeneous, (1.4.36) simplifies to

(1.4.37)¶\[R_{ij}(\boldsymbol{x},\tau) = R_{ij} [f*f]\left(\frac{\tau U(\boldsymbol{x})}{\sigma(\boldsymbol{x})}\right)\]

Thus the two-time correlation of the signal at time \(\tau\) is simply the autocorrelation function of \(f\) at separation distance \(\tau U /\sigma\). By integrating the above equation it can be proved that the integral time scale of the signal writes \(T = \sigma/U C_f\) where \(C_f\) is a coefficient only depends on the choice of \(f\). Since the synthetic velocity is a stationary process, the information the two-time cross-correlation tensor \(R_{ij}(\boldsymbol{x},\tau)\) contains can be re-expressed in terms of the wave number velocity spectrum tensor which writes

(1.4.38)¶\[\phi_{ij}(\boldsymbol{x},\omega) = \mathcal{F}_{\omega}\{R_{ij}(\boldsymbol{x},\tau)\}\]

Using again the convolution theorem as expressed, the above expression simplifies to

(1.4.39)¶\[\phi_{ij}(\boldsymbol{x},\omega) = R_{ij}\frac{\sigma}{|U|} |\mathcal{F}_{\omega\sigma / |U|}\{f\}|^2\]

Commonly used velocity shape functions¶

We list three commonly used velocity shape functions \(f\) below for reference. There are the tent function, the step function and the truncated Gaussian function.

Tent function

(1.4.40)¶\[\begin{split}f(x) = \begin{cases} \sqrt{\frac{3}{2}}(1-|x|), & 0 \leq |x| < 1 \\ 0, & |x| \geq 1 \end{cases}\end{split}\]

(1.4.41)¶\[\begin{split}[f*f](r) = \begin{cases} 1-\frac{3}{2}r^2+\frac{3}{4}|r|^3, & 0 \leq |r| < 1 \\ 2-3|r|+\frac{3}{2}r^2-\frac{1}{4}|r|^3, & 1 \leq |r| <2 \\ 0, & |r|\geq 2 \end{cases}\end{split}\]
Step function

(1.4.42)¶\[\begin{split}f(x) = \begin{cases} \frac{1}{\sqrt{2}}, & 0 \leq |x| < 1 \\ 0, & |x| \geq 1 \end{cases}\end{split}\]

(1.4.43)¶\[\begin{split}[f*f](r) = \begin{cases} 1-\frac{|r|}{2}, & 0 \leq |r| < 2 \\ 0, & |r|\geq 2 \end{cases}\end{split}\]
Truncated Gaussian function

(1.4.44)¶\[\begin{split}f(x) = \begin{cases} Ce^{-9x^2/2}, & 0 \leq |x| < 1 \\ 0, & |x| \geq 1 \end{cases}\end{split}\]

(1.4.45)¶\[\begin{split}[f*f](r) = \begin{cases} e^{-9r^2/2} & \leq |r| < 2 \\ 0, & |r|\geq 2 \end{cases}\end{split}\]

where \(C\) is a constant that ensures \(f\) satisfies the normalization (1.4.4).

Digital filtering method¶

The synthetic digital filtering method (SDFM) initiated by [KSJ03] attempts to model the spatial and temporal coherence of turbulent inflow through the digital filtering uncorrelated random data, and account for inhomogeneity and anisotropy using the method proposed by [LWS98]. It is relatively easy to implement and is able to reproduce the first and second order one-point statistics as well as autocorrelation function. However, the synthetic turbulence generated by SDFM does not satisfy the continuity equation. [KCX13] offered a promising approach to enforce the divergence-free constraint in the SDFM by inserting the synthetic turbulence on a transverse plane near the inlet and relying on pressure-velocity coupling to do the correction. From a computational wind engineering point of view, the ability of SDFM to impose a two-point spatial correlation directly is very attractive.

We now briefly introduce the filtering method by [KSJ03]. In order to create two-point correlations, let \(r_m\) be a series of random data with zero mean and unity variance, then

(1.4.46)¶\[u_m = \sum_{n=-N}^N b_n r_{m+n}\]

defines a convolution or a digital linear non-recursive filter. The \(b_n\) are filter coefficients and \(N\) is related to the length of the filter. The independence between two different random numbers \(r_m\) and \(r_n\) implies that \(\langle r_m r_n \rangle = 0\) for \(m \neq n\) and consequently the two-point correlation between \(u_{m}\) and \(u_{m+k}\) writes

(1.4.47)¶\[R_{uu}(k\varDelta x) = \frac{\langle u_{m} u_{m+k} \rangle}{\langle u_{m} u_{m} \rangle} = \sum_{j=-N+k}^N b_j b_{j-k} / \sum_{-j=-N}^N b_j^2\]

where \(\varDelta x\) is the grid spacing. Note that \(u_{m}\) and \(u_{m+k}\) can be interpolated as the values of a random variable field (e.g., velocity) at two distinct grid points with a distance \(k\varDelta x\) defined on a one dimensional axis. It is straightforward to tell (1.4.47) defines a relation between the filter coefficients and the correlation function of \(u_m\) (denoted by \(R_{uu}\) hereafter). This suggests that a prescribed correlation function can be reproduced through a careful determination of the filter coefficients. Also note that the coefficients should be determined such that the resulting correlation function fulfil some basic properties like \(R_{uu}(0)=1\), \(R_{uu}(\infty) = 0\) and the prescribed integral length scales.

For a general target correlation function, the filter coefficients \(b_n\) can be computed by solving a system of non-linear equations in the form of (1.4.47) with a multidimensional Newton method. The procedure can be taken from a standard textbook and needs no further comment. However, for a Gaussian or an exponential type of correlation function, there exists a simple but approximate prescribed solution. More specifically, for a Gaussian correlation function in the form of

(1.4.48)¶\[R(r) = \mathrm{exp}\left(-\frac{\pi r^2}{4L^2}\right)\]

where \(r\) is the distance and \(L\) is the length scale. It is possible to approximately reproduce (1.4.48) by computing the filter coefficients as

(1.4.49)¶\[b_k = \tilde{b}_k / \left( \sum_{j=-N}^N \tilde{b}_j^2 \right)^{1/2}\]

where

(1.4.50)¶\[\tilde{b}_k = e^{-\frac{\pi k^2}{4n^2}}\]

The width \(N\) of the filter should be chosen such that \(N\geq 2n\) (where \(n=L\varDelta x\)) to ensure the accuracy of the approximation. On the other hand, for an exponential correlation function

(1.4.51)¶\[R(r) = \mathrm{exp}\left(-\frac{\pi |r|}{2L}\right)\]

It is suggested [XC08] to evaluate the filter coefficients using

(1.4.52)¶\[b_k = \tilde{b}_k / \left( \sum_{j=-N}^N \tilde{b}_j^2 \right)^{1/2}\]

where

(1.4.53)¶\[\tilde{b}_k = e^{-\frac{\pi|k|}{n}}\]

Again, the width \(N\) of the filter should be chosen such that \(N\geq 2n\) (where \(n=L\varDelta x\)) to ensure the accuracy of the approximation. Now we have finished the discussion of the digital filtering method for one-dimensional case. Such a technique of generating spatially (or temporally) correlated data from general random numbers can be easily extended to three dimensional case by introducing multi-index filter coefficients \(b_{ijk}\) defined as

(1.4.54)¶\[b(i,j,k) = b_{ijk} = b_i \cdot b_j \cdot b_k\]

An algorithm for generating inflow data may look like this (alternatively one can generate a large volume of data, store it and convect it through the inflow plane by applying Taylor’s hypothesis):

Choose for each coordinate direction corresponding to the inflow plane a length scale \(L_y = n_y\varDelta y\), \(L_z = n_z\varDelta z\), a time scale \(T\) and determine the filter width \(N_{\alpha}\) (\(\alpha =x,y,z\)) accordingly.
Initialize and store three random fields \(R_{\alpha}\) (again \(\alpha =x,y,z\)) of dimensions \([-N_x:N_x,-N_y+1:M_y+N_y,-N_z+1:M_z+N_z]\) where \(M_y \times M_z\) denotes the dimensions of computational gird of the inflow plane.

(c) Compute the filter coefficients \(b(i,j,k)\) with a prescribed function or by a multidimensional Newton method such that the resulting correlation function (1.4.47) meets the target one.

Applying the following filter operation for \(j=1,\ldots,M_y\), \(k=1,\ldots,M_z\)

(1.4.55)¶\[\Psi_{\alpha}(j,k) = \sum_{i'=-N_x}^{N_x}\sum_{j'=-N_y}^{N_y}\sum_{k'=-N_z}^{N_z}b(i',j',k')R_{\alpha}(i',j+j',k+k')\]

which yields the two-dimensional arrays of spatially correlated data \(\Psi_{\alpha}\), \(\alpha =x,y,z\).
Output velocity data with the transformation

(1.4.56)¶\[u_i(j,k) = U_i + a_{ij}\Psi_j(j,k)\]

where the coefficients \(a_{ij}\) are given by (1.4.8). This step ensures the synthetic velocity reproduces the target mean velocity and Reynolds stress tensor.
Discard the first \((y,z)\)-plane of \(\Psi_{\alpha}\) and shift the whole data: \(\Psi_{\alpha}(i,j,k) := R_{\alpha}(i+1,j,k)\). Fill the plane \(R_{\alpha}(N_x,j,k)\) with new random numbers.
Repeat the steps (d) \(\sim\) (g) for each time step.

If the target correlation function is an exponential function, an alternative approach by [XC08] can be adopted for generating inflow turbulence which turns out to be much more efficient than the method of [KSJ03]. Instead of using the filtering operation discussed above, Xie and Castro’s method obtain the temporal correlation with the expression

(1.4.57)¶\[\Psi_{\alpha}(t+\varDelta t,j,k) = \Psi_{\alpha}(t,j,k)\mathrm{exp}\left(-\frac{\pi \varDelta t}{2T} \right)+\varPsi_{\alpha}(t,j,k)\left[1-\mathrm{exp}\left(-\frac{\pi \varDelta t}{2T} \right)\right]^{0.5}\]

where \(\Psi_{\alpha}(t,j,k)\) and \(\varPsi_{\alpha}(t,j,k)\) are two set of spatially-correlated random data resulting from a two dimensional filtering operation. For simplicity, we write \(\Psi_{\alpha,0}\), \(\Psi_{\alpha,k}\), \(\varPsi_{\alpha,0}\) and \(\varPsi_{\alpha,k}\) for \(\Psi_{\alpha}(t,j,k)\), \(\Psi_{\alpha}(t+k\varDelta t,j,k)\), \(\varPsi_{\alpha}(t,j,k)\) and \(\varPsi_{\alpha}(t+k\varDelta t,j,k)\), respectively. One easily verifies that

(1.4.58)¶\[\begin{split}\begin{split} \left\langle \Psi_{\alpha,0}\Psi_{\alpha,k} \right\rangle &= \left\langle \Psi_{\alpha,0}\left\{\Psi_{\alpha,k-1}\left(-\frac{\pi \varDelta t}{2T} \right)+ \varPsi_{\alpha,k-1}\left[1-\mathrm{exp}\left(-\frac{\pi \varDelta t}{2T} \right)\right]^{0.5}\right\}\right\rangle \\ & = \left\langle \Psi_{\alpha,0} \Psi_{\alpha,k-1} \right\rangle \mathrm{exp}\left(-\frac{\pi \varDelta t}{2T}\right) \\ & \cdots \\ & = \mathrm{exp}\left(-\frac{k\pi \varDelta t}{2T}\right) \end{split}\end{split}\]

which reproduces an exponential function. An overall algorithm for generating the inflow velocity supported by the method of [XC08] can be stated as follows

Choose for each coordinate direction corresponding to the inflow plane a length scale \(L_y = n_y\varDelta y\), \(L_z = n_z\varDelta z\), a time scale \(T\) and determine the filter width \(N_{\alpha}(\alpha =x,y,z)\) accordingly.
Initialize and store three random fields \(R_{\alpha}\) (again \(\alpha =x,y,z\)) of dimensions \([-N_y+1:M_y+N_y,-N_z+1:M_z+N_z]\) where \(M_y \times M_z\) denotes the dimensions of computational gird in the inflow plane.
Compute the filter coefficients \(b(j,k)\) with a prescribed function or by a multidimensional Newton method such that the resulting correlation function meet the target one.
Applying the following filter operations for \(j=1,\ldots,M_y\), \(k=1,\ldots,M_z\)

(1.4.59)¶\[\varPsi_{\alpha}(j,k) = \sum_{j'=-N_y}^{N_y}\sum_{k'=-N_z}^{N_z}b(j',k')R_{\alpha}(j+j',k+k')\]

which yields the two-dimensional arrays of spatially correlated data \(\varPsi_{\alpha}\), \(\alpha =x,y,z\).
Compute \(\Psi_{\alpha}(j,k)\) with (1.4.57) and output the velocity signal with the transformation

(1.4.60)¶\[u_i(j,k) = U_i + a_{ij}\Psi_j(j,k)\]

where the coefficients \(a_{ij}\) are given by (1.4.8). Again, this step ensures the synthetic velocity reproduces the target mean velocity and Reynolds stress tensor.
Repeat the steps (d) \(\sim\) (f) for each time step.

References¶

1: Xingquan Guan, Henry Burton, and Thomas Sabol. Python-based computational platform to automate seismic design, nonlinear structural model construction and analysis of steel moment resisting frames. Engineering Structures, 224:111199, 2020.
2: Sejun Park and DongHun Yeo. Database-assisted design and equivalent static wind loads for mid- and high-rise structures: concepts, software, and user's manual. Technical Report NIST TN 2000, National Institute of Standards and Technology, Gaithersburg, MD, June 2018. doi:10.6028/NIST.TN.2000.
3: Xiaohua Wu. Inflow turbulence generation methods. Annual Review of Fluid Mechanics, 49:23–49, 2017.
4(1,2): Robert H Kraichnan. Diffusion by a random velocity field. The physics of fluids, 13(1):22–31, 1970.
5(1,2): Masaru Hoshiya. Simulation of multi-correlated random processes and application to structural vibration problems. In Proceedings of the Japan Society of Civil Engineers, number 204, 121–128. Japan Society of Civil Engineers, 1972.
6(1,2): Abiy F Melaku and Girma T Bitsuamlak. A divergence-free inflow turbulence generator using spectral representation method for large-eddy simulation of abl flows. Journal of Wind Engineering and Industrial Aerodynamics, 212:104580, 2021. URL: https://doi.org/10.1016/j.jweia.2021.104580.
7(1,2,3,4): M. Klein, A. Sadiki, and J. Janicka. A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations. Journal of Computational Physics, 186(2):652–665, 2003.
8(1,2,3): Nicolas Jarrin, Sofiane Benhamadouche, Dominique Laurence, and Robert Prosser. A synthetic-eddy-method for generating inflow conditions for large-eddy simulations. International Journal of Heat and Fluid Flow, 27(4):585–593, 2006.
9(1,2): Masanobu Shinozuka and C-M Jan. Digital simulation of random processes and its applications. Journal of sound and vibration, 25(1):111–128, 1972.
10: A Smirnov, S Shi, and I Celik. Random flow generation technique for large eddy simulations and particle-dynamics modeling. Journal of fluids engineering, 123(2):359–371, 2001.
11: Rixin Yu and Xue-Song Bai. A fully divergence-free method for generation of inhomogeneous and anisotropic turbulence with large spatial variation. Journal of Computational Physics, 256:234–253, 2014.
12: SH Huang, QS Li, and JR Wu. A general inflow turbulence generator for large eddy simulation. Journal of Wind Engineering and Industrial Aerodynamics, 98(10-11):600–617, 2010.
13: Hugo G Castro, Rodrigo R Paz, Javier L Mroginski, and Mario A Storti. Evaluation of the proper coherence representation in random flow generation based methods. Journal of Wind Engineering and Industrial Aerodynamics, 168:211–227, 2017.
14: Haitham Aboshosha, Ahmed Elshaer, Girma T Bitsuamlak, and Ashraf El Damatty. Consistent inflow turbulence generator for les evaluation of wind-induced responses for tall buildings. Journal of Wind Engineering and Industrial Aerodynamics, 142:198–216, 2015.
15: R Poletto, T Craft, and A Revell. A new divergence free synthetic eddy method for the reproduction of inlet flow conditions for les. Flow, turbulence and combustion, 91(3):519–539, 2013.
16: Thomas S Lund, Xiaohua Wu, and Kyle D Squires. Generation of turbulent inflow data for spatially-developing boundary layer simulations. Journal of computational physics, 140(2):233–258, 1998.
17: Yusik Kim, Ian P Castro, and Zheng-Tong Xie. Divergence-free turbulence inflow conditions for large-eddy simulations with incompressible flow solvers. Computers & Fluids, 84:56–68, 2013.
18(1,2,3): Zheng-Tong Xie and Ian P Castro. Efficient generation of inflow conditions for large eddy simulation of street-scale flows. Flow, turbulence and combustion, 81(3):449–470, 2008.
19: L. E. Wittig and A. K. Sinha. Simulation of multicorrelated random processes using the FFT algorithm. The Journal of the Acoustical Society of America, 58(3):630–634, 1975. URL: https://doi.org/10.1121/1.380702, doi:10.1121/1.380702.
20: J. C. Kaimal, J. C. Wyngaard, Y. Izumi, and O. R. Coté. Spectral characteristics of surface-layer turbulence. Quarterly Journal of the Royal Meteorological Society, 98(417):563–589, 1972. URL: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.49709841707, doi:10.1002/qj.49709841707.
21: Emil Simiu and Robert H Scanlan. Wind effects on structures: fundamentals and application to design. Book published by John Willey & Sons Inc, 1996.
22: AG Davenport. The dependence of wind loading on meteorological parameters. In Proc. of Int. Res. Seminar, Wind Effects On Buildings & Structures, NRC, Ottawa. 1967.