**3.3. Model performance in atmospheric boundary layer flows**

The new algebraic nonlinear closure (Eqn. (9)) is assessed through a systematic comparison with well-established empirical formulations and theoretical predictions of a variety of flow statistics in a neutral atmospheric boundary layer (see 58 for more details) and a a stable boundary layer (see 11, 60 for more details). Obtained from different resolution simulations of the neutral ABL case, figure 2 shows the values of the nondimensional vertical gradients of the resolved streamwise velocity as a function of vertical position, Φ*<sup>M</sup>* = *<sup>κ</sup><sup>z</sup> u*∗ *∂*<*u*> *<sup>∂</sup><sup>z</sup>* , and the scalar counterpart, the values of the nondimensional vertical gradients of the mean resolved scalar concentration as a function of vertical position, Φ*<sup>θ</sup>* = *<sup>κ</sup><sup>z</sup> θ*∗ *∂*<*θ* > *<sup>∂</sup><sup>z</sup>* . On the basis of experimental results and dimensional analysis [e.g., 15, 88, 93], it has been found in neutral cases that Φ*<sup>M</sup>* = 1, and Φ*<sup>θ</sup>* = 0.74 for all *z* in the surface layer. The new closure yields the value of Φ*<sup>M</sup>* that remains close to 1 (indicative of the expected logarithmic velocity profile), and the value of Φ*<sup>θ</sup>* that remains close to 0.74. Further, figure 3 shows the nondimensional one-dimensional spectra of the simulated streamwise velocity and the nondimensional one-dimensional power spectra of the resolved scalar concentration obtained from the 1283 simulation using the new closure, computed at different heights. In the inertial subrange (*k*1*<sup>z</sup>* 1) all the normalized spectra are in good agreement with the −5/3 power-law scaling.

10.5772/57051

199

http://dx.doi.org/10.5772/57051

θ

**Figure 3.** Averaged non-dimensional 1-D spectra of (a) the streamwise velocity and (b) the resolved scalar concentration

Φ<sup>θ</sup>

**Figure 4.** Nondimensional (a) velocity gradient and (b) temperature gradient at different resolution simulations of a stable ABL case. The solid and dashed lines correspond to the formulations according to equations (12) and (13). Figure is modified from

, <sup>Φ</sup>*<sup>θ</sup>* <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>z</sup>*

where the coefficients are *a* = 1, *b* = 2/3, *c* = 5 and *d* = 0.35. These formulations are plotted along with the LES Φ*<sup>M</sup>* and Φ*<sup>θ</sup>* results as functions of *z*/*L* in figure 4. It is evident that all of our simulation results agree quite well with the empirical relations. In the bulk of the surface

*<sup>L</sup>* , <sup>Φ</sup>*<sup>θ</sup>* <sup>=</sup> 0.74 <sup>+</sup> 4.7 *<sup>z</sup>*

 *a* 1 + 2 3 *az <sup>L</sup>* <sup>+</sup> *be*<sup>−</sup> *dz L* 

*L*

*<sup>L</sup>* , (12)

<sup>1</sup> <sup>+</sup> *<sup>c</sup>* <sup>−</sup> *dz L*

, (13)

obtained from the 128<sup>3</sup> simulation of a neutral ABL case. Figure is modified from Lu and Porté-Agel [58, 60].

(a) (b)

(a) (b)

<sup>Φ</sup>*<sup>M</sup>* <sup>=</sup> <sup>1</sup> <sup>+</sup> 4.7 *<sup>z</sup>*

layer the results have better agreement with equation (12) than equation (13).

and nonlinear relations derived from Beljaars and Holtslag [12]

<sup>1</sup> <sup>+</sup> *<sup>c</sup>* <sup>−</sup> *dz L*

Φ

Lu and Porté-Agel [60].

<sup>Φ</sup>*<sup>M</sup>* <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>z</sup>*

*L*

*<sup>a</sup>* + *be*<sup>−</sup> *dz L*  θ

Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research

**Figure 2.** Nondimensional vertical gradients of (a) the mean resolved streamwise velocity and (b) the mean resolved scalar concentration obtained from the simulations of a neutral ABL case. Figure is modified from Lu and Porté-Agel [58, 60].

Stable boundary layers (SBLs) are relatively shallow and are characterized by strong vertical shear and a relatively high wind near the top of the boundary layer, thus are particularly challenging to simulate accurately. Φ*<sup>M</sup>* and Φ*<sup>θ</sup>* are key parameters for surface parameterizations in large-scale models and assessment of a SGS model. Owing to the existence of non-zero mean spanwise velocity component, the definition equation is modified as Φ*<sup>M</sup>* = *<sup>κ</sup><sup>z</sup> u*∗ *∂*<*u*> *∂z* 2 + *∂*<*v*> *∂z* 2 . In the surface layer, Φ*<sup>M</sup>* and Φ*<sup>θ</sup>* are usually parameterized as functions of *z*/*L*, where *L* is the Obukhov length. For instance, the well-known linear relations [15, 88]

<sup>198</sup> Computational and Numerical Simulations Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research 9 10.5772/57051 Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research http://dx.doi.org/10.5772/57051 199

8 Computational and Numerical Simulations

power-law scaling.

modified as Φ*<sup>M</sup>* = *<sup>κ</sup><sup>z</sup>*

*u*∗

well-known linear relations [15, 88]

*∂*<*u*> *∂z*

2 + *∂*<*v*> *∂z*

(a) (b)

**3.3. Model performance in atmospheric boundary layer flows**

The new algebraic nonlinear closure (Eqn. (9)) is assessed through a systematic comparison with well-established empirical formulations and theoretical predictions of a variety of flow statistics in a neutral atmospheric boundary layer (see 58 for more details) and a a stable boundary layer (see 11, 60 for more details). Obtained from different resolution simulations of the neutral ABL case, figure 2 shows the values of the nondimensional vertical gradients

the scalar counterpart, the values of the nondimensional vertical gradients of the mean

basis of experimental results and dimensional analysis [e.g., 15, 88, 93], it has been found in neutral cases that Φ*<sup>M</sup>* = 1, and Φ*<sup>θ</sup>* = 0.74 for all *z* in the surface layer. The new closure yields the value of Φ*<sup>M</sup>* that remains close to 1 (indicative of the expected logarithmic velocity profile), and the value of Φ*<sup>θ</sup>* that remains close to 0.74. Further, figure 3 shows the nondimensional one-dimensional spectra of the simulated streamwise velocity and the nondimensional one-dimensional power spectra of the resolved scalar concentration obtained from the 1283 simulation using the new closure, computed at different heights. In the inertial subrange (*k*1*<sup>z</sup>* 1) all the normalized spectra are in good agreement with the −5/3

**Figure 2.** Nondimensional vertical gradients of (a) the mean resolved streamwise velocity and (b) the mean resolved scalar concentration obtained from the simulations of a neutral ABL case. Figure is modified from Lu and Porté-Agel [58, 60].

Stable boundary layers (SBLs) are relatively shallow and are characterized by strong vertical shear and a relatively high wind near the top of the boundary layer, thus are particularly challenging to simulate accurately. Φ*<sup>M</sup>* and Φ*<sup>θ</sup>* are key parameters for surface parameterizations in large-scale models and assessment of a SGS model. Owing to the existence of non-zero mean spanwise velocity component, the definition equation is

2

parameterized as functions of *z*/*L*, where *L* is the Obukhov length. For instance, the

*u*∗ *∂*<*u*> *<sup>∂</sup><sup>z</sup>* , and

*<sup>∂</sup><sup>z</sup>* . On the

θ

*θ*∗ *∂*<*θ* >

θ

. In the surface layer, Φ*<sup>M</sup>* and Φ*<sup>θ</sup>* are usually

of the resolved streamwise velocity as a function of vertical position, Φ*<sup>M</sup>* = *<sup>κ</sup><sup>z</sup>*

resolved scalar concentration as a function of vertical position, Φ*<sup>θ</sup>* = *<sup>κ</sup><sup>z</sup>*

**Figure 3.** Averaged non-dimensional 1-D spectra of (a) the streamwise velocity and (b) the resolved scalar concentration obtained from the 128<sup>3</sup> simulation of a neutral ABL case. Figure is modified from Lu and Porté-Agel [58, 60].

**Figure 4.** Nondimensional (a) velocity gradient and (b) temperature gradient at different resolution simulations of a stable ABL case. The solid and dashed lines correspond to the formulations according to equations (12) and (13). Figure is modified from Lu and Porté-Agel [60].

$$
\Phi\_M = 1 + 4.7 \frac{z}{L}, \qquad \Phi\_\theta = 0.74 + 4.7 \frac{z}{L} \, , \tag{12}
$$

and nonlinear relations derived from Beljaars and Holtslag [12]

$$\Phi\_M = 1 + \frac{z}{L} \left( a + b e^{-\frac{dz}{L}} \left( 1 + c - \frac{dz}{L} \right) \right), \\ \Phi\_\theta = 1 + \frac{z}{L} \left( a \sqrt{1 + \frac{2}{3} \frac{az}{L}} + b e^{-\frac{dz}{L}} \left( 1 + c - \frac{dz}{L} \right) \right), \tag{13}$$

where the coefficients are *a* = 1, *b* = 2/3, *c* = 5 and *d* = 0.35. These formulations are plotted along with the LES Φ*<sup>M</sup>* and Φ*<sup>θ</sup>* results as functions of *z*/*L* in figure 4. It is evident that all of our simulation results agree quite well with the empirical relations. In the bulk of the surface layer the results have better agreement with equation (12) than equation (13).

10.5772/57051

201

http://dx.doi.org/10.5772/57051

**Figure 6.** Cyclonic two-dimensional coherent structures appearing in rotating turbulence as indicated by iso-surfaces of vorticity, planar contours of kinetic energy and velocity vectors: (a) statistically steady state of small scale forced rotating

Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research

**Figure 7.** Assessment of the 643 LES of large scale forced rotating turbulence at statistically steady state: (a) 3D & 2D kinetic energy spectra obtained from the new nonlinear mixed model; (b) PDF of *ω*<sup>3</sup> obtained from the filtered DNS (gray line), the Smagorinsky model (△), the mixed scale-similarity model () and the new model (). Figure is modified from Lu et al. [62].

In order to examine the cyclonic/anti-cyclonic asymmetry, vortex regions rather than the whole domain are concentrated. The vortex regions can be identified using the criterion *λ*<sup>2</sup> < 0<sup>1</sup> for a vortex region [35], and points with *λ*<sup>2</sup> < (1/6) min(*λ*2) < 0 are sampled to obtain the probability density function *PDF*(*ω*3). Figure 7(b) compares the Smagorinsky model, the mixed scale-similarity model and the new model with respect to the (PDF) of *ω*3.

<sup>2</sup> (*∂ui*/*∂xj* <sup>−</sup> *<sup>∂</sup>uj*/*∂xi*)

case. (b) statistically steady state of large scale forced rotating case. Figure is adapted from Lu et al. [61].

(a) (b)

<sup>1</sup> *λ*<sup>2</sup> is the second large eigenvalue of tensor *SimSmj* + Ω*im*Ω*mj*, where Ω*ij* = <sup>1</sup>

**Figure 5.** Nondimensional skewness of vertical velocity in an unstable ABL flow.

In unstable ABL flows, the land surface is warmer than the surrounding air, in response to solar heating. The warmer land surface leads to a positive (upward) heat flux, which creates a thermal instability and generates turbulence. Vertical-velocity skewness in an unstable ABL flow is indicative of the structure of the motion, when a positive value means that updrafts are narrower than surrounding downdrafts. A number of puzzling features of the vertical-velocity skewness are found in observations and LES of the unstable ABL. Figure 5 includes the vertical-velocity skewness derived from observations (the Minnesota experiment, 98, and AMTEX, 51) and from LESs. There is fairly good agreement over the lower portion of the convective ABL; however in the upper portion, previous LESs [63, 70, 71] indicate a further increase of the vertical-velocity skewness while the observations show a nearly constant value. The new closure, evidently, predicts much more accurate vertical-velocity skewness.

### **3.4. Model performance in rotating and isotropic turbulences**

It has been noted that rotating turbulence with a moderate Rossby number below an *O*(1) critical value accompanies asymmetry in favor of large scale cyclonic vortical columns [e.g. 83, 84]. This cyclone/anti-cyclone asymmetry is widely observable in atmospheric science. The direction of rotation of large scale wind flow (e.g., hurricanes & typhoons) is counterclockwise in the northern hemispheres (clockwise in the south). It is interesting to note that small scale rotating turbulent flows in the atmosphere, in which the direct influence of planetary Coriolis effect is inconsequential, are usually triggered by large scale cyclonical storms, and approximately more than 90% of tornadoes rotate in a cyclonic direction. Figure 6 shows that rotating turbulence captures this asymmetry, and a quasi two-dimensional flow - in other words, reduced variations along the rotation axis [e.g. 17, 18, 61, 62, 83].

Figure 7(a) shows the 3D & 2D kinetic energy spectra in a statistically steady state obtained from the new mixed nonlinear model (Eqn. (7)) at the resolution of 643. The flow is forced at large scales, and *Roω*<sup>3</sup> = 0.12. The new model facilitates the two-dimensionalization process resulting in very close 3D & 2D energy spectra at large scales, *k* <∼ 7. All other examined models fail to deliver this quasi 2D structure at large scales.

<sup>200</sup> Computational and Numerical Simulations Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research 11 10.5772/57051 Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research http://dx.doi.org/10.5772/57051 201

10 Computational and Numerical Simulations

vertical-velocity skewness.

**Figure 5.** Nondimensional skewness of vertical velocity in an unstable ABL flow.

**3.4. Model performance in rotating and isotropic turbulences**

models fail to deliver this quasi 2D structure at large scales.

In unstable ABL flows, the land surface is warmer than the surrounding air, in response to solar heating. The warmer land surface leads to a positive (upward) heat flux, which creates a thermal instability and generates turbulence. Vertical-velocity skewness in an unstable ABL flow is indicative of the structure of the motion, when a positive value means that updrafts are narrower than surrounding downdrafts. A number of puzzling features of the vertical-velocity skewness are found in observations and LES of the unstable ABL. Figure 5 includes the vertical-velocity skewness derived from observations (the Minnesota experiment, 98, and AMTEX, 51) and from LESs. There is fairly good agreement over the lower portion of the convective ABL; however in the upper portion, previous LESs [63, 70, 71] indicate a further increase of the vertical-velocity skewness while the observations show a nearly constant value. The new closure, evidently, predicts much more accurate

It has been noted that rotating turbulence with a moderate Rossby number below an *O*(1) critical value accompanies asymmetry in favor of large scale cyclonic vortical columns [e.g. 83, 84]. This cyclone/anti-cyclone asymmetry is widely observable in atmospheric science. The direction of rotation of large scale wind flow (e.g., hurricanes & typhoons) is counterclockwise in the northern hemispheres (clockwise in the south). It is interesting to note that small scale rotating turbulent flows in the atmosphere, in which the direct influence of planetary Coriolis effect is inconsequential, are usually triggered by large scale cyclonical storms, and approximately more than 90% of tornadoes rotate in a cyclonic direction. Figure 6 shows that rotating turbulence captures this asymmetry, and a quasi two-dimensional flow

Figure 7(a) shows the 3D & 2D kinetic energy spectra in a statistically steady state obtained from the new mixed nonlinear model (Eqn. (7)) at the resolution of 643. The flow is forced at large scales, and *Roω*<sup>3</sup> = 0.12. The new model facilitates the two-dimensionalization process resulting in very close 3D & 2D energy spectra at large scales, *k* <∼ 7. All other examined


**Figure 6.** Cyclonic two-dimensional coherent structures appearing in rotating turbulence as indicated by iso-surfaces of vorticity, planar contours of kinetic energy and velocity vectors: (a) statistically steady state of small scale forced rotating case. (b) statistically steady state of large scale forced rotating case. Figure is adapted from Lu et al. [61].

**Figure 7.** Assessment of the 643 LES of large scale forced rotating turbulence at statistically steady state: (a) 3D & 2D kinetic energy spectra obtained from the new nonlinear mixed model; (b) PDF of *ω*<sup>3</sup> obtained from the filtered DNS (gray line), the Smagorinsky model (△), the mixed scale-similarity model () and the new model (). Figure is modified from Lu et al. [62].

In order to examine the cyclonic/anti-cyclonic asymmetry, vortex regions rather than the whole domain are concentrated. The vortex regions can be identified using the criterion *λ*<sup>2</sup> < 0<sup>1</sup> for a vortex region [35], and points with *λ*<sup>2</sup> < (1/6) min(*λ*2) < 0 are sampled to obtain the probability density function *PDF*(*ω*3). Figure 7(b) compares the Smagorinsky model, the mixed scale-similarity model and the new model with respect to the (PDF) of *ω*3.

<sup>1</sup> *λ*<sup>2</sup> is the second large eigenvalue of tensor *SimSmj* + Ω*im*Ω*mj*, where Ω*ij* = <sup>1</sup> <sup>2</sup> (*∂ui*/*∂xj* <sup>−</sup> *<sup>∂</sup>uj*/*∂xi*)

The new model successfully delivers the cyclonic/anti-cyclonic asymmetry (the dominance of positive *ω*<sup>3</sup> vortices); others fail to do so.

10.5772/57051

203

http://dx.doi.org/10.5772/57051

(a) (b)

(a) (b)

the standard dynamic mixed nonlinear model.

energy at small scales.

model, and the � represents the results from the new model. Figure is modified from Lu [57].

**Figure 9.** (a) Evolution of resolved kinetic energy; (b) comparison of energy spectra at *t* = 0.4 (normalized using initial eddy turnover time). The gray line represents the results from the filtered DNS, the × represents the results from simulation without a model, the △ represents the results from the Smagorinsky model, the ▽ represents the results from the standard gradient

Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research

For further testing of the new algebraic model, two high-Reynolds-number cases were conducted. Regarding the first case, the original simulation was performed at the Johns Hopkins University using 1024 grid points in each direction [52]. The database contains a 1024<sup>4</sup> space-time history of an incompressible isotropic turbulent flow in 3D. The initial condition for decaying LES runs, downloaded from "turbulence.pha.jhu.edu" at time equals 2 without space interpolation, bears *Re<sup>λ</sup>* = 430. A direct numerical simulation of decay was lacking; thus comparisons could be performed only against statistical theories of turbulence. Figure 10(a) shows the kinetic energy spectra obtained from the 128<sup>3</sup> LES at different times. They follow the −5/3 power-law behavior until the dissipation range starts to be captured through LES at a late period, and importantly, there are no improper accumulations of kinetic

**Figure 10.** Energy spectra obtained from the LES (a) starting from high-Reynolds-number DNS data; Figure is modified from [57] (b) starting from high-Reynolds-number wind-tunnel measured spectrum, and the gray lines represent the results from measurement, the solid lines represent the results from the new model, and the dash-dotted lines represent the results from

Figure 8 examines the 64<sup>3</sup> intermediate scale forced large eddy simulations with *Roω*<sup>3</sup> = 0.16. As shown in figure 8(a), the new model dissipates the kinetic energy at small scales and captures the reverse energy transfer to large scales. Figure 8(b) compares SGS models with respect to *PDF*(*ω*3), and shows that the new model is capable of delivering the cyclonic/anti-cyclonic asymmetry.

**Figure 8.** Assessment of the 643 LES of small scale forced rotating turbulence at statistically steady state: (a) 3D & 2D kinetic energy spectra obtained from the new nonlinear mixed model; (b) PDF of *ω*<sup>3</sup> obtained from the filtered DNS (gray line), the Smagorinsky model (△), the mixed scale-similarity model () and the new model (). Figure is modified from Lu et al. [62].

A direct numerical simulation case has been simulated at the University of Wisconsin-Madison by means of decaying, and has been adopted in model assessment studies [57, 61, 62]. The initial Taylor micro-scale Reynolds number (*Reλ*) is approximately 85. Thus, the 128<sup>3</sup> DNS has resolved the flow of all scales, and the DNS results can be used to verify the accuracy of SGS models and identify their problems. Figure 9(a) shows the evolution of the resolved kinetic energy (normalized by its initial value) obtained from DNS and LESs. The decaying case starts at *Re<sup>λ</sup>* = 85, and is beyond the capability of 64<sup>3</sup> simulation (typically *Re<sup>λ</sup>* ∼ 50). The total resolved kinetic energy obtained from the 64<sup>3</sup> simulation without a model, which simply omits *τij*, yields the worst prediction. In order to obtain proper kinetic energy decay rates, turbulence modeling is needed for coarser grids. It is clear that the results obtained using the new algebraic nonlinear model are in the excellent agreement with the filtered DNS results.

The Smagorinsky model yields a higher kinetic energy decay rate than that is found in the filtered DNS results. Figure 9(b) shows that using the Smagorinsky model, kinetic energy at small scales is dissipated excessively. By construction, the standard gradient model allows for energy "backscatter." However, over a period of time, it yields a lower kinetic energy decay rate, and kinetic energy at small scales is accumulating (cannot be dissipated effectively). As its revision, the new algebraic model has shown significant improvement in energy-spectrum accuracy.

<sup>202</sup> Computational and Numerical Simulations Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research 13 10.5772/57051 Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research http://dx.doi.org/10.5772/57051 203

12 Computational and Numerical Simulations

cyclonic/anti-cyclonic asymmetry.

of positive *ω*<sup>3</sup> vortices); others fail to do so.

(a) (b)

agreement with the filtered DNS results.

accuracy.

The new model successfully delivers the cyclonic/anti-cyclonic asymmetry (the dominance

Figure 8 examines the 64<sup>3</sup> intermediate scale forced large eddy simulations with *Roω*<sup>3</sup> = 0.16. As shown in figure 8(a), the new model dissipates the kinetic energy at small scales and captures the reverse energy transfer to large scales. Figure 8(b) compares SGS models with respect to *PDF*(*ω*3), and shows that the new model is capable of delivering the

**Figure 8.** Assessment of the 643 LES of small scale forced rotating turbulence at statistically steady state: (a) 3D & 2D kinetic energy spectra obtained from the new nonlinear mixed model; (b) PDF of *ω*<sup>3</sup> obtained from the filtered DNS (gray line), the Smagorinsky model (△), the mixed scale-similarity model () and the new model (). Figure is modified from Lu et al. [62].

A direct numerical simulation case has been simulated at the University of Wisconsin-Madison by means of decaying, and has been adopted in model assessment studies [57, 61, 62]. The initial Taylor micro-scale Reynolds number (*Reλ*) is approximately 85. Thus, the 128<sup>3</sup> DNS has resolved the flow of all scales, and the DNS results can be used to verify the accuracy of SGS models and identify their problems. Figure 9(a) shows the evolution of the resolved kinetic energy (normalized by its initial value) obtained from DNS and LESs. The decaying case starts at *Re<sup>λ</sup>* = 85, and is beyond the capability of 64<sup>3</sup> simulation (typically *Re<sup>λ</sup>* ∼ 50). The total resolved kinetic energy obtained from the 64<sup>3</sup> simulation without a model, which simply omits *τij*, yields the worst prediction. In order to obtain proper kinetic energy decay rates, turbulence modeling is needed for coarser grids. It is clear that the results obtained using the new algebraic nonlinear model are in the excellent

The Smagorinsky model yields a higher kinetic energy decay rate than that is found in the filtered DNS results. Figure 9(b) shows that using the Smagorinsky model, kinetic energy at small scales is dissipated excessively. By construction, the standard gradient model allows for energy "backscatter." However, over a period of time, it yields a lower kinetic energy decay rate, and kinetic energy at small scales is accumulating (cannot be dissipated effectively). As its revision, the new algebraic model has shown significant improvement in energy-spectrum

**Figure 9.** (a) Evolution of resolved kinetic energy; (b) comparison of energy spectra at *t* = 0.4 (normalized using initial eddy turnover time). The gray line represents the results from the filtered DNS, the × represents the results from simulation without a model, the △ represents the results from the Smagorinsky model, the ▽ represents the results from the standard gradient model, and the � represents the results from the new model. Figure is modified from Lu [57].

For further testing of the new algebraic model, two high-Reynolds-number cases were conducted. Regarding the first case, the original simulation was performed at the Johns Hopkins University using 1024 grid points in each direction [52]. The database contains a 1024<sup>4</sup> space-time history of an incompressible isotropic turbulent flow in 3D. The initial condition for decaying LES runs, downloaded from "turbulence.pha.jhu.edu" at time equals 2 without space interpolation, bears *Re<sup>λ</sup>* = 430. A direct numerical simulation of decay was lacking; thus comparisons could be performed only against statistical theories of turbulence. Figure 10(a) shows the kinetic energy spectra obtained from the 128<sup>3</sup> LES at different times. They follow the −5/3 power-law behavior until the dissipation range starts to be captured through LES at a late period, and importantly, there are no improper accumulations of kinetic energy at small scales.

**Figure 10.** Energy spectra obtained from the LES (a) starting from high-Reynolds-number DNS data; Figure is modified from [57] (b) starting from high-Reynolds-number wind-tunnel measured spectrum, and the gray lines represent the results from measurement, the solid lines represent the results from the new model, and the dash-dotted lines represent the results from the standard dynamic mixed nonlinear model.

The second high-Reynolds-number is based on measurements of nearly isotropic turbulence downstream of an active grid [41]. The initial Taylor micro-scale Reynolds number (*Reλ*) is approximately 720. The new model is implemented in LES of decaying isotropic turbulence with initial conditions that match the measured energy spectra at *x*1/*M* = 20. Figure 10(b) shows energy spectra at different times, for comparison, the results obtained using the standard dynamic nonlinear mixed model [41] were also included in the plot. The new model clearly gives more accurate results at small scales that *<sup>k</sup>* > 6 m<sup>−</sup>1.

10.5772/57051

205

http://dx.doi.org/10.5772/57051

*relc*(*CL***e***<sup>L</sup>* <sup>+</sup> *CD***e***D*) , (14)

Figure 11(b) shows one discretized element of a blade, and figure 11(c) shows a cross-sectional element at radius *r* defining the airfoil in the (*θ*, *x*) plane, where *x* is the axial direction. With the tangential and axial velocities of the incident flow denoted as *V<sup>θ</sup>* and *Va*, respectively, the local velocity relative to the rotating blade is given as **<sup>V</sup>***rel* = (*V<sup>θ</sup>* − <sup>Ω</sup>*r*, *Va*). The angle of attack is defined as *<sup>α</sup>* <sup>=</sup> *<sup>ϕ</sup>* <sup>−</sup> *<sup>γ</sup>*, where *<sup>ϕ</sup>* <sup>=</sup> *tan*−1(*Va*/(Ω*<sup>r</sup>* <sup>−</sup> *<sup>V</sup><sup>θ</sup>* )) is the angle between **V***rel* and the rotor plane, and *γ* is the local pitch angle. The turbine-induced force

Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research

where *CL* = *CL*(*α*, *Re*) and *CD* = *CD*(*α*, *Re*) are the lift coefficient and the drag coefficient, respectively, *ρ* is the air density, *c* is the chord length, and **e***<sup>L</sup>* and **e***<sup>D</sup>* denote the unit vectors in the directions of the lift and the drag, respectively. The flow of interest is essentially inviscid, and viscous effects from the boundary layer of blade are introduced only as integrated quantities through the use of airfoil data. The airfoil data and subsequent loading are determined by computing local angles of attack from the movement of the blades and the local flow field. The model enables us to study in detail the dynamics of the wake and the tip vortices and their influence on the induced velocities in the rotor plane. The applied blade forces are distributed smoothly to avoid singular behavior and numerical instability. In practice, the aerodynamic blade forces are distributed along and away from the actuator lines in a three-dimensional Gaussian manner through the convolution of the computed local

*<sup>ǫ</sup>*3*π*3/2 exp

where *d* is the distance between grid points and points at the actuator line, and *ǫ* is a parameter that serves to adjust the concentration of the regularized load. The value of this regularization parameter, *ǫ*, is typically on the order of 1 − 3 grid sizes and should be as small as possible so that the turbine-induced force is distributed over an area representing the chord distribution and does not affect the wake structure. However, small values will deliver significant numerical oscillations; thus, as a trade-off between numerical stability and

The main advantage of representing the blades by airfoil data is that many fewer grid points are needed to capture the influence of the blades than would be needed for simulating the actual geometry of the blades. Therefore, the ALM is well suited for wake studies since grid points can be concentrated in a larger part of the wake while keeping the computing costs at a reasonable level. Moreover, the ALM encompasses blade motions as well as their mixing mechanism, which is crucial for simulating more realistic wind-turbine wakes and achieving improved predictions with respect to the standard actuator disk model [7, 78, 97]. However, the ALM does not resolve detailed flows on blade surfaces and relies on tabulated two-dimensional airfoil data for providing lift and drag coefficients. This makes it dependent on both the quality of airfoil data and the method used to model the influence of dynamically

−*d*<sup>2</sup> *ǫ*2 , (15)

per radial unit length is given by the following equation

load, **f**, and a regularization kernel *ηǫ* as shown below

changing angles of attack and stall.

**<sup>f</sup>** <sup>=</sup> **dF**

*dr* <sup>=</sup> <sup>1</sup> 2 *ρV*<sup>2</sup>

**<sup>f</sup>***<sup>ǫ</sup>* <sup>=</sup> **<sup>f</sup>** <sup>⊗</sup> *ηǫ* , *ηǫ* <sup>=</sup> <sup>1</sup>

accuracy [34, 90], it is set to be 1.5 times of the grid size in the rotor plane.
