On the Feasibility of Using Large-Eddy Simulations for Real-Time Turbulent-Flow Forecasting in the Atmospheric Boundary Layer

Abstract

We investigate the feasibility of using large-eddy simulation (LES) for real-time forecasting of instantaneous turbulent velocity fluctuations in the atmospheric boundary layer. Although LES is generally considered computationally too expensive for real-time use, wall-clock time can be significantly reduced by using very coarse meshes. Here, we focus on forecasting errors arising on such coarse grids, and investigate the trade-off between computational speed and accuracy. We omit any aspects related to state estimation or model bias, but rather look at the size and evolution of restriction errors, subgrid-scale errors, and chaotic divergence, to obtain a first idea of the feasibility of LES as a forecasting tool. To this end, we set-up an idealized test scenario in which the forecasting error in a neutral atmospheric boundary layer is investigated based on a fine reference simulation, and a series of coarser LES grids. We find that errors only slowly increase with grid coarsening, related to restriction errors that increase. Unexpectedly, modelling errors slightly decrease with grid coarsening, as both chaotic divergence and subgrid-scale error sources decrease. A practical example, inspired by wind-energy applications, reveals that there is a range of forecasting horizons for which the variance of the forecasting error is significantly reduced compared to the turbulent background variance, while at the same time, associated LES wall times are up to 300 times smaller than simulated time.

Appendices

A Restriction and Interpolation

We provide further details on the interpolation and restriction operators introduced in Sect. 3.1. First of all, formally, we define $\varvec{u}^{i} = [u_1^i,u_2^i,u_3^i] \in {\mathbb {R}}^{N^i}$, with $N^i = 3N_x^iN_y^iN_z^i-N_x^iN_y^i$ (cf. staggered arrangement of variables discussed in Sect. 3.1). Similarly, $\varvec{u}^{j}\in {\mathbb {R}}^{N^j}$, further using the convention that $i<j$ (so that j is the coarser grid). Consequently, for the interpolation and restriction operators in Eqs. 5 and 6, we have ${\mathcal {I}}_{j}^{i} \in {\mathbb {R}}^{N^i\times N^j}$, and ${\mathcal {R}}_{i}^{j} \in {\mathbb {R}}^{N^j\times N^i}$.

Since we use a Cartesian mesh, we split the interpolation and restriction operators in three consecutive one-dimensional operators, so that ${\mathcal {I}}_{j}^{i} = I_{j,z}^i I_{j,y}^i I_{j,x}^i$, and ${\mathcal {R}}_{i}^{j}=R_{i,z}^j R_{i,y}^j R_{i,x}^j$. The matrix $I_{j,z}^i$ has dimensions $N^i \times N^{iij}$ with $N^{iij}= 3N_x^i N_y^i N_z^j - N_x^i N_y^i$. The dimensions of $I_{j,y}^i$ are $N^{iij} \times N^{ijj}$, with $N^{ijj}= 3N_x^i N_y^j N_z^j - N_x^i N_y^j$, and the dimensions of $I_{j,x}^i$ correspond to $N^{ijj} \times N^j$. Similar dimensions follow straightforwardly for $R_{i,x}^j$, $R_{i,y}^j$, and $R_{i,z}^j$.

The rows of $I_{j,x}^i$, $I_{j,y}^i$, $I_{j,z}^i$ contain one-dimensional interpolation stencils (and similar for the restriction matrices). Therefore, below, we provide the stencils that we use based on a simple scalar function $\phi ^i$ and $\phi ^j$ along one-dimensional grids ${\varvec{r}}^i$ and ${\varvec{r}}^j$. The allocation of the different coefficients in these stencils to elements in the different rows of $I_{j,x}^i$, $I_{j,y}^i$, etc., is straightforward, and not further detailed for sake of brevity.

For the interpolation in the x- and y-directions, spectral interpolation is used, simply leading to

$$\begin{aligned} \phi ^i_k = \frac{1}{N^j}\sum _{m=0}^{N^j-1}\sum _{l=0}^{N^j-1}\mathrm {exp}\left( {-2\uppi {\mathrm{i}} m\frac{r^i_l-r^j_k}{L}}\right) \phi _l^j. \end{aligned}$$

(20)

where $\phi ^i_k$ and $\phi _l^j$ correspond to fine- and coarse-grid values on locations $r^i_k$ and $r^j_l$ respectively. In practice, we do not implement the interpolation in real space, but instead perform the operation in Fourier space.

For the interpolation in the z-direction, we use a polynomial interpolation of order p, where we take $p=4$, in analogy with our vertical discretization scheme. First to simplify notation, we define the operator $min_c(a,{\varvec{b}})$, which returns a set of the $c\in {\mathbb {N}}$ closest points in set ${\varvec{b}}\in {\mathbb {R}}^N$ to a scalar $a\in {\mathbb {R}}$. This simply gives for the interpolation operator

$$\begin{aligned} \phi ^i_k = \sum _{l\in {\mathcal {P}}}\left( \prod _{\begin{array}{c} m\in {\mathcal {P}}\\ m\ne l \end{array}}\frac{r_k^i -r_m^j}{r_l^j-r_m^j}\right) \phi ^j_l, \text { with } {\mathcal {P}} = min_p\left( r_k^i,{\varvec{r}}^j \right) , \end{aligned}$$

(21)

In analogy, the rows of $R_{i,x}^j$, $R_{i,y}^j$, $R_{i,z}^j$, contain the one-dimensional restriction stencils. For the restriction in x- and y-directions, a spectral cut-off filter is used in combination with simple injection to the coarse grid, leading to

$$\begin{aligned} \phi _l^j = \frac{1}{N^j}\sum _{m=0}^{N^j-1}\sum _{k=0}^{N^i-1}\mathrm {exp}\left( {-2\uppi {\mathrm{i}} m\frac{r^j_k-r^i_l}{L}}\right) \phi ^i_k. \end{aligned}$$

(22)

For the restriction in the z-direction, we use a combination of a box filter and an injection. For the box filter we use a width of $\varDelta _z^j$, which comes down to $s=\varDelta _z^j/\varDelta _z^i$ cells on the fine grid. It is easily shown that the following relation holds to filter a field $\phi ^i$, which is assumed to have been filtered with a width $\varDelta _z^i$, to a field $\phi ^j$ with a width $\varDelta _z^j$

$$\begin{aligned} {\widetilde{\phi }}^i_l = \frac{1}{s}\sum _{k=-s/2}^{s/2-1}\phi ^i_{l+k+1/2}, \end{aligned}$$

(23)

where the interpolation of $\phi ^i_{l+k+1/2}$ happens with the same fourth-order interpolation as is described above.

For the refinement experiment we use a Gaussian filter where the standard deviation is chosen as $\sigma ^2=(s^2-1)/12$, and where the factor 1 / 12 is determined such that the second moments of the Gaussian and box filter are equal [see Leonard (1975) for a derivation], and the factor $-1$ appears due to the successive filtering [see e.g. Pope (2000)], such that choosing $s=1$ leaves the field unaltered. This leads to the following relation

$$\begin{aligned} {\widetilde{\phi }}^i_l = \frac{\sum _{k=-\infty }^{\infty }\mathrm {exp}\left( {-\frac{1}{2}\left( \frac{k}{\sigma }\right) ^2}\right) \phi ^i_{l+k}}{\sum _{k=-\infty }^{\infty }\mathrm {exp}\left( -\frac{1}{2}\left( \frac{k}{\sigma }\right) ^2\right) }. \end{aligned}$$

(24)

In a further step the field is restricted to the coarser grid. Due to mismatching cell locations for the u and v velocity components an additional interpolation is needed. For this we again use the fourth-order polynomial interpolation, which leads to the following expression

$$\begin{aligned} \phi ^j_l = \sum _{k\in {\mathcal {P}}}\left( \prod _{\begin{array}{c} m\in {\mathcal {P}}\\ m\ne k \end{array}}\frac{r_l^j -r_m^i}{r_k^i-r_m^i}\right) {\widetilde{\phi }}^i_k, \text { with } {\mathcal {P}} = \mathrm {min}_p\left( r_l^j,{\varvec{r}}^i \right) . \end{aligned}$$

(25)

B Comparison of Time-Averaged Mean Fields

For the sake of completeness, we provide a comparison of time-averaged velocity and turbulent kinetic energy fields obtained on the different grids, which is the standard basis for comparing LES results (using different grids, models, or codes). In contrast to the error analysis in the main text, we present long time averages that omit the initial transient that occurs when initializing with a turbulent field that is not in statistical equilibrium on the simulation grid. To this end the simulations on the different grids are spun up until a statistical steady state is reached. Afterwards, averaging is performed over a period of ${8000}\,\hbox {s}$, ensuring sufficient statistical convergence.

Results are shown in Fig. 9, and overall, it is appreciated that profiles of the mean flow match closely. Profiles of turbulent kinetic energy show a more pronounced grid dependency close to the wall. This is quite standard, as the integral length scale decreases proportional with the distance to the wall, so that less large-scale motions are resolved in this region on coarser grids.

Fig. 9

Time-averaged equilibrium streamwise velocity component, $U_1^i$ (left) and turbulent kinetic energy, $E^{i}$ (right) for the different grids. Grid numbers: 0 ( ), 1 ( ), 2 ( ) and 3 ( )