First, the isotropy of small-scale fluctuations in the generated steady turbulent field is examined. Figure 2(a) illustrates how the isotropy of small-scale fluctuations varies with the bulk Reynolds number, using the following definition:where 〈 〉 denotes ensemble average for both temporal and spatial directions, and Isotropy in the small-scale eddies corresponds to I_s = 1. The figure confirms that, with isotropic forcing, this index stays fixed at unity for all investigated Reynolds numbers. This indicates that when the isotropic forcing field is applied, the small-scale fluctuations are isotropic, validating the isotropic forcing field used in this study. When the anisotropic forcing field is used, I_s has smaller values at low bulk Reynolds numbers. However, as the Reynolds number increases, I_s also increases, reaching and maintaining values around I_s = 1 for Re ≥ 100. This result is qualitatively consistent for the two different constant values C used in the anisotropic forcing field. These results suggest that for bulk Reynolds numbers around Re ≥ 100, local turbulence can be considered isotropic. Conversely, for smaller Reynolds numbers, it is confirmed that the anisotropy of the forcing field causes the small-scale turbulence to also exhibit anisotropy.

The forthcoming results detail the turbulent Reynolds number for the analyzed flow. The turbulent Reynolds number, based on the Taylor microscale, Re_λ, is defined as follows:where K represents the turbulent kinetic energy; λ_x, λ_y, and λ_z denote the Taylor microscale in the x, y, and z directions, respectively. As shown in Fig. 2(b), for steady turbulence generated by isotropic forcing, the turbulent Reynolds number can be approximated as a power function of the bulk Reynolds number as follows:where m is the exponent, which is determined from the figure to be m = 2/3. Figure 2 also includes results for turbulence generated by anisotropic forcing. As shown, the distribution of turbulent Reynolds numbers for C = 2 does not coincide with the distribution for isotropic forcing, while the distribution for C = 1 is observed to coincide with the distribution for isotropic forcing. This finding indicates that setting the anisotropic forcing coefficient C to C = 1 ensures that the turbulent Reynolds numbers generated by anisotropic and isotropic forcing fields are consistent. Therefore, in this study, subsequent analyses of anisotropic forcing will use the results obtained with C = 1.

Figure 3 presents the relationship between the bulk Reynolds number and the turbulent kinetic energy K. Here, the turbulent kinetic energy is a statistical quantity obtained as the spatially averaged value over the entire computational domain, as follows K = (1/2)(〈u²〉 + 〈v²〉 + 〈w²〉). In the case of isotropic forcing, as shown in the figure, the turbulent kinetic energy increases with the bulk Reynolds number. In this analysis, the magnitude of the forcing is predefined, and the spatial length scale of the large-scale flow is also determined by the forcing field. Meanwhile, the dissipation coefficient is observed to decrease with increasing Reynolds number in the low turbulence Reynolds number range. These factors are considered to contribute to the monotonic increase in the mean turbulent kinetic energy with the Reynolds number observed in this study. Figure 3 also includes the turbulent kinetic energy of turbulence maintained by an anisotropic forcing field with C = 1. As shown, unlike the case of isotropic forcing, the turbulent kinetic energy of turbulence under anisotropic forcing exhibits a minimum at intermediate Reynolds numbers. For Reynolds numbers higher than this minimum, the turbulent kinetic energy of turbulence driven by anisotropic forcing increases with increasing Reynolds number.

Figure 4 (a) illustrates the variation of large-scale isotropy with the bulk Reynolds number, defined as follows:where I_L = 1 corresponds to the case where the turbulence is isotropic. When an anisotropic forcing field is applied, the nature of the field ensures that 〈v²〉 = 〈u²〉. As shown in the figure, in the case of turbulence maintained by isotropic forcing, the large-scale isotropy remains at I_L =1 regardless of the bulk Reynolds number. This result validates the turbulence generated by the isotropic forcing field used in this study. The figure also shows the large-scale isotropy results for turbulence driven by an anisotropic forcing field. As shown, in this case the isotropy I_L increases with Reynolds number at lower bulk Reynolds numbers. For Re ≥ 100, the I_L becomes approximately constant, independent of further increases in Reynolds number. The minimum Reynolds number at which the I_L reaches constant isotropy corresponds to the minimum Re at which the small-scale turbulence becomes isotropic, as shown in the figure.

To investigate the degree of anisotropy in large-scale turbulence, the anisotropy tensor [1] is used. The anisotropy tensor b_ij is given as follows:where i = 1 – –3, j = 1 – –3, and the summation is taken over k. Using the parameters η and ξ, defined as follows:The characteristics of the turbulence anisotropy can be classified using the above parameters. Specifically: if η = ξ = 0, the turbulence is isotropic. When ξ = –1/6 and η = 1/6, the turbulence is two-component anisotropic. If ξ = 1/3 and η = 1/3, the turbulence is one-component anisotropic. When η = ξ, the turbulence has axisymmetric anisotropy with a large eigenvalue. When η = – ξ, the turbulence has axisymmetric anisotropy with a small eigenvalue. When , the turbulence has two-component anisotropy. These relationships are organized using the Lumley triangle [1] shown in Fig. 4(b). In the figure, the triangle is represented by dashed lines. The turbulence values generated by isotropic and anisotropic forcing fields are plotted on this triangle. As shown in the figure, in the case of turbulence generated by isotropic forcing, the plots consistently lie at η = ξ = 0, regardless of the bulk Reynolds number. This again confirms that the turbulence generated by isotropic forcing is isotropic. Figure 4(b) similarly presents the Reynolds-number trend for turbulence produced under anisotropic forcing. As shown, the plots lie along the line η = –ξ regardless of the Reynolds number. This indicates that the anisotropic turbulence has an axisymmetric anisotropy with a small eigenvalue. As shown in the figure, for Re = 10, the plot is closest to the point ξ = –1/6 and η = 1/6. As the Reynolds number increases, the plots move further away from this point and approach a convergence point. This behavior corresponds to the earlier results for large-scale turbulence isotropy, where I_L = 〈w²〉/〈u²〉 increases with the Reynolds number and eventually converges to a certain value.

Figure 5 presents a comparison of the plots of turbulence generated by anisotropic forcing obtained from LES with those from DNS using the Lumley triangle. In addition to the results from the three SGS models, the results from implicit LES are also included. As demonstrated in the figure, the result observed in DNS, where the plots lie along the line η = –ξ, is similarly observed in all LES results. For the range of Reynolds numbers from Re = 10 to 50, the LES plots generally align with the DNS results. This range corresponds to the region where local turbulence is not isotropic, as indicated in the earlier figure showing the anisotropy of small-scale turbulence. Conversely, for Reynolds numbers ranging from 100 to 300, the LES plots appear to deviate from those of DNS. This discrepancy is observed not only in the results from the three LESs using SGS models but also in the results from implicit LES without SGS models. As indicated in the preceding figure, the range of Reynolds numbers exhibiting this discrepancy matches the interval in which the small-scale motions are effectively isotropic.

Figure 6 compares the Reynolds number dependence of turbulent kinetic energy obtained using LES with that obtained from DNS. The results are presented for both anisotropic and isotropic forcing cases. The DNS results presented in the figure are a reproduction of those shown earlier. For turbulence generated by anisotropic forcing, as shown in Fig. 6(a), the LES results generally align with the DNS results within the range of Re = 10 to Re = 40. This observation is consistently observed across all three SGS models. However, within the range of Re = 50 to Re = 200, the LES results deviate from the DNS results, exhibiting higher values. The Vreman model and coherent structure model, despite their capacity to automatically reduce SGS eddy viscosity to zero in laminar regions such as near wall surfaces, also fail to match the DNS results. Their output is generally consistent with the plots from the Smagorinsky model. Figure 6(b) compares LES results with DNS results for turbulence generated by isotropic forcing. As shown in the figure, in the case of isotropic forcing, the LES results demonstrate a high degree of agreement with the DNS results across all Reynolds numbers. This agreement is consistently observed across all three SGS models. In sharp contrast to the findings for isotropic forcing, discrepancies are observed in the LES results for anisotropic forcing.

Focusing on turbulence produced under anisotropic forcing, the study assesses Reynolds-number effects on turbulent kinetic energy through the transport equation for relevant statistical quantities [1]. Specifically, the balance of the transport equation for velocity fluctuation intensities, derived from the decomposition of turbulent kinetic energy, is calculated. The transport equations for velocity fluctuation intensities corresponds to the case where i = j in the transport equations for Reynolds stress 〈u_iu_j〉_t, where 〈〉_t denotes time-averaging operation. The transport equations for Reynolds stress 〈u_iu_j〉_t is expressed as follows:where U_i is temporal mean velocity for i direction. The transport equation consists of the following terms: the convection term, the production term, the viscous dissipation term, the pressure-strain correlation term, the turbulent diffusion term, the pressure diffusion term, the viscous diffusion term, and a contribution term from forcing. It is expected that energy is injected into the system per unit time via the contribution term from forcing, considering the entire computational domain. The viscous dissipation term quantifies the dissipation of velocity fluctuation intensities as an energy quantity due to viscosity, while the pressure-strain correlation term redistributes energy among the components of velocity fluctuation intensities. The three diffusion terms mitigate spatial inhomogeneities in the distribution of velocity fluctuation intensities. In this study, the analysis of the transport equation for Reynolds stress components is conducted exclusively using DNS results. This decision is based on the fundamental difference between DNS and LES formulations. In LES, spatial filtering introduces additional SGS stress terms that are modeled by eddy-viscosity-based assumptions, resulting in governing equations that are structurally and physically distinct from those of DNS. Due to these differences, the transport equation terms in LES, such as SGS dissipation and modeled pressure-strain effects, cannot be directly compared with their DNS counterparts. A term-wise comparison between DNS and LES would therefore be inconsistent and potentially misleading. To ensure physical interpretability without the influence of SGS model dependencies, the present analysis focuses solely on DNS-based equations.

As demonstrated in Fig. 7, the balance of velocity fluctuation intensities is contingent upon the Reynolds number. In the context of turbulence driven by anisotropic forcing, the results for both intensities of the velocity fluctuations denoted as 〈u²〉 and 〈w²〉 are presented. Here, the magnitude of each term shown is the spatially averaged value. As demonstrated in Fig. 7(a), for turbulence driven by anisotropic forcing, the energy introduced via the forcing term is balanced by the sum of the redistribution of energy to the 〈w²〉 component through the pressurestrain correlation terms and viscous dissipation. The magnitudes of the other terms are spatially averaged to zero over the entire domain. The energy budget for velocity-fluctuation intensities shows that contributions from the pressure-strain correlation terms to the componentwise kinetic energy are exactly counteracted by viscous dissipation. In contrast, as shown in Fig. 7(b), for turbulence driven by isotropic forcing, the energy input via the forcing term balances only with viscous dissipation. The spatially averaged values of each term remain approximately constant with respect to Reynolds number, except for the low Reynolds number condition. This result, where the spatially averaged values of each term exhibit no significant Reynolds number dependence, can also be observed in the case of turbulence driven by anisotropic forcing. However, it differs in distribution shape from the Reynolds number dependence of turbulent kinetic energy for anisotropic forcing, which exhibits a minimum value, as shown in the earlier figure. This discrepancy suggests that the spatially averaged balance of velocity fluctuation intensities may not adequately explain the turbulent kinetic energy distribution with a minimum value observed in the case of anisotropic forcing.

The present study focuses on the spatial inhomogeneity of the terms in the transport equations, in light of the preceding results. Given the spatial variations of the forcing field defined in this study, the statistical distributions of the turbulence stabilized by this field may be spatially inhomogeneous. As shown in Fig. 8(a), the root-mean-square (RMS) deviations of the terms in the transport equation from their spatially averaged values are presented for turbulence generated by the anisotropic forcing. As the figure illustrates, the turbulent diffusion terms, pressure diffusion terms, and pressure-strain correlation terms demonstrate higher spatial inhomogeneity compared to other terms. Furthermore, this inhomogeneity is observed to depend on the bulk Reynolds number, with the spatial inhomogeneity of these three terms reaching a minimum around Re = 50. This observation is consistent with the findings in Fig. 6(a), which depicts the bulk Reynolds number dependence of turbulent kinetic energy for turbulence generated by anisotropic forcing, also exhibiting a minimum around Re = 50. The distribution shapes of the spatial inhomogeneity of these three terms, as indicated by this result, qualitatively resemble the distribution of turbulent kinetic energy. Figure 8(b) presents the RMS deviations of the terms in the transport equation for turbulence generated by isotropic forcing as a function of Re. As demonstrated in the figure, even for turbulence generated by isotropic forcing, the RMS deviations of the turbulent diffusion term, pressure diffusion term, and pressure-strain correlation term are larger than those of the other terms. However, the Reynolds number dependence of the RMS deviations for these three terms is smaller in the case of isotropic forcing than in the case of anisotropic forcing.

Figure 9 illustrates the spatial RMS deviation of velocity fluctuation intensities, 〈u²〉 and 〈w²〉, from their spatial averages as a function of the bulk Reynolds number. In the context of spatially homogeneous velocity fluctuation intensities, the spatial RMS deviation is equivalent to zero. As demonstrated in Fig. 9, for turbulence driven by the anisotropic forcing, the spatial RMS deviation of 〈u²〉 varies with the Reynolds number. This deviation attains a minimum around Re = 100. At lower Reynolds numbers, the RMS deviation decreases with increasing Re, approaching this minimum. At higher Reynolds numbers, the RMS deviation increases as Re rises beyond the minimum. The anisotropic forcing field applied in this study remains constant across all Reynolds numbers. Consequently, this result suggests that the spatial homogenity of the generated turbulence changes, reaching its highest level around Re = 100. However, this Reynolds number dependence of the RMS deviation is not observed for the second intensity 〈w²〉. In contrast, for turbulence driven by isotropic forcing, the spatial RMS deviation remains approximately constant across all Reynolds numbers and is smaller than that observed in the case of anisotropic forcing. This indicates that the spatial uniformity of turbulence statistics under isotropic forcing remains nearly constant with respect to the Reynolds number and is higher than that under the anisotropic forcing.

Figure 10 now offers a four-way comparison of flow visualizations obtained from DNS and the three eddy-viscosity SGS models (Smagorinsky, Vreman, and the coherent-structure model). All frames display green isosurfaces of negative static-pressure fluctuation because this quantity remains robust even on the coarser LES grids and therefore provides a scale-consistent indicator of the large-scale tube-like vortices generated by the anisotropic forcing. Although the Q-criterion (the second invariant of the velocity-gradient tensor) is also plotted as white isosurfaces in the DNS panels, it is intentionally omitted from the LES panels: the enhanced numerical and subgrid dissipation in LES, combined with the larger grid spacing, suppresses the small-scale vortical signal and yields sparsely distributed Q isosurfaces that obscure, rather than clarify, the intended large-scale comparison. This unified visualization therefore isolates the capability of each SGS model to reproduce (or over-predict) the deterministic two-dimensional tubular structures while avoiding artefacts introduced by under-resolved Q-criterion fields. For the DNS cases shown in Fig. 10(a1)–(c1), the white surfaces signify the isosurfaces of the second invariant of the velocity gradient tensor, thus emphasizing small-scale turbulence structures. The results are presented for Re = 10, 100, 300 based on the Reynolds-number-dependent distribution of turbulent kinetic energy. As demonstrated in the figures, for the DNS results at both Re = 10 and Re = 300, two-dimensional tubular flow structures, rendered by the green isosurfaces, are discernible. Conversely, at Re = 100, such tubular structures are not observed, and appear more fragmented. The stable presence of tubular flow structures is considered to enhance the spatial inhomogeneity of the turbulence, while the absence of such structures at a Reynolds number of 100 is thought to reduce the spatial inhomogeneity of the turbulence. This observation is qualitatively consistent with the earlier findings of the three terms in the transport equation and the low spatial inhomogeneity of velocity fluctuations observed at a Reynolds number between 50 and 100. Figure 10(a2)–(c2), (a3)–(c3), and (a4)–(c4) illustrate the visualization outcomes of anisotropic turbulence, as predicted by the Smagorinsky model, vreman model, and the coherent structure model, respectively, as functions of the Reynolds number. As demonstrated, for specific values of the Reynolds number, namely, 10 and 300, the flow fields obtained from both LES analyses exhibit two-dimensional tubular structures that are analogous to those observed in the DNS results.

Conversely, at Re = 100, such tubular structures are not clearly observed and appear more fragmented. This result was confirmed to be robust across a range of negative pressure thresholds used for the visualization, as well as at multiple time instances under statistically steady conditions. It should be noted that the formation of the tubular structures is induced by the anisotropic spatial structure of the forcing field, and their appearance at high Reynolds numbers is a consequence of the periodicity and symmetry embedded in the forcing—thus, such structures are expected to emerge deterministically. However, at Re = 100, the breakdown of these tubular structures was found to occur intermittently over time, suggesting that their absence is not static but rather temporally varying. This behavior is clearly illustrated in Fig. 11, where largescale vortex structures characterized by negative static pressure fluctuations are visualized at several time instances. According to the figure, well-defined tubular formations are less apparent in the DNS results for Re = 100 than for the higher Reynolds numbers. and the overall flow field is dominated by fragmented and unstable patterns. Nevertheless, temporal observations confirm that the disappearance of these structures occurs in an intermittent and transient manner. In other words, the dynamics at Re = 100 suggest the existence of periodic or quasi-periodic collapse and reconstruction of the tubular structures. Such breakdown phenomena are not prominently observed in the LES results, where twodimensional vortex tubes are stably sustained at many time instances. The mechanism of this intermittent disappearance and reappearance of the tubular structures remains an open issue and is expected to be addressed in future studies through time-resolved statistical and spectral analyses.

This finding is in qualitative disagreement with the DNS visualization results. The Reynolds-number-dependent distribution of turbulent energy obtained from DNS demonstrates that the increase in turbulence inhomogeneity corresponds qualitatively to the increase in turbulent energy. These visualization results suggest that the LES analyses reproduce two-dimensional tubular structures that are not observed in the DNS results, indicating that the inhomogeneity of turbulence is enhanced, thereby increasing turbulent kinetic energy. At Re = 100, the turbulent kinetic energy in LES was significantly overestimated compared to DNS, particularly with the Smagorinsky model. This discrepancy is attributed to the inability of the eddyviscosity SGS models to capture the transient collapse of coherent structures observed in the DNS. While LES predicts persistent twodimensional tubular structures, DNS results show intermittent fragmentation, which leads to lower spatial inhomogeneity and reduced energy. These differences in structure and transport terms provide a physical basis for the divergence in flow predictions between LES and DNS.