1. Theory
1.1. Introduction
In computational wind engineering (CWE) applications, 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. 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 generating inflow turbulence. Research activities on synthetic turbulence generation have been vigorous over the past decades. They have branched out into several categories of techniques, including the synthetic random Fourier method, the synthetic digital filtering method, and the synthetic eddy methods. Other schemes that promise to enhance the inflow simulation process involve the use of GPUs [5]. This is another area that can benefit from machine learning as one could tap on a database of simulated data matching prescribed inflow kinematics and dynamics.
In the following, an inflow toolbox TInF at the SimCenter offers an off-the-shelf simulation of inflow conditions using a host of stochastic simulation schemes and is presented with the theoretical background. The task of simulating inflow conditions involves the generation of a turbulent velocity field \(\mathbf{u}(\mathbf{x},t)\) given by
where \(\mathbf{U}(\mathbf{x},t)\) and \(\mathbf{u}'(\mathbf{x},t)\) are the mean and fluctuating velocities at the position \(\mathbf{x}\). The turbulent velocity field \(\mathbf{u}\) and its fluctuation \(\mathbf{u}'\) need to satisfy a number of properties which are list below:
\(\mathbf{u}'\) should be spatially and temporally correlated.
\(\mathbf{u}'\) needs to have prescribed Reynolds stresses tensor \(R_{ij}(\mathbf{x}) = \langle u_i'u_j'\rangle(\mathbf{x})\) where \(u_i'\) (\(i=1,2,3\)) is the \(i\)-th component of \(\mathbf{u}'\) and the angles denote the time average.
\(\mathbf{u}'\) needs to have prescribed integral length scales, \(L_{ij}^{\mathbf{e}}(\mathbf{x})\), with respect to a direction \(\mathbf{e}\) as
where \(\rho_{ij}(\mathbf{x},\mathbf{e})\) is the correlation function given by
\(\mathbf{u}\) should fulfil the divergence free constraint \(\nabla \cdot \mathbf{u} = 0\).
\(\mathbf{u}'\) should have prescribed correlation functions \(\rho_{ij}(\mathbf{x},\mathbf{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, the 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 [6], including the synthetic random Fourier method [7, 8], the synthetic digital filtering method [9] and the synthetic eddy methods [2]. 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 the artificial turbulence to evolve into real turbulence. Fulfilling the properties above with the synthetic turbulence is important to minimize the adaptation time or length.
1.2. 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 [8] and [10] 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 differs from the target values.
The second group of the SRFM was initiated by the work in [7] on divergence-free homogeneous isotropic turbulence synthesis through the superposition of random harmonic functions. [11] took a step forward by combing 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 constraint for inhomogeneous turbulence was discussed in [12]. A major drawback of RFG technique is that the power-spectra of the generated turbulence only follows Gaussian’s spectra model, so it is not suitable for simulating flows in atmospheric boundary layer. [13] revisited Kraichnan’s method and proposed a technique called DSRFG (for discretizing and synthesizing random flow generation) which allows to generate turbulence 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 modelled by a parameter which is not a function of frequency but a constant value. Inspired by the DSRFG method, [14] 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 statistic quantities of synthetic turbulence on spectra discretization resolution. [15] 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.
1.3. Digital Filtering Method
The digital filtering method (DFM) initiated by [9] attempts to model the spatial and temporal coherence of inflow turbulence through the digital filtering uncorrelated random data, and account for inhomogeneity and anisotropy using the method proposed by [16]. It is relatively easy to implement and is able to reproduce the first and second order one-point statistics as well as auto-correlation function. However, the synthetic turbulence generated by DFM does not satisfy the continuity equation. [17] offered a promising approach to enforce the divergence-free constraint in the DFM 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 DFM to impose a two-point spatial correlation directly is very attractive.
A brief introduction on the filtering method by [9] is stated as follows. In order to create two-point correlations, let \(r_m\) be a series of random data with zero mean and unity variance, then
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
where \(\Delta 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\Delta x\) defined on a one dimensional axis. It is straightforward to tell (1.5) 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.5) with a multidimensional Newton iteration 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
where \(r\) is the distance and \(L\) is the length scale. It is possible to approximately reproduce (1.6) by computing the filter coefficients as
where
The width \(N\) of the filter should be chosen such that \(N\geq 2n\) (where \(n=L / \Delta x_1\)) to ensure the accuracy of the approximation. On the other hand, for an exponential correlation function
It is suggested by [1] to evaluate the filter coefficients using
where
Again, the width \(N\) of the filter should be chosen such that \(N\geq 2n\) (where \(n=L\Delta 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
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):
Step (a). Choose for each coordinate direction corresponding to the inflow plane a length scale \(L_{22} = n_2\Delta x_2\), \(L_{33} = n_3\Delta x_3\), a time scale \(T\) and determine the filter width \(N_{\alpha}\) (\(\alpha =1,2,3\)) accordingly.
Step (b). Initialize and store three random fields \(R_{\alpha}\) (again \(\alpha =1,2,3\)) of dimensions \([-N_1:N_1,-N_2+1:M_2+N_2,-N_3+1:M_3+N_3]\) where \(M_2 \times M_3\) denotes the dimensions of computational gird of the inflow plane.
Step (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.5) meets the target one.
Step (d). Applying the following filter operation for \(j=1,\ldots,M_2\), \(k=1,\ldots,M_3\)
which yields the two-dimensional arrays of spatially correlated data \(\Psi_{\alpha}\), \(\alpha =1,2,3\).
Step (e). Output velocity data with the transformation
where the coefficients \(a_{ij}\) are computed from the Cholesky decomposition of the user-defined Reynolds stress tensor. The detailed information can be found in (1.23). This step ensures the synthetic velocity reproduces the target mean velocity and Reynolds stress tensor.
Step (f). Discard the first \((x_2,x_3)\)-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_1,j,k)\) with new random numbers.
Repeat the steps (d):math:sim (g) for each time step.
If the target correlation function is an exponential function, an alternative approach by [1] can be adopted for generating inflow turbulence which turns out to be much more efficient than the method of [9]. Instead of using the filtering operation discussed above, the method of [1] obtain the temporal correlation with the expression
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\Delta t,j,k)\), \(\varPsi_{\alpha}(t,j,k)\) and \(\varPsi_{\alpha}(t+k\Delta t,j,k)\), respectively. One easily verifies that
which reproduces an exponential function. An overall algorithm for generating the inflow velocity supported by the method of [1] can be stated as follows
Step (a). Choose for each coordinate direction corresponding to the inflow plane a length scale \(L_{22} = n_2\Delta x_2\), \(L_{33} = n_3\Delta x_3\), a time scale \(T\) and determine the filter width \(N_{\alpha}\) (\(\alpha =1,2,3\)) accordingly.
Step (b). Initialize and store three random fields \(R_{\alpha}\) (again \(\alpha =1,2,3\)) of dimensions \([-N_2+1:M_2+N_2,-N_3+1:M_3+N_3]\) where \(M_2 \times M_3\) denotes the dimensions of computational gird in the inflow plane.
Step (c). 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.
Step (d). Applying the following filter operations for \(j=1,\ldots,M_2\), \(k=1,\ldots,M_3\)
which yields the two-dimensional arrays of spatially correlated data \(\varPsi_{\alpha}\), \(\alpha =1,2,3\).
Step (e). Compute \(\Psi_{\alpha}(j,k)\) with (1.15) and output the velocity signal with the transformation
where the coefficients \(a_{ij}\) are given by (1.23). Again, this step ensures the synthetic velocity reproduces the target mean velocity and Reynolds stress tensor.
Step (f). Repeat the Steps (d) - (f) for each time step.
The digital filtering approach implemented in the TInF provides three filter types: Gaussian, exponential and bessel to define the correlation function. Users can either select one filter type through the TInF application which will automatically modify the corresponding OpenFOAM files, or directly specify it in the U file by following the instruction in Code Implementation.
1.4. Synthetic Eddy Method
The synthetic eddy method (SEM) initiated by [2] 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 modelled 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 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 auto-correlation function, and consequently the power spectrum of the synthetic turbulence. Enforcement of the continuity condition in the SEM was discussed in [3] which will be introduced later.
A brief introduction on the SEM presented by [2] 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
The inflow plane on which we want to generate the synthetic turbulence with the SEM is basically a finite set of points \(S = \{\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{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
where
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 of
where \(\mathbf{x}\) represents the coordinates of computational points and \(\mathbf{x}^k\) represents the coordinates of eddies. The number of eddies noted by \(N\) can be set to be \((\text{max}(x_i)-\text{min}(x_i))/\sigma\). The coefficient \(a_{ij}\) results from the Cholesky decomposition of a prescribed Reynolds stress tensor \(R_{ij}\). Specifically, the calculation of \(a_{ij}\) can be expressed in the following
The coefficient \(\epsilon_j^k\) (\(j=1,2,3\)) is the uniformly random intensity factor of values \(+1\) or \(-1\), and \(f_{\mathbf{\sigma}(\mathbf{x})} (\mathbf{x}-\mathbf{x}^k)\) is the velocity distribution at \(\mathbf{x}\) of the eddy located at \(\mathbf{x}^k\) defined as follows:
where \(\mathbf{\sigma}=(\sigma_1,\sigma_2,\sigma_3)^T\). The position of the eddies \(\mathbf{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 \(\mathbf{U}(\mathbf{x})\). At each time step, the new position of eddy \(k\) is given by
where \(\Delta 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
with a new random intensity vector \(\epsilon_j^k\). \(B_{\Delta t}\) denotes the region in which regenerated eddy \(\mathbf{x}^k(t) \in B_{\Delta t}\) dose not affect 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.22) can be expressed as
where the angles denote the mean operator. The independence between the random variables \(\mathbf{x}^k\) and \(\varepsilon_j^k\) in the mean operator implies that
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
The Reynolds stresses \(\langle u_i u_j \rangle\) of the synthesized can be written as
Using again the independence between the random variables \(\mathbf{x}^k\) and \(\varepsilon_j^k\), the above equation is reduced to
The term
verifies that \(\mathbf{x}^k\) follows a uniform distribution over \(B\), i.e.
Finally, we arrive at
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
where \(\mathbf{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.22) and the linearity of the statistical mean, we obtain
Using again the independence between the positions \(\mathbf{x}^k\) and the intensities \(\varepsilon^k\) of the eddies, this yields
By (1.33), the term in the mean operator is written as
Inserting (1.38) back to (1.37) and using (1.24), this yields
where \(*\) denotes the convolution product. For homogeneous turbulence where integral length scales \(\mathbf{\sigma}(\mathbf{x}) = \mathbf{\sigma}(\mathbf{x}+\mathbf{r}) =(\sigma,\sigma,\sigma)^T\), the two-point cross-correlation tensor \(R_{ij}(\mathbf{x},\mathbf{r})\) only depends on \(\mathbf{r}\) and consequently (1.39) is simplified to
Recall the integral length scale \(L_{ij}\) is defined as the integral of the two-point correlation \(R_{ij}(\mathbf{x},\mathbf{r})\) in a particular direction and is thus proportional to \(\sigma\). By integrating (1.40), 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
Recall the convolution theorem for cross-correlation states that
Hence the spatial velocity spectrum tensor can be expressed as
where \(\mathbf{k} = (k_1,k_2,k_3)\). More specifically for instance, the one-dimensional spectra in the \(x\) direction is
Two-time Correlation
The two-time correlation tensor of the velocity, denoted by \(R_{ij}(\mathbf{x},\tau)\), is the correlation between \(u_i(\mathbf{x},t)\) and \(u_j(\mathbf{x},t+\tau)\) at times \(t\) and \(t+\tau\) respectively, i.e.,
By (1.22) and the linearity of the statistical mean, we have
The independence between the position \(\mathbf{x}^k\) and intensity \(\varepsilon_m^k\) of different eddies implies that, for \(k \neq l\), the statistical mean in (1.46) can be split as follows
Consequently (1.46) is reduced to
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\)
If \(\mathbf{x}^k(t)\in B_{\tau}\), then it is going to be recycled between time \(t\) and \(t+\tau\) and hence both \(\mathbf{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 \(\mathbf{x}^k(t) \in B_{\tau}\) to the term in the angles of (1.48) is thus zero. On the contrary if \(\mathbf{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 \(\mathbf{x}^k(t+\tau) =\mathbf{x}^k(t)+\tau\mathbf{U}(\mathbf{x}^k)\). Thus both \(\varepsilon_n^k(t+\tau) = \varepsilon_n^k(t)\) and \(\mathbf{x}^k(t+\tau)\) depend on the previous position \(\mathbf{x}^k(t)\) of eddy \(k\) relative to \(B_{\tau}\). By (1.34) and the definition of \(B_{\tau}\), (1.35) can then be replaced by
The condition \(\mathbf{y}\in B_{\tau}\) leads to \(f_{\mathbf{\sigma}(\mathbf{x})}(\mathbf{x}-(\mathbf{y}+\tau\mathbf{U}))=0\). Thus, the integral over \(B/B_{\tau}\) in (1.39) can be extended to an integral over \(B\). Besides \(\mathbf{y}\in B\) suggests \(f_{\mathbf{\sigma}(\mathbf{x})}(\mathbf{x}-\mathbf{y})=0\) as previously demonstrated, therefore the integral in (1.50) can be further extended to an integral over \(\mathbb{R}^3\). Using (1.24), this yields
In the case where the mean velocity is in the \(x_1\)-direction only \(\mathbf{U} = (U,0,0)\) and the target turbulence is homogeneous, (1.51) is simplified to
Thus, the two-time correlation of the signal at time \(\tau\) is simply the auto-correlation 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}(\mathbf{x},\tau)\) contains can be re-expressed in terms of the wave number velocity spectrum tensor which is written as
Using again the convolution theorem, the above expression is simplified to
Commonly Used Velocity Shape Functions
The three commonly used velocity shape functions for \(f\) which are the tent function, the step function and the truncated Gaussian function are implemented in the TInF. Their definitions are given below for references.
Tent function
which yields
Step function
which yields
Truncated Gaussian function
which yields
where \(C\) is a constant that ensures \(f\) satisfies the normalization (1.19).
1.5. Divergence Free Synthetic Eddy Method
One drawback of the synthetic eddy method is that the inlet velocity field is not necessarily divergent-free, which may result in large pressure fluctuations near the inlet. One route to obtain a divergence free method, as suggested by [3], is to apply the original SEM methodology to the vorticity field, which is then transformed back to the velocity field by taking the curl of it. One easily verifies that vorticity and velocity fields are linked by the following:
Because of the hypothesis of incompressible flow, the first term on the right hand side of the above equation vanishes, leading to a Poisson equation for the velocity field. The solution of this Poisson equation, achieved by using the Biot–Savart kernel, leads to the fluctuating velocity field expressed as:
where \(\mathbf{r}^k=(\mathbf{x}-\mathbf{x}^k)/\sigma^k\), \(q_{\sigma}(|\mathbf{r}^k|)\) is a suitable shape function and \(\alpha_i^k\) are random numbers with zero average which represent the eddy intensities.
Despite the similarities between (1.22) and (1.62), i.e. both including a user defined shape function and a random level eddy intensity, the lack of Lund coefficients in the second formulation poses a significant problem. Their role in the original SEM was crucial in order to allow any given turbulence state to be generated, however, they cannot be re-introduced to (1.62) without forgoing the divergence-free condition.
In order to increase the turbulence anisotropy reproduction capabilities, the method presented by [3] employs the formulation of (1.62), but with an anisotropic length-scale, \(\sigma_i\), employed in each of the coordinate directions, and allows a different shape function to be associated with each direction. However, such a form no longer automatically satisfies the divergence-free condition ensured by (1.62), and further constraints on the shape functions need to be considered in order to retain a divergence-free field. [3] redefine the shape functions to be of the form \(q_{sigma}=q|\mathbf{r}^k|^3\) where \(q\) is a function which depends on the locations of \(\mathbf{x}\) and \(\mathbf{x}^k\), and \(\mathbf{r}^k\) differs slightly from its previous definition, as it now takes into account the length-scale anisotropy \(\mathbf{r}_{\beta}^k=(\mathbf{x}_{\beta}-\mathbf{x}_{\beta}^k)/\sigma_{\beta}^k\). The new general formulation for the velocity fluctuations thus becomes:
In (1.63), the cross product presented in (1.62) has been rewritten using the index notation for tensors, where \(\epsilon_{ijl}\) is the Levi–Civita symbol, and no summation is implied over Greek subscripts. As noted above, with the redefined shape functions, the form of (1.63) no longer automatically satisfies the divergence-free condition. However, on substituting it into the condition that \(\nabla\cdot\mathbf{u}'=0\), a sufficient condition for ensuring a divergence-free velocity field can be found as:
A simple analytically function for \(q_i\) that satisfies the above restrictions is:
where \(d^k = \sqrt{(r_j^k)^2}\).
The function \(q_i\) with the form (1.65) is continuous everywhere, but its derivative is not strictly defined for \(d^k=1\), where it is only possible to define a right or left sided derivative. The above formulation thus defines a divergence-free velocity field everywhere except at the eddy surface (\(d^k=1\)), although this formal omission is not believed to result in serious problems. The expression for the velocity obtained by substituting (1.65) into (1.63) can be written as:
Time-averaging the product of (1.66) with itself leads to an expression for the Reynolds stresses, from which one can examine how the prescription of the length scales, \(\sigma_i^k\), and intensities, \(\alpha_i^k\), affect the stress anisotropy associated with the synthetically generated field given by (1.66):
On examining (1.67), it is clear that the eddies are independent of each other, and that their intensities are uncorrelated (so \(\left<\alpha_l^k\alpha_n^k\right>\) for \(l \neq m\)); as such, the predicted shear stresses (\(\left<u_{\beta}'u_{\gamma}'\right>\) for \(\beta\neq\gamma\)) will be zero. In order to overcome this problem, fluctuations in the global coordinate system are computed via a standard rotation transformation of the eddies generated in the local principal axes coordinate system (where the Reynolds stress tensor is diagonal):
where \(R_{im}^{P\rightarrow G}\) is the rotation transformation matrix from the principal to the global coordinate system, \(u_i^{P}\) and \(u_i^{G}\) are the velocity fluctuations in the principal axes and global systems respectively, and \(C_1\) is a normalization coefficient required in order to have \(\left<(u_i')^2\right>=1\) when \(\left<(\alpha_l^k)^2\right>=1\):
where \(V_0\) is the eddy box volume. For the normal stresses, the contribution from the \(k\)-th eddy in (1.67) thus gives
where all the terms not explicitly reported are represented by \(C_2\). The remaining issue is to choose appropriate length-scales and eddy intensities to ensure the above will return the desired Reynolds stress statistics, over a wide range of stress anisotropy levels. For any choice of length-scale ratios (\(\sigma_1/\sigma_2\) and \(\sigma_1/\sigma_3\)), varying the intensity \(\alpha_l^k\) allows one to reproduce possible turbulence anisotropy states over a particular region of the Lumley triangle. For given length scales, the Reynolds stresses are reproduced by defining the following intensities:
where \(\lambda\) are the normal stresses in the local principal reference system. Since the right hand side of (1.71) must be positive, for any value of \(\Gamma\) it is only possible to reproduce a part of the Lumley triangle. By defining a series of ratios \(\Gamma = \frac{\sigma_1}{\sigma_2} = \frac{\sigma_1}{\sigma_3}\), each allowing one of the regions of the Lumley triangle to be mapped, the Reynolds stress statistics with a wide range of stress anisotropy levels can be reproduced.
1.6. Turbulent Spot Method with Anisotropic Vorton
Earlier versions of the synthetic eddy method and the turbulent spot method were only able to produce turbulence which does not obey the continuity constraint in general. This was only possible for special cases (isotropic or near-isotropic turbulence). The work of [4] utilizes another approach for introducing the anisotropy into the turbulent spots, which obeys continuity and allows to reproduce strong levels of anisotropy at the same time. Their approach is basically a continuation of the vorton formulation described in [18]. The generation is performed in the coordinate system \((x_1,x_2,x_3)\) determined by principle axes of the Reynolds stresses. The Reynolds stresses in any other system \((x_1',x_2',x_3')\) are calculated by \(R_{ij}' = E_{pi}R_{pq}E_{qj}\), where \(E_{ij}\) is the rotation matrix describing coordinate transformation between \((x_1,x_2,x_3)\) and \((x_1',x_2',x_3')\) axes system. Integral lengths in different coordinate systems can be found from the relation:
In (1.72), the vector the vector potential is scaled by a function with spherical symmetry which in the case of the spectrum of decaying turbulence gives an analytic expression:
Also note that other spectra \(E(k)\) could be used in principle. This would result in different shapes of the inner velocity distribution. For the sake of simplicity and because it yields reasonably simple formulas, the work of [4] is restricted to this spectrum.
The spherical symmetry of \(\mathbf{A}(x_1,x_2,x_3)\) is the reason for the isotropy of turbulence generated using these vortons. At this level, anisotropy can be introduced by stretching the coordinates individually, i.e.,
With this, the vector potential and velocity induced by the vorton are now written as:
Note that the multiplication with the coordinates is introduced to make the resulting Reynolds stress tensor a diagonal tensor which is identified with the diagonal matrix of eigenvalues from a principal component analysis of the prescribed Reynolds stress tensor. By aligning the \(x_1\), \(x_2\), \(x_3\)-directions of the vorton with the principal directions of the Reynolds stress tensor, arbitrary anisotropic Reynolds stresses can be reproduced. The vorton sizes \(\sigma_1\), \(\sigma_2\) and \(\sigma_3\) and strength vector components \(\gamma_1\), \(\gamma_2\) and \(\gamma_3\) are free parameters of the vorton and can be used to match the prescribed Reynolds stresses and integral length scales.
Statistical Properties of Anisotropic Vortons
Statistical properties can be analytically derived for homogeneous turbulence. Consider a set of fully uncorrelated vortons, i.e. \(\langle\gamma_{ik}\gamma_{jm}\rangle\) for each pair of \(k\)-th and \(m\)-th vortons with strength components \(i\) and \(j\), respectively. Then the Reynolds stress \(R_{ij}\) of the total field is equal to the sum of Reynolds stresses produced by each vorton
Without loss of generality, the magnitude of the strength is set as unity, i.e. \(|\pm\mathbf{\gamma}|=1\). Then, the expected value of the Reynolds stress \(R_{ij}\) at the point \((0,0,0)\) is
where \(P(\mathbf{x})\) is the probability density function in the event that the vorton is placed at the point \(\mathbf{x}\). For the uniform distribution \(P(\mathbf{x}) = N/V\) is the vorton density. If the computational domain becomes infinite \(N\) should increase, so that the vorton density remains constant:
Substitution of velocity induced by an anisotropic vorton (1.76) in (1.79) results in a simple formula for the principal Reynolds stresses:
Unfortunately, there is no possibility to change the auto-correlation function since the velocities in the form of (1.76) uniquely predetermine them
The calculation gives:
with
Also the one dimensional spectra can be calculated analytically:
As seen, the dependence of the one-dimensional spectra on wave number \(e^{k^2}\) is the same as that in isotropic decaying turbulence.
Determination of Anisotropic Vorton Parameters
Integration of the auto-correlation functions (1.82) reveals a simple and clear interpretation of stretching parameters \(\sigma_i\):
Therefore, the parameters \(\sigma_i\) are uniquely determined from the last formulas \(\sigma_i =L_{ii}/\sqrt{\pi}\). The vorton strength vector \(\gamma\) is found from the condition for Reynolds stresses:
Since the determinant of (1.86) is zero, a solution of the system (1.86) is only possible if the following condition is satisfied:
or
The signs before different terms are independent of each other. Therefore, the integral lengths can not be arbitrary. If two length scales \(L_{11}\) and \(L_{22}\) are prescribed the remaining length should satisfy the conditions above. Particularly, this solution is wrong for the isotropic turbulence since, if \(R_{11} = R_{22} = R_{33}\) and \(L_{11} = L_{33} = L\), the third length is \(L_{22} = L/2\) although all integral lengths based on longitudinal auto-correlation functions should be equal. One possible remedy of the length scale restriction is to superpose two statistically independent anisotropic vortons with different parameters at the same position. It is then possible to prescribe arbitrary combinations of length scales and Reynolds stresses. Two approaches for determination of the free parameters are proposed by [4]:
(a) Analytic Determination (Type R)
In this formulation, it is intended to fulfill all prescribed Reynolds stresses at once. As becomes obvious from (1.87), then at least one length scale has to be constrained. A natural choice is to constrain \(L_{33}\), which is associated with minimum principal Reynolds stress direction \(R_{33}\). Thus, \(L_{11}\) and \(L_{22}\) remain unchanged while \(L_{33}\) follows from:
Finally, the vorton parameter \(\gamma_i\) follows from
or
with \(\sigma_i\) given by (1.85). By determining parameters this way, Reynolds stresses are fulfilled and length scales are only approximately fulfilled. The anisotropic turbulent spot method with this kind of parameter determination is labeled by “Type R” in the work of [4].
(b) Pseudo Inverse (Type L)
In this approach, it is intended to fulfill all prescribed length scales. Although it is not possible to account for all Reynolds stresses at the same time, it is possible to formulate a minimization problem and solve for the closest possible Reynolds stress state. First, the \(\sigma_i\) are computed from the prescribed length scales (1.85). With fixed length scales, the Reynolds stresses follow from the following system of equations, which is unsolvable because of the zero diagonal:
The minimization problem is solved by using the Moore-Penrose pseudo-inverse \(\mathbf{M}_L^+\) of \(\mathbf{M}_L\):
This way, all length scales are fulfilled, but the Reynolds stresses only approximately. The anisotropic turbulent spot method with this kind of parameter determination is labeled by “Type L”. The TInF provides both types of turbulent spot for generating inflow turbulences.
2. Code Implementation
The turbulence inflow tool provides a simple and efficient solution for generating a spatially and temporally correlated turbulent velocity field at the inflow of the computational domain via the open-source code, OpenFOAM. It is developed based on the digital filtering method [1, 9], the synthetic eddy method [2, 3] and the turbulent spot method [4]. The back-end of the turbulence inflow tool is a custom-designed C++ library which contains a list of turbulent velocity boundary conditions developed within OpenFOAM. These boundary conditions can be used to feed the inflow plane with time-varying turbulent velocity signals. Note that the inflow tool provides the source code for two latest standard public versions of the OpenFOAM, i.e., version 6.0 and 7.0, respectively.
The tool can be mainly in two approaches. For users who are quite familiar with OpenFOAM, one can download the source code of the presented boundary conditions, compile the code locally on personal computers (or servers) and then carry out simulations as if those boundary conditions are originally available in the standard OpenFOAM. Similar to other boundary conditions available in OpenFOAM, the use of the presented boundary conditions requires adding some specific entries to the related files in an OpenFOAM project. For users who are not familiar with OpenFOAM, a user-interface (i.e., the front-end of the turbulence inflow tool), which help users to customize the related entries in the OpenFOAM project, is provided.
Before an introduction on the input entries required, it is good to have a preliminary understanding of the file system of a general OpenFOAM project first. To start with, there are usually three folders, i.e., the 0, constant and system folders, inside a valid project. The 0 folder contains the files which specify the initial conditions and boundary conditions for different variable fields, e.g., velocity, pressure and turbulent kinetic energy. The constant folder includes the files storing the geometry information of the mesh and some constant parameters related to the viscosity of the fluid and the turbulence model (if required). The system folder contains files which concern the specification of the time-length (or steps) of simulations, the numerical schemes employed for spatially and temporal discretization, the non-linear equation solution method and other related numerical settings.
As mentioned earlier, the source code of the turbulence inflow tool provides some custom-developed boundary conditions for velocities. To implement those boundary conditions during a simulation, there are several files needs to be either modified or created accordingly. These files include: the U file in the 0 folder, the inflowProperties file in the constant folder and the controlDict file in the system folder. Compared to the U and controlDict files which are required by default in a standard OpenFOAM project, the inflowProperties file is a special item required (when necessary) by the presented boundary conditions. It is designed to store some statistical properties of the target turbulence which include mean velocities, Reynolds stresses and integral length scales.
2.1. Modifications regarding the controlDict file
We now introduce the entries or scripts need to be defined in or added to the above mentioned three files respectively. To start with, we consider the controlDict file (inside the system folder) and an example of this file is given as follows.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
libs ("libturbulentInflow.so");
application pisoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 1e-3;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
The only entry needs to be added to this file is the line libs(“libturbulentInflow.so”). This line tells the fluid solver to include the complied library file libturbulentInflow.so prior to the simulation so that the presented boundary conditions contained in this library become available for the current simulation.
2.2. Modifications regarding the U file
Subsequently, we focus on the U file located inside the 0 folder. It contains the information of the discrete internal field and boundary fields for the velocity at a specific time instant (see the example given below).
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (10 0 0);
boundaryField
{
inlet
{
type turbulentDFMInlet;
filterType exponential;
gridFactor 1;
filterFactor 4;
periodicInY true;
periodicInZ false;
cleanRestart false;
value $internalField;
}
inletOld
{
type fixedValue;
value uniform (10 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //
The boundary conditions for velocities are all specified in the boundaryField dictionary, and there are several sub-dictionaries in the boundaryField dictionary with the names such as inlet, outlet, wall, etc. The name of each sub-dictionary corresponds to the name of a particular boundary patch of the mesh, and the entries contained in each sub-dictionary are the related to the boundary condition for the velocity field at the corresponding boundary patch. Let’s focus on the sub-dictionary associated with the inflow patch, i.e., the inlet sub-dictionary. The commonly used boundary condition for the velocity at the inflow is the fixedValue condition (see the sub-dictionary inletOld). For a fixedValue velocity boundary, the velocities at the inlet patch are constant and fixed to the vector specified inside the baskets (coming after the the values entry) during the simulation. In our case, we would like to the velocities at the inlet patch to be stochastic and time-varying. For this purpose, the developed turbulence inflow package currently provides four boundary conditions in total, i.e. turbulentDFMInlet, turbulentSEMInlet, turbulentDFSEMInlet and turbulentATSMInlet. These four boundary conditions correspond to the synthetic turbulent method by [9], [2], [3] and [4], respectively.
Also note that each presented boundary condition is, as a matter of fact, a derived class of the fixedValue boundary condition. Since the initialization of the fixedValue boundary condition requires the specification of the value entry. Therefore, this entry should be defined for the four presented boundary conditions as well. However, when one of the four presented boundary conditions is employed, the velocities at the boundary will be later overwritten by the generated turbulence velocity field and consequently the value entry does not have any effect here. Thus, this entry can be just specified as $internalField; Once the type entry takes the turbulentDFMInlet, turbulentSEMInlet, turbulentDFSEMInlet or turbulentATSMInlet boundary condition, there are some unique entries related to each boundary condition accordingly.
Apart from the type and value entries, there are also some identical entries shared by the four presented boundary conditions. The entries periodicInY and periodicInZ determine whether the synthetic turbulence is periodic in \(x_2\) and \(x_3\)-direction or not, respectively. Before we discuss the cleanRestart entry, it should be mentioned that the four presented boundary conditions are capable of restarting the generation of the synthetic turbulence from an old time-step. This feature is achieved by storing the necessary information and numerical quantities (obtained at the current time-step) which will affect the generation of the turbulence at the next-step. The entry cleanRestart determines whether to disregard those information and generate a new turbulent velocity field at the beginning of a new simulation.
entry name |
type |
descriptions |
default values |
|---|---|---|---|
periodicInY |
bool |
if the synthetic turbulent is periodic in \(x_2\)-direction |
false |
periodicInZ |
bool |
if the synthetic turbulent is periodic in \(x_3\)-direction |
false |
cleanRestart |
bool |
whether to disregard old turbulence or not |
false |
The turbulentDFMInlet boundary condition
For the turbulentDFMInlet boundary condition, the unique entries to be specified are filterType, filterFactor and gridFactor. The filterType entry refers to the type of the prescribed function for computing the filter coefficients, and requires a string input. When this entry is taken as gaussian, the coefficients are computed with (1.6); when this entry is taken as exponential, the coefficients follow from (1.10). As discussed earlier, the digital filtering method should be employed on uniform spacing Cartesian grids. To make this method to be applied to more general cases, a virtual uniform spacing Cartesian grid system will be generated during the implementation. The velocity fluctuations will be first generated on this virtual grid and then interpolated to the discrete points on the inflow patch. The gridFactor entry defines ratio between the virtual grid spacing \(\Delta\) with the square root of the area of the smallest face element on the inflow plane. Higher gridFactor leads to larger grid spacing and vice versa. Finally, the filterFactor entry denotes the values of \(N/n\) where \(N\) and \(n\) follow from the notation in Section 1.3, and consequently it requires an integer input. The entries related to the turbulentDFMInlet boundary, their input variable types and limitations are listed in Table 2.2.
entry name |
type |
descriptions |
default values |
|---|---|---|---|
filterType |
string |
gaussian, exponential or bessel |
exponential |
gridFactor |
float |
\(\geq 1\) |
1 |
filterFactor |
int |
\(\geq 4\) |
4 |
The turbulentSEMInlet boundary condition
For the turbulentSEMInlet boundary condition, the unique entries to be specified are eddyType and eddyDensity. The eddyType entry refers to the type of the velocity shape function \(f\). It takes an string input and currently the available options are tent, step and gaussian which correspond to (1.55), (1.57) and (1.59), respectively. The eddyDensity refers to the ratio between \(V_B\) and the sum of the eddy volumes (see Section 1.4). Higher eddyDensity leads to a larger amount of synthetic eddies. It takes an scalar input with a value no less than one. The lower bound of eddyDensity (taken as 1 by default) ensures the eddy box (introduced in Section 1.4), from a statistical point of view, can be covered up by the eddies. The unique entries for the turbulentSEMInlet condition are summarized in Table 2.3.
entry name |
type |
descriptions |
default values |
|---|---|---|---|
eddyType |
string |
tent, step or gaussian |
gaussian |
eddyDensity |
float |
\(\geq 1\) |
1 |
The turbulentDFSEMInlet boundary condition
For the turbulentDFSEMInlet boundary condition, the only unique entry to be specified is eddyDensity whose meaning is identical to that of the turbulentSEMInlet boundary condition.
The turbulentATSMInlet boundary condition
For the turbulentATSMInlet boundary condition, the unique entries to be specified are vortonType and vortonDensity. The vortonType entry refers the approach parameter determination for anisotropic vortons discussed in Section 1.6. It takes an string input and the options currently available are typeR and typeL which correspond to (1.57) and (1.59), respectively. The vortonDensity entry is similar to the eddyDensity of the turbulentSEMInlet boundary condition. We summarize the unique entries for the turbulentATSMInlet condition in Table 2.4 for reference.
entry name |
type |
descriptions |
default values |
|---|---|---|---|
vortonType |
string |
typeR or typeL |
typeR |
vortonDensity |
float |
\(\geq 1\) |
1 |
A common feature shared by the entries filterFactor, eddyDensity and vortonDensity mentioned above is that when they are assigned with higher values, the resulting velocity fluctuations will usually have better a quality. However higher values also means a larger consumption of mathematical calculations and computational memory. Therefore, users are suggested to balance the needs between accuracy and efficient during the implementation of those boundary conditions.
2.3. Specification of the statistics of the target turbulence
The aforementioned entries defined in the boundary conditions for velocities only concern the selection of the method and the corresponding parameters associated with the selected method. The generation of the synthetic turbulence also requires the statistical information of the target turbulence to be reproduced which include mean velocities, Reynolds stresses and integral length scales. There are mainly three approaches for the input of those information, i.e., direct specification, interpolation and prescribed function.
Direct specification
The direct specification approach is suitable for the case in which the mean velocity, Reynolds stress and length scales on each face element of the inflow plane are known in prior. In this approach, the values of those quantities can be directly specified in the inlet dictionary. If they are uniformly distributed on the inflow plane, the related entries are in the form of
boundaryField
{
inlet
{
type turbulentDFMInlet;
filterType exponential;
gridFactor 1;
filterFactor 4;
periodicInY true;
periodicInZ false;
cleanRestart false;
value $internalField;
U uniform 10;
R uniform (2.0 0.5 0.5 1.5 -0.5 1.0);
L uniform (0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1);
}
}
The three entries displayed in the inlet dictionary are U, R and L which represent mean velocity magnitude, Reynolds stress and length scale. If U, R and L are not uniformly distributed on the inlet patch, the related entries should be defined in the form of
boundaryField
{
inlet
{
type turbulentDFMInlet;
filterType exponential;
gridFactor 1;
filterFactor 4;
periodicInY true;
periodicInZ false;
cleanRestart false;
value `internalField;
U
{
1.0
2.0
...
6.0
};
R
{
(2.0 0.5 0.5 1.5 -0.5 1.0)
(2.1 0.6 0.6 1.6 -0.6 1.1)
...
(2.5 1.0 1.0 2.0 -1.0 1.5)
};
L
{
(0.30 0.30 0.30 0.20 0.20 0.20 0.10 0.10 0.10)
(0.31 0.31 0.31 0.21 0.21 0.21 0.11 0.11 0.11)
...
(0.35 0.35 0.35 0.25 0.25 0.25 0.15 0.15 0.15)
};
}
}
It is noted that each element in R defines a six-component symmetric tensor of the form \((R_{11} \ R_{21} \ R_{31} \ R_{22} \ R_{32} \ R_{33})\). For the turbulentDFMInlet and turbulentSEMInlet boundary conditions, each element in L defines a nine-component tensor of the form (\(L_{11}^{x_1}\) \(L_{11}^{x_2}\) \(L_{11}^{x_3}\) \(L_{22}^{x_1}\) \(L_{22}^{x_2}\) \(L_{22}^{x_3}\) \(L_{33}^{x_1}\) \(L_{33}^{x_2}\) \(L_{33}^{x_3})\). For the turbulentATSMInlet boundary condition, each element in L defines a three-component vector of the form \((L_{11}^{x_1}\) \(L_{22}^{x_2}\) \(L_{33}^{x_3})\). For the turbulentDFSEMInlet boundary condition, each element in L is a scalar. The main difficultly in specifying the entries U, R and L directly is to make sure that the sequence of the elements in each entry is properly sorted coping with the corresponding faces on the inflow plane.
Specification via interpolation
The mean velocities, Reynolds stresses and integral length scales on the cell-faces of the inflow plane can also be specified through face interpolation. The detailed procedures are as follows:
Create a folder with the name boundaryData inside the constant folder of a standard OpenFOAM project.
Fig. 2.1 The boundaryData folder
Create another folder using the name of the inflow patch inside the above created folder.
Fig. 2.2 The inlet folder
Add four empty files with the names U, R, L and points, respectively inside the inlet folder. These four files will be used to store the given values of the mean velocity magnitude, Reynolds stress, integral scale at some specific points and the locations of those points in the global coordinate system (of the computational domain).
The U file, which stores the information of mean velocity magnitudes, should be written in the format of
(
0.000000
0.029538
0.118110
...
0.029538
0
)
The entries within the braces are the magnitudes of mean velocity at the corresponding points.
The R file, which stores the information of Reynolds stresses, should be written in the format of
(
(0.000000 0.000000 0.000000 0.000000 0.000000 0.000000)
(0.000137 0.000000 0.000000 0.000000 0.000000 0.000053)
(0.002183 0.000002 0.000001 0.000000 0.000000 0.000827)
...
(0.000137 0.000000 0.000000 0.000000 0.000000 0.000053)
(0.000000 0.000000 0.000000 0.000000 0.000000 0.000000)
)
The entries within the braces are the Reynolds stress tensors at the corresponding points.
The L file, which stores the information of turbulence length scales, should be written in the format of
(
(0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000)
(0.000024 0.000024 0.000024 0.000024 0.000024 0.000024 0.000024 0.000024 0.000024)
(0.000096 0.000096 0.000096 0.000096 0.000096 0.000096 0.000096 0.000096 0.000096)
...
(0.000024 0.000024 0.000024 0.000024 0.000024 0.000024 0.000024 0.000024 0.000024)
(0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000)
)
The entries within the braces are the components of the turbulence length scales at the corresponding points. Again, both the turbulentDFMInlet and turbulentSEMInlet boundary conditions employ a nine-component length scale. The turbulentATSMInlet and turbulentDFSEMInlet boundary conditions employ a three-component and one-component length scale, respectively.
The points file, which stores the locations of the points in the global coordinate system, should be written in the format of
(
(0.000000 0.000000 0.000000)
(0.000000 0.000075 0.000000)
(0.000000 0.000301 0.000000)
...
(0.000000 2.000000 0.000000)
(0.000000 2.000000 3.000000)
)
It is noted that the points listed in this file should be able to define a single plane. That is to say the points should not locate on a single line or different planes. In summary, specification of statistics of target turbulent via interpolation is suitable for the case in which the mean velocities, Reynolds stresses and integral length scales at a group of selected points are available.
Specification via prescribed functions
Finally, the specification of the mean velocities, Reynolds stresses and integral length scales can also be done with the prescribed functions embedded in the presented boundary conditions. The entries related to this approach are all defined in the inflowProperties file in the constant folder. To employ this approach, the first step is to define a local coordinate system and its origin. The local \((x_1',x_2',x_3')\) coordinate system is constructed with its \(x_1'\)-axis being parallel to normal of the inflow plane and pointing towards the interior of the computational domain.
Fig. 2.3 A sketch of the Euler angles and the \(N\)-axis
To determine the orientation of the \(x_2'\)- and \(x_3'\)-axis, an entry named as Naxis is defined which takes a vector input. This vector represents the direction of the \(N\)-axis (see Fig. 2.3) where \(N\) is the intersection line of the \((x_1,x_2)\)-plane and \((x_1',x_2')\)-plane. The \(x_2'\)- and \(x_3'\)-axis are then determined by \(\mathbf{x}_2' = \mathbf{x}_3'\times\mathbf{x}_1'\) and \(\mathbf{x}_3' = \mathbf{x}_1'\times\mathbf{N}\). The reason why we specify the direction of the \(N\)-axis instead of the \(x_2'\)- or \(x_3'\)-axis is because, once the direction of the \(x_1'\)-axis is determined, the \(x_2\)- and \(x_3\)-axis cannot be specified arbitrarily since they must be parallel to the plane orthogonal to the \(x_1\)-axis. In contrast, the direction of the \(N\)-axis can be specified arbitrarily as long as the \(x_3\)-components of \(N\) is zero. The entry Naxis is set as \((0 \ 0 \ 0)\) by default. In this case, the \((x_1',x_2')\)-plane and \((x_1,x_2)\)-plane are parallel to each other and consequently the \(x_3\)- and \(x_3'\)-axis are align.
Now, the remaining issue is the determination of the origin of the local coordinate system. Consider a rectangular defined on the \((x_2,x_3)\)-plane bounding all the nodes of the inflow patch. The origin of the local coordinate system is chosen as the lower left corner of this rectangular by default. We also allow users to offset the origin through the use of the entries yOffset and zOffset which take scalars for input.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 6
\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object inflowProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Naxis (0 0 0);
offset (0 0 0);
UDict
{
referenceValue 1;
profile uniform;
}
RDict
{
referenceValue (0.1 0 0 0.1 0 0.1);
profile uniform;
}
LDict
{
referenceValue (0.3 0.2 0.1);
profile uniform;
Gamma (1 1 1);
}
// ************************************************************************* //
In addition to the entries mentioned above, there are a total of three sub-dictionaries need to be specified in the inflowProperties file, i.e. UDict, RDict and LDict, which are responsible for the specifications of the mean velocity magnitude, the Reynolds stress tensor \((R_{11} \ R_{21} \ R_{31} \ R_{22} \ R_{32} \ R_{33})\) and the integral length scales. For the turbulentDFMInlet and turbulentSEMInlet boundary conditions, the integral length scales should be defined in the form of (\(L_{11}^{x_1}\) \(L_{11}^{x_2}\) \(L_{11}^{x_3}\) \(L_{22}^{x_1}\) \(L_{22}^{x_2}\) \(L_{22}^{x_3}\) \(L_{33}^{x_1}\) \(L_{33}^{x_2}\) \(L_{33}^{x_3})\). For the turbulentATSMInlet boundary condition, the integral length scales should be defined in the form of \((L_{11}^{x_1} \ L_{22}^{x_2} \ L_{33}^{x_3})\). Note that \(L_{11}^{x_1}\) refers to the integral length scale obtained by integrating the two-point correlation function \(R_{11}(\mathbf{x})\) with respect to the local \(x_1\)-direction (i.e., the stream-wise direction) and similarly for \(L_{22}^{x_2}\) and \(L_{33}^{x_3}\).
Each sub-dictionary contains the entries required for the computation of the values of the corresponding variable on the inlet boundary which are almost identical to each other. In all sub-dictionaries mentioned above, the two basic entries need to be specified are the referenceValue and the profile. For the UDict sub-dictionary, the referenceValue entry requires a scalar value for input, while a symmetric tensor is required for the RDict sub-dictionary. Finally, the referenceValue entry of the LDict sub-dictionary demands a vector.
On the other hand, the profile entry, as its name implies, defines the profile function of the corresponding variable on the inlet boundary. It takes a string variable for input, and there are only two the valid options are uniform and exponential at the present stage (expecting more in the coming future). When the profile entry takes uniform, the corresponding variable will be set universally to the reference value on the entire inflow patch; When the profile entry takes exponential, the corresponding variable will be therefore computed with an exponential function and additional entries entitled as alpha, referenceAngl and referenceDist should be defined as a supplement, see examples below.
UDict
{
referenceValue 1;
profile exponential;
referenceAngl 0;
referenceDist 1.0;
alpha 0.1;
}
RDict
{
referenceValue (0.1 0 0 0.1 0 0.1);
profile exponential;
referenceAngl 0;
referenceDist 1.0;
alpha (0.3 0.2 0.1);
}
While the entries required for the same type of profile function are almost the same, the form of the exponential function differs in scalar variables (e.g., the mean velocity magnitude) and symmetric tensor variables (e.g., the turbulent intensity). For a scalar variable namely \(\phi\), the corresponding exponential function has the form
where \(\bar{\phi}\) is the reference value (defined in the referenceValue entry), \(\alpha\) is the exponential coefficient (defined by the alpha entry), and \(d_0\) is a reference distance (defined by the referenceDist entry). The symbol \(\mathbf{n}\) is a direction vector (located on the inflow plane) to which the spatial coordinate \(\mathbf{x}\) (defined in the local \((x_1,x_2,x_3)\) coordinate system) is projected. The direction of \(\mathbf{n}\) is specified indirectly by defining the relative angle between the local \(x_3\)-axis and \(\mathbf{n}\), and this angle (ranging from 0 to 180 degrees) is specified in the referenceAngl entry. This reason behind such an approach is obvious, the \(\mathbf{n}\) can not be arbitrarily specified for it is located on the local \(x_1x_2\) (i.e., the inflow) plane. Also for obvious reasons, the reference values of the mean velocity magnitude and integral length scales should be larger than zero.
Finally, a symmetric tensor variable (e.g., the turbulent intensity) denoted by \(\mathbf{\phi}\) for demonstration, the corresponding exponential function writes
where \(\bar{\phi}_{\gamma}\) (\(\gamma=1,2,3\)) are principal values of the reference symmetric tensor \(\bar{\mathbf{\phi}}\) (i.e., the one defined in the referenceValue entry). \(\mathbf{M}_{\gamma}\) (\(\gamma=1,2,3\)) are tensors defined by
where \(\mathbf{a}_{\gamma}\) are the eigenvectors of \(\bar{\mathbf{\phi}}\). Other symbols in (2.2) follow from the notation in (2.1). Also note that in (2.2), the alpha entry requires a three-component vector, i.e., \((\alpha_1 \ \alpha_2 \ \alpha_3)\), for input. The referenceValue requires a six-component vector, i.e., \((\bar{\phi}_{11} \ \bar{\phi}_{21} \ \bar{\phi}_{31}\ \bar{\phi}_{22} \ \bar{\phi}_{32}\ \bar{\phi}_{33})\), for input, and for obvious reasons the \(\bar{\phi}_{11}\), \(\bar{\phi}_{22}\) and \(\bar{\phi}_{33}\) components should be larger than zero.
entry name |
type |
limits |
default values |
|
|---|---|---|---|---|
referenceValue |
UDict |
float |
>0 |
|
LDict |
a vector of floats |
>0 |
||
RDict |
a vector of floats |
positive definite |
||
profle |
string |
uniform, linear or exponential |
uniform |
|
referenceDist |
float |
>0 |
1 |
|
referenceAngl |
float |
(0 180) |
0 |
|
Naxiss |
float |
(0 180) |
0 |
|
2.4. References
- 1(1,2,3,4,5)
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.
- 2(1,2,3,4,5)
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.
- 3(1,2,3,4,5,6)
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.
- 4(1,2,3,4,5,6)
Hannes Kröger and Nikolai Kornev. Generation of divergence free synthetic inflow turbulence with arbitrary anisotropy. Computers and Fluids, 165:78–88, 2018.
- 5
Fei Ding and Ahsan Kareem. Generation of inflow velocity field for CFD analyses using GPUs. In Proceedings of the 15th International Conference on Wind Engineering. 2019.
- 6
Xiaohua Wu. Inflow turbulence generation methods. Annual Review of Fluid Mechanics, 49:23–49, 2017.
- 7(1,2)
Robert H Kraichnan. Diffusion by a random velocity field. The physics of fluids, 13(1):22–31, 1970.
- 8(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.
- 9(1,2,3,4,5,6)
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.
- 10
Masanobu Shinozuka and C-M Jan. Digital simulation of random processes and its applications. Journal of sound and vibration, 25(1):111–128, 1972.
- 11
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.
- 12
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.
- 13
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.
- 14
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.
- 15
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.
- 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
Kornev Nikolai and Hassel Egon. Synthesis of homogeneous anisotropic divergence-free turbulent fields with prescribed second-order statistics by vortex dipoles. Physics of Fluids, 19(6):68–101, 2007.