## 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 [35]. 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.

Conceptually, the STF method proposed here seeks to induce a synthetic turbulence field through the addition of a time-dependent, nonstochastic forcing term in the momentum equation. The forcing term is constructed to drive the instantaneous, local velocity toward a time-dependent target velocity that satisfies the user-specified first- and second-order one-point statistics for the turbulence, i.e., the mean velocity vector and Reynolds stress tensor. To introduce the method we consider the general form of the continuity, momentum, and energy equations for single-phase, single-species, compressible flow:
$∂ρ∂t+∂∂xj(ρuj) =0$
(1)
$∂∂t(ρui)+∂∂xj(ρuiuj) =−∂p∂xi+∂∂xj(σij)$
(2)
$∂∂t(ρE)+∂∂xj(ρujH)=∂∂xj(qj+uiσij)$
(3)
In the above, the flow variables $ρ$,$u$,$p$, $E$, and $H$ may represent local instantaneous (DNS) or filtered (LES) quantities. Likewise, the viscous stress tensor $σij$ and heat flux vector $qj$ include both molecular and, for LES, subfilter contributions. The STF method is implemented by adding a forcing term, $fi$, to the momentum and energy equations
$∂∂t(ρui)+∂∂xj(ρuiuj) =−∂p∂xi+∂∂xj(σij)+fi$
(4)
$∂∂t(ρE)+∂∂xj(ρujH)=∂∂xj(qj+uiσij)+ujfj$
(5)
The source term is constructed such that during each time-step of a simulation, the resolved velocity vector is forced toward a target velocity vector that would yield a desired target statistical distribution for the time-varying velocity field. The general form of the forcing term is
$fi= ρτf(ui*−ui)$
(6)
Here, $ui*$is the target local, instantaneous velocity, and $τf$ is a characteristic time scale for the forcing term. Inputs to the model include prescription of a local target mean velocity, $u¯i*$, and turbulent stress tensor
$ui′uj′¯*=(uiuj¯−u¯iu¯j)*$
(7)

where the overbar denotes either Reynolds or Favre (mass-weighted) averaging.

The key aspect of the method is the calculation of the target velocity vector$ui*$. It is first noted that the transformation proposed by Lund et al. [2] can be used to map an ensemble of isotropic velocity fluctuations $v′$ to an ensemble of fluctuations that satisfy a target statistical distribution $Tij=ui′uj′¯*$ as follows:
$ui′*=Bijvj′$
(8)
$Bij= [T1100T21/B11T22−B2120T31/B11(T32−B21B31)/B22T33−B312−B322]$
(9)
Similarly, an ensemble of resolved fluctuations satisfying a particular statistical distribution $Rij=ui′uj′¯$ can be mapped to an isotropic distribution $v′$ using the inverse of the Lund coefficient matrix
$vi′=Aij−1uj′$
(10)
$Aij−1= [1/A1100−A21/(A11A22)1/A220(A21A32−A31A22)/(A11A22A33)−A32/(A22A33)1/A33]$
(11)
$Aij= [R1100R21/A11R22−A2120R31/A11(R32−A21A31)/A22R33−A312−A322]$
(12)
It is therefore possible to define a mapping from a distribution of resolved velocity fluctuations $ui′$ with known statistical second moment tensor (turbulent stress) $Rij$ to a distribution $ui′*$ with target turbulent stress $Tij$ as
$ui′*=Cijuj′$
(13)
$Cij= BikAkj−1$
(14)
The instantaneous target velocity used in the forcing function (Eq. (6)) includes the target fluctuating velocity as well as the target mean velocity
$ui*=u¯i* +Cijuj′$
(15)
where the fluctuating velocity is defined relative to the mean:
$ui′= ui−u¯i$
(16)

In practice, the method is implemented as follows: First, a target statistical velocity distribution is specified prior to the simulation in terms of $u¯i*$ and$ui′uj′¯*$. As the simulation proceeds, the resolved statistics $u¯i$ and $ui′uj′¯$ are obtained using an appropriate averaging technique. At each time-step, the transformation tensor $Cij$ is computed at each point in the domain based on $ui′uj′¯*$ and $ui′uj′¯$. During each iteration, the fluctuating velocity $ui′$ is computed, and the target instantaneous velocity $ui*$ is found using Eq. (15). The forcing term $fi$ is then computed using Eq. (6) and included as an additional source term in the momentum and energy equations.

### 2.2 Ensemble Averaging.

In theory, any appropriate approximation to the Reynolds-averaging operation can be used in the CFD simulations. In this study, two different methods are implemented and tested for simulation of statistically stationary, homogeneous turbulence. The first is volume averaging, for which the Reynolds-averaged value of any arbitrary variable $φ$ is defined as
$φ¯(t)=1V∫φ(xi,t) dV$
(17)
and the integral is performed over the entire simulation domain with volume $V$. The second is time averaging, defined as
$φ¯(xi,t)=1t∫0tφ(xi,τ) dτ$
(18)
where $t$ is the physical simulation time. In practice this is achieved using a discrete running time average, for which the averaged value at each point in the domain can be determined by
$φ¯(xi,t)=n−1nφ¯(xi,t−Δt)+1nφ(xi,t)$
(19)

Here, $Δt$ is the time-step size and $n$ is the current number of time steps in the simulation. In the limit $n→∞$, the averaged value becomes constant for statistically stationary flow.

### 2.3 Spatial Filtering.

Spatial filtering is implemented in the STF method to allow the user to have some measure of control of the turbulent length scale. All simulations adopting spatial filtering in this study were performed using a second-order differential elliptic filter [18]. In this method, for an isotropic filter, the filtered resolved velocity ($ûi$) is obtained by solution of Eq. (20), where $Δ$ is the filter width or size
$∂∂xj(Δ2∂ûi∂xj)=ûi− ui$
(20)
Spatial filtering is implemented in the STF method by first redefining fluctuating velocity ($ui′)$ in Eq. (16) as filtered fluctuating velocity ($ui′′)$
$ui′′=(ui−ûi)−(u¯i−ûi¯)$
(21)
The instantaneous target velocity used in the forcing function includes the target filtered fluctuating velocity ($Cijuj′′)$ as well as the target mean velocity $(u¯i*)$
$ui*=u¯i*+(ûi−u¯i)+Cijuj′′$
(22)
At every iteration, the filtered transformation tensor $Cij$ is computed at each point in the domain based on a modified target turbulent stress, $Tij$, and resolved turbulent stress, $Rij$, where $Mij$ is a term that represents the contribution to the turbulent stress by the interaction between the filtered and instantaneous velocity
$Tij=ui′uj′¯*−Mij$
(23)
$Rij=ui′uj′¯−Mij$
(24)
$Mij=12[(ûiuj¯+uiûj¯)−(ûi¯u¯j+u¯iûj¯)]$
(25)
Equation (20) represents an isotropic spatial filtering operation. In practice, the STF method adopts a generalized anisotropic filter defined by
$∂∂xj(τT2 uj′uk′¯*∂ûi∂xk)=ûi− ui$
(26)
Anisotropic filtering is a relatively simple way to incorporate the fact that turbulence can have different length scales in different directions. The method assumes that regardless of their length scale, large eddies (i.e., velocity fluctuations) will share the same time scale. The tensorial filter width is defined to depend on a characteristic turbulent time scale,$τT$, and the target turbulent stress tensor,$ui′uj′¯*$. For homogeneous anisotropic turbulence, spatial variations and turbulent fluctuations are not statistically uniform in all directions and the filtered resolved velocity ($ûi$) is obtained by solution of Eq. (26). For isotropic turbulence,the anisotropic filtering operation is formally similar to the isotropic filter defined in Eq. (20). For this case, spatial variations and turbulent fluctuations are statistically uniform in all directions, hence the anisotropic filter is equivalent to a scalar isotropic filter for which the filter size is defined as
$Δ2=13τT2 (uk′uk′¯*)$
(27)

### 2.4 Prescription of the Forcing Time Scale.

The forcing time scale $τf$ is arbitrary. In principle, a smaller value will increase the magnitude of the forcing term and drive the velocity more rapidly toward its target value, but too small a value may result in instability or may constrain the flow from developing naturally once resolved turbulence has been introduced. A user can select a relevant time scale based on the known flow physics of the problem under investigation, numerical and stability considerations, and/or trial and error. It would be advantageous, however, to incorporate a universal time scale that takes into account both the physical and numerical aspects of the simulation and simplifies the user input requirements. It is proposed, therefore, that an appropriate time scale is of the form
$τf ∼ τT$
(28)

where $τT$ 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 $k/ε$ can be used to specify $τT$.

For the simulations presented in this paper that use spatial filtering, the form of the forcing term is
$fi=ρCfτT(ui*−ui)$
(29)
Turbulent time scale, $τT$, is specified by the user, and the effect that different values have on the simulation is investigated. The coefficient $Cf$ dictates the overall strength of the forcing term. For cases with spatial filtering, $τT$ and $Cf$ must be independently selected. For cases without spatial filtering, it is only the ratio $Cf/τT$ that impacts the simulation by controlling the strength of the forcing term by dictating the effective forcing time scale $τf$.

## 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.

The majority of the simulations presented here use the Smagorinsky eddy-viscosity-based subgrid stress model [21], for which the deviatoric part of the subgrid stress tensor is expressed as
$τijSGS=2νTSij$
(30)
The eddy viscosity is formulated as
$νT=(Cs Δg)22SijSij$
(31)

where $Δg$ is the characteristic mesh size, equal to the cube root of cell volume in the current simulations, and the coefficient $Cs$= 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, $τT$, and forcing coefficient, $Cf$

• 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.

Fig. 1
Fig. 1
Close modal

## 5 Results and Discussion

### 5.1 Homogeneous Isotropic Turbulence.

Homogeneous turbulence simulations were first performed for an isotropic target field corresponding to
$u¯1*=u¯2*=u¯3*=0$
(32)
$u1′u1′¯*=u2′u2′¯*=u3′u3′¯*=v′2$
(33a)
$u1′u2′¯*=u1′u3′¯*=u2′u3′¯*=0$
(33b)
The target turbulent kinetic energy is therefore $k*=32v′2$. The characteristic velocity scale for the target flow field is $v′$, and a characteristic simulation time scale is defined as
$τs= L/v′$
(34)

#### 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 $τf/τs$ = 0.45, 0.045, 0.0045, and 0.00033. Forcing was applied during the time interval$0. 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 $t=2τs$, at which point the fluctuating velocity field had reached an apparently stationary state. The results are shown for a forcing time scale $τf/τs$ = 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.

Fig. 2
Fig. 2
Close modal

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 $t=2τs$. 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.

Fig. 3
Fig. 3
Close modal

Figure 6 shows the time history of the ratio of resolved turbulent kinetic energy,$k=12ui′uj′¯$, to the target value, $k*$. 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 $t=3τs$.

Fig. 4
Fig. 4
Close modal

Figure 7 highlights the resolved turbulence level reached in the simulations during the forcing time interval, expressed as the ratio $k/k*$, 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, $τf$, 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 $τf=0.45$, the ratio increases from 0.84 to 0.87 as the mesh is refined from 323 to 1923 cells.

Fig. 5
Fig. 5
Close modal

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 $κ≤EI$, the dissipation range is indicated as$κ≥DI$, and the inertial range is indicated as $EI≤κ≤DI.$ 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.

Fig. 6
Fig. 6
Close modal

For the 1283 grid results shown in Fig. 8(b), energy spectra obtained using different forcing time scales ($τf=0.45, 0.045, 0.0045, and 0.00033$) are all in good quantitative agreement, with only slight differences evident in the low wavenumbers ($κ≤EI$). 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 $τf$. 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.

Fig. 7
Fig. 7
Close modal

#### 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$τf=$ 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 $τf$ 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 $τf$ 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, $τs$. 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.

Fig. 8
Fig. 8
Close modal
Fig. 9
Fig. 9
Close modal
Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal
Fig. 12
Fig. 12
Close modal
Table 1

Percentage ratio of resolved-to-target turbulent kinetic energy for different averaging techniques and forcing time scales at $t=3τs$

τfs = 4.5τfs = 0.45τfs = 0.045
TA57%74%102%
VA37%88%99%
τfs = 4.5τfs = 0.45τfs = 0.045
TA57%74%102%
VA37%88%99%
a

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 $Δ =v′τT$. 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$τT/ τs=0.45,$ and$Cf=100$. This corresponds to a forcing time scale of $τf=0.0045$, and for the spatially filtered case corresponds to a filter size $ΔL=0.45$. 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.

Fig. 13
Fig. 13
Close modal
Fig. 14
Fig. 14
Close modal

Figure 15 shows the normalized energy spectra for the simulations with and without spatial filtering, for $τT/ τs=0.45$ and $τf=0.0045$. The effect of spatial filtering is evident, with the maximum energy density occurring at $κ≈4/L$, which corresponds to $κ≈1.8/Δ$. The remainder of the spectrum, including the inertial and dissipation ranges, is qualitatively similar for the two cases.

Fig. 15
Fig. 15
Close modal

To further highlight the effect of spatial filtering, simulations were run using four different values of the characteristic large-eddy timescale ($τT$), and a constant value of the forcing coefficient $Cf=10$. 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 $Cf$ is held constant, as $τT$ is reduced the effective forcing timescale $τf$ is correspondingly reduced and the magnitude of the forcing term increases. As a consequence, as $τT$ 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 $Cf$. For comparison purposes, the simulation with $τTτs=0.45$ and $Cf=100$ 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 $t/τs=3$, the turbulent energy decay rate increases (i.e., turbulence decays faster) as $τT$ is reduced.

Fig. 16
Fig. 16
Close modal
Fig. 17
Fig. 17
Close modal

Figure 18(a) shows normalized energy spectra for the cases with varying characteristic large eddy time scale ($τT$), and a constant value of the forcing coefficient $Cf=10$. The effect of spatial filtering is apparent, as the peak energy density moves toward higher wavenumbers as $τT$ 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 $τT$, there is an evident inertial range that approximately follows Kolmogorov −5/3 scaling.

Fig. 18
Fig. 18
Close modal

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 ($τf/τs$) 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.

Table 2

Summary of turbulence statistics for homogeneous isotropic STF cases (green color implies target and over 90% of target statistics)

Target Reynolds stress
$u1′u1′¯/k*$$u2′u2′¯/k*$$u3′u3′¯/k*$
0.66660.66660.6666
Averaging methodSpatially filtered$τT/τs$$τf/τs$Grid$u1′u1′¯/k*$$u2′u2′¯/k*$$u3′u3′¯/k*$Max % difference(k/k*) ratio %
VolumeNo0.0453230.66010.65600.66651.5999.13
VolumeNo0.453230.57250.54330.568418.4984.21
VolumeNo0.0456430.65620.65220.64952.5697.90
VolumeNo0.456430.57250.55050.570117.4184.66
VolumeNo0.04512830.65670.65320.65482.0198.23
VolumeNo0.4512830.58110.54330.577618.4985.10
VolumeNo4.512830.22110.23590.215167.7333.60
VolumeYes0.450.04512830.62110.62030.61917.1393.03
VolumeYes0.450.4512830.48810.50090.499626.7874.43
VolumeYes0.454.512830.20980.22380.204169.3931.88
VolumeNo0.4519230.57780.58870.580013.3387.32
VolumeNo0.0003312830.66650.66650.66650.0199.98
VolumeNo0.003312830.66550.66580.66490.2699.81
VolumeNo0.004512830.66510.66510.66550.2399.78
VolumeNo0.04512830.65630.65150.65392.2798.08
VolumeNo0.06712830.63790.65010.64564.3196.68
VolumeYes0.450.4512830.48810.50090.499626.7874.43
VolumeYes0.2250.22512830.48750.49470.491326.8773.68
VolumeYes0.1120.11212830.49730.49690.498525.4574.64
VolumeYes0.0560.05612830.49120.49090.490226.4673.62
TimeYes0.450.04512830.61280.59680.588011.8089.88
TimeNo0.04512830.65440.63780.63634.5496.42
TimeNo0.4512830.54640.51080.492326.1477.47
TimeNo4.512830.34670.28710.316156.9347.50
Target Reynolds stress
$u1′u1′¯/k*$$u2′u2′¯/k*$$u3′u3′¯/k*$
0.66660.66660.6666
Averaging methodSpatially filtered$τT/τs$$τf/τs$Grid$u1′u1′¯/k*$$u2′u2′¯/k*$$u3′u3′¯/k*$Max % difference(k/k*) ratio %
VolumeNo0.0453230.66010.65600.66651.5999.13
VolumeNo0.453230.57250.54330.568418.4984.21
VolumeNo0.0456430.65620.65220.64952.5697.90
VolumeNo0.456430.57250.55050.570117.4184.66
VolumeNo0.04512830.65670.65320.65482.0198.23
VolumeNo0.4512830.58110.54330.577618.4985.10
VolumeNo4.512830.22110.23590.215167.7333.60
VolumeYes0.450.04512830.62110.62030.61917.1393.03
VolumeYes0.450.4512830.48810.50090.499626.7874.43
VolumeYes0.454.512830.20980.22380.204169.3931.88
VolumeNo0.4519230.57780.58870.580013.3387.32
VolumeNo0.0003312830.66650.66650.66650.0199.98
VolumeNo0.003312830.66550.66580.66490.2699.81
VolumeNo0.004512830.66510.66510.66550.2399.78
VolumeNo0.04512830.65630.65150.65392.2798.08
VolumeNo0.06712830.63790.65010.64564.3196.68
VolumeYes0.450.4512830.48810.50090.499626.7874.43
VolumeYes0.2250.22512830.48750.49470.491326.8773.68
VolumeYes0.1120.11212830.49730.49690.498525.4574.64
VolumeYes0.0560.05612830.49120.49090.490226.4673.62
TimeYes0.450.04512830.61280.59680.588011.8089.88
TimeNo0.04512830.65440.63780.63634.5496.42
TimeNo0.4512830.54640.51080.492326.1477.47
TimeNo4.512830.34670.28710.316156.9347.50

### 5.2 Homogeneous Anisotropic Turbulence.

In principle, there is no restriction of the STF method to an isotropic target Reynolds stress. To demonstrate the ability of the method to synthetically generate an arbitrary anisotropic turbulent velocity field, simulations were run on the 1283 cell grid with $τf/τs$ = 4.5, 0.45, and 0.045. Results were obtained using volume averaging and no spatial filtering. A first case was run for unequal normal stresses but zero shear stress (denoted ZSS). The target statistics for the ZSS case were prescribed as
$u¯1*=u¯2*=u¯3*=0$
(35)
$u1′u1′¯*=0.5 k*, u2′u2′¯*=k*, u3′u3′¯*=1.5 k*$
(36)
$u1′u2′¯*=u1′u3′¯*=u2′u3′¯*=0$
(37)
An additional case was run with the same mesh and simulation conditions, but with finite shear stress components (denoted FSS). The target statistics for this case were prescribed as
$u¯1*=u¯2*=u¯3*=0$
(38)
$u1′u1′¯*=0.5 k*, u2′u2′¯*=k*, u3′u3′¯*=1.5 k*$
(39)
$u1′u2′¯*=0.36 k*, u1′u3′¯*=0.44 k*, u2′u3′¯*=0.62 k*$
(40)

Figure 19 shows the contours of instantaneous velocity magnitude for both anisotropic test cases with forcing time scale $τf/τs$ = 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.

Fig. 19
Fig. 19
Close modal

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 $τf$, the resolved stress components agree more closely with their target values. For the ZSS case at $τf=0.045$, 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.

Table 3

Summary of turbulence statistics for homogeneous anisotropic STF cases (green color implies target and over 90% of target statistics)

Target Reynolds stress
$u1′u1′¯/k*$$u1′u2′¯/k*$$u2′u2′¯/k*$$u1′u3′¯/k*$$u3′u3′¯/k*$$u2′u3′¯/k*$
0.33330.24000.66660.29330.99990.4133
Resolved Reynolds stress
Averaging methodFSS$τf/τs$Grid$u1′u1′¯/k*$$u1′u2′¯/k*$$u2′u2′¯/k*$$u1′u3′¯/k*$$u3′u3′¯/k*$$u2′u3′¯/k*$Max % difference(k/k*) ratio %
VolumeNo (ZSS)0.04512830.33540.65200.98492.1898.61
VolumeNo (ZSS)0.4512830.33680.52140.829121.7984.36
VolumeNo (ZSS)4.512830.19640.23360.299370.0736.46
VolumeYes (FSS)0.4512830.39400.17950.56180.19720.79060.278032.7787.32
VolumeYes (FSS)4.512830.09310.02200.11800.03340.16780.040590.8418.95
Target Reynolds stress
$u1′u1′¯/k*$$u1′u2′¯/k*$$u2′u2′¯/k*$$u1′u3′¯/k*$$u3′u3′¯/k*$$u2′u3′¯/k*$
0.33330.24000.66660.29330.99990.4133
Resolved Reynolds stress
Averaging methodFSS$τf/τs$Grid$u1′u1′¯/k*$$u1′u2′¯/k*$$u2′u2′¯/k*$$u1′u3′¯/k*$$u3′u3′¯/k*$$u2′u3′¯/k*$Max % difference(k/k*) ratio %
VolumeNo (ZSS)0.04512830.33540.65200.98492.1898.61
VolumeNo (ZSS)0.4512830.33680.52140.829121.7984.36
VolumeNo (ZSS)4.512830.19640.23360.299370.0736.46
VolumeYes (FSS)0.4512830.39400.17950.56180.19720.79060.278032.7787.32
VolumeYes (FSS)4.512830.09310.02200.11800.03340.16780.040590.8418.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).

## References

1.
Spalart
,
P. R.
, and
Watmuff
,
J. H.
,
1993
, “
Experimental and Numerical Study of a Turbulent Boundary Layer With Pressure Gradients
,”
J. Fluid Mech.
,
249
(
1
), pp.
337
371
.10.1017/S002211209300120X
2.
Lund
,
T. S.
,
Wu
,
X.
, and
Squires
,
K. D.
,
1998
, “
Generation of Turbulent Inflow Data for Spatially-Developing Boundary Layer Simulations
,”
J. Comp. Phys.
,
140
(
2
), pp.
233
258
.10.1006/jcph.1998.5882
3.
Spalart
,
P. R.
,
Strelets
,
M.
, and
Travin
,
A.
,
2006
, “
Direct Numerical Simulation of Large-Eddy-Break-Up Devices in a Boundary Layer
,”
Int. J. Heat Fluid Flow
,
27
(
5
), pp.
902
910
.10.1016/j.ijheatfluidflow.2006.03.014
4.
Shur
,
M. L.
,
Spalart
,
P. R.
,
Strelets
,
M. K.
, and
Travin
,
A. K.
,
2011
, “
A Rapid and Accurate Switch From RANS to LES in Boundary Layers Using an Overlap Region
,”
Flow Turbul. Combust.
,
86
(
2
), pp.
179
206
.10.1007/s10494-010-9309-9
5.
Araya
,
G.
,
Castillo
,
L.
,
Meneveau
,
C.
, and
Jansen
,
K.
,
2011
, “
A Dynamic Multi-Scale Approach for Turbulent Inflow Boundary Conditions in Spatially Developing Flows
,”
J. Fluid Mech.
,
670
, pp.
581
605
.10.1017/S0022112010005616
6.
Schlüter
,
J. U.
,
Pitsch
,
H.
, and
Moin
,
P.
,
2005
, “
Outflow Conditions for Integrated Large Eddy Simulation/Reynolds-Averaged Navier–Stokes Simulations
,”
AIAA J.
,
43
(
1
), pp.
156
164
.10.2514/1.11007
7.
Kraichnan
,
R. H.
,
1970
, “
Diffusion by a Random Velocity Field
,”
Phys. Fluids
,
13
(
1
), pp.
22
31
.10.1063/1.1692799
8.
Rogallo
,
R.
,
1981
, “
Numerical Experiments in Homogeneous Turbulence
,” NASA, NASA Ames Research Center Moffett Field, CA, Report No. 81315.
9.
Orszag
,
S.
, and
Patterson
,
G.
,
1972
, “
Numerical Simulation of Three-Dimensional Homogeneous Isotropic Turbulence
,”
Phys. Rev. Lett.
,
28
(
2
), pp.
76
79
.10.1103/PhysRevLett.28.76
10.
Lee
,
S.
,
Lele
,
S.
, and
Moin
,
P.
,
1992
, “
Simulation of Spatially Evolving Compressible Turbulence and the Application of Taylor's Hypothesis
,”
Phys. Fluids
,
4
(
7
), pp.
1521
1530
.10.1063/1.858425
11.
Lundgren
,
T. S.
,
2003
, “
Linearly Forced Isotropic Turbulence
,”
Annual Research Briefs
,
Center for Turbulence Research
,
Stanford, CA
, p.
461
.
12.
Rosales
,
C.
, and
Meneveau
,
C.
,
2005
, “
Linear Forcing in Numerical Simulations of Isotropic Turbulence: Physical Space Implementations and Convergence Properties
,”
Phys. Fluids
,
17
(
9
), p.
095106
.10.1063/1.2047568
13.
Jarrin
,
N.
,
Prosser
,
R.
,
Uribe
,
J. C.
,
,
S.
, and
Laurence
,
D.
,
2006
, “
A Synthetic-Eddy-Method for Generating Inflow Conditions for Large Eddy Simulations
,”
Int. J. Heat Fluid Flow
,
27
(
4
), pp.
585
593
.10.1016/j.ijheatfluidflow.2006.02.006
14.
Jarrin
,
N.
,
Prosser
,
R.
,
Uribe
,
J. C.
,
,
S.
, and
Laurence
,
D.
,
2009
, “
Reconstruction of Turbulent Fluctuations for Hybrid RANS/LES Simulations Using a Synthetic-Eddy Method
,”
Int. J. Heat Fluid Flow
,
30
(
3
), pp.
435
442
.10.1016/j.ijheatfluidflow.2009.02.016
15.
Keating
,
A.
,
Piomelli
,
U.
,
Balaras
,
E.
, and
Kaltenbach
,
H.-J.
,
2004
, “
A Priori and a Posteriori Tests of Inflow Conditions for Large Eddy Simulation
,”
Phys. Fluids (1994-Present)
,
16
(
12
), pp.
4696
4712
.10.1063/1.1811672
16.
Spille-Kohoff
,
A.
, and
Kaltenbach
,
H.-J.
,
2001
, “
Generation of Turbulent Inflow Data With a Prescribed Shear-Stress Profile
,”
Proceedings of the Third AFOSR International Conference on DNS/LES, DNS/LES Progress and Challenges
,
C.
,
Liu
,
L.
Sakell
, and
T.
Beutner
, eds.,
Greyden Press
,
Columbus
, Aug. 5–8, pp.
319
326.
17.
De Laage de Meux
,
B.
,
Audebert
,
B.
,
Manceau
,
R.
, and
Perrin
,
R.
,
2015
, “
Anisotropic Linear Forcing for Synthetic Turbulence Generation in Large Eddy Simulation and Hybrid RANS/LES Modeling
,”
Phys. Fluids
,
27
(
3
), p.
035115
.10.1063/1.4916019
18.
Luke
,
E. A.
,
Tong
,
X.
, and
Chamberlain
,
R.
,
2016
, “
FlowPsi: An Ideal Gas Flow Solver-The User Guide
,” sourceforge.net, Starkville, pp.
29
40, accessed Feb. 2, 2022,
https://github.com/libm3l/FlowPsi
19.
Subbareddy
,
P. K.
, and
Candler
,
G. V.
,
2009
, “
A Fully Discrete, Kinetic Energy Consistent Finite Volume Scheme for Compressible Flows
,”
J. Comp. Phys.
,
228
(
5
), pp.
1347
1364
.10.1016/j.jcp.2008.10.026
20.
Germano
,
M.
,
1986
, “
Differential Filters for the Large Eddy Numerical Simulation of Turbulent Flows
,”
Phys. Fluids
,
29
(
6
), pp.
1755
1766
.10.1063/1.865649
21.
Smagorinsky
,
J.
,
1963
, “
General Circulation Experiments With the Primitive Equations
,”
Mon. Weath. Rev.
,
91
(
3
), pp.
99
164
.10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2
22.
Fureby
,
C.
, and
Grinstein
,
F. F.
,
1999
, “
Monotonically Integrated Large Eddy Simulation of Free Shear Flows
,”
AIAA J.
,
37
(
5
), pp.
544
556
.10.2514/2.772