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 Reference Garnier, Adams and Sagaut2009). 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 Reference Tabor and Baba-Ahmadi2010; Dhamankar, Blaisdell & Lyrintzis Reference Dhamankar, Blaisdell and Lyrintzis2017).
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. Reference Dhamankar, Blaisdell and Lyrintzis2017). 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 Reference Lund, Wu and Squires1998; El-Askary, Schroeder & Meinke Reference El-Askary, Schroeder and Meinke2003). 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 Reference Yousif, Yu and Lim2022); 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 Reference Kraichnan1970; Smirnov, Shi & Celik Reference Smirnov, Shi and Celik2001), digital-filter-based method (Klein, Sadiki & Janicka Reference Klein, Sadiki and Janicka2003; Kempf, Klein & Janicka Reference Kempf, Klein and Janicka2005; Kim, Castro & Xie Reference Kim, Castro and Xie2013), synthetic-eddy method (Benhamadouche et al. Reference Benhamadouche, Jarrin, Addad and Laurence2006; Jarrin et al. Reference Jarrin, Benhamadouche, Laurence and Prosser2006; Mathey et al. Reference Mathey, Cokljat, Bertoglio and Sergent2006; Subbareddy et al. Reference Subbareddy, Peterson, Candler and Marusic2006), turbulent spot method (Kornev & Hassel Reference Kornev and Hassel2007; Kröger & Kornev Reference Kröger and Kornev2018) and diffusion-based method (Kempf et al. Reference Kempf, Klein and Janicka2005). 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 Reference Poletto, Craft and Revell2013). 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. Reference Dhamankar, Blaisdell and Lyrintzis2017). Many types of synthetic-turbulence methods, such as digital filters (Klein et al. Reference Klein, Sadiki and Janicka2003) and synthetic-eddy methods (Jarrin et al. Reference Jarrin, Benhamadouche, Laurence and Prosser2006), 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. Reference Kim, Castro and Xie2013; Poletto et al. Reference Poletto, Craft and Revell2013), these improvements are typically limited by the methodological framework. For example, the improved method of Poletto et al. (Reference Poletto, Craft and Revell2013) 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. Reference Dhamankar, Blaisdell and Lyrintzis2017). 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. Reference Kim, Castro and Xie2013). In the method proposed by Kempf et al. (Reference Kempf, Klein and Janicka2005), if the projection step (Lee, Lele & Moin Reference Lee, Lele and Moin1992) 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. Reference Mathey, Cokljat, Bertoglio and Sergent2006).
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. (Reference Smirnov, Shi and Celik2001) 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 Reference Yu and Bai2014). 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 Reference Yu and Bai2014). 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 (Reference Patruno and Ricci2018). Their improved spectrum-based method combines the ideas of Poletto et al. (Reference Poletto, Craft and Revell2013) and Kröger & Kornev (Reference Kröger and Kornev2018) 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. Reference Bervida, Patruno, Stanič and de Miranda2020).
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. (Reference Bechara, Bailly, Lafon and Candel1994), the velocity vector in Cartesian space for a single sampling can be expressed in a Fourier series as follows:
\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 Reference Comte-Bellot and Corrsin1971), 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,
\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
\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. Reference Smirnov, Shi and Celik2001; Yu & Bai Reference Yu and Bai2014). 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:
$$\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}$$
$$\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,
\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. Reference Shur, Spalart, Strelets and Travin2014; Patterson, Balin & Jansen Reference Patterson, Balin and Jansen2021).
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
\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
\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
\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
\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 (Reference Kim, Moin and Moser1987).
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
\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 Reference Yu and Bai2014). 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. Reference Saad, Cline, Stoll and Sutherland2016) 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. Reference Saad, Cline, Stoll and Sutherland2016) 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:
\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 Reference Patruno and Ricci2018; Bervida et al. Reference Bervida, Patruno, Stanič and de Miranda2020). 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 (Reference Patruno and Ricci2018), 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. (Reference Bervida, Patruno, Stanič and de Miranda2020), 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. (Reference Bervida, Patruno, Stanič and de Miranda2020), 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 Reference Comte-Bellot and Corrsin1971), many existing studies have employed various energy-spectrum models to ensure the use of complex flows (von Kármán Reference von Kármán1948; Pao Reference Pao1965; Kraichnan Reference Kraichnan1970; Bailly & Juve Reference Bailly and Juve1999; Yu & Bai Reference Yu and Bai2014). 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 (Reference Kraichnan1970), almost all high-Reynolds-number models can be expressed in their general form, as follows:
\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 Reference Comte-Bellot and Corrsin1971), 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 Reference Pope2001). The final spectrum model that accounts for high-Reynolds-number-flow is defined as follows:
\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.
(i) The correlation reconstruction tensor field
$L_ij$ is calculated from (2.6) based on the given Reynolds stress distribution.(ii) The unit direction vector
$\sigma _i^{(n)}$ and random unit vector $\xi _i^{(n)}$ are generated for each Fourier mode.(iii) The intermediate direction vector
$\tilde {\sigma }_i^{(n)}$ is calculated from (2.12).(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)}$.(v) The intermediate velocity field
$v_i$ is calculated through (2.4).(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 Reference Nicoud and Ducros1999). 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 Reference Comte-Bellot and Corrsin1971).
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 Reference Yu and Bai2014). 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:
\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:
\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 Reference Kraichnan1970; Smirnov et al. Reference Smirnov, Shi and Celik2001). 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:
$$\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}$$
$$\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:
\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. Reference Kim, Moin and Moser1987), and further information is provided by Mansour, Kim & Moin (Reference Mansour, Kim and Moin1988). 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.
Funding
This project was supported by the National Natural Science Foundation of China (51922060) and the National Science and Technology Major Project (2017-III-0003-0027).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Details of shifter version divergence-free correction method
The direction vector 
$\sigma _i^{(n)}$ exists in the same form in 
$Er_4$ and 
$Er_5$ such that it satisfies
\begin{equation} \left(\cos\alpha^{(n)} L_{ij,i} - \sin\alpha^{(n)} \kappa_i^{(n)} L_{ij}\right) \sigma_j^{(n)} =\tilde{\kappa}_j^{(n)} \sigma_j^{(n)} = 0, \end{equation}
that is, 
$\sigma _j^{(n)}$ and 
$\kappa _i^{(n)}$ are perpendicular to each other, and the divergence of the generated turbulence is strictly zero. Therefore, by improving the construction strategy of the direction vector, we obtain a shifter version of the divergence-free spectrum-based method suitable for the generation of inhomogeneous and anisotropic turbulence. First, the intermediate wavevector is calculated as follows:
\begin{equation} \tilde{\kappa}_j^{(n)} =\cos\alpha^{(n)} L_{ij,i} - \sin\alpha^{(n)} \kappa_i^{(n)} L_{ij}. \end{equation}A mode direction vector that is perpendicular to the intermediate wavevector is then constructed. The shifter version of this divergence-free method simultaneously corrects for the correlation reconstruction errors caused by inhomogeneity and anisotropy, and it guarantees the fully solenoidal characteristics of the generated turbulence. However, in practical applications, we noticed that some problems in the shifter version require further improvement.
A.1. Correlation deviation problem
Equation (A2) includes trigonometric functions, and the intermediate wave vector shifts in the plane based on the divergence of the reconstruction matrix 
$L_{ij,i}$ and vector 
$\kappa _i^{(n)} L_{ij}$ with different spatial positions; hence, this version is referred to as the ‘shifter version’. When the inhomogeneity of the turbulent field is significantly small (
$L_{ij,i}$ is close to or strictly equal to zero), the existence of the 
$\sin \alpha ^{(n)}$ function leads to the frequent flipping of the supposed constant vector 
$\sigma _i^{(n)}$, if the perpendicular vector is constructed using (2.2). This limitation leads to a considerable deviation in the synthetic turbulence from the given energy spectrum. Moreover, the method is unable to correctly recover to the original version presented in § 2.1 when considering homogeneous turbulence. The flipping occurs because the restriction introduced by the theory only provides a perpendicular line; however, the direction vector has two options. The final direction depends on the algorithm that numerically employs the vector product operation. However, the one-time vector product directly introduces the positive and negative effects of trigonometric functions. This results in a direction vector jumping between the two possibilities according to spatial positions. To overcome this limitation, we introduce two vector-product operations in the direction computation step, that is, we substitute (2.2) using the following procedure:
$$\begin{gather} \zeta^{(n)}_i =\left(\varepsilon_{ijk} \xi^{(n)}_j \tilde{\kappa}^{(n)}_k\right)_{normalize}, \end{gather}$$
$$\begin{gather}\sigma^{(n)}_i =\left(\varepsilon_{ijk} \tilde{\kappa}^{(n)}_j \zeta_k^{(n)} \right)_{normalize}. \end{gather}$$
This treatment consists of two vector products, and the projection of the random unit vector 
$\xi _i^{(n)}$ on a plane perpendicular to 
$\kappa _i^{n}$ is calculated. Moreover, this treatment is more reasonable in physics than the one-time vector product version because the quantity that contains the permutation symbol 
$\varepsilon _{ijk}$ is not a tensor in essence (it does not satisfy Galilean invariance in a mirror coordinate transformation).
Although the correction method of the shifter version strictly ensures the solenoidal condition in a large application, it still requires solving the correlation deviation problem as mentioned in § 2.6. Consequently, the covariance matrix of the intermediate velocity 
$v_i$ deviates from the identity matrix, leading to a departure from the desired Reynolds stress. However, by adopting a special algorithm of (A3) and (A4), we can transform the correlation deviation problem into two physical and predictable problems, namely the anisotropic degeneration and TKE reduction. The TKE reduction can be almost eliminated by scaling it according to the two types of anisotropy intensities, whereas the solution is not simple for anisotropy degeneration. The primary reason for this deviation is that the random wavevector originally follows a spherical distribution; however, it cannot maintain this characteristic after a series of specific transformations. A further discussion concerning anisotropic degeneration and TKE reduction issues is presented in detail in § 4.2.1.
In summary, the shifter version of the correction method strictly ensures divergence-free properties of inhomogeneous and anisotropic turbulence. Although the generated turbulence field will exhibit a deviation in the correlation function, we observed that the shifter version, which preserves a solid spatial correlation, can generate the turbulent boundary and initial field. However, it requires an increased recovery length (boundary condition) or recovery time (initial condition), which is still less than that of other methods that include decay in turbulence. Owing to the effectiveness of the spatial correlation, the initial deviation of the synthetic turbulence produced by the shifter version method can rapidly develop and recover to the correct state without decay.
A.2. Correction for TKE reduction
As shown in § 4.2.1 and figure 7, the ratio of reduced TKE over the desired one can be estimated as
$$\begin{gather} \frac{k_{{cal}}}{k_{{ref}}} = \sqrt{1 - \frac{A_d}{2}}, \end{gather}$$
$$\begin{gather}\frac{k_{{cal}}}{k_{{ref}}} = 1 - \sqrt{3}A_n. \end{gather}$$Therefore, a scaling operation can be employed after (2.5) to fix the reduction as
\begin{equation} [u_i]_{{scaled}}=\frac{ u_i\left(x_{(k)}, t\right) }{ \left( \sqrt{1 - \dfrac{A_d}{2}} \right) \left( 1 - \sqrt{3}A_n \right) + \epsilon }, \end{equation}
where 
$\epsilon$ is a small number to prevent the denominator from being zero.
A.3. Algorithm of the shifter-version correction method
The algorithm for the shifter version is as follows.
(i) The correlation reconstruction tensor field
$L_ij$ is calculated from (2.6) based on the given Reynolds stress distribution.(ii) Random unit vectors
$\hat {\kappa }_i^{(n)}$ and $\xi _i^{(n)}$ are generated for each Fourier mode.(iii) The temporary or intermediate wavevector
$\tilde {\kappa }_i^{(n)}$ is calculated from (A2).(iv) The direction vector
$\sigma _i^{(n)}$ is obtained based on $\xi _i^{(n)}$ via two vector products using (A3) and (A4).(v) The intermediate velocity field
$v_i$ is calculated using (2.4).(vi) The final velocity field
$u_i$ is calculated from (2.5) using the correlation reconstruction tensor $L_{ij}$.(vii) The scaling operation using (A7) is conducted to correct the TKE.
Evidently, the difference between the inverter version and shifter version is in steps 2–4 and step 7, particularly in step 3. Because the trigonometric functions in the latter are eventually used in (2.4), the construction of the intermediate wavevector in step 3 adds only a few multiplication operations. Similarly, in the shifter version of the correction method, the implementation of the correlation reconstruction for inhomogeneity and anisotropy extension does not involve any coordinate transformation or eigenvalue computation algorithm; this significantly improves the computation speed, particularly in HPC. Similar to the inverter version, the shifter version is concise and efficient, and the codes of other primitive spectrum-based methods require only minimal changes to realize this divergence-free extension.
Further, similar to the inverter version, the shifter version of the correction method only adds a small amount of memory usage to the original method, and these runtime storages comprise field data that can be arranged in distributed machines with corresponding meshes. The algorithm only requires local information; therefore, it is easy to parallelize in HPC, which is suitable for providing both the initial and boundary conditions for high-cost LESs or DNSs.
