Abstract
Computational fluid dynamics (CFD) results are presented for synthetic turbulence generation by a proposed statistically targeted forcing (STF) method. The new method seeks to introduce a fluctuating velocity field with a distribution of first and second moments that approximate a user-specified target mean velocity and Reynolds stress tensor, by incorporating deterministic time-dependent forcing terms into the momentum equation for the resolved flow. The STF method is formulated to extend the applicability of previously documented methods and provide flexibility in regions where synthetic turbulence needs to be generated or sustained, for use in engineering level large eddy and hybrid large eddy/Reynolds-averaged Navier–Stokes CFD simulations. The objective of this study is to evaluate the performance of the proposed STF method in large-eddy simulation (LES) simulations of isotropic and anisotropic homogeneous turbulent flow test cases. Results are interrogated and compared to target statistical velocity and turbulent stress distributions and evaluated in terms of energy spectra. Analysis of the influence of STF model parameters, mesh resolution, and LES subgrid stress model on the results is investigated. Results show that the new method can successfully reproduce desired statistical distributions in a homogeneous turbulent flow.
1 Introduction
Computational fluid dynamics (CFD) methods for effectively reproducing time-dependent turbulence boundary and/or initial (B/I) conditions have significant potential value for improving simulation of many engineering systems. Because advances in computational resources have made turbulent scale-resolving simulations such as large eddy simulation (LES) feasible for some industrially relevant flows, accurate and efficient methods for prescribing these complex conditions are increasingly needed.
Research on generation of turbulence B/I conditions has been active over the past two decades, ranging from library-based methods to recycling/rescaling to synthetic turbulence generation (STG) with controlled forcing. These methods become increasingly important when the application of steady-state B/I conditions causes development of resolved turbulence fluctuating velocities in the simulation that is either delayed, inaccurate, or nonexistent. It is often not feasible to include the source of turbulent B/I conditions within a simulation. A simple example is simulation of the interaction of an aerodynamic vehicle with a turbulent freestream flow. The source of the freestream turbulence is in fact due to the interaction of the atmospheric boundary layer with the ground, but available resources would typically preclude simulating the atmospheric boundary layer flow in addition to the local vehicle aerodynamics. The goal of B/I condition methods such as STG is to replace turbulent content obtained from fully resolved simulations with a reasonable approximation of turbulence for a substantially lower computational cost. In addition, with STG methods, turbulent content can be selectively located in specific regions of the computational domain or on the boundaries, providing flexibility such that turbulence B/I conditions are only used in regions where they are needed.
One well-known method for prescribing turbulence boundary conditions is recycling/rescaling. For recycled turbulent content, streamwise periodic boundary conditions are imposed on the domain or a portion of the domain such that the turbulent flow leaving the outlet is reintroduced at the inlet. Rescaling of the velocity field can be performed to ensure that the turbulent statistics remain appropriately spatially developing. This method was used, for example, by Spalart et al. [1] to perform LES of a turbulent boundary layer. Lund et al. [2] used the recycling/rescaling method to perform an auxiliary simulation of a turbulent boundary layer, and then extracted planes of time-dependent velocity data to be mapped to the inlet of a simulation with a more complex geometry. Several other studies have extended the recycling/rescaling approach to simulate complex wall-bounded flows [3–5]. Schlüter et al. [6] used the recycling/rescaling method to impose fluctuating velocities at the outlet to an LES region of a simulation to impose the statistics obtained from a Reynolds-averaged Navier–Stokes (RANS) solution in the downstream region.
A general class of methods that represent an alternative to recycling/rescaling is STG. For applications of practical engineering interest on complex geometries, STG methods have the potential to reproduce turbulent fluctuations at desired locations and with desired statistical distributions, without the need to run an auxiliary simulation. STG methods can be used to specify inflow boundary conditions or as initial conditions for a simulation.
Kraichnan [7] proposed one of the first STG methods for isotropic turbulence, by utilizing a spectral approach to artificially produce an isotropic turbulent velocity field from random Fourier modes. This approach of generating isotropic velocity fields with a specified energy spectrum has been used, for example, to generate initial conditions for DNS of isotropic turbulence [8,9]. Lee et al. [10] similarly proposed a Fourier transform-based STG method to generate inflow boundary conditions, however, one limitation of this method is that it is not applicable to wall-bounded flows due to statistical inhomogeneity in the wall-normal direction.
Using a similar approach for isotropic turbulence, Lundgren [11] defined a forcing term in the momentum equations that is proportional to the fluctuating velocity component. This isotropic linear forcing (ILF) term imitates the natural production mechanism in the turbulent kinetic energy equation. This ILF forcing can be restricted to low wave number modes when using spectral numerical methods. Rosales et al. [12] extended the method in Ref. [11] by formulating the forcing term in physical space.
A different algorithmic approach proposed by Jarrin et al. [13,14] is the synthetic eddy method (SEM), which is used to generate realistic synthetic eddies at the inflow of an LES simulation. Results have shown that the synthetic eddy field can evolve to physically realistic turbulent flow after a relatively short distance downstream of the inlet. Some limitations exist in the SEM method such as depletion of the smaller scales of turbulence. This has motivated modification of the SEM to include momentum source terms that energize the velocity fluctuations for some distance downstream.
Keating et al. [15] explain how inappropriate modeling of the scale and structure of synthetic turbulence can lead to a rapid dissipation of velocity fluctuations and an increase in the time/distance required for the flow to recover into a fully turbulent state. Therefore, most recent STG methods attempt to introduce some degree of spatial and temporal coherence through artificial control forcing techniques. Spille-Kohoff and Kaltenbach [16] proposed an inflow STG method based on an added forcing source term in the wall-normal momentum equation. This forcing term enhances the velocity fluctuations in that direction, to match a desired target profile of Reynolds shear stress. This technique enhances the wall-normal fluctuations at discrete locations, with amplitude proportional to the difference between the calculated Reynolds shear stress and a provided target profile. This is an example of an STG method with a controlled forcing feedback loop to achieve a target statistical distribution. They documented that the method reduces the error in the Reynolds shear stress to acceptable values within five channel half heights, although the coefficient of friction and the turbulent kinetic energy required longer downstream distances to reach their fully developed values.
More recently, de Laage de Meux et al. [17] proposed a method to impose target statistics of the flow in terms of mean velocity and resolved turbulent stress, using a method denoted as anisotropic linear forcing. The time-dependent forcing function is proportional to the instantaneous velocity via a tensor transformation. The method was found to provide accurate results for isotropic, anisotropic and spatially developing turbulence test cases for LES and hybrid RANS-LES simulations.
The objective of this study is to investigate a newly proposed statistically targeted forcing (STF) method for synthetic turbulence generation. The STF method is a variant of STG with controlled forcing within the simulation domain, implemented via added source terms in the momentum and energy equations. The method can be viewed to act as a restoring force toward a target statistical state within either a time-averaging or volume-averaging framework. Simulation results are presented for homogeneous isotropic and anisotropic turbulent flow, and results are evaluated in terms of one-point statistics and spectral characteristics.
2 Statistically Targeted Forcing Method
2.1 General Description of the Method.
where the overbar denotes either Reynolds or Favre (mass-weighted) averaging.
In practice, the method is implemented as follows: First, a target statistical velocity distribution is specified prior to the simulation in terms of and. As the simulation proceeds, the resolved statistics and are obtained using an appropriate averaging technique. At each time-step, the transformation tensor is computed at each point in the domain based on and . During each iteration, the fluctuating velocity is computed, and the target instantaneous velocity is found using Eq. (15). The forcing term is then computed using Eq. (6) and included as an additional source term in the momentum and energy equations.
2.2 Ensemble Averaging.
Here, is the time-step size and is the current number of time steps in the simulation. In the limit , the averaged value becomes constant for statistically stationary flow.
2.3 Spatial Filtering.
2.4 Prescription of the Forcing Time Scale.
where is the characteristic large-eddy turbulent time scale used in the spatial filter defined by Eq. (26). Depending on the type of simulation performed, the turbulent time scale can be approximated by a characteristic imposed length scale, or from the source of the target statistical distribution. For example, if a precursor RANS simulation is used to specify the target mean velocity and Reynolds stress tensor, the dissipation time scale can be used to specify .
3 Simulation Details
All simulations in this study were performed using the open source CFD code flowpsi [19], a finite volume density-based solver constructed in C++ using the Loci framework. FlowPsi uses high-resolution approximate Riemann solvers and implicit numerical methods. For this study, all simulations were run with sufficiently low Mach number to simulate incompressible flow conditions. Inviscid fluxes are reconstructed using a modified skew symmetric flux (SSF) scheme. The SSF scheme is a generalization of the kinetic energy consistent central difference scheme of Subbareddy and Candler [20], blended with a small second-order upwind flux contribution. For all simulations in this paper, the blending distribution was 95% central difference and 5% upwind. The SSF scheme has been shown to provide low numerical dissipation and effective resolution of high wavenumber velocity and pressure modes in unsteady turbulent flow simulations.
where is the characteristic mesh size, equal to the cube root of cell volume in the current simulations, and the coefficient = 0.1.
One test simulation was run using an implicit LES approach similar to monotonically integrated large eddy simulation (MILES) [22], for which the subgrid dissipation was assumed to be modeled by the numerical dissipation inherent in the blended upwind portion of the SSF inviscid flux formulation. As shown in the results section, this led to a build up of energy at high wavenumbers, and as a consequence all of the other simulations used the Smagorinsky model. We note that the implicit LES results could perhaps be improved by increasing the upwind contribution to the inviscid flux.
4 Test Cases
This study consists of two test cases, isotropic and anisotropic homogeneous turbulence. The STF method is used to produce and sustain time-dependent, stationary flow in the domain. The resulting flow field is interrogated to evaluate the accuracy of the method for producing realistic turbulent flow conditions with the prescribed first- and second-order one-point statistics. For the anisotropic turbulence test case, unequal normal stress components and a nonzero shear stress component were imposed as the target statistical flow field. Key issues investigated include:
Effect of mesh resolution
Capability of the method to accurately reproduce mean velocity and Reynolds stresses in isotropic and anisotropic turbulent flows
Effect of averaging techniques and method parameters such as turbulence time scale, , and forcing coefficient,
Spectral characteristics of the turbulence generated by the forcing method
The domain for the homogeneous isotropic and anisotropic test cases (“turbulence in a box”) is a cube with side length L, with all boundaries periodic. Flow was initialized for all cases to zero velocity and gage pressure. Target Reynolds stress distributions were selected such that the Mach number based on maximum instantaneous velocity was approximately 0.1, which effectively reproduces incompressible flow conditions. For the isotropic turbulence cases, the Reynolds number based on the Taylor microscale,, and kinematic viscosity, varied from 471 for the finest grid to 1154 for the coarsest grid. Four structured, uniform, Cartesian meshes were used for the simulations, corresponding to 323, 643, 1283, and 1923 cells, and they are illustrated in Fig. 3.

Planar view of the computational grids for homogeneous turbulence cases showing four mesh resolution levels (a) 323, (b) 643, (c) 1283, and (d) 1923
5 Results and Discussion
5.1 Homogeneous Isotropic Turbulence.
5.1.1 Baseline Results.
Initial simulations employing the STF model were run using volume averaging for turbulence statistics, and no spatial filtering. Four different values of the forcing time scale were investigated, corresponding to = 0.45, 0.045, 0.0045, and 0.00033. Forcing was applied during the time interval. Simulations were run for an additional interval after forcing was removed in order to observe the behavior of the velocity field during turbulence decay after being initialized by the STF method.
Figure 4 shows contours of velocity magnitude on the periodic bounding surfaces of the domain, corresponding to a simulation time of , at which point the fluctuating velocity field had reached an apparently stationary state. The results are shown for a forcing time scale = 0.045. The velocity field exhibits qualitative features of turbulent flow, including a visible range of spatial scales that resemble turbulent eddies. As the mesh resolution is increased from 323 to 1923, the size of the smallest resolved eddies is reduced. The effect of mesh resolution on the overall qualitative flow structure, however, appears to be minimal.

Contours of instantaneous velocity magnitude for forcing simulation of homogeneous isotropic turbulence with forcing time scale , taken at a simulation time . Four mesh resolution levels are shown: (a) 323, (b) 643, (c) 1283, and (d) 1923.
Figure 5 shows contours of velocity magnitude for the 1283 mesh for simulations using three different values of the forcing time scale. Results are shown on one bounding surface of the domain, at a simulation time of . It is apparent that the qualitative flow structure is similar for all three cases. The maximum resolved velocity increases as the forcing time scale is reduced. This is expected since a smaller time scale yields an overall larger forcing term, in effect driving the local velocity more strongly toward the target velocity vector computed by Eq. (15) at each time-step.

Planar view contours of instantaneous velocity magnitude for forcing simulation of homogeneous isotropic turbulence on 1283 mesh, taken at a simulation time . Results are shown for three different values of the forcing time scale.
Figure 6 shows the time history of the ratio of resolved turbulent kinetic energy,, to the target value, . Figure 6(a) compares different mesh resolution levels and Fig. 6(b) compares different values of the forcing time scale on the 1283 mesh. The turbulence level can be seen to increase relatively rapidly, and on all grids the ratio levels off at an apparently stationary state near the target value in under one characteristic simulation time. As seen in Fig. 6(b), the time required to reach the stationary state is reduced as the forcing time scale is reduced. The energy level remains nearly constant throughout the time when forcing is applied, and once forcing is removed a power law decay is apparent. The resolved turbulence level during stationary forcing appears to be insensitive to the mesh resolution level. As seen in Fig. 6(b), the decay rate once forcing is removed changes as the mesh resolution is increased but appears to be nearly grid independent as the mesh size is increased from 1283 to 1923. The decay rate is relatively insensitive to the forcing time scale, which is expected since the forcing term has no direct effect on the flow dynamics after .

Time evolution of turbulent kinetic energy (k/k*) for forcing simulation of homogeneous isotropic turbulence with volume averaging and no spatial filtering: (a) effect of mesh resolution and (b) effect of forcing time scale
Figure 7 highlights the resolved turbulence level reached in the simulations during the forcing time interval, expressed as the ratio , for different mesh resolution levels and forcing time scales. As shown previously, the level of resolved turbulence is strongly dependent on the forcing time scale. On the 1283 mesh, the turbulence level ratio increases from 0.85 to 0.999 as the time scale, , is reduced from 0.45 to 0.00033. The level of resolved turbulence is only weakly dependent on the mesh resolution level. For example, for a forcing time scale of , the ratio increases from 0.84 to 0.87 as the mesh is refined from 323 to 1923 cells.

Ratio of resolved-to-target turbulent kinetic energy, k/k* for (a) different mesh sizes and (b) on 1283 mesh at different values of the forcing time scale parameter
To understand the spectral characteristics of the turbulence generated by the synthetic forcing method, the fast Fourier transform was applied to the velocity field and integrated over spherical shells in wavenumber ( space to obtain the energy density as a function of wavenumber. The results are shown in Fig. 8. Also shown for reference is the Kolmogorov −5/3 scaling for the inertial subrange. In Fig. 8(b), the energy containing range is indicated as , the dissipation range is indicated as, and the inertial range is indicated as The spectra obtained from all four grids look similar at the lowest wavenumbers while significant difference is apparent at the highest wavenumbers. The two most refined grids (1283 and 1923) appear to show an inertial range for which the behavior qualitatively matches the −5/3 law, while the 323 and 643 grids do not reproduce this behavior. This may explain the decay behavior seen in Fig. 6(a). Once the synthetic generated turbulence field is sufficiently refined to include an inertial range and proper energy cascade dynamics, the decay rate once forcing is removed is relatively insensitive to further refinement. On all meshes, the energy is damped rapidly for wavenumbers larger than 1/N, where N is the number of cells in each of the three coordinate directions.

Normalized energy density spectrum for forcing simulation of homogeneous isotropic turbulence with volume averaging on (a) different meshes at and (b) 1283 mesh,
For the 1283 grid results shown in Fig. 8(b), energy spectra obtained using different forcing time scales () are all in good quantitative agreement, with only slight differences evident in the low wavenumbers (). Note that the plots are scaled by the maximum value of energy density. Since the total resolved energy increases as forcing time scale is reduced, the unscaled plots would show a vertical shift for different values of . Overall, the figure shows spectral characteristics indicative of three-dimensional turbulent flow for all values of forcing time scale and for meshes that are sufficiently refined to yield an inertial scaling range.
To confirm that the STF method works independently of the specific LES model used, one set of simulations was run using no explicit subgrid stress term, for which the numerical dissipation inherent in the blended upwind inviscid flux term was used to dissipate energy in the small resolved scales. This approach is an informal implementation of the MILES modeling methodology. Results are shown in Fig. 9 for simulations on the 1283 grid, using three different values of the forcing time scale. For all three cases, the overall level of turbulent kinetic energy was comparable to the cases using the Smagorinsky SGS model. Again, it is evident that the level of resolved turbulence increases as the forcing time scale is reduced, similar to results with the Smagorinsky model shown in Fig. 5. The spectral behavior was also similar for the low wavenumber portion, displaying a clearly identifiable region with −5/3 inertial scaling. In the higher wavenumber region near the filter cutoff, there is evidence of energy pile up indicating that there is insufficient dissipation to effectively represent the forward scatter of energy to the subfilter velocity scales. It is likely that increasing the upwind contribution to the inviscid flux term would improve the result. The differences between results in Figs. 8 and 9 arise due to the SGS model (or lack thereof) rather than the STF method itself, and the results indicate that the STF method is agnostic to the details of LES model, as expected.

(a) Contours of instantaneous velocity magnitude and (b) normalized energy spectra for implicit LES simulation with no subgrid stress model, on the 1283 grid, using three different values of
5.1.2 Effect of Averaging Technique.
Figure 10 compares contours of instantaneous velocity magnitude for equivalent simulations using volume averaging (VA) versus time-averaging (TA) to compute turbulence statistics on the 1283 grid at 4.5, 0.45, and 0.045. The results show no obvious qualitative differences with regard to the flow structure, and regardless of averaging method used, the maximum velocity magnitude increases as the forcing time scale is decreased. Figure 11 shows the time evolution of turbulent kinetic energy for different values of forcing time scale. Regardless of averaging method, as the value of is decreased, the resolved turbulence level more quickly approaches the target value, as shown previously. However, cases using volume averaging show a nearly monotonic increase to the final value, while the cases with time averaging show oscillatory behavior, with the magnitude of oscillation increasing as is decreased. This effectively represents a temporal lag in the statistics that are used to compute the forcing term. It is expected that the results will converge to a nearly constant value given a sufficiently long run time, but as shown in the figures the time required to reach this is longer than a single characteristic simulation time, . Figure 12 compares the energy spectra obtained for three different forcing time scales, using volume and time averaging. The overall shape of the spectra is similar regardless of which averaging method is used. Table 1 compares the quasi-stationary turbulent kinetic energy level achieved for each case. The key point is that consistent with results shown previously, a lower value for the forcing time scale produces more accurate results, regardless of which averaging method is used.

Contours of instantaneous velocity magnitude for forcing simulation of homogeneous isotropic turbulence using VA (above) and TA (below) to compute turbulence statistics on 1283 grid at

Time evolution of turbulent kinetic energy, k/k* for forcing simulation of homogeneous isotropic turbulence using VA (solid line) and TA (dash line) to compute turbulence statistics on 1283 grid, at

Normalized energy spectra for forcing simulation of homogeneous isotropic turbulence using VA (solid line) and TA (dash line) to compute turbulence statistics on 1283 grid, at

Contours of instantaneous velocity magnitude for forcing simulation of homogeneous isotropic turbulence, comparing results of (a) spatially filtered with (b) nonspatially filtered simulations, on 1283 grid, ( and)

Time evolution of turbulent kinetic energy, k/k* for forcing simulation of homogeneous isotropic turbulence, comparingresults of SF with Non_SF simulations, on 1283 grid, ( and)
Percentage ratio of resolved-to-target turbulent kinetic energy for different averaging techniques and forcing time scales at
τf/τs = 4.5 | τf/τs = 0.45 | τf/τs = 0.045 | |
---|---|---|---|
TA | 57% | 74% | 102% |
VA | 37% | 88% | 99% |
τf/τs = 4.5 | τf/τs = 0.45 | τf/τs = 0.045 | |
---|---|---|---|
TA | 57% | 74% | 102% |
VA | 37% | 88% | 99% |
Note: Accuracy of the result increases with decrease in dimensionless forcing time scale.
5.1.3 Effect of Spatial Filtering.
For all of the results shown in Secs. 5.1.1 and 5.1.2, no spatial filtering was used in the forcing term and the characteristic large-eddy length scale was effectively imposed by the size of the domain. This is apparent in the previous energy spectra plots, which show peak energy density at a wavenumber slightly less than 2/L. As discussed in Sec. 2.3, a spatial filter is applied by specifying a characteristic turbulent time scale as well as the target Reynolds stress tensor. For the isotropic turbulence simulations, the characteristic (isotropic) spatial filter size, as indicated in Eq. (27), is . Several simulations were performed on the 1283 grid to demonstrate the effect of spatial filtering.
Figure 13 shows contours of instantaneous velocity magnitude comparing a spatially filtered (SF) case with the nonspatially filtered (Non_SF) case on the 1283 grid, with and. This corresponds to a forcing time scale of , and for the spatially filtered case corresponds to a filter size . Volume averaging was used to compute velocity statistics. The figure shows that the maximum instantaneous velocity magnitude is approximately equal for both cases, but the turbulent length scales are visibly smaller for the spatially filtered case. The corresponding plots of evolution of turbulent kinetic energy are shown in Fig. 14. It is apparent that the spatially filtered case reaches a stationary state faster than the nonfiltered case, and that the overall level of resolved turbulence is lower. This is due to the fact that the effective turbulence production introduced by the forcing term is balanced in the simulation by the SGS (and to a lesser extent numerical) dissipation, and dissipation is greater when the resolved turbulent structures are smaller. It is also interesting to note that, once forcing is removed, the spatially filtered case shows more rapid decay of turbulent energy due to the fact that the initialized flow field obtained from the STF has a smaller integral length scale.

Normalized energy spectra for forcing simulation of homogeneous isotropic turbulence, comparing results of SF with Non_SF simulations, on 1283 grid, using volume-averaging results at and

Contours of instantaneous velocity magnitude for forcing simulation of homogeneous isotropic turbulence, with different time-scale targets imposing different turbulent length scales, on 1283 grid, using volume-averaging results with , at (a) (b), (c) , (d) , and (e)
Figure 15 shows the normalized energy spectra for the simulations with and without spatial filtering, for and . The effect of spatial filtering is evident, with the maximum energy density occurring at , which corresponds to . The remainder of the spectrum, including the inertial and dissipation ranges, is qualitatively similar for the two cases.

Time evolution of turbulent kinetic energy, k/k* for forcing simulation of homogeneous isotropic turbulence, with different time-scale targets imposing different turbulent length scales, on 1283 grid, using volume-averaging results with , and 100, at ,, and
To further highlight the effect of spatial filtering, simulations were run using four different values of the characteristic large-eddy timescale (), and a constant value of the forcing coefficient . Figure 16 shows instantaneous velocity magnitude contours for the four different cases. It is apparent that the scale of the visible resolved turbulent flow structures decreases as the filter timescale is reduced, as expected. The temporal evolution of the resolved turbulent energy for the four cases is shown in Fig. 17. It should be noted that, since is held constant, as is reduced the effective forcing timescale is correspondingly reduced and the magnitude of the forcing term increases. As a consequence, as is decreased the simulations more rapidly reach a quasi-stationary state. Interestingly, the level of turbulence reached is apparently insensitive to the imposed time (and length) scale, and the values reached are all approximately equal for a given value of the forcing coefficient . For comparison purposes, the simulation with and is also shown. For that case, the resulting turbulence level is much closer to the target value. Note also that the once forcing is removed after a simulation time , the turbulent energy decay rate increases (i.e., turbulence decays faster) as is reduced.

Normalized energy spectra for forcing simulation of homogeneous isotropic turbulence, with different time-scale targets imposing different turbulent length scales, on 1283 grid, using volume-averaging results with: , at ,, and

Contours of instantaneous velocity magnitude for forcing simulation of homogeneous anisotropic turbulence with volume averaging, on 1283 grid, for ZSS results (left) with = 4.5 and 0.45 and imposed FSS results (right) with = 4.5 and 0.45
Figure 18(a) shows normalized energy spectra for the cases with varying characteristic large eddy time scale (), and a constant value of the forcing coefficient . The effect of spatial filtering is apparent, as the peak energy density moves toward higher wavenumbers as is decreased. The shift in peak wavenumber scales inversely with the effective spatial filter size, , suggesting that filtering is an effective means of controlling the large eddy length scale in the synthetically generated turbulent velocity field. For all cases except that with the smallest value of , there is an evident inertial range that approximately follows Kolmogorov −5/3 scaling.

Time evolution of turbulent kinetic energy, k/k* for forcing simulation homogeneous anisotropic turbulence with volume averaging, on 1283 grid, (a) zero shear stress (ZSS) results, at = 0.45 and 0.045, (b) ZSS and FSS results, at = 0.45, and (c) ZSS and FSS results, at = 4.5
Results for isotropic turbulence are summarized in Table 2. For all cases the statistically targeted forcing method reproduces some level of resolved turbulence that approximates an isotropic flow. As seen in the previous results, for an appropriately refined mesh in which sufficient separation exists between the largest resolved length scale and the mesh size, a portion of the turbulent energy spectrum arises that approximates Kolmogorov inertial scaling. The degree to which the prescribed target Reynolds stress tensor is reproduced depends directly on the magnitude of the forcing time scale, which controls the overall strength of the forcing term. As seen in Table 2, all of the cases that produced at least 90% of the target turbulent kinetic energy had forcing time scale ratios () less than 0.05. It is worth noting that none of the cases studied exhibited any numerical instability, which suggests that the STF method can be effectively implemented with a judicious choice of forcing time scale that allows both accurate and stable solutions for synthetic turbulence generation.
Summary of turbulence statistics for homogeneous isotropic STF cases (green color implies target and over 90% of target statistics)
Target Reynolds stress | |||||||||
---|---|---|---|---|---|---|---|---|---|
0.6666 | 0.6666 | 0.6666 | |||||||
Averaging method | Spatially filtered | Grid | Max % difference | (k/k*) ratio % | |||||
Volume | No | 0.045 | 323 | 0.6601 | 0.6560 | 0.6665 | 1.59 | 99.13 | |
Volume | No | 0.45 | 323 | 0.5725 | 0.5433 | 0.5684 | 18.49 | 84.21 | |
Volume | No | 0.045 | 643 | 0.6562 | 0.6522 | 0.6495 | 2.56 | 97.90 | |
Volume | No | 0.45 | 643 | 0.5725 | 0.5505 | 0.5701 | 17.41 | 84.66 | |
Volume | No | 0.045 | 1283 | 0.6567 | 0.6532 | 0.6548 | 2.01 | 98.23 | |
Volume | No | 0.45 | 1283 | 0.5811 | 0.5433 | 0.5776 | 18.49 | 85.10 | |
Volume | No | 4.5 | 1283 | 0.2211 | 0.2359 | 0.2151 | 67.73 | 33.60 | |
Volume | Yes | 0.45 | 0.045 | 1283 | 0.6211 | 0.6203 | 0.6191 | 7.13 | 93.03 |
Volume | Yes | 0.45 | 0.45 | 1283 | 0.4881 | 0.5009 | 0.4996 | 26.78 | 74.43 |
Volume | Yes | 0.45 | 4.5 | 1283 | 0.2098 | 0.2238 | 0.2041 | 69.39 | 31.88 |
Volume | No | 0.45 | 1923 | 0.5778 | 0.5887 | 0.5800 | 13.33 | 87.32 | |
Volume | No | 0.00033 | 1283 | 0.6665 | 0.6665 | 0.6665 | 0.01 | 99.98 | |
Volume | No | 0.0033 | 1283 | 0.6655 | 0.6658 | 0.6649 | 0.26 | 99.81 | |
Volume | No | 0.0045 | 1283 | 0.6651 | 0.6651 | 0.6655 | 0.23 | 99.78 | |
Volume | No | 0.045 | 1283 | 0.6563 | 0.6515 | 0.6539 | 2.27 | 98.08 | |
Volume | No | 0.067 | 1283 | 0.6379 | 0.6501 | 0.6456 | 4.31 | 96.68 | |
Volume | Yes | 0.45 | 0.45 | 1283 | 0.4881 | 0.5009 | 0.4996 | 26.78 | 74.43 |
Volume | Yes | 0.225 | 0.225 | 1283 | 0.4875 | 0.4947 | 0.4913 | 26.87 | 73.68 |
Volume | Yes | 0.112 | 0.112 | 1283 | 0.4973 | 0.4969 | 0.4985 | 25.45 | 74.64 |
Volume | Yes | 0.056 | 0.056 | 1283 | 0.4912 | 0.4909 | 0.4902 | 26.46 | 73.62 |
Time | Yes | 0.45 | 0.045 | 1283 | 0.6128 | 0.5968 | 0.5880 | 11.80 | 89.88 |
Time | No | 0.045 | 1283 | 0.6544 | 0.6378 | 0.6363 | 4.54 | 96.42 | |
Time | No | 0.45 | 1283 | 0.5464 | 0.5108 | 0.4923 | 26.14 | 77.47 | |
Time | No | 4.5 | 1283 | 0.3467 | 0.2871 | 0.3161 | 56.93 | 47.50 |
Target Reynolds stress | |||||||||
---|---|---|---|---|---|---|---|---|---|
0.6666 | 0.6666 | 0.6666 | |||||||
Averaging method | Spatially filtered | Grid | Max % difference | (k/k*) ratio % | |||||
Volume | No | 0.045 | 323 | 0.6601 | 0.6560 | 0.6665 | 1.59 | 99.13 | |
Volume | No | 0.45 | 323 | 0.5725 | 0.5433 | 0.5684 | 18.49 | 84.21 | |
Volume | No | 0.045 | 643 | 0.6562 | 0.6522 | 0.6495 | 2.56 | 97.90 | |
Volume | No | 0.45 | 643 | 0.5725 | 0.5505 | 0.5701 | 17.41 | 84.66 | |
Volume | No | 0.045 | 1283 | 0.6567 | 0.6532 | 0.6548 | 2.01 | 98.23 | |
Volume | No | 0.45 | 1283 | 0.5811 | 0.5433 | 0.5776 | 18.49 | 85.10 | |
Volume | No | 4.5 | 1283 | 0.2211 | 0.2359 | 0.2151 | 67.73 | 33.60 | |
Volume | Yes | 0.45 | 0.045 | 1283 | 0.6211 | 0.6203 | 0.6191 | 7.13 | 93.03 |
Volume | Yes | 0.45 | 0.45 | 1283 | 0.4881 | 0.5009 | 0.4996 | 26.78 | 74.43 |
Volume | Yes | 0.45 | 4.5 | 1283 | 0.2098 | 0.2238 | 0.2041 | 69.39 | 31.88 |
Volume | No | 0.45 | 1923 | 0.5778 | 0.5887 | 0.5800 | 13.33 | 87.32 | |
Volume | No | 0.00033 | 1283 | 0.6665 | 0.6665 | 0.6665 | 0.01 | 99.98 | |
Volume | No | 0.0033 | 1283 | 0.6655 | 0.6658 | 0.6649 | 0.26 | 99.81 | |
Volume | No | 0.0045 | 1283 | 0.6651 | 0.6651 | 0.6655 | 0.23 | 99.78 | |
Volume | No | 0.045 | 1283 | 0.6563 | 0.6515 | 0.6539 | 2.27 | 98.08 | |
Volume | No | 0.067 | 1283 | 0.6379 | 0.6501 | 0.6456 | 4.31 | 96.68 | |
Volume | Yes | 0.45 | 0.45 | 1283 | 0.4881 | 0.5009 | 0.4996 | 26.78 | 74.43 |
Volume | Yes | 0.225 | 0.225 | 1283 | 0.4875 | 0.4947 | 0.4913 | 26.87 | 73.68 |
Volume | Yes | 0.112 | 0.112 | 1283 | 0.4973 | 0.4969 | 0.4985 | 25.45 | 74.64 |
Volume | Yes | 0.056 | 0.056 | 1283 | 0.4912 | 0.4909 | 0.4902 | 26.46 | 73.62 |
Time | Yes | 0.45 | 0.045 | 1283 | 0.6128 | 0.5968 | 0.5880 | 11.80 | 89.88 |
Time | No | 0.045 | 1283 | 0.6544 | 0.6378 | 0.6363 | 4.54 | 96.42 | |
Time | No | 0.45 | 1283 | 0.5464 | 0.5108 | 0.4923 | 26.14 | 77.47 | |
Time | No | 4.5 | 1283 | 0.3467 | 0.2871 | 0.3161 | 56.93 | 47.50 |
5.2 Homogeneous Anisotropic Turbulence.
Figure 19 shows the contours of instantaneous velocity magnitude for both anisotropic test cases with forcing time scale = 4.5 and 0.45. Similar to the isotropic turbulence results, the maximum velocity magnitude increases with a decrease in the forcing time scale. Likewise, the structure of the velocity field is qualitatively similar to the isotropic results shown in Fig. 5. The time history of the turbulent kinetic energy for the anisotropic case is shown in Fig. 20. The behavior again resembles that of the isotropic test case. It is interesting to note in Fig. 20(b) that the decay of turbulent kinetic energy differs for the two different anisotropic cases, with more rapid decay occurring for the case with finite shear stress components. Figure 21 shows the energy spectra for the anisotropic cases. Consistent with the other results, the spectra are similar to those for the isotropic test cases, with the peak wavenumber determined by the domain size and the dissipation cutoff wavenumber corresponding to the characteristic mesh size. As for the isotropic cases, a portion of the spectrum is apparent that approximates −5/3 inertial range scaling.

Normalized energy spectra for forcing simulation of homogeneous anisotropic turbulence with volume averaging, on 1283 grid, (a) zero shear stress (ZSS) at = 0.45 and 0.045, (b) ZSS and finite shear stress (FSS) results, at = 0.45 , and (c) ZSS and FSS results, at = 4.5
Table 3 summarizes the results in terms of the resolved turbulence statistics. The maximum difference between the target and resolved Reynolds stress components for each case is shown next to the far-right column. As the strength of the forcing term is increased by decreasing the value of the forcing time scale , the resolved stress components agree more closely with their target values. For the ZSS case at , the maximum difference is 2.18%. For both isotropic (Table 2) and anisotropic (Table 3) cases, the results indicate that selection of coefficients to provide a sufficiently large forcing term can successfully be used to reproduce the desired target statistics.
Summary of turbulence statistics for homogeneous anisotropic STF cases (green color implies target and over 90% of target statistics)
Target Reynolds stress | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|
0.3333 | 0.2400 | 0.6666 | 0.2933 | 0.9999 | 0.4133 | ||||||
Resolved Reynolds stress | |||||||||||
Averaging method | FSS | Grid | Max % difference | (k/k*) ratio % | |||||||
Volume | No (ZSS) | 0.045 | 1283 | 0.3354 | 0.6520 | 0.9849 | 2.18 | 98.61 | |||
Volume | No (ZSS) | 0.45 | 1283 | 0.3368 | 0.5214 | 0.8291 | 21.79 | 84.36 | |||
Volume | No (ZSS) | 4.5 | 1283 | 0.1964 | 0.2336 | 0.2993 | 70.07 | 36.46 | |||
Volume | Yes (FSS) | 0.45 | 1283 | 0.3940 | 0.1795 | 0.5618 | 0.1972 | 0.7906 | 0.2780 | 32.77 | 87.32 |
Volume | Yes (FSS) | 4.5 | 1283 | 0.0931 | 0.0220 | 0.1180 | 0.0334 | 0.1678 | 0.0405 | 90.84 | 18.95 |
Target Reynolds stress | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|
0.3333 | 0.2400 | 0.6666 | 0.2933 | 0.9999 | 0.4133 | ||||||
Resolved Reynolds stress | |||||||||||
Averaging method | FSS | Grid | Max % difference | (k/k*) ratio % | |||||||
Volume | No (ZSS) | 0.045 | 1283 | 0.3354 | 0.6520 | 0.9849 | 2.18 | 98.61 | |||
Volume | No (ZSS) | 0.45 | 1283 | 0.3368 | 0.5214 | 0.8291 | 21.79 | 84.36 | |||
Volume | No (ZSS) | 4.5 | 1283 | 0.1964 | 0.2336 | 0.2993 | 70.07 | 36.46 | |||
Volume | Yes (FSS) | 0.45 | 1283 | 0.3940 | 0.1795 | 0.5618 | 0.1972 | 0.7906 | 0.2780 | 32.77 | 87.32 |
Volume | Yes (FSS) | 4.5 | 1283 | 0.0931 | 0.0220 | 0.1180 | 0.0334 | 0.1678 | 0.0405 | 90.84 | 18.95 |
6 Summary and Conclusion
A new method for synthetic turbulence generation in scale-resolving turbulent flow CFD simulations is presented. The new method, denoted STF, incorporates a source term in the momentum equation to drive the local, instantaneous velocity vector toward a target value. The target value is computed at each instant and location in the simulation based on a mapping of the resolved first- and second-order statistics to the desired, target single-point statistics. The resolved statistics can be computed during the simulation using either volume averaging for homogeneous turbulence or time averaging for stationary turbulence.
The method was evaluated by performing several test simulations of homogeneous turbulence, using LES with the Smagorinsky subgrid stress model and or the MILES (implicit LES) modeling approach. Qualitative and quantitative results were presented which showed that the method was able to reproduce the target Reynolds stress components after a relatively short period of simulation time. The simulations using volume averaging responded more quickly than those using time averaging since there is no inherent lag in statistical calculations. It was demonstrated that the forcing time scale can be varied to control the strength of the forcing term, and reducing the time scale was found to be an effective means of reducing the time required for evolution to a quasi-stationary turbulent state, and for improving the accuracy of the forced simulation.
Energy spectra showed that the simulations reproduced the characteristic inertial range scaling when a sufficiently fine mesh was used, although energy pile up occurred near the cutoff wavenumber when an appropriate dissipative subgrid stress model was not employed. It was also demonstrated that spatial filtering can be used to control the effective large-eddy length scale. In sum, the results indicate that the STF method is capable of reproducing a synthetic homogeneous turbulence field with prescribed first- and second-order statistics and appropriate spectral content, which can be used to specify initial and/or boundary conditions for LES simulations. The method is relatively simple to implement, nonstochastic, stable, and computationally efficient. The STF method may therefore offer an attractive alternative for synthetic turbulence generation in three-dimensional Navier–Stokes CFD codes. Future work will investigate extension of the method for generation of inlet freestream and boundary layer conditions for LES of wall-bounded turbulent flows.
Acknowledgment
The computing for this project was performed at the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU).