An efficient and low-divergence method for generating inhomogeneous and anisotropic turbulence with arbitrary spectra

Abstract

This paper proposes a low-divergence method to generate inhomogeneous and anisotropic turbulence. Based on the idea of correlation reconstruction, the proposed method uses the Cholesky decomposition matrix to re-establish turbulence correlation functions. This is in contrast with the time-consuming procedure used in conventional methods to solve eigenvalues and eigenvectors in each location required by coordinate transformation, thereby reducing the computational complexity and improving the efficiency of synthetic turbulence generation. The proposed method is based on the classical spectrum-based method that is widely used to generate homogeneous and isotropic turbulence. By adjusting the generation strategy of specific random vectors, inhomogeneous and anisotropic turbulence with a relatively low divergence level can be obtained in practice, with almost no additional computational burden. Two versions of the proposed method are presented: the inverter and shifter versions. Both the versions are highly efficient, easy to implement, and compatible with high-performance computing. This method is also suitable for providing high-quality initial or boundary conditions for scale-resolving turbulence simulations with large grid numbers (such as direct numerical simulations or large eddy simulations); it can be rapidly implemented using various open-source computational fluid dynamics codes or common commercial software.

1. Introduction

With the rapid development of high-performance computing (HPC) and related hardware capabilities, scale-resolving turbulence simulations, such as direct numerical simulations (DNSs) and large eddy simulations (LESs), have been gradually applied to increasingly complex geometries and working conditions (Garnier, Adams & Sagaut 2009). Despite the considerable progress in this field, acquiring the initial or boundary conditions of general turbulence remains an unavoidable problem and has not yet been perfectly solved (Tabor & Baba-Ahmadi 2010; Dhamankar, Blaisdell & Lyrintzis 2017).

The methods used to acquire the initial and boundary conditions for highly resolved turbulence simulations can be divided into two general categories: the database methods and synthetic turbulence methods. Detailed descriptions of these categories can be found in the literature (Dhamankar et al. 2017). The database methods can obtain field data that are consistent with or close to real turbulence with high precision and accuracy by directly using experimental data, storing the results of a DNS, or constructing real-time mapping from additional computational domains (Lund, Wu & Squires 1998; El-Askary, Schroeder & Meinke 2003). However, the scope of application for these methods is very limited. The data storage and access burden of these methods are typically large to preserve the original information of real turbulence as much as possible, and sometimes, the corresponding computational cost is of the same order as the main computational domain of focused turbulent flows. Most methods are not well-adapted to parallel programming. Moreover, a single method is typically only applicable to a specific type of turbulence; therefore, meeting the requirements of generating more complex turbulence in engineering may be difficult. Some previous studies have attempted to build turbulence generators based on databases through deep learning (Yousif, Yu & Lim 2022); however, the synthetic turbulence methods are undoubtedly more direct and effective at solving this general problem. Synthetic turbulence methods are used to synthesize turbulence based on particular information to provide the necessary initial and boundary conditions. They cover a wide range of specific method types, including the spectral-representation-based method (Kraichnan 1970; Smirnov, Shi & Celik 2001), digital-filter-based method (Klein, Sadiki & Janicka 2003; Kempf, Klein & Janicka 2005; Kim, Castro & Xie 2013), synthetic-eddy method (Benhamadouche et al. 2006; Jarrin et al. 2006; Mathey et al. 2006; Subbareddy et al. 2006), turbulent spot method (Kornev & Hassel 2007; Kröger & Kornev 2018) and diffusion-based method (Kempf et al. 2005). The spectrum-based method is widely applied to many in-house codes and general computational fluid dynamics (CFD) software, owing to its valid spatial correlation and scale-distribution property, which can realize arbitrary spectra. For example, both the open-source code OpenFOAM (ESI-OpenCFD 2021) and the commercial software Ansys Fluent (2016) have implemented versions of the spectrum-based method.

One of the most important requirements for synthetic turbulence is ensuring a divergence-free or solenoidal property condition in an incompressible flow. Synthetic incompressible turbulence that does not fulfil this condition may rapidly decay or cause excessive pressure fluctuations (Poletto, Craft & Revell 2013). However, for compressible flows, the artificially introduced velocity divergence produces spurious noise, which may interfere with other shocks or expansion waves in high-speed flows (Dhamankar et al. 2017). Many types of synthetic-turbulence methods, such as digital filters (Klein et al. 2003) and synthetic-eddy methods (Jarrin et al. 2006), in their widely used classical forms do not naturally generate divergence-free results because they primarily focus on producing precise spatial/temporal correlations or physical statistics. Although some studies have attempted to improve the digital-filter or synthetic-eddy methods to meet the solenoidal requirements (Kim et al. 2013; Poletto et al. 2013), these improvements are typically limited by the methodological framework. For example, the improved method of Poletto et al. (2013) was derived only based on homogeneous anisotropic turbulence and inherits some difficulties of the synthetic-eddy method by nature in specifying the length scale and shape of turbulent structures (Dhamankar et al. 2017). Most improved versions of these methods either lose some of the advantages that the classical versions provide or attach more limitations to the methods to ensure the divergence-free property. For example, the improved method proposed by Kim et al. cannot simultaneously meet the requirements of maintaining a constant mass flux and ensuring a divergence-free property (Kim et al. 2013). In the method proposed by Kempf et al. (2005), if the projection step (Lee, Lele & Moin 1992) is employed to correct the divergence error, problems arise with the boundary conditions, and desired correlations may deviate under certain circumstances. The vortex method using the Biot–Savart law can guarantee divergence-free conditions only if there are no streamwise fluctuations (Mathey et al. 2006).

In contrast to other synthetic methods, the classical form of the spectrum-based method generates homogeneous and isotropic turbulence and produces only a divergence-free result. However, subsequent improvements that extend it to cover inhomogeneous and anisotropic turbulence have eliminated this solenoidal property. The method proposed by Smirnov et al. (2001) is one of the most widely used spectrum-based methods for handling inhomogeneity and anisotropy. This method uses a coordinate transformation and scaling operation to achieve an inhomogeneity-and-anisotropy extension; however, the divergence of the resulting field can no longer be ensured if the target turbulence is not homogeneous or isotropic. A recent study showed that even the approximate zero divergence claimed by Smirnov et al. does not hold in many cases (Yu & Bai 2014). Moreover, to obtain the orthogonal matrix for coordinate transformation, all eigenvalues and eigenvectors of the Reynolds stress tensor must be calculated for each grid point, which significantly increases the computational cost of this method type. Yu et al. also introduced coordinate transformation based on the extended version by Smirnov et al., and they improved the method to ensure divergence-free results using the curl of a vector potential field (Yu & Bai 2014). Although their method fulfils the solenoidal requirement, the time-consuming eigenvalue and eigenvector calculations are also inherited from Smirnov's method because coordinate-transform matrices are still required for every grid point. In addition, the introduction of the vector potential makes their method fundamentally different from the conventional spectrum-based method, and the computational complexity is significantly increased. Therefore, most of the existing spectrum-based-method code frameworks cannot be sufficiently modified to implement this new version, which considerably limits its practical application. A more recent progress was made by Patruno & Ricci (2018). Their improved spectrum-based method combines the ideas of Poletto et al. (2013) and Kröger & Kornev (2018) and innovatively corrects the divergence-free error caused by directional anisotropy. However, the original $N$ modes superposition is transformed to an $M \times N$ modes superposition in the improved method, which clearly increases the computational complexity. Moreover, the divergence-free error caused by strong inhomogeneity may be amplified by the superposition of $M$ modes compared with other previous methods (Bervida et al. 2020).

This study proposes a new divergence-free spectrum-based method, which is based on Cholesky decomposition for correlation reconstruction instead of coordinate transformation, to solve the problem of inhomogeneity and anisotropy extension. This method acquires divergence-free turbulence by correcting the calculation strategy of specific defective steps according to divergence-free error terms that result from correlation reconstruction; thus, it completely inherits the framework and advantages of the classical spectrum-based method. The algorithm of this new method is quite concise, requiring the modification of only two steps of the original method to ensure the divergence-free property in inhomogeneous and anisotropic flow. Because the conventional coordinate transformation is eliminated, no eigenvalue and eigenvector solving algorithms are involved. Thus, a faster computation speed can be achieved when implementing a large number of grid points, particularly in high-resolution DNSs/LESs. The effectiveness of the proposed method was verified using four types of test cases: homogeneous and isotropic turbulence, anisotropic turbulence, inhomogeneous turbulence, and typical inhomogeneous and anisotropic turbulence of the channel flow.

The remainder of this paper is organized as follows. Section 2 explains the theory and derivation of the proposed method in detail. Section 3 analyses the algorithms and computational complexity. In § 4, four representative cases are employed to verify the proposed method and thoroughly test its practical performance. Section 5 presents the concluding remarks of the study.

2. Theory and methodology

2.1. Methods for generating homogeneous and isotropic turbulence

Adopting the notation of Bechara et al. (1994), the velocity vector in Cartesian space for a single sampling can be expressed in a Fourier series as follows:

(2.1)\begin{equation} u_i\left(x_{(k)},t\right) = \sum_n 2 p^{(n)}\cos\left(\kappa^{(n)}_j x_j + \omega^{(n)} t + \varphi^{(n)}\right)\sigma^{(n)}_i, \end{equation}

where the superscript $^{(n)}$ represents the $n$th Fourier mode and does not participate in the tensor summation. Mode phase $\varphi ^{(n)} \sim U(0,2{\rm \pi} )$ and mode amplitude $p^{(n)} = \sqrt {E_{k}(\kappa ^{(n)}) \Delta \kappa ^{(n)}}$. The energy spectrum can be obtained using experimental data such as CBC data (Comte-Bellot & Corrsin 1971), which are widely used in the calculations of homogeneous and isotropic turbulence decay. When computing each mode in the algorithm used to generate the fluctuation velocity, the unit wavevector $\hat {\kappa }_i^{(n)}$ of the mode is first obtained by generating a random unit vector. The symbol $\hat {}$ represents the corresponding unit vector after normalization; for example, $\kappa _i^{(n)} = \kappa ^{(n)} \hat {\kappa }_i^{(n)}$. The vector product or cross product is then used to obtain the mode direction perpendicular to the wave vector,

(2.2)\begin{equation} \sigma^{(n)}_i =\left(\varepsilon_{ijk} \xi^{(n)}_j \hat{\kappa}^{(n)}_k\right)_{normalize}, \end{equation}

where $\xi ^{(n)}$ is a random unit vector that follows the spherical distribution as well. If we let $\alpha ^{(n)}(x_{(k)},t) = \kappa ^{(n)}_j x_j + \omega ^{(n)} t + \varphi ^{(n)}$, we can obtain $\alpha_{i}^{(n)} = \kappa _i^{(n)}$ and $\partial _t \alpha ^{(n)} = \omega ^{(n)}$.

For a given time instant, the divergence of the turbulent velocity field can be expressed as

(2.3)\begin{equation} u_{i,i} =\sum_n 2 p^{(n)} \sin \alpha^{(n)} \kappa^{(n)}_i \sigma^{(n)}_i. \end{equation}

As long as the wavevector and direction vector are strictly perpendicular, the generated turbulence satisfies the divergence-free condition.

2.2. Extended method for inhomogeneous and anisotropic turbulence

The generation method represented by (2.1) only accounts for homogeneous and isotropic turbulence and requires an extension technique to include inhomogeneous turbulence. A typical approach is to introduce coordinate transformations and scaling operations (Smirnov et al. 2001; Yu & Bai 2014). However, the orthogonal matrix computation for coordinate transformation requires the calculation of all eigenvalues and eigenvectors of the Reynolds stress tensor at each spatial point, where synthetic turbulence is required. Consequently, the computational cost increases significantly. To ensure the efficiency of the computation, this study adopts a correlation reconstruction transformation based on Cholesky decomposition matrices. The velocity is generated using the following formula:

(2.4)$$\begin{gather} v_i\left(x_{(k)},t\right) = \sum_n 2 q^{(n)}\cos\left(\kappa^{(n)}_j x_j + \omega^{(n)} t + \varphi^{(n)}\right) \sigma^{(n)}_i, \end{gather}$$

(2.5)$$\begin{gather}u_i\left(x_{(k)}, t\right) = L_{ij}\left(x_{(k)}, t\right) v_j, \end{gather}$$

where $q^{(n)} = p^{(n)}/u_t$. $u_t$ is the characteristic velocity of turbulence corresponding to the integration scale, and it can be calculated using $\sqrt { \frac 13{\langle u_i u_i \rangle } }$.

The correlation reconstruction matrix $L_{ij}$ is obtained based on the Cholesky decomposition of the Reynolds stress tensor $R_{ij}$, that is, $R_{ij} = L_{ik} L_{jk}$. Therefore,

(2.6)\begin{equation} \left.\begin{gathered} L_{11} = \sqrt{R_{11}},\quad L_{21} = \frac{R_{21}}{L_{11}},\quad L_{22} = \sqrt{R_{22} - L_{21}^2}, \\ L_{31} = \frac{R_{31}}{L_{11}},\quad L_{32} = \frac{R_{32} - L_{31} L_{21}}{L_{22}},\quad L_{33} = \sqrt{R_{33} - L_{31}^2 - L_{32}^2}, \\ L_{12} = L_{13} = L_{23} = 0. \end{gathered}\right\} \end{equation}

The correlation reconstruction operation presented in (2.4)–(2.6) realizes the extension of the original spectrum-based method to inhomogeneous and anisotropic turbulence. It has been adopted in some previous studies for different types of generation methods to achieve similar purposes (Shur et al. 2014; Patterson, Balin & Jansen 2021).

2.3. Divergence-free errors caused by correlation reconstruction

Although the correlation reconstruction matrix is relatively concise, this extended method is limited because the divergence is no longer strictly zero. After the derivation, the velocity divergence at this point is

(2.7)\begin{align} u^{(n)}_{i,i} &= 2 \sum_n\left[\underbrace{\cos\alpha^{(n)} q^{(n)}_{,i} L_{ij} \sigma^{(n)}_j }_{Er_1} +\underbrace{\cos\alpha^{(n)} q^{(n)} L_{ij} \sigma^{(n)}_{j,i} }_{Er_2} -\underbrace{\sin\alpha^{(n)} q^{(n)} L_{ij} \sigma^{(n)}_{j} \kappa^{(n)}_{k,i} x_{k} }_{Er_3} \right. \nonumber\\ &\quad \left.+\underbrace{\cos\alpha^{(n)} q^{(n)} L_{ij,i} \sigma^{(n)}_j }_{Er_4} -\underbrace{\sin\alpha^{(n)} q^{(n)} \kappa_i^{(n)} L_{ij} \sigma^{(n)}_{j} }_{Er_5}\right]. \end{align}

There are five error terms in (2.7), which are responsible for the non-zero divergence in inhomogeneous and anisotropic turbulence. The first term $Er_1$ arises from the inhomogeneity of the amplitude $q^{(n)}$ and is present in almost all improved spectrum-based methods for inhomogeneous turbulence. Because $Er_1$ is calculated from the normalized energy spectrum, it is primarily affected by the change in the turbulence integral scale in space. However, its inhomogeneity effect is limited. As shown in § 4.4, the results of a spectrum using a constant scale and inhomogeneous scale are similar to each other. Only in a statistical way through the entire computational domain could a small difference be observed. This implies that the effect of this term is minimal compared with that of other terms and can be considered negligible. In the previous method, the second and third terms are strictly zero. After modifications, the corresponding errors are introduced only if the unit direction vector $\sigma _i^{(n)}$ or wavevector $\kappa _i^{(n)}$ exhibits strong variation. Here $Er_2$ and $Er_3$ are essentially bound to additional adjustments; therefore, it is difficult to design the correction in advance to ensure that these two terms naturally become zero. Fortunately, most of the improvements result in both of these errors being negligibly small in practical calculations. Furthermore, our verification cases presented in §§ 4.2 and 4.3 support this conclusion. In these test cases, the divergence-free errors caused by these two items were always present and not processed. If $Er_2$ and $Er_3$ contribute significantly to the divergence level of the results, the limited effect of the current divergence-free correction ignoring these two terms would be observed. However, the correction method results are very effective (nearly 100 %); this indicates that the errors caused by $Er_2$ and $Er_3$ are not significant.

As a result, the main contribution of the divergence-free error is obtained as

(2.8)\begin{equation} \left[ u^{(n)}_{i,i}\right]_{{main}} = 2 \sum_n\left[\underbrace{\cos\alpha^{(n)} q^{(n)} L_{ij,i} \sigma^{(n)}_j }_{Er_4} -\underbrace{\sin\alpha^{(n)} q^{(n)} \kappa_i^{(n)} L_{ij} \sigma^{(n)}_{j} }_{Er_5}\right]. \end{equation}

The increase in the divergence magnitude produced by the correlation reconstruction operation is mainly because of $Er_4$ and $Er_5$; therefore, the correction method to strictly ensure the divergence-free property of the resulting turbulence starts with these two error terms. The fourth term $Er_4$ represents the inhomogeneity of the correlation reconstruction matrix, which is caused by the spatial distribution of the Reynolds stress tensor field $R_{ij}$. When $R_{ij}$ is uniform in the computational domain, $Er_4 = 0$. The fifth term $Er_5$ is the influence of anisotropy of the correlation reconstruction, which is caused by the anisotropy of the local $R_{ij}$. When $R_{ij}$ is isotropic, the fifth term degenerates to the form $\sin \alpha ^{(n)} q^{(n)} \frac {1}{3} L_{kk} \kappa _i^{(n)} \delta _{ij} \sigma_j^{(n)} = 0$.

2.4. Weak inhomogeneity assumption

In (2.8), only $Er_5$ contains the wavenumber magnitude $\kappa ^{n}$. This indicates that a large wavevector $\kappa _i^{(n)}$ or small length scale tends to cause a larger anisotropy error term $Er_5$ of a specific Fourier mode, while the inhomogeneity error term $Er_4$ remains unaffected. When the wavenumber $\kappa ^{(n)} = \|\kappa _i^{(n)}\|$ is sufficiently large to exceed a critical value, the overall divergence-free error caused by the correlation reconstruction is expected to be dominated by $Er_5$. In this scenario, accounting only for $Er_5$ can also lead to a significantly low divergence level. This phenomenon is termed weak inhomogeneity, and the aforementioned modes are referred to as the scale satisfying weak inhomogeneity assumption.

It is observed that, when

(2.9)\begin{equation} \frac{ \|L_{ij,i}\|_{{2\text{-}norm}} }{ \kappa^{(n)} \| L_{ij} \|_{{2\text{-}norm}}}\ll 1, \end{equation}

$Er_4$ is negligible compared with $Er_5$. If (2.9) is satisfied for a certain wavenumber, the corresponding scale can be considered as fulfilling the weak inhomogeneity assumption. The physical meaning of (2.9) is that the spectrum-based method essentially superimposes Fourier modes/series to construct/reproduce turbulent fluctuations. However, the largest scale that can be identified using the discrete Fourier transform is the largest spatial scale corresponding to the computational domain. The zero-wavenumber mode is a constant without spatial distribution. Therefore, spatial inhomogeneity is a signal with wave numbers beyond the recognition range of the discrete Fourier series belonging to $[0,2{\rm \pi} /L_d]$. As the wavenumber increases, the microscale turbulence tends to be homogeneous, and the influence of this small-wavenumber signal becomes weaker. Therefore, large-scale turbulence is mainly affected by macroscopic inhomogeneity, and the divergence-free errors due to correlation reconstruction are also concentrated in large-scale modes close to the scale of the computational domain. Thus, there exists a critical mode $n_c$, where the effect of $Er_4$ is negligible for arbitrary Fourier modes fulfilling $\kappa ^{(n)} > \kappa ^{(n_c)}$, that is, the error term associated with the spatial inhomogeneity generated by these modes can be considered to be zero.

However, for realistic turbulence, its integral scale/energy-containing scale is typically in a large-scale range as well. Therefore, it is not sufficient for evaluating the divergence-free error of flows to only consider whether the scale satisfies the weak inhomogeneity assumption. The relative size of the energy-containing scale of the specific turbulent flow should be taken into account as well. This is referred to as determining whether the turbulence satisfies the weak inhomogeneity assumption. Whether $Er_4$ is negligible for flows depends on the relative size of the critical mode and integral scale of the turbulence. If it satisfies

(2.10)\begin{equation} \frac{\kappa_c}{\kappa_l} =\frac{ \|L_{ij,i}\|_{{2\text{-}norm}} }{ 2{\rm \pi} \sqrt{ \rho \left(R_{ij}\right) } } l \ll 1, \end{equation}

the reconstruction divergence-free error owing to inhomogeneity is negligible. In (2.10), $l$ is the integral scale of the local turbulence, and $\rho (R_{ij})$ represents the spectral radius of the local $R_{ij}$. As discussed in §§ 4.3 and 4.4, the condition of (2.10) is typically satisfied for practical applications.

As regards whether the practical applications satisfy the weak inhomogeneity assumption (i.e. whether (2.10) is satisfied), we conducted a theoretical calculation and evaluation. A dimensionless variable $\kappa _c/\kappa _l$, which indicates the influence of inhomogeneity, is calculated based on channel-flow data and (2.10). The distribution of this wavenumber ratio is indicated by the red line shown in figure 1. Evidently, the critical wavenumber of inhomogeneity $\kappa _c$ in the entire domain is less than 8 % of the energy-containing wavenumber $\kappa _l$, and this ratio remains below 4 % in most areas. This illustrates that, in the channel flow, $Er_4$ related to inhomogeneity is considerably less than $Er_5$. Therefore, the turbulent flow satisfies the weak inhomogeneity assumption.

Figure 1. The proportion of influence of inhomogeneity and the spatial distribution of the turbulence integral scale in channel flow. Results are based on DNS data of Kim, Moin & Moser (1987).

It should be emphasized that the weak inhomogeneity assumption does not assume that turbulence is homogeneous. Using the correlation reconstruction extension in § 2.2, the method can still obtain accurate, inhomogeneous and anisotropic turbulence spatial correlations. The concept of the weak inhomogeneity assumption is that the divergence-free error $Er_4$ caused by the inhomogeneity is sufficiently small enough; hence, only the anisotropy error of $Er_5$ requires correction to obtain the inhomogeneous and anisotropic turbulent field with considerably low divergence level in practice.

2.5. Divergence-free correction: inverter version

The inverter version divergence-free correction method is based on the weak inhomogeneity assumption. It accounts only for the anisotropy error term $Er_5$; this implies that

(2.11)\begin{equation} \kappa_i^{(n)} L_{ij} \sigma^{(n)}_{j} = 0 \end{equation}

should be satisfied for all modes.

However, it is not feasible to use the conventional operation directly because it would result in the deviation of the turbulence correlation. This is common in previous studies that sought to enforce a divergence-free property for various primitive methods. For example, the use of the projection step as mentioned in § 1 will result in this type of deviation from desired correlations (Yu & Bai 2014). Furthermore, this deviation often depends on specific algorithms, in particular random unit vector generation and manipulation. Therefore, deviations resulting from methods discussed in previous studies are often unpredictable and extremely severe in certain instances; for example, a minimum principal stress may be obtained in the direction of the expected maximum principal stress. Therefore, the conventional operation in previous research (Saad et al. 2016) should not be adopted. In contrast, if the weak inhomogeneity assumption is satisfied, this problem can be avoided using a special computational strategy. The possible correlation deviation caused by modifying the direction vectors $\sigma _i^{(n)}$ was observed. Thus, in the inverter version, we neglect the previously proposed conventional idea (Saad et al. 2016) to obtain $\sigma _i^{(n)}$ by first calculating the inhomogeneous temporary wavevector field according to the uniform wavevector. Instead, the uniform $\sigma _i^{(n)}$ of each mode is first obtained, and the temporary direction vector is calculated as follows:

(2.12)\begin{equation} \tilde{\sigma}_i^{(n)} = L_{ij} \sigma_j^{(n)}. \end{equation}

The final wavevector direction $\hat {\kappa }_i^{(n)}$, which is perpendicular to $\tilde {\sigma }_i^{(n)}$, is then obtained. There is no trigonometric function in (2.12); therefore, the vector product can be used once or twice (detailed in Appendix A) to calculate the perpendicular vectors. The inverter version of the divergence-free spectrum-based method completely avoids the correlation deviation, and it transfers all possible effects to the wavevector of each mode, which has exhibited a negligible effect on the energy spectrum during our practical computation test. This is also supported by our verification cases discussed in §§ 4.2 and 4.3, for which no deviation from the desired correlations is observed in the results of the inverter method.

For practical turbulent flow that fulfils the weak inhomogeneity assumption, the inverter version correction method exhibits good performance not only in the low-divergence property but also in the accuracy of the desired correlation tensors. For homogeneous anisotropic turbulence in particular, this version of the generation method strictly corrects the divergence-free error and can be termed a divergence-free method, whereas for inhomogeneous anisotropy, the inverter version can still correct $Er_5$ for anisotropy. The inverter version naturally returns to the primitive extended method only in the synthesis of inhomogeneous and isotropic turbulence. Even under this circumstance, the turbulence correlation tensor produced by the inverter method is accurate and does not show the aforementioned deviation owing to the employed strategy.

Notably, the inverter correction method handles the component (polarization) anisotropy, and there is no additional processing regarding the different length scales in different directions (directional anisotropy) as in certain prior studies (Patruno & Ricci 2018; Bervida et al. 2020). Therefore, the generated turbulence relies on subsequent calculations to restore this directional anisotropy. However, compared with a total computation requirement of $M \times N$ modes presented by Patruno & Ricci (2018), the inverter method uses a concise operation of (2.12) and maintains a total number of $N$ modes in the current framework, which has advantages in terms of computational complexity. Moreover, as mentioned by Bervida et al. (2020), the divergence-free error resulting from strong inhomogeneity may increase through the superposition of $M$ modes. In contrast, although the divergence-free error of the inverter method may increase in processing strongly inhomogeneous flow, the method can ensure that the generated turbulence spatial correlations are accurate without deteriorating.

2.6. Divergence-free correction: shifter version

If the inhomogeneity of the turbulence is significantly strong, the assumption of weak inhomogeneity of the flow is no longer valid. As presented by Bervida et al. (2020), the divergence-free error caused by inhomogeneity will contribute significantly to the results under this circumstance. At this time, although the inverter method can still provide precise inhomogeneous turbulence statistics, the resulting divergence level will increase significantly. Therefore, we propose the shifter version correction method to correct the two terms in (2.8) for both inhomogeneous and anisotropic turbulence. However, the correlation deviation can no longer be simply eliminated as in the inverter method. Although this deviation problem cannot be solved perfectly, we noted that the actual behaviour of deviation depends on the random unit vector operations. Therefore, using special algorithms that generate and manipulate unit direction vectors, we successfully transformed the correlation deviation which is a mathematical or statistical problem into two more physical and predictable problems, namely the anisotropic degeneration and TKE reduction. Moreover, the latter can be effectively compensated for the decrease by adding a scaling operation to the algorithm (a more detailed explanation of the final shifter version correction method is discussed in Appendix A). In general, although the shifter method still has the problem of anisotropic degeneration when dealing with inhomogeneous and anisotropic flows, such errors can be estimated, and are more consistent with the physical laws. Hence, the generated turbulence is less likely to undergo abnormal decay in subsequent calculations. In addition, the results generated by the shifter method are strictly divergence-free regardless of whether the correlation matrix is deviated or not.

2.7. Energy spectrum models

A distinct advantage of the spectrum-based method is that it can produce turbulence based on arbitrary forms of energy spectra. In addition to some specific flows that can be acquired by interpolating experimental data (Comte-Bellot & Corrsin 1971), many existing studies have employed various energy-spectrum models to ensure the use of complex flows (von Kármán 1948; Pao 1965; Kraichnan 1970; Bailly & Juve 1999; Yu & Bai 2014). For the spectrum-based method, the adopted model determines whether the method can be applied for high-Reynolds-number turbulence. To ensure that the method applies to more extensive flows, a high-Reynolds-number model with high accuracy is necessary. Except for the low-Reynolds-number energy spectrum model used by Kraichnan (1970), almost all high-Reynolds-number models can be expressed in their general form, as follows:

(2.13)\begin{equation} E_k(\kappa) = C \varepsilon^{2/3} \kappa^{-5/3} f_l(\kappa l) f_\eta(\kappa \eta), \end{equation}

where the functions $f_l$ and $f_\eta$ describe the distribution form of the spectrum in the energy-containing and dissipation ranges, respectively. We compared various energy-spectrum models, including the low-Reynolds-number model, with the CBC spectrum data (Comte-Bellot & Corrsin 1971), and the results are illustrated in figure 2.

Figure 2. Comparison of different energy spectrum models with experimental data.

Figure 2 shows that the high-Reynolds-number model that adopts the energy-containing range model of $p_0 = 4$ and the dissipation range model of $\beta = 5.2$ exhibits the best performance. Therefore, this combination of the energy-spectrum model is implemented in the code. Detailed information concerning particular parameters and the calculation of model coefficients can be found in Pope's book (Pope 2001). The final spectrum model that accounts for high-Reynolds-number-flow is defined as follows:

(2.14)\begin{equation} E_k(\kappa) =1.5 \varepsilon^{2/3} \kappa^{-5/3} \left[\frac{\kappa l}{\sqrt{{(\kappa l)}^2 + c_L}}\right]^{17/3} \exp\left({-5.2\left[\sqrt[4]{(\kappa \eta)^4 + c_\eta^4} - c_\eta \right]}\right), \end{equation}

where $c_L = 1.10075$ and $c_\eta = 0.401685$.

2.8. Remarks on the inverter and shifter methods

This section presents a summary of the features and application range of the proposed inverter version and shifter version divergence-free correction method. The inverter method can ensure accurate spatial correlations under general conditions regarding inhomogeneous and anisotropic turbulence. However, the weak inhomogeneity assumption must be satisfied to keep the divergence level of generated turbulence sufficiently low. If the weak inhomogeneity assumption is not valid, the divergence-free error will gradually increase with the enhancement of inhomogeneity. In several applications of complex flow, strict verification of weak inhomogeneity is difficult to implement. Hence, we recommend using the inverter version correction method first. If the generated turbulence has an unusually high divergence level or an abnormal decay is observed in the follow-up calculation, the shifter version should be adopted to regenerate the turbulence. The shifter method can obtain divergence-free results for any inhomogeneous and anisotropic turbulence. Furthermore, the method can transform the unpredictable correlation distortion into predictable anisotropic degradation. Owing to the advantages of the spectrum-based method in spatial correlation, such degradation generally does not cause abnormal decay in subsequent calculations and the correct anisotropy can be recovered; however, its possible drawback in practical use must be noted.

An advantage of the proposed method (both the inverter and shifter versions) is that the restrictions on the specific flow types applicable are minimal. A variety of internal and external flows can be applied, and the requirements for input turbulence information are flexible. The quality of the results generated using this method is not related to the specific flow type, but to the detail of turbulence information input to the method in practice. If the energy spectrum data, Reynolds stress tensor distribution, and dissipation rate distribution corresponding to a unique flow can be provided, the turbulence generated by the proposed method would be close to the real condition. However, even if the turbulence input parameters are not sufficiently detailed (for example, several external flow problems can only provide a turbulence intensity of the upstream), the proposed method can still adopt the default energy spectrum model (2.14). In addition, it can estimate the missing parameters (e.g. integral scale and correlation tensor) according to the default settings and then obtain a result that satisfies the conventional energy spectrum distribution while maintaining a strong spatial correlation. The generated turbulence is less likely to decay owing to the spatial correlation associated with the spectrum-based method. If there is a deviation from the real turbulence owing to a lack of input parameters, it will be rapidly recovered in the subsequent calculation. To ensure a strong correlation and prevent decay, some other more special turbulence characteristics (such as directional anisotropy) are not considered in the method at present. Therefore, these characteristics can only be gradually recovered using subsequent calculations; this is a major limitation of the current method.

3. Algorithm

In this section, specific implementation algorithm of the inverter-version correction methods is described, and the corresponding computational complexity and runtime storage requirements are analysed. For detailed descriptions and analysis of the shifter version algorithm, please refer to § A.3.

3.1. Implementation recipe

The algorithm for the inverter version is as follows.

1. (i) The correlation reconstruction tensor field $L_ij$ is calculated from (2.6) based on the given Reynolds stress distribution.

2. (ii) The unit direction vector $\sigma _i^{(n)}$ and random unit vector $\xi _i^{(n)}$ are generated for each Fourier mode.

3. (iii) The intermediate direction vector $\tilde {\sigma }_i^{(n)}$ is calculated from (2.12).

4. (iv) The wavevector $\kappa _i^{(n)} = \kappa ^{(n)} \hat {\kappa }_i^{(n)}$ is calculated using a one-time or two-time cross product operation using $\xi _i^{(n)}$.

5. (v) The intermediate velocity field $v_i$ is calculated through (2.4).

6. (vi) The final velocity field $u_i$ is calculated from (2.5) using the correlation reconstruction tensor $L_{ij}$.

3.2. Computational complexity

Compared with earlier spectrum-based methods that are applicable to the generation of homogeneous and isotropic turbulence, the additional computational complexity according to the above steps mainly consists of the Cholesky decomposition of the Reynolds stress field and the construction of intermediate wavevectors. The computational cost for the construction of the intermediate direction vector $\tilde {\sigma }_i^{(n)}$ using (2.12) is rather small. Cholesky decomposition needs to operate on the entire domain; however, it only needs to be performed once, regardless of mode, and it does not involve any algorithm for computing eigenvalues or eigenvectors. Compared with the algorithm that requires coordinate transformation and scaling operations, the proposed method is significantly more efficient. The computational procedure above is concise and easy to implement, and it requires only a slight modification of the spectrum-based method code that generates homogeneous and isotropic turbulence to realize this new improvement.

3.3. Runtime storage requirements

The additional runtime memory usage in the inverter version of the proposed divergence-free method is only the storage of the correlation reconstruction tensor field $L_{ij}$, and it is almost negligible. Moreover, because $L_{ij}$ is a lower triangular matrix, only six components must be stored for each grid point in practice, rather than the nine components required of a typical orthogonal matrix. All the above steps use only local storage and local information; therefore, the proposed method is easy to parallelize naturally and is suitable for HPC for the scale-resolving turbulence simulation.

4. Results and discussion

Four types of verification calculation are performed to verify the effectiveness of the new method. The code implementation and corresponding CFD computations for each method are realized using OpenFOAM (ESI-OpenCFD 2021). The number of modes in all cases is 5000, and the interpolation method of logarithmically distributed modes is employed. Case 1 is a classical benchmark of homogeneous and isotropic turbulence decay, which is mainly used to illustrate the correct realization of the classical spectrum-based method and the effectiveness of the CFD codes and to verify that the new method accurately reverts to the original method in the application of homogeneous and isotropic turbulence. Another important purpose of Case 1 is to obtain the benchmark value of the divergence level at which the box-geometry-type case can be considered divergence-free in practical calculations under the same number of grids. Cases 2 and 3 used the box geometry of Case 1 in a similar manner. Case 2 focuses on the performance of the method for handling anisotropic turbulence. Therefore, homogeneous synthetic turbulence with different anisotropy types is examined in detail, and the performances in producing the desired statistics between the shifter version and inverter version are compared. The anisotropy degeneration and TKE decay issues produced by the shifter version are also studied in Case 2. Case 3 verifies the treatment of inhomogeneity. An inhomogeneous turbulence with distribution only in the $x$ direction and a wavelength of one computational domain size is assumed. The purpose of this artificial design is to maintain the $y$–$z$ plane for spatial averaging in a single realization/sampling as well as ensemble averaging for multiple realizations/samplings. In Case 3, which mainly adopts the inverter version of the divergence-free method, the magnitudes of the reconstructed divergence-free errors caused by inhomogeneity and anisotropy are compared. Case 4 considers a typical application of inhomogeneous and anisotropic turbulence in engineering practice: channel flow. First, we analyse whether inhomogeneous divergence-free errors of the actual turbulence can be ignored. Case 4 is primarily carried out using the inverter version as well. Turbulence generation using a uniform integral scale and a spatial-distributed integral scale are compared in order to investigate the effect of the inhomogeneity of the normalized energy spectrum. Finally, the proposed method is applied to generate the initial field of the channel flow to verify its practical performance in a scale-resolving turbulence simulation.

In Cases 1 and 4, LESs are carried out with second-order temporal and spatial discretizations. The subgrid model of the LES in Case 1 is the WALE model (Nicoud & Ducros 1999). In Case 4, the subgrid model is not employed because a mesh with a DNS resolution was used.

4.1. Homogeneous and isotropic turbulence

Case 1 covers a benchmark of the decay of homogeneous and isotropic turbulence, where the size of the computational domain is $2{\rm \pi} \times 9$ cm, and the number of grids is $256^3$. When generating homogeneous and isotropic turbulence, the shifter version, inverter version, and extended methods without divergence-free correction are equivalent to the original spectrum-based method, and there is no difference in the generated initial field. There LESs are performed based on OpenFOAM, a finite-volume method open-source CFD code. A second-order Euler scheme is used for time advancement, and the spatial discretization employs the central difference scheme. The energy spectra of the computational results at different times are compared with the experimental data, as shown in figure 3. The energy spectrum of the simulation is in good agreement with the CBC data, which indicates that both the method implementation and CFD code are valid.

Figure 3. Energy spectrum validation. (a) The smallest turbulence scale resolved by different cell sizes. Dash line represents high-$Re$ spectrum model shown in (2.14). (b) Energy spectrum of homogeneous and isotropic turbulence at different times during decay. Experimental data were obtained from the literature (Comte-Bellot & Corrsin 1971).

Although the divergence is zero in theory, the numerical divergence of the practical computation after discretization has a non-zero value, owing to the influence of a series of factors, such as the floating-point number accuracy, numerical error, flow type and divergence-calculation method. Previous studies have shown that when the number of grids is sufficient, high-order interpolation methods exhibit limited improvement in the accuracy of divergence calculations (Yu & Bai 2014). Therefore, the second-order central difference method is used to calculate the divergence in all the test cases. Undoubtedly, the solenoidal property of the original spectral method that addressed inhomogeneity and isotropy is recognized. Therefore, we consider the realistic divergence value of the proposed method, which produces homogeneous and isotropic turbulence as the reference or benchmark, thereby indicating that it is divergence-free. This treatment eliminates other irrelevant factors in the subsequent discussion of the divergence-free error of the same type of box geometry. Figure 4 shows the actual divergence magnitude level $|u_{i,i}|$ for homogeneous and isotropic turbulence at different grid resolutions. The blue line in the figure represents the spatial average of the divergence magnitude $E$ (sample expectation), and the grey area represents the range of $[\max (0,E-D),D]$, where $D$ is the sample standard deviation of the divergence magnitude distribution that is obtained based on spatial averaging.

Figure 4. Reference value indicating a divergence-free property under different mesh resolutions (obtained by homogeneous and anisotropic turbulence-generation computations).

4.2. Anisotropic turbulence

Case 2 adopts the same box geometry as Case 1, and four types of anisotropic turbulence are tested while the TKE remains the same. Detailed component information is presented in table 1. In the table, $R_0 = k/6 = \langle u_i u_i \rangle / 12$. The off-diagonal elements of the Reynolds stress of types A, B and C are all zero, and the three diagonal elements are the three principal stresses. The type-D tensor contains non-zero off-diagonal components. Type A represents the anisotropy caused by one principal stress being significantly large; type B represents the anisotropy caused by one principal stress being equal to the average stress, and the other two principal stresses being larger or smaller; and type C represents the anisotropy caused by one principal stress being smaller than the other two. Because type D contains non-zero off-diagonal elements, the coordinate axis is not the same as the stress principal axis, which is a more common situation. The magnitude relationship between each component is designed to satisfy the stress-magnitude relationship in a typical channel flow (i.e. Case 4): $\langle u_1 u_1 \rangle$ is the maximum, $\langle u_3 u_3 \rangle$ is slightly greater than $\langle u_2 u_2 \rangle$, $\langle u_1 u_2 \rangle$ is negative, and the remaining components are all zero.

Table 1. Four types of tested anisotropy.

The shifter and inverter versions were used to generate these four types of anisotropic turbulence, and the corresponding statistical results based on spatial averaging are shown in table 2. The data in the table exhibit the turbulence-field results obtained by only one realization (one sample), where the seeds of the random-number engine in the two versions are guaranteed to be the same. It is evident from the table that the inverter-version results are in good agreement with the four given anisotropies, which ensures the accurate reproduction of the desired second-order correlation function. However, the Reynolds stress components depicted by the shifter version exhibit an evident deviation under all four anisotropies, which leads to tensor deformation. As mentioned above, the reason for this result is that the unit direction vector cannot maintain the original spherical distribution; thus, the covariance matrix deviates from the unit matrix.

Table 2. The accuracy of desired anisotropic Reynolds stress generated by different versions of the proposed method.

Figure 5 shows the increase in the divergence-free error if the divergence correction is not applied (i.e. the extended method in § 2.2). The error-growth factor in the figure is defined as follows:

(4.1)\begin{equation} R_X = \frac{X_{{uncorrected}}-X_{{standard}}}{X_{{standard}}}, \end{equation}

where $X$ is the statistical indicator used; in figure 5(a), this is the mean value $E$, and in figure 5(b), this is the deviation $D$. Here $X_{{standard}}$ represents the zero-divergence benchmark under the same mesh resolution, which is the result of the calculation in § 4.1 (figure 4). Four grid numbers ($64^3$, $128^3$, $192^3$, and $256^3$) are adopted. It is clear from the figure that the more anisotropic the turbulence, the larger the divergence-free error, particularly for types A and D. The behaviour of the anisotropy-intensity relationship is consistent with the relative magnitudes of $A_d$ and $A_n$ calculated using (4.3) and (4.4), respectively. In addition, the relative growth rate of the error $R_X$ also exhibited a mild dependence on the mesh resolution, that is, the larger the number of grids, the higher the relative growth rate of the error. The behaviour of the mean value $E$ in the figure is the same as that of the deviation value $D$, which means that the observation is independent of the chosen indicator.

Figure 5. Variation of the divergence-free error due to different anisotropy types with respect to increasing mesh resolution: (a) mean value/sample expectation; (b) deviation value/sample standard deviation.

Figure 6 illustrates the effectiveness of the inverter-version spectrum-based method for correcting divergence growth. The correction effectiveness in the figure is calculated using the following formula:

(4.2)\begin{equation} R_X = \frac{X_{{corrected}} - X_{{standard}}}{X_{{uncorrected}} - X_{{standard}}}, \end{equation}

where $X_{{corrected}}$ represents the calculation result of the divergence-free method. When $\eta _X$ is close to zero, the method does not correct the error. As $\eta _X$ approaches 1, the method corrects the 100 % increase in divergence caused by anisotropy. Owing to the existence of factors such as numerical errors, the correction effectiveness calculated in this manner may exceed unity. Figures 6(a) and 6(b) show the effectiveness represented by the mean and deviation values, respectively. The figure shows that the inverter version corrects the four types of anisotropy close to 100 %, and the minimum value in figure 6(a) also reaches unity. Owing to the characteristics of the standard deviation itself, the behaviour of the correction performance in figure 6(b) is different from that in figure 6(a). The minimum value of the effectiveness is also above 90 % of the total effectiveness, which fully illustrates the qualification of the inverter version in correcting the divergence-free error caused by anisotropy. Moreover, although a mesh dependence is observed in the uncorrected errors shown in figure 5, the behaviour in figure 6 is independent of grid size, which indicates that the inverter method is grid-independent for the correction of a divergence-free error.

Figure 6. Effectiveness of the inverter version on the correction of a divergence-free error for different anisotropy types: (a) mean value/sample expectation, (b) deviation value/sample standard deviation.

4.2.1. Correlation deviation issue of the shifter method

Further investigation of the results in table 2 shows that the correlation deformation of the shifter version integrating the special technique exhibits a common pattern, which we call anisotropy degeneration: any anisotropic tensor tends to degenerate towards an isotropic tensor $\delta _{ij}$. In the deformed tensors of types A, B and C, this pattern is expressed as the three principal stresses being more uniform: the larger principal stress decreases, and the smaller principal stress increases. In the results for type D, this isotropy trend indicates that the off-diagonal elements are closer to zero. This anisotropic degradation pattern is acceptable because it does not introduce errors or distortions that violate physics. The anisotropy intensity of the resultant turbulent field is weakened; however, it can be quickly restored to the real anisotropic state after development (calculation). Moreover, it is important to note that the anisotropic degradation presented in the shifter version described above does not necessarily apply to other spectrum-based methods that adopt similar ideas but different vector normalization algorithms (Kraichnan 1970; Smirnov et al. 2001). This is because the change in the statistical distribution owing to normalization at different steps is typically different. For example, when three directions of a vector follow three independent Gaussian distributions, the covariance matrix is the identity matrix. However, after normalizing it to a unit vector, each component has a normal distribution of range $[-1,1]$ and a 1/3 variance, and the three components are no longer independent. Although three times the covariance matrix of the latter is still the identity matrix, the two random vectors are already different from a statistical perspective. After an identical tensor operation, for example, tensor multiplication, the distributions of these two vectors will be significantly different (even if the covariance matrix is no longer the same). We attempted to use only the Gaussian distribution to generate vector components and did not perform unit-vector normalization in the intermediate process. We only performed an operation similar to dividing by three to ensure that the covariance matrix of the vector remained unchanged. As a result, the Reynolds-stress distortion obtained by this treatment was even worse than that of the shifter version, and there was no similar pattern of anisotropy degeneration. Unreasonable oversizing and undersizing were observed. However, this treatment did not perform normalization operations; therefore, the analytical formula of the statistics could be obtained through theoretical analysis. Among them, we obtained two physical quantities that characterize the degree of anisotropy:

(4.3)$$\begin{gather} A_d =\frac{ \left(R_{11}-R_{22}\right)^2 + \left(R_{11}-R_{33}\right)^2 + \left(R_{22}-R_{33}\right)^2 }{ \left(R_{11} + R_{22} + R_{33}\right)^2 }, \end{gather}$$

(4.4)$$\begin{gather}A_n = \frac{ R_{12}^2 + R_{13}^2 + R_{23}^2 }{ \left(R_{11} + R_{22} + R_{33}\right)^2}. \end{gather}$$

In these equations, $A_d$ denotes the anisotropy of the diagonal elements of the matrix, and $A_n$ denotes that of the off-diagonal elements.

In addition to the anisotropic degradation, the shifter version leads to a decay in the TKE, which is also not present in the inverter version. This TKE reduction is also caused by the change in the statistical distribution of the unit direction vector $\sigma _i^{(n)}$. As mentioned above, theoretically deriving the exact distribution form and covariance matrix can be difficult and tedious because the intermediate step of the shifter version correction method introduces vector-normalization operations. Therefore, we used numerical experiments to test the behaviour of this TKE decay with respect to anisotropy intensity $A_d$ and $A_n$, and the relevant results are shown in figure 7. Each data point in figure 7 is calculated by averaging 20 000 samples. Figures 7(a) and 7(b) show how the reduction of TKE varies with the diagonal anisotropy intensity $A_d$ and off-diagonal anisotropy intensity $A_n$, respectively. The off-diagonal components of the input stress tensor in figure 7(a) are zero, and each diagonal component varies in range from 1 to 100 by two orders of magnitude. The figure clearly shows that at this time, different types of anisotropy do not completely fall into a single curve, but present a clear distribution area. The range of this distribution area is large, at $A_d = 0.5\unicode{x2013} 1.5$, and it tends to an overlapping curve when the anisotropy is large and small. The blue dot set in figure 7(a) shows the influence of type A anisotropy, in which a particular principal stress is significantly larger, on the TKE decay. Additionally, in this case, the diagonal anisotropy intensity is close to its upper limit of two. It can be observed that this type basically determines the upper limit of the TKE decay, which coincides with the analytical curve of $k_{{cal}} = k_{{ref}} \sqrt {1 - 0.5 A_d}$. The anisotropy of the principal stress is observed to have a strong influence on TKE reduction. When $A_d$ tends to two, the TKE will be close to zero. In figure 7(b), the diagonal components are all maintained at unity, and the off-diagonal components are varied from 0 to 1/3 to ensure that the matrix is always positive definite. The influence of off-diagonal anisotropy has a law that is similar to that of diagonal anisotropy, i.e. it does not coincide with a single curve, but its influence is significantly smaller than that of the diagonal anisotropy. Even when the three off-diagonal elements reach the maximum value, the calculated TKE reduction is less than 20 %. The red dot in the figure shows the case in which the three off-diagonal elements are equal. The upper limit of the decay is determined again, and it is satisfied by the analytical curve $k_{{cal}} = k_{{ref}} ( 1 - \sqrt {3} A_n )$.

Figure 7. The reduction in TKE owing to the mapped direction vector: (a) diagonal element anisotropy, and (b) off-diagonal element anisotropy.

Therefore, the shifter version exhibits two undesirable characteristics: anisotropic degradation and TKE reduction. The latter can be scaled to recover according to $A_d$ and $A_n$ using the analytical function, whereas the former has no effective correction strategy if the framework of the shifter version is retained. However, the degenerated anisotropy can be quickly recovered to the real state in the actual calculation, as revealed by the applications using isotropic turbulence as the initial condition for anisotropic computations. However, the inverter version does not encounter the above two issues and performs well in reproducing the desired correlations and statistics.

4.3. Inhomogeneous turbulence

Case 3 also uses box geometry; however, to study the performance of processing inhomogeneous turbulence, an artificial distribution of correlations along the $x$ direction is set as follows:

(4.5)\begin{equation} \langle u_i u_j \rangle =\langle u_i u_j \rangle_{{mean}} \left[ 1 - A_0 \cos \left( 2{\rm \pi}\frac{x}{x_d} \right) \right], \end{equation}

where $x_d$ is the total length of the box domain in the $x$ direction. The $y$ and $z$ directions are still set to maintain homogeneity. Thus, in addition to ensemble averaging for multiple realizations, the planar averaging operation can still be performed to study the performance of the new method in a single realization. According to the Fourier series in (2.4), the minimum wavenumber of the generated velocity $u_i$ is $2{\rm \pi} /x_d$. Therefore, the minimum wave number of $u_i u_j$ is ${\rm \pi} /x_d$, that is, the maximum wavelength identified is $0.5x_d$. The Fourier mode of wavelength $x_d$ in (4.5) is clearly not in the recognizable range. Therefore, as mentioned above, it is a good choice for characterizing macroscopic inhomogeneity. By controlling the parameter $A_0$ ($0 \leq A_0 \leq 1$), distributed turbulence with different inhomogeneities can be obtained. The calculation results in this section are all acquired using the inverter version.

Isotropic but inhomogeneous turbulence and type A anisotropic turbulence satisfying (4.5) were generated under the condition that the spatial-averaged TKE is the same. The spatial distribution of the stresses after 100 consecutive generations and ensemble averaging are presented in figure 8. The $z$ cross-section in the figure shows the principal stress in the $x$ direction, and the two $y$ cross-sections show the distribution of the other two principal stresses. The figure shows that both isotropic and anisotropic turbulence accurately restore the inhomogeneous distribution of the cosine modes. Compared with the situation in which the three principal stress distributions in figure 8(a) are approximately the same, a significantly larger value of components $\langle u_1 u_1 \rangle$ is observed in figure 8(b), which indicates that the method also provides good accuracy in reproducing inhomogeneous anisotropy.

Figure 8. Spatial distribution of statistical results of Reynolds stress components (100 realizations of turbulence generation): (a) isotropic turbulence; (b) anisotropic turbulence (type A).

Distributions of the TKE and Reynolds stress components in figure 9 are obtained by a single realization and cross-plane averaging to further verify the reduction effect of the proposed method on the inhomogeneity. Figure 9(a) shows the TKE distribution for different inhomogeneity intensities reflected by the $A_0$ magnitude. It is evident that the new method has good accuracy for all distributions of $k$, even for a single realization. Figure 9(b) shows the verification results of the anisotropy of type D. The reproduction characterizations of all distributions are in good agreement with either the three normal stresses or the one non-zero shear stress. The distributions of $\langle u_1 u_3 \rangle$ and $\langle u_2 u_3 \rangle$ remain at zero with mild fluctuations. The larger fluctuations in the components $\langle u_1 u_3 \rangle$ are due to larger velocities in the two associated directions.

Figure 9. The reproduction of the target inhomogeneous distribution of turbulence correlations using the divergence-free method (single realization, $y$–$z$ plane averaging): (a) TKE $k$ for different inhomogeneities, and (b) Reynolds stress component of the type D anisotropy for $A_0 = 0.8$.

The special-averaged divergence levels of isotropy, type-A anisotropy and type-D anisotropy at $A_0 = 0$ (homogeneity), 0.2, 0.5 and 0.8, were computed to observe the absolute divergence increase caused by anisotropy and inhomogeneity in both anisotropic and inhomogeneous turbulence. The results are presented in table 3. The turbulence intensity varies at different spatial locations owing to the inhomogeneity. Therefore, to eliminate the influence of this factor on statistics, the averaged results in the table are dimensionless normalized divergence $|u_{i,i}|/u_t$.

Table 3. Divergence-free errors (uncorrected) of synthesized turbulence with different inhomogeneity intensities.

The calculated results of the mean and deviation of the divergence level in the table show that the increase in inhomogeneity intensity has almost no effect on the increase in divergence level. When $A_0$ increases from 0 to 0.8, the divergence level of both isotropic and anisotropic turbulence remains nearly the same. Taking isotropic turbulence as an example, its mean value remains at approximately 40, whereas the standard deviation remains almost unchanged at approximately 51 for all inhomogeneities. However, when the flow changes from isotropic to anisotropic, the divergence-free error increased significantly. For the type A anisotropy, the mean value increased from 44 to 106, and the deviation increased from 51 to 85. For the type D anisotropic turbulence, the mean value increased from 44 to 108, and the deviation increased from 51 to 87. The results in table 3 show that under an inhomogeneous distribution of this form (cosine distribution of wavelength $x_d$), the increase in the error caused by anisotropy exceeds that caused by inhomogeneity. When $A_0 = 0.8$, which is where the most intensive inhomogeneity of this wavenumber is approached, the increase in divergence-free errors brought about by inhomogeneity is still not clearly observed. Clearly, increasing the wavenumber (frequency) of the mode can indeed induce a larger local gradient and inhomogeneity, but this part of the mode can already be identified and constructed using spectrum-based methods. It is doubtful whether these high-frequency signals can be regarded as macroscopic inhomogeneity rather than as small-scale turbulence.

Table 4 shows the inverter version's correction performance for the two types of anisotropy in table 3 (isotropic data directly follows the results of table 3). The table shows that both the mean and deviation values of the two types of anisotropic turbulence generated by the correction method reach the actual divergence level of isotropic turbulence under the same conditions; therefore, the divergence-free requirements can be considered satisfied in practice. The inverter-version correction does not account for inhomogeneity errors; however, the correlation reconstruction matrix varies spatially, which leads to an effective correction for inhomogeneous anisotropy as well. This is reflected in table 4, which shows that the results of different intensities of $A_0$ have similar correction effectiveness.

Table 4. Divergence-free errors (after the inverter-version correction) of synthesized turbulence with different inhomogeneity intensities.

4.4. Typical inhomogeneous and anisotropic turbulent flow

The fourth case is a typical inhomogeneous and anisotropic turbulent flow that is common in practical engineering: a channel flow with fully developed boundary layers. The first three cases have verified the method in terms of the correlation accuracy and divergence-free property; however, Case 4 is mainly used to test the performance of the proposed method (inverter version) in practical applications used for practical turbulent flow. Additionally, the influence of inhomogeneity on divergence errors is checked analytically for the channel flow. The periodic channel-flow parameters and grids selected in this case are consistent with those in the DNS literature (Kim et al. 1987), and further information is provided by Mansour, Kim & Moin (1988). In this case, the LES does not use any subgrid stress model owing to the DNS-resolution mesh. (Complete LES calculations have been performed, and the validity of the CFD code is not shown here owing to content constraints and limited relevance.) The reference data used for comparison and the data used for the theoretical analysis in this section are the DNS calculation results in the literature.

In channel flows, the turbulent integral length scale (energy-containing scale) is typically not constant; it has a specific spatial distribution, as shown by the lilac line in figure 1. The integral scale of the region near the core is larger than that of the region near the wall. However, the non-uniform integral scale directly leads to different energy spectra at different spatial locations (multiple items in (2.13) are directly or indirectly dependent on $l$ or the turbulent dissipative rate). As shown in (2.7), this will cause the divergence-free error of $Er_1$ to occur simultaneously. To analyse the influence of the distributed integral scale, Case 4 generates inhomogeneous and anisotropic turbulence of the channel flow from the uniform integral scale of a single constant (the mean value represented by the grey dashed lines in figure 1) and the distributed integral scale profile, respectively. The differences in the statistical accuracy and divergence-free characteristics of the two are compared.

Figure 10 shows the probability density distributions of the divergence-free error of turbulence generated by the distributed and constant integral scale in a single realization. The solid red line represents an approximate fitting line. First, it should be noted that no common distribution can accurately describe the probability distribution function of the two, and the relative frequencies of both distributions approach zero rapidly, which reflects the effectiveness of the divergence-free correction. The black dot is the relative cumulative frequency of the divergence level of the probability distribution function, and the red dotted line is the fitting curve in an exponential form. The figure also shows the specific divergence level at cumulative frequencies of 20 %, 50 % and 90 %. The horizontal axis in the figure represents the dimensionless result normalized by $u_\tau / \nu$. The figure shows that there are differences in the statistical distribution of the divergence under the two integral scale choices. The results of the non-uniform integral scale in figure 10(a) have a larger frequency near zero, and the frequency in the range of 0–0.5 is smaller than that of the constant integral scale in figure 10(b). However, both the mean and deviation were very close to each other. From the perspective of cumulative frequency, the divergence level of the varying scale at 20 % is slightly lower than that of the constant scale, the divergence level at 50 % is close to that of the constant scale, and the divergence level at 90 % is slightly higher than that of the constant scale. Although the figure shows a difference in the absolute divergence between the two cases, the effect is minimal, and it is acceptable to ignore $Er_1$ when the distributed integral scale is adopted.

Figure 10. Statistical distribution of divergence-free errors in the entire computation domain: (a) results of the spectrum adopting the non-uniform integral scale profile, and (b) results of the spectrum adopting a constant integral length scale.

The generated correlation between the two energy spectra using different integral scale models was averaged using 100 samples, as shown in figure 11. In the figure, the $z$ cross-section shows the shear stress contour, whereas the other three planes indicate the normal stress contour. The figure illustrates that the given correlation distribution is well produced by both scale models. Notably, it is difficult to distinguish between figures 11(a) and 11(b), which indicates that the difference between the two is negligible and that the constant integral scale model is sufficiently accurate. The results of figure 11(a) are further spatially averaged in figure 12 to obtain normal stress profiles in figure 12(a) and the shear stress profile in figure 12(b). The figures show that all profiles are in good agreement with the benchmark data.

Figure 11. Spatial distribution of ensemble-averaged Reynolds stress components in the channel flow (100 samples): (a) results of the spectrum adopting the non-uniform integral scale profile, and (b) results of the spectrum adopting a constant integral length scale.

Figure 12. The reproduction performance of the stress profile by the proposed method (100 samples, spatial averaging in $x$–$z$ plane): (a) normal stress; (b) shear stress.

Finally, the new method is applied to the LES calculation of channel flow. To illustrate the turbulence decay problem, a comparison is made with a classical white-noise generator. The calculation results for a coarser mesh are added to discuss the performance of the method under different mesh resolutions. Although figure 13(a) shows only the instantaneous velocity at a probe near the core area, the instantaneous velocity and surface average Reynolds stress at other locations are also monitored, and their behaviour is consistent. Figure 13(a) clearly shows that the turbulence initial field generated by the white-noise method has a very clear decay phenomenon in the CFD calculation, and this decay exhibits a strong dependence on the grid size. The initial fluctuation under the coarse grid is strong, the decay speed is slow, and the correct state can be regenerated by the grid itself. However, the initial fluctuation in the finer grid is weaker, and the decay process occurs considerably faster. Generating turbulence again after the initial turbulence becomes totally laminar is difficult (the figure only shows the calculation results for 300 s, but there is still no fluctuation when the actual calculation time reaches 5000 s). This is caused by the lack of spatial correlation and distortion of the energy-spectrum distribution in the generated white noise. Previous literature typically refers to this generator as the white-noise method because the energy of the signal is the same at each frequency; however, this conclusion only addresses the temporal characteristics of boundary turbulence. When analysing the energy spectrum, we observed that the energy spectrum of the white-noise method exhibits the behaviour presented in figure 13(b), owing to the spatial characteristics of the turbulence. The energy of the large-scale turbulence was completely missing. Moreover, the smaller the grid size, the smaller the energy-containing scale. Thus, for a grid with DNS resolution, the kinetic energy of the white-noise method mostly falls into the dissipation scale. This feature leads to an unphysically large dissipation of the turbulence. The proposed spectrum-based method in figure 13(a) is based on the inverter version; barely any turbulence-decay phenomenon is observed, which shows the effectiveness of the method as an initial-field generator. The actual performance of the shifter version is very similar to that of the inverter version in figure 13(a) (not shown here), which indicates that the shifter version is also efficient as an initial-field-generation method, even if anisotropy degeneration occurs.

Figure 13. The effectiveness of different types of synthetic turbulence methods when generating initial fields for LES computation. (a) Comparison of turbulence decay between the proposed spectrum-based method and white-noise generator in coarse and fine meshes, and (b) energy spectra of white-noise generator in different grid sizes with respect to the length scale.

5. Conclusions

This study proposed a highly efficient and low-divergence inhomogeneous and anisotropic turbulence-generation method that can use arbitrary energy spectra. This method was verified by thorough test cases, and the results show that this method can accurately restore the distribution of turbulence statistics and effectively correct divergence-free errors. The main conclusions are as follows.

(i) A high-Reynolds-number spectrum model is adopted in the proposed method. The verification showed that the non-uniform spectrum has almost no effect on the correlation function and divergence level of the generated turbulence. Therefore, it is acceptable to use either a distributed integral scale or constant scale in practical applications.

(ii) If the weak inhomogeneity assumption is satisfied, an inverter version of the proposed method is recommended. Its effectiveness in anisotropy error correction was well illustrated in several test cases. Although the correction for an inhomogeneous divergence-free error was not presented, this study demonstrated through several verification cases that an increase in the divergence level in common applications is primarily owing to anisotropy, and the error caused by inhomogeneity is typically negligible. Therefore, for typical inhomogeneous and anisotropic turbulent flows in engineering, the results obtained using the inverter version are both accurate and effectively solenoidal.

(iii) If the target turbulence exhibits strong inhomogeneity, the shifter version of the method is recommended. This method ensured a strict divergence-free feature in inhomogeneous and anisotropic turbulence; however, issues of anisotropy degeneration and kinetic energy reduction may occur. The latter can be effectively corrected using scaling functions, whereas the former may require a longer recovery time or length. However, using the proposed method of this version as an initial field generator is still effective because the degenerated anisotropy can rapidly recover to its real state without unphysical turbulence decay.

(iv) The proposed method uses the correlation reconstruction technique to extend the method to the application of inhomogeneous and anisotropic turbulence without an algorithm that involves coordinate transformation or eigenvalue computations. Compared with the spectrum-based methods proposed in previous studies, our method provides the advantages of high computational efficiency, low runtime memory cost and easy implementation. Moreover, the required information is stored locally. Therefore, this method can undergo multiprocess parallelization in HPC, which is suitable for the practical computing requirements of large-scale scale-resolving turbulence simulations, such as DNSs and LESs.