On the suitability of second-order accurate finite-volume solvers for the simulation of atmospheric boundary layer flow

The present work analyzes the quality and reliability of an important class of general-purpose, second-order accurate finite-volume (FV) solvers for the large-eddy simulation of a neutrally stratified atmospheric boundary layer (ABL) flow. The analysis is carried out within the OpenFOAM^® framework, which is based on a colocated grid arrangement. A series of open-channel flow simulations are carried out using a static Smagorinsky model for subgrid scale momentum fluxes in combination with an algebraic equilibrium wall-layer model. The sensitivity of the solution to variations in numerical parameters such as grid resolution (up to 160³ control volumes), numerical solvers, and interpolation schemes for the discretization of nonlinear terms is evaluated and results are contrasted against those from a well-established mixed pseudospectral–finite-difference code. Considered flow statistics include mean streamwise velocity, resolved Reynolds stresses, velocity skewness and kurtosis, velocity spectra, and two-point autocorrelations. A quadrant analysis along with the examination of the conditionally averaged flow field are performed to investigate the mechanisms responsible for momentum transfer in the flow. It is found that at the selected grid resolutions, the considered class of FV-based solvers yields a poorly correlated flow field and is not able to accurately capture the dominant mechanisms responsible for momentum transport in the ABL. Specifically, the predicted flow field lacks the well-known sweep and ejection pairs organized side by side along the cross-stream direction, which are representative of a streamwise roll mode. This is especially true when using linear interpolation schemes for the discretization of nonlinear terms. This shortcoming leads to a misprediction of flow statistics that are relevant for ABL flow applications and to an enhanced sensitivity of the solution to variations in grid resolution, thus calling for future research aimed at reducing the impact of modeling and discretization errors.

1 Introduction

An accurate prediction of atmospheric boundary layer (ABL) flows is of paramount importance across a wide range of fields and applications, including weather forecasting, complex terrain meteorology, agriculture, air quality modeling, and wind energy (Whiteman, 2000; Fernando, 2010; Calaf et al., 2010; Oke et al., 2017; Shaw et al., 2019).

Since the early work of Deardorff (1970), the large-eddy simulation (LES) technique has spurred considerable insight on the fundamental dynamics of ABL flow over rough surfaces (Anderson and Meneveau, 2010; Salesky et al., 2017; Momen et al., 2018), over and within plant and urban canopies (Yue et al., 2007b; Bailey and Stoll, 2013; Pan et al., 2014; Tseng et al., 2006; Bou-Zeid et al., 2009; Giometto et al., 2017; Li and Bou-Zeid, 2019), and for wind energy applications (Calaf et al., 2010; Abkar and Porté-Agel, 2013; Stevens and Meneveau, 2017), amongst others.

The majority of the past work has relied on fully or partially dealiased mixed pseudospectral–finite-difference (PSFD) solvers – the go-to approach for LES studies since the works of Moin et al. (1978) and Moeng (1984). Such solvers are known to yield accurate flow fields up to the LES cutoff frequency and to produce good results when used in conjunction with dynamic subgrid scale (SGS) models (Germano et al., 1991; Lilly, 1992). However, single domain PSFD-based solvers are limited to regular domains, are not suitable for the simulation of nonperiodic flows, have sharp variations in the flow field such as shocks or fluid–solid interfaces in boundary layer flows, and are typically difficult to parallelize owing to the global support of their spatial representation (see, e.g., Margairaz et al., 2018). With the increasing need to account for complex geometries and multiphysics, several efforts have been devoted to the mitigation of the aforementioned limitations (Fang et al., 2011; Li et al., 2016; Chester et al., 2007). However, the solutions are often ad hoc or validated only for specific applications, thus introducing a degree of uncertainty in model results that is hard to quantify and generalize.

There is hence a growing interest from the ABL community in LES solvers based on compact spatial schemes via structured or unstructured meshes (Orlandi, 2000; Ferziger and Peric, 2002). The parallelized large-eddy simulation model (Raasch and Schröter, 2001; Maronga et al., 2015) and the weather research and forecasting model (Skamarock and Klemp, 2008; Powers et al., 2017) are prominent examples of said efforts. Both the approaches are based on a high-order finite-difference discretization, with nonlinear terms approximated by using high-order upwind biased differencing schemes. The latter are suitable for LES in complex geometries with arbitrary grid stretching factors and outflow boundary conditions (Beaudan and Moin, 1994; Mittal and Moin, 1997) but are dissipative and do not strictly conserve energy. On the other hand, if central schemes are used instead for the evaluation of nonlinear terms, no numerical dissipation is introduced, but truncation errors can have an overwhelming impact on the computed flow field (Ghosal, 1996; Kravchenko and Moin, 1997). These limitations typically result in a strong sensitivity of the solution to properties of the spatial discretization and numerical scheme (Meyers et al., 2006, 2007; Meyers and Sagaut, 2007; Vuorinen et al., 2014; Rezaeiravesh and Liefvendahl, 2018; Breuer, 1998; Montecchia et al., 2019). Further, truncation errors corrupt the high wavenumber range of the solution, restricting the ability to adopt dynamic LES closure models that make use of information from the smallest resolved scales of motion to evaluate the SGS diffusion (Germano et al., 1991). Notwithstanding these limitations, central schemes have been heavily employed in the past in both the geophysical and engineering flow communities and are the de facto standard in the wind engineering community, where most of the numerical simulations are carried out using second-order accurate finite-volume (FV)-based solvers (Stovall et al., 2010; Churchfield et al., 2010; Balogh et al., 2012; Churchfield et al., 2013; Shi and Yeo, 2016, 2017; García-Sánchez et al., 2017; García-Sánchez and Gorlé, 2018).

Motivated by the aforementioned needs, the present study aims at characterizing the quality and reliability of an important class of second-order accurate FV solvers for the LES of neutrally stratified ABL flows. The analysis is conducted in the open-channel flow setup (no Coriolis acceleration) via the OpenFOAM^® framework (Weller et al., 1998; De Villiers, 2006; Jasak et al., 2007). A suite of simulations is carried out varying physical and numerical parameters, including grid resolution (up to 160³ control volumes), the numerical solver, and interpolation schemes for the discretization of the nonlinear term. Predictions from the FV solvers are contrasted against the results from the Albertson and Parlange (1999) PSFD code in terms of flow statistics, including mean streamwise velocity, resolved Reynolds stresses, two-point velocity autocorrelations, and mechanisms supporting momentum transport. The end goal is to provide a more nuanced understanding of the capabilities of general-purpose, second-order, FV-based solvers in predicting ABL flow.

The work is organized as follows. Section 2 summarizes the setup of the problem, the simulation database, and the postprocessing procedure. Results are shown in Sect. 3 and conclusions are drawn in Sect. 4. A further discussion on the sensitivity of the solution to model constants, interpolation schemes, and numerical solvers is provided in the Appendix.

2 Methodology

2.1 Governing equations and numerical schemes

We use index notation in a Cartesian reference system. The spatially filtered Navier–Stokes equations are considered,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial u_i}{\partial t}}+u_j{\displaystyle \frac{\partial u_i}{\partial x_j}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{p}}{\partial x_i}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}^{\mathrm{SGS},\mathrm{dev}}}{\partial x_j}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial P}{\partial x_i}},\text{(2)}& {\displaystyle \frac{\partial u_i}{\partial x_i}}=\mathrm{0},\end{array}$$

where $u_i=(u,v,w)$ is the spatially filtered velocity field along the streamwise (x), cross-stream (y), and vertical (z) coordinate directions, respectively, t is the time, ρ is the constant fluid density (Boussinesq approximation), $\stackrel{\mathrm{̃}}{p}\equiv p+\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\tau }}_{kk}^{\mathrm{SGS}}$ is a modified pressure term, τ_ij is the filtered viscous stress tensor, and ${\mathit{\tau }}_{ij}^{\mathrm{SGS},\mathrm{dev}}$ is the deviatoric part of the SGS stress tensor. In addition, the term $-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial P}{\partial x_i}$ is an imposed constant pressure gradient driving the flow. The spatially filtered viscous tensor is ${\mathit{\tau }}_{ij}=-\mathrm{2}\mathit{\nu }{S}_{ij}$, where ν=const is the kinematic viscosity of the Newtonian fluid and S_ij is the resolved (in the LES sense) rate of strain tensor. For the SGS stress tensor, the static Smagorinsky model is used,

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{ij}^{\mathrm{SGS},\mathrm{dev}}=-\mathrm{2}{\mathit{\nu }}^{\mathrm{SGS}}{S}_{ij}=-\mathrm{2}\left(C_{\mathrm{S}}\mathrm{\Delta }{)}^{\mathrm{2}}\right|S|S_{ij},\end{array}$$

where ν^SGS is the SGS eddy viscosity, C_S is the Smagorinsky coefficient (Smagorinsky, 1963), $\mathrm{\Delta }=(\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z{)}^{\mathrm{1}/\mathrm{3}}$ is a local length scale based on the volume of the computational cell (Scotti et al., 1993), and $\left|S\right|=\sqrt{\mathrm{2}{S}_{ij}{S}_{ij}}$ quantifies the magnitude of the rate of strain. In the present work, C_S=0.1, unless otherwise specified. Note that dynamic Smagorinsky models are preferred to the static one for the LES of ABL flows (Germano et al., 1991; Lilly, 1992; Meneveau and Lund, 1997; Porté-Agel, 2004; Bou-Zeid et al., 2005). Dynamic models evaluate SGS stresses via first-principles-based constraints, feature improved dissipation properties when compared to the static Smagorinsky model (especially in the vicinity of solid boundaries), and are free of explicit modeling parameters. The choice made in the present study is motivated by problematics encountered when using the available dynamic Lagrangian model in preliminary tests. However, while SGS dissipation plays a crucial role in PSFD solvers, truncation errors may overshadow SGS stress contributions in the second-order FV-based ones (Kravchenko and Moin, 1997). The static Smagorinsky SGS model used herein might hence perform similarly to dynamic SGS models for the considered flow setup. This conjecture is supported by the results of Majander and Siikonen (2002).

The large scale separation between near-surface and outer-layer energy-containing ABL motions poses stringent resolution requirements to numerical modelers, if all the energy containing motions have to be resolved. To reduce the computational cost of such simulations, the near-surface region is typically bypassed and a phenomenological wall-layer model is leveraged instead to account for the impact of near-wall (inner-layer) dynamics on the outer-layer flow (Piomelli, 2008; Bose and Park, 2018). This approach is referred to as wall-modeled large-eddy simulation (WMLES) and is used herein. An algebraic wall-layer model for surfaces in a fully rough aerodynamic regime was implemented based on the logarithmic equilibrium assumption, i.e.,

$$\begin{array}{}\text{(4)}& \left|\stackrel{\mathrm{̃}}{\mathit{u}}\right|={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right),\end{array}$$

where $\left|\stackrel{\mathrm{̃}}{\mathit{u}}\right|\equiv \sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}$ is the norm of the velocity at a certain distance from the ground level, u_* is the friction velocity (see Sect. 2.2 for details), κ is the von Kármán constant, z is the distance from the ground level, and z₀ is the so-called aerodynamic roughness length, a length scale used to quantify the drag of the underlying surface. In this work, the values κ=0.41 and z₀=0.1 m are set. The kinematic wall shear stress is assumed to be proportional to the local velocity gradient (Boussinesq hypothesis),

$$\begin{array}{}\text{(5)}& {\mathit{\tau }}_{iz,\mathrm{w}}=(\mathit{\nu }+{\mathit{\nu }}_{\mathrm{t}}){\displaystyle \frac{\partial u_i}{\partial z}}{|}_{\mathrm{w}},i=x,y,\end{array}$$

where ν_t is the total eddy viscosity. Employing the no-slip condition for the velocity field, the standard FV approximation of the shear stress at the wall gives (Mukha et al., 2019)

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{iz,\mathrm{w}}=(\mathit{\nu }+{\mathit{\nu }}_{\mathrm{t}}{)}_{\mathrm{f}}{\displaystyle \frac{u_{i,\mathrm{c}}}{\mathrm{\Delta }z}},i=x,y,\end{array}$$

where the subscript f is used to denote the evaluation at the center of the wall face, the subscript c denotes the evaluation at the center of the wall-adjacent cell, and Δz is the distance from the wall. From the logarithmic law (Eq. 4) evaluated at the first cell center, one can write $u_{*}=\mathit{\kappa }|\stackrel{\mathrm{̃}}{\mathit{u}}{|}_{\mathrm{c}}/\ln (\frac{\mathrm{\Delta }z}{z_{\mathrm{0}}})$. Using the definition of friction velocity $u_{*}=\sqrt{{\mathit{\tau }}_{\mathrm{w}}^{\mathrm{2}}}$, where τ_w is the magnitude of the kinematic wall shear stress vector, along with Eq. (5), and rearranging, the total eddy viscosity at the wall can be written as

$$\begin{array}{}\text{(7)}& {\mathit{\nu }}_{\mathrm{t},\mathrm{f}}={\left({\displaystyle \frac{\mathit{\kappa }|\stackrel{\mathrm{̃}}{\mathit{u}}{|}_{\mathrm{c}}}{{\displaystyle }\ln \left({\displaystyle \frac{\mathrm{\Delta }z}{z_{\mathrm{0}}}}\right)}}\right)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{\Delta }z}{|\stackrel{\mathrm{̃}}{\mathit{u}}{|}_{\mathrm{c}}}}-\mathit{\nu },\end{array}$$

which is the formulation implemented herein. Note that $\mathit{\nu }+{\mathit{\nu }}_{\mathrm{t}}\approx {\mathit{\nu }}_{\mathrm{t}}$ in the boundary layer flows in the fully rough aerodynamic regime, so ν could be neglected without loss of accuracy.

In the present work, the computational grid is colocated, being the only colocated grid arrangement available within the OpenFOAM^® framework. Note that although advantageous in complex domains when compared to staggered grids (Ferziger and Peric, 2002), the colocated arrangement is known to cause difficulties with pressure–velocity coupling, hence requiring specific procedures to avoid oscillations in the solution. OpenFOAM^® offers the standard Rhie–Chow correction (Rhie and Chow, 1983), which is known to negatively affect the energy-conservation properties of central schemes (Ferziger and Peric, 2002). In addition, when approximating the integrals over the surfaces bounding each control volume (as a consequence of the Gauss divergence theorem), the unknowns are evaluated at face centers and are assumed to be constant at each face, yielding an overall second-order spatial accuracy (Churchfield et al., 2010). Since the divergence form of the convective term is used in combination with a low-order scheme over a nonstaggered grid, the solution is inherently unstable (Kravchenko and Moin, 1997). The present work makes use of the linear and QUICK interpolation schemes (Ferziger and Peric, 2002) to evaluate the unknowns at face centers (more details are provided in Sect. 2.2). The numerical solver is based on the PISO algorithm (Issa, 1985) for the pressure–velocity calculation and on an implicit Adams–Moulton scheme for time integration (Ferziger and Peric, 2002). In Appendix A2, the performances of an alternative solver with a Runge–Kutta time-advancement scheme and a projection method for the pressure–velocity coupling (Vuorinen et al., 2014) are analyzed.

2.2 Problem setup

A series of WMLES of ABL flow (open-channel flow setup) is performed. Tests are carried out in the domain $[\mathrm{0},L_x]\times [\mathrm{0},L_y]\times [\mathrm{0},L_z]$ with L_x=2πh, $L_y=\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }h$, L_z=h, where h=1000 m denotes the width of the open channel. Symmetry is imposed at the top of the computational domain, no-slip applies at the lower surface, and periodic boundary conditions are enforced along each side. A kinematic pressure gradient term $-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial P}{\partial x}=\mathrm{1}\mathrm{m}/{\mathrm{s}}^{\mathrm{2}}$ drives the flow along the x coordinate direction, yielding $u_{*}=\mathrm{1}\mathrm{m}/\mathrm{s}$. The kinematic viscosity is set to a nominal value of ${\mathrm{10}}^{-\mathrm{7}}{\mathrm{m}}^{\mathrm{2}}/\mathrm{s}$, which results in an essentially inviscid flow.

The computational mesh is Cartesian, with a uniform stencil along each direction. Three simulations are run over 64³, 128³, and 160³ control volumes, with the linear interpolation scheme for the evaluation of the unknowns at the face centers (simulations FV64, FV128, and FV160, respectively). Three additional simulations are run, at the same grid resolutions, with the linear scheme for the approximation of every term except for the nonlinear one, for which the QUICK scheme is used instead (simulations FV64*, FV128*, and FV160*). The cases span different grid resolutions at the same aspect ratio $\mathrm{\Delta }x/\mathrm{\Delta }z=\mathrm{2}\mathit{\pi }$. Note that the chosen grid resolutions are in line with those typically used in studies of ABL flow with the pseudospectral approach (see, e.g., Salesky et al., 2017). All the calculations satisfy the Courant–Friedrichs–Lewy (CFL) condition C≲0.1, where C is the Courant number. Runs are initialized from a fully developed open-channel flow simulation in statistically steady state (dynamic equilibrium), and time integration is carried out for 100 eddy turnover times, where the eddy turnover time is defined as $h/u_{*}$. Flow statistics are the result of an averaging procedure over the horizontal plane of statistical homogeneity of turbulence (xy) and in time over the last 60 eddy turnover times. The procedure yields well-converged statistics throughout the considered cases. In the following, the horizontal and temporal averaging operation is denoted by 〈⋅〉. The results from the present study are contrasted against the corresponding ones from the Albertson and Parlange (1999) mixed PSFD code (simulations PSFD64, PSFD128, and PSFD160). The code is based on an explicit second-order accurate Adams–Bashforth scheme for time integration and on a fractional-step method for solving the system of equations. Simulations from the PSFD solver are carried out using a static Smagorinsky SGS model with C_S=0.1, a rough wall-layer model with z₀=0.1 m, and C≲0.1. A summary of the runs is given in Table 1 along with the acronyms used in this study.

Table 1 Tabulated list of cases.

3 Results

This section is devoted to the analysis of velocity central moments (Sect. 3.1), spectra and spatial autocorrelations (Sect. 3.2), and momentum transfer mechanisms (Sect. 3.3).

Figure 1 Vertical structure of mean streamwise velocity $\langle u\rangle /u_{*}$ (a), streamwise velocity RMS $u_{\mathrm{RMS}}^{\prime }/u_{*}$ (b), cross-stream velocity RMS $v_{\mathrm{RMS}}^{\prime }/u_{*}$ (c), and vertical velocity RMS $w_{\mathrm{RMS}}^{\prime }/u_{*}$ (d). The red line in (a) denotes the reference logarithmic profile and the red line in (b) is a reference profile from Hultmark et al. (2013).

3.1 Mean velocity, Reynolds stresses, and higher-order statistics

Figure 1 shows first- and second-order statistics for all the considered cases. The mean streamwise velocity is shown in Fig. 1a in a comparison with the phenomenological logarithmic-layer profile. The velocity at the first two cell centers off the wall is consistently underpredicted, whereas a positive log-layer mismatch (LLM) is observed in the bulk of the flow (Kawai and Larsson, 2012). The LLM is particularly pronounced for the cases using the QUICK interpolation scheme. This behavior could have been anticipated, as the wall shear stress is evaluated using the instantaneous horizontal velocity at the first cell center off the wall. A number of procedures has been proposed to alleviate the LLM, including modifying the SGS stress model in the near-wall region (Sullivan et al., 1994; Porté-Agel et al., 2000; Chow et al., 2005; Wu and Meyers, 2013), shifting the matching location further away from the wall (Kawai and Larsson, 2012), and carrying out a local horizontal/temporal filtering operation (Bou-Zeid et al., 2005; Xiang et al., 2017). In preliminary runs, the approach of Kawai and Larsson (2012) was implemented in an attempt to alleviate the LLM. However, no apparent improvement was observed and the solution became very sensitive to grid resolution and matching location. This finding suggests that alternative procedures might need to be devised to overcome the LLM in ABL flow simulations when using the considered class of FV solvers. Note that profiles from the PSFD solver also feature a positive LLM in spite of a spatial, low-pass filtering operation that is carried out on the horizontal velocity field before evaluating the surface shear stress (Bou-Zeid et al., 2005).

The vertical structure of turbulence intensities is also shown in Fig. 1, where $(\cdot {)}_{\mathrm{RMS}}^{\prime }$ denotes the root mean square (RMS) of the fluctuations. Profiles from the FV-based solver start off relatively slow at the wall when compared to those from the PSFD-based solver and to the reference profile from Hultmark et al. (2013). This behavior is due to a combination of SGS and discretization errors, which damp the energy of high-wavenumber modes and whose accurate quantification remains an open challenge in LES (see, e.g., Meyers et al., 2006; Meyers and Sagaut, 2007; Meyers et al., 2007). Further aloft, $u_{\mathrm{RMS}}^{\prime }$ ($w_{\mathrm{RMS}}^{\prime }$) features relatively stronger (weaker) peak values when compared to the corresponding PSFD profile, the overprediction (underprediction) being more apparent in the simulations with the QUICK scheme. The overshoot in the peak of $u_{\mathrm{RMS}}^{\prime }$ is a well-known problem of FV-based WMLES (Bae et al., 2018). Lack of energy redistribution via pressure fluctuation from shear generated $u_{\mathrm{RMS}}^{\prime }$ to $v_{\mathrm{RMS}}^{\prime }$ and $w_{\mathrm{RMS}}^{\prime }$ is the root cause of said behavior, and possible mitigation strategies include allowing for wall transpiration (Bose and Moin, 2014). Grid refinement shifts the velocity RMS peaks closer to the surface and increases the magnitude of the velocity RMS therein, but leads to no improvement in the $max\left(u_{\mathrm{RMS}}^{\prime }\right)$ and only marginally improves the estimation of the $max\left(w_{\mathrm{RMS}}^{\prime }\right)$. A quantitative measure of the relative error on $u_{\mathrm{RMS}}^{\prime }$ with respect to the reference profile $u_{\mathrm{RMS},\mathrm{ref}}^{\prime }$ from Hultmark et al. (2013) in the $z_{\mathrm{0}}/h\le z/h\le \mathrm{0.4}$ interval is shown in Table 2. The FV-based solver performs worse than the PSFD-based one and the convergence is not monotonic. Note that nonmonotonic convergence is relatively common in LES at relatively coarse resolutions and is due to the interaction between discretization and modeling errors, whose impact on the solution cannot be a priori quantified (Meyers et al., 2007).

Table 2 Relative error on the turbulence intensities $\left|\right|u_{\mathrm{RMS}}^{\prime }-u_{\mathrm{RMS},\mathrm{ref}}^{\prime }|{|}_{L^{\mathrm{2}}}/|\left|u_{\mathrm{RMS},\mathrm{ref}}^{\prime }\right|{|}_{L^{\mathrm{2}}}$ w.r.t. the reference profile from Hultmark et al. (2013) in the interval $z_{\mathrm{0}}/h\le z/h\le \mathrm{0.4}$.

Skewness and kurtosis of the streamwise velocity (S_uu and K_uu, respectively) are shown in Fig. 2. The profiles of S_uu obtained with the FV-based solver and the QUICK scheme as well as those obtained with the PSFD-based solver are in good agreement with experimental results from Monty et al. (2009), here taken as a reference. On the contrary, the FV-based solver overpredicts S_uu when the linear interpolation scheme is used, with the skewness remaining positive throughout the whole extent of the surface layer. Note that a positive skewness of streamwise velocity represents a flow field where negative fluctuations are more likely to happen than the corresponding positive ones. The kurtosis obtained with the FV-based solver is consistently overpredicted, representing a flow field populated by a greater number of extreme events. Again, profiles from all cases feature a nonmonotonic convergence to the reference ones, as shown in Table 3, where the relative error on skewness and kurtosis with respect to the measurements from Monty et al. (2009) is reported in the interval $z_{\mathrm{0}}/h\le z/h\le \mathrm{0.4}$.

Figure 2 Vertical structure of skewness of streamwise velocity (a) and kurtosis of streamwise velocity (b). Lines are defined in Fig. 1. The red x-marks denote the measurements from Monty et al. (2009), digitalized by the authors.

Table 3 Relative error on skewness $\left|\right|S_{uu}-S_{uu,\mathrm{meas}}|{|}_{L^{\mathrm{2}}}/|\left|S_{uu,\mathrm{meas}}\right|{|}_{L^{\mathrm{2}}}$ and kurtosis $\left|\right|K_{uu}-K_{uu,\mathrm{meas}}|{|}_{L^{\mathrm{2}}}/|\left|K_{uu,\mathrm{meas}}\right|{|}_{L^{\mathrm{2}}}$ w.r.t. the measurements from Monty et al. (2009) in the interval $z_{\mathrm{0}}/h\le z/h\le \mathrm{0.4}$.

Figure 3 (a) Normalized one-dimensional spectra of streamwise velocity at $z/h\approx \mathrm{0.1}$. The solid red line depicts the (k_xz)⁻¹ production range and $(k_{x}z{)}^{-\mathrm{5}/\mathrm{3}}$ inertial subrange scaling. All other lines as in Fig. 1. (b) Premultiplied one-dimensional spectra of streamwise velocity at $z/h\approx \mathrm{0.1}$.

3.2 Spectra and autocorrelations

One-dimensional spectra of streamwise velocity fluctuations (E_uu) are shown in Fig. 3a. The profiles are contrasted against the phenomenological production range and inertial subrange power-law profiles (k⁻¹ and $k^{-\mathrm{5}/\mathrm{3}}$, respectively). Predictions from the PSFD-based solver feature a relatively good agreement with the phenomenological power-law profile, especially at high grid resolution. For example, the cases PSFD128 and PSFD160 exhibit a slope of −1.2 in the production range (here defined as k_xz<1). Profiles from the FV-based solver, on the contrary, exhibit strong sensitivity to grid resolution and are unable to capture the expected power-law behavior. In the production range, velocity spectra from the FV solver start off relatively shallow at small wavenumber, especially when using the linear scheme. A narrow band can be identified where $E_{uu}\sim (k_{x}z{)}^{-\mathrm{1}}$, followed by a rapid decay in energy density – the decay being particularly pronounced when using the QUICK interpolation scheme because of the associated numerical dissipation. Overall, the energy density in the production range and in the inertial subrange is not well captured by the FV-based solver and grid refinement does not help circumvent this limitation, at least at the considered resolutions. The authors note that this fact might limit the use of dynamic procedures based on the Germano et al. (1991) identity. A further characterization of the energy distribution in the wavenumber space is given in Fig. 3b, where premultiplied velocity spectra $k_x{E}_{uu}{u}_{*}^{-\mathrm{2}}$ are shown. The usual reason for considering these quantities is to create a plot in semi-log scale where equal areas under the profiles correspond to equal energy. In addition, premultiplied spectra provide information on the coherence of the flow, in particular on the so-called large and very large scale motions (LSMs and VLSMs, respectively). These structures are responsible for carrying more than half of the kinetic energy and Reynolds shear stress and are a persistent feature of the surface and outer layers of both aerodynamically smooth and rough walls (Kim and Adrian, 1999; Balakumar and Adrian, 2007; Monty et al., 2007; Hutchins and Marusic, 2007; Fang and Porté-Agel, 2015). The current domain is of modest dimensions and is able to accommodate only LSMs (Lozano-Durán and Jiménez, 2014), which are identified in premultiplied spectra by a local maximum at the streamwise wavenumber $k_x/h\approx \mathrm{1}$. The location of the peaks from the FV-based solver with linear interpolation scheme shifts toward higher wavenumber with grid refinement, with a maximum at $k_x/h\approx \mathrm{4}$ for the FV160 case. This fact signals a flow field where the streamwise extent of energetic modes (a.k.a., coherent structures) reduces as the grid is refined. On the contrary, the FV-based solver in combination with the QUICK scheme predicts the peak in premultiplied energy density at the expected wavenumber, hence suggesting that this approach is able to capture LSMs. The PSFD-based solver features a peak at the expected wavenumber (k_x=1) only at the lowest resolution (PSFD64). Profiles from the higher-resolution cases feature high energy densities at the lowest wavenumber, highlighting an artificial “periodization” of energy-containing structures in the streamwise direction. This behavior is linked to the limited horizontal extent of the computational domain. The authors have indeed verified that a larger domain (twice as large along each horizontal direction) enables one to capture LSMs with the PSFD solver at resolutions matching the one of the PSFD128 case (not shown). A corresponding single run was carried out with the FV solver over the said larger domain and premultiplied spectra were found to be in good agreement with those presented herein, supporting the conjecture that the proposed domain size suffices to capture the range of variability of FV solvers for the problem under consideration.

Figure 4 Contours of two-dimensional spatial autocorrelation of streamwise velocity at height $z/h\approx \mathrm{0.1}$ from the simulations FV64 (a), FV64* (b), PSFD64 (c), FV160 (d), FV160* (e), and PSFD160 (f). Contour levels from 0.1 to 0.9 with increments of 0.1.

To gain better insight on the spatial coherence of the flow field, the contour lines of the two-dimensional autocorrelation of the streamwise velocity $R_{uu}^{\mathrm{2}\mathrm{D}}$ in the xy plane are shown in Fig. 4. The $R_{uu}^{\mathrm{2}\mathrm{D}}=\mathrm{0.1}$ contour is often used to identify the boundaries of coherent structures populating the flow field. The contours from the FV-based solver with linear scheme (Figs. 4a, d) are representative of a poorly correlated flow field with a streamwise extent of the $R_{uu}^{\mathrm{2}\mathrm{D}}=\mathrm{0.1}$ contour of 0.5h and 0.1h along the streamwise and cross-stream directions, respectively. On the contrary, the contours from the FV-based solver with the QUICK scheme (Fig. 4b, e) depict a flow field characterized by larger spatial autocorrelation, in line with results from the PSFD-based solver. Note that the flow statistics presented above should not be impacted by the fact that the current domain size prevents some of the contour lines (simulations FV64*, FV160*, PSFD64, and PSFD160) from closing, as discussed in Lozano-Durán and Jiménez (2014).

Figure 5 One-dimensional spatial autocorrelation of streamwise velocity at height $z/h\approx \mathrm{0.1}$ along the streamwise direction (a) and along the cross-stream direction (b). Lines as in Fig. 1.

The one-dimensional spatial autocorrelation (R_uu), shown in Fig. 5 along the streamwise and cross-stream directions, further corroborates the above findings. From Fig. 5a it is apparent that the extension of the selected domain does not enable the flow to become completely uncorrelated in the streamwise direction for the PSFD solver and for the FV solver using QUICK; R_uu remains finite in the available $r_x/h$ range across resolutions. On the other hand, profiles from the FV-based solver using the linear interpolation rapidly decay towards zero. Along the cross-stream direction (Fig. 5b), profiles from the PSFD-based solver feature the expected negative lobes, highlighting the presence of high- and low-momentum streamwise-elongated streaks flanking each other in the said direction. This behavior is in line with findings from previous studies on the coherence of wall-bounded turbulence and with standard turbulence theory. Profiles from the FV-based solver exhibit a similar profile, albeit featuring a more rapid decay and less prominent negative lobes, especially for the high-resolution cases using the linear interpolation scheme. A quantitative measure of the coherence of the flow field is provided in Table 4, where the integral lengths ${\mathrm{\Lambda }}_{r_x,u}$ and ${\mathrm{\Lambda }}_{r_y,u}$ are reported for all the considered cases and compared against direct numerical simulations of a channel flow at Re_τ=2000 from Sillero et al. (2014). The integral lengths in Table 4 are evaluated at $z/h\approx \mathrm{0.15}$ since the data from Sillero et al. (2014) are available at this height. Although ${\mathrm{\Lambda }}_{r_x,u}$ might not be meaningful across the considered cases, owing to the lack of a zero crossing of the autocorrelation function, it is apparent that the FV-based solver underestimates the integral lengths when compared to the PSFD cases and the reference DNS values, especially when the linear interpolation scheme is used.

Sillero et al. (2014)

Table 4 Integral lengths at height $z/h\approx \mathrm{0.15}$.

Figure 6 Instantaneous snapshots of normalized streamwise velocity fluctuations at $z/h\approx \mathrm{0.1}$ from the simulations FV64 (a), FV64* (b), PSFD64 (c), FV160 (d), FV160* (e), and PSFD160 (f). The normalized velocity fluctuation is defined as $(u-\langle u{\rangle }_{xy})/u_{*}$, where averages (and fluctuations therefrom) are evaluated in space over the selected horizontal plane.

Instantaneous snapshots of streamwise velocity fluctuations over a horizontal plane support the above findings (see Fig. 6). Artificially periodized, streamwise-elongated bulges of uniform high and low momentum are indeed apparent in the snapshots from the PSFD-based solver (Fig. 6c, f). On the contrary, the instantaneous streamwise velocity field from the FV solver is populated by smaller regions of uniform momentum, especially when using the linear scheme, and the size of energetic structures diminishes with increasing grid resolution (see, e.g., Fig. 6a, d).

3.3 Momentum transfer mechanisms

This section is devoted to the analysis of momentum transfer mechanisms in the ABL with a focus on quadrant analysis (Lu and Willmarth, 1973) and on statistics of conditionally averaged flow fields.

The quadrant hole analysis is a technique based on the decomposition of the velocity fluctuations into four quadrants: the first and third quadrants, outward interactions ($u^{\prime }>\mathrm{0}$, $w^{\prime }>\mathrm{0}$) and inward interactions ($u^{\prime }<\mathrm{0}$, $w^{\prime }<\mathrm{0}$), respectively, are negative contributions to the momentum flux, whereas the second and fourth quadrants, a.k.a. ejections of low-speed fluid outward from the wall ($u^{\prime }<\mathrm{0}$, $w^{\prime }>\mathrm{0}$) and sweeps of high-speed fluid toward the wall ($u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$), represent positive contributions. A range of flow statistics can be defined based on this decomposition and used to provide insight on the mechanisms supporting momentum transfer in the ABL.

Figure 7 Stress fractions at $z/h\approx \mathrm{0.1}$. The profiles are normalized so that the sum of the stress fractions for H=0 is unity across the cases. Lines are defined in Fig. 1.

Figure 7 features the quadrant-hole analysis, where the notation is the same as in Yue et al. (2007a), with H being the hole size, S_i,H the resolved Reynolds shear stress contribution to the ith quadrant at hole size H, and $S_{i,H}^f$ is the corresponding quadrant fraction. Stress fractions are presented for values of the hole size H ranging from 0 to 8, where larger hole sizes correspond to contributions to the resolved Reynolds shear stress from more extreme events. Clearly, the FV-based solver with the linear scheme underpredicts ejections (Fig. 7a), outward interactions (Fig. 7b), and inward interactions (Fig. 7c), and overpredicts sweeps at large hole size H (Fig. 7d). On the contrary, the FV solver with the QUICK scheme underpredicts all the profiles except for the ejections, which are captured fairly well instead (see Fig. 7a). Note that ejections are violent events, concentrated over a very thin region in the cross-stream direction of the ABL (Fang and Porté-Agel, 2015).

Figure 8 Vertical structure of event ratios: (a) ratio of negative to positive contributions to the momentum flux; (b) ratio of sweeps to ejections. Lines as in Fig. 1.

To gain insight on the vertical structure of momentum transfer mechanisms, the exuberance ratio and the ratio of sweeps to ejections are analyzed in the following. Figure 8a shows the exuberance ratio, defined as the ratio of negative to positive contributions to the momentum flux, $(S_{\mathrm{1},\mathrm{0}}+S_{\mathrm{3},\mathrm{0}})/(S_{\mathrm{2},\mathrm{0}}+S_{\mathrm{4},\mathrm{0}})$ (Shaw et al., 1983). The exuberance ratios from the PSFD-based solver are larger in absolute value than the correspondent ones from the FV-based solver across the whole surface layer except very close to the surface. Profiles highlight that outward and inward interactions have a significant impact on the resolved Reynolds stress in the PSFD-based solver, whereas the flow simulated with the FV-based solver is characterized by a predominance of sweeps and ejections. This behavior is consistent throughout the ABL. Figure 8b shows the ratio of sweeps to ejections at the lowest portion of the ABL ($z/h\le \mathrm{0.4}$). Profiles obtained with the QUICK scheme are in line with predictions from the PSFD-based solver and with findings from measurements of surface-layer flow over rough surfaces, where ejections are identified as the dominant momentum transport mechanism in the ABL (Raupach et al., 1991). On the contrary, the FV-based solver with a linear scheme tends to favor sweeps over ejections as the mechanisms for momentum transfer in the surface layer.

Figure 9 Visualization of the conditionally averaged velocity field in the cross-stream vertical plane at $\mathrm{\Delta }x/h=\mathrm{0}$ from simulations FV64 (a), FV64* (b), PSFD64 (c), FV160 (d), FV160* (e), and PSFD160 (f). The conditional event is a positive streamwise velocity fluctuation at $\mathrm{\Delta }x/h=\mathrm{0}$, $\mathrm{\Delta }y/h=\mathrm{0}$, and $z/h=\mathrm{0.5}$. Colors are used to represent the magnitude of the streamwise component and vectors denote the cross-stream and vertical components.

Figure 10 Conditionally averaged flow field from simulations FV64 (a), FV64* (b), PSFD64 (c), FV160 (d), FV160* (e), and PSFD160 (f). The conditional average is computed as in Fig. 9. Red isosurfaces show positive fluctuations (>0.7, top; >0.65, bottom); blue isosurfaces show negative fluctuations ($<-\mathrm{0.55}$, top; $<-\mathrm{0.5}$, bottom).

To conclude the analysis on the mechanisms responsible for momentum transfer, velocity statistics from a conditionally averaged flow field are discussed next. The approach of Fang and Porté-Agel (2015) is adopted to compute the conditionally averaged flow field, where the conditional event is a positive streamwise velocity fluctuation at $\mathrm{\Delta }x/h=\mathrm{0}$, $\mathrm{\Delta }y/h=\mathrm{0}$, and $z/h=\mathrm{0.5}$. Figure 9 features a pseudocolor and vector plot of the conditionally averaged velocity field in a cross-stream vertical plane for selected cases, whereas Fig. 10 displays a three-dimensional isosurface thereof. The flow structure in the equilibrium surface layer is known to be characterized by counter-rotating rolls and low- and high-momentum streamwise-elongated streaks flanking each other in the cross-stream direction. Rolls and streaks are indeed the dominant flow mechanism responsible for tangential Reynolds stress (Ganapatisubramani et al., 2003; Lozano-Durán et al., 2012). As apparent from Fig. 9, the PSFD conditionally averaged velocity field exhibits counter-rotating patterns associated with positive and negative streamwise velocity fluctuations (corresponding to the aforementioned streaks). Throughout the ABL, the roll modes feature a diameter that is consistent with findings from the literature (d≈h). Moreover, positive and negative velocity fluctuations are approximately of the same magnitude $(\approx u_{*})$. From Fig. 10, it is apparent that the considered isosurfaces extend about 4h along the streamwise direction. Quite surprisingly, the FV-based solver is not able to predict the roll modes, irrespective of the interpolation scheme and grid resolution, and severely underpredicts the magnitude of the low-momentum streaks. Further, Figs. 9 and 10 both depict a FV conditionally averaged flow field that is poorly correlated along the cross-stream and streamwise directions, resulting in significantly smaller momentum-carrying structures. This fact supports previous findings from the two-dimensional spatial autocorrelation (Fig. 4). The lack of roll modes implies that the FV-based solvers used here are not able to capture the fundamental mechanism supporting momentum transfer in the ABL, at least at the considered grid resolutions. This limitation is likely to be the root cause of several of the observed problematics associated with the FV-solver solution, including the relatively high (low) streamwise-velocity skewness when using linear (QUICK) schemes (see Fig. 2,a) and the observed imbalance between sweeps and ejections (Figs. 1 and 8).

4 Conclusions

The present work provides insight on the quality and reliability of an important class of general-purpose, second-order accurate FV-based solvers for the wall-modeled LES of neutrally stratified ABL flow. The considered FV-based solvers are part of the OpenFOAM^® framework, make use of the divergence form for the nonlinear term, and are based on a colocated grid arrangement.

A suite of simulations was carried out in an open-channel flow setup, varying the grid resolution up to 160³ control volumes, the interpolation schemes for the discretization of the nonlinear term, the value of the Smagorinsky coefficient, the pressure-velocity coupling method, and the time-advancement scheme. Several flow statistics were contrasted against profiles from a well-established PSFD-based solver and against experimental measurements when these were available. Considered flow statistics include mean velocity, turbulence intensities, velocity skewness and kurtosis, velocity spectra, and spatial autocorrelations. An analysis of mechanisms supporting momentum transfer in the flow field was also proposed. The main findings are summarized below.

With the exception of the FV solver with the projection method and the Runge–Kutta time-advancement scheme, mean velocity profiles from the PSFD and FV solvers all feature a positive LLM. Existing techniques to alleviate this limitation led to no apparent improvement, thus calling for alternative approaches.

Near-surface streamwise velocity fluctuations are consistently overpredicted by both the PSFD and FV solvers, irrespective of the grid resolution. The overshoot is particularly pronounced for the cases based on the QUICK interpolation scheme. This behavior can be related to a deficit of pressure redistribution in the budget equations for the velocity variances, which results in a pile-up of shear-generated streamwise velocity fluctuations and deficit in the vertical and cross-stream velocity fluctuation components.

The interpolation scheme used for the discretization of the nonlinear term plays a role in determining the remaining flow statistics. Specifically, FV solvers with a linear interpolation scheme lead to

• a positive streamwise velocity skewness throughout the surface layer, which is at odds with experimental findings;

• a severe overprediction of the streamwise velocity kurtosis;

• a poorly correlated streamwise velocity field in the horizontal directions, especially at high grid resolutions;

• a severe underprediction of outward and inward interactions and ejection events;

• a lack of organized high- and low-momentum streaks and associated roll modes in the conditionally averaged flow field.

Grid resolution either does not affect the above quantities or leads to larger departures from the expected behavior. The QUICK scheme, on the other hand, leads to

• an improved prediction of the streamwise velocity skewness and kurtosis, especially as the grid stencil is reduced;

• a streamwise velocity field that is more correlated along the horizontal directions, but integral length scales remain only a fraction of those from the PSFD and reference DNS results;

• an underprediction of inward and outward interactions;

• a lack of organized high- and low-momentum streaks and associated roll modes in the conditionally averaged flow field.

To summarize, the considered class of FV-based solvers predicts a flow field that is less correlated than the one obtained with the PSFD solver and does not capture the salient mechanisms responsible for momentum transfer in the ABL, at least at the considered grid resolutions. These limitations appear to be the root cause of many of the observed discrepancies between FV flow statistics and the corresponding PSFD or experimental ones, including the mispredicted streamwise-velocity skewness (Fig. 2a), the imbalance between sweeps and ejections (Figs. 1 and 8), and the overall sensitivity of flow statistics to variations in the grid resolution. Higher grid resolutions might help alleviate some of these shortcomings, but given that grid resolutions used herein are state-of-the-art for general-purpose FV-based solvers and that computing power increases relatively slowly with time (Moore, 1965), the aforementioned limitations are likely to persist for years to come, thus introducing a degree of uncertainty in model results that needs to be addressed. These limitations call for research aimed at reducing the impact of discretization errors in this class of solvers, or for alternative approaches such as using discretizations based on staggered grid arrangements and higher-order spatial discretization schemes.

Appendix A

A1 Smagorinsky constant

We here test the sensitivity of selected flow statistics to variations in the Smagorinsky constant C_S. The values C_S=0.1, C_S=0.12, C_S=0.14, C_S=0.16, and C_S=0.1678 (the default value in OpenFOAM^®) are considered, and all tests are carried out at 64³ control volumes.

As shown in Fig. A1a, the Smagorinsky constant has a relatively important and nonmonotonic impact on the mean velocity profile. The case at C_S=0.1 results in the largest positive LLM, in agreement with the predictions from the PSFD-based solver, whereas the cases at larger C_S exhibit a smaller, albeit still positive, LLM. The Smagorinsky coefficient also has a discernible impact on the velocity RMSs. Specifically, as C_S is increased, the magnitude of the near-surface maximum for both $u_{\mathrm{RMS}}^{\prime }$ (Fig. A1b) and $w_{\mathrm{RMS}}^{\prime }$ (Fig. A1c) is reduced, and the location of the maximum is shifted away from the surface – possibly the result of a higher near-surface energy dissipation. In addition, larger values of C_S yield a more apparent departure from the corresponding profiles obtained with the PSFD-based solver.

The one-dimensional spectra (Fig. A2a) show that larger values of the Smagorinsky coefficient result in a more rapid decay of energy density and in a shift of profiles toward the inertial subrange. No value of the Smagorinsky coefficient seems suitable for capturing the k⁻¹ power law in the production range of turbulence. Increasing C_S leads to a modest improvement in the two-point autocorrelation profiles (Fig. A2b, c).

A2 Solvers

The performance of an alternative solver within the OpenFOAM^® framework is considered here, and the results are contrasted against those previously shown (obtained with the PISO algorithm in combination with an Adams–Moulton time-advancement scheme). The solver is based on a projection method coupled with the Runge–Kutta 4 time-advancement scheme (Ferziger and Peric, 2002). Details on the implementation can be found in Vuorinen et al. (2015). The performances of the two solvers are compared at moderate Reynolds number in Vuorinen et al. (2014), where it is pointed out that the projection method coupled with the Runge–Kutta 4 time-advancement scheme provides similar results at lower computational cost. In the following, the performances of the solver are tested for the considered ABL flow. Two grid resolutions are considered based on 64³ (case FV64RKp) and 128³ (case FV128RKp) control volumes.

The vertical profile of the mean streamwise velocity is shown in Fig. A3a. The use of the projection Runge–Kutta 4 solver leads to an underprediction of the velocity at the wall as for the simulations FV64 and FV128, but no apparent LLM in the surface layer. $u_{\mathrm{RMS}}^{\prime }$ exhibits the previously observed near-surface peaks (Fig. A3b) whereas $w_{\mathrm{RMS}}^{\prime }$ is overpredicted above $z/h=\mathrm{0.15}$ (Fig. A3c).

Figure A1 Vertical structure of streamwise velocity (a), streamwise velocity RMS (b), and vertical velocity RMS (c). Red lines denote the phenomenological logarithmic-layer profile (a) and analytical expressions from similarity theory (Stull, 1988) (b, c).

Figure A2 Normalized one-dimensional spectra of streamwise velocity at $z/h\approx \mathrm{0.1}$ (a); one-dimensional spatial autocorrelation of streamwise velocity at $z/h\approx \mathrm{0.1}$ along the streamwise direction (b) and along the cross-stream direction (c). Lines as in Fig. A1. The red line denotes (k_xz)⁻¹.

Figure A3 Vertical structure of streamwise velocity (a), streamwise velocity RMS (b), and vertical velocity RMS (c). Red lines denote the phenomenological logarithmic-layer profile (a) and the analytical expressions from similarity theory (Stull, 1988) (b, c).