2 Description of the numerical database

The set-up of the LES by Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014) corresponds to a flat plate where a statistically two-dimensional ZPG TBL, with density $\unicode[STIX]{x1D70C}$ , viscosity $\unicode[STIX]{x1D708}$ and free-stream velocity $U_{\infty }$ develops along the streamwise direction. The dimensions of the domain are $L_{x}\times L_{y}\times L_{z}=13\,500\times 400\times 540$ in the streamwise, wall-normal and spanwise directions, where non-dimensionalization with the displacement thickness at the inlet $\unicode[STIX]{x1D6FF}_{0}^{\ast }$ is considered. Here we employ a Reynolds decomposition of the velocity field into a spanwise and temporal average $U(x,y)$ and the fluctuating velocities, $u(x,y,z,t)$ , $v(x,y,z,t)$ and $w(x,y,z,t)$ . The outer velocity scale is $U_{\infty }$ , and the outer length scale is $\unicode[STIX]{x1D6FF}_{99}$ , i.e. the wall-normal position where the mean velocity $U(x,y)$ reaches 99 % of $U_{\infty }$ . The Reynolds number $Re_{\unicode[STIX]{x1D703}}=U_{\infty }\unicode[STIX]{x1D703}/\unicode[STIX]{x1D708}$ is then defined in terms of the momentum thickness $\unicode[STIX]{x1D703}$ . The inner scale is based on the friction velocity $u_{\unicode[STIX]{x1D70F}}=\sqrt{\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}}$ , which is obtained from the mean shear stress at the wall, $\unicode[STIX]{x1D70F}_{w}$ and the inner length scale is the viscous length $l_{\ast }=\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}}$ . Inner scaling will be denoted by the superscript $+$ . The friction Reynolds number is then defined as $Re_{\unicode[STIX]{x1D70F}}=u_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6FF}_{99}/\unicode[STIX]{x1D708}$ .

Following a procedure similar to that in wind-tunnel experiments, laminar flow enters the domain and is then forced to transition via a tripping (see Schlatter & Örlü (Reference Schlatter and Örlü2012) for a complete description of the method). The flow is tripped close to the leading edge, at $Re_{\unicode[STIX]{x1D703}}=180$ , where the tripping is more efficient. Transition to turbulence is considered to be complete at $Re_{\unicode[STIX]{x1D703}}\geqslant 600$ , as shown in Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014). The highest Reynolds number to be considered in this work is of $Re_{\unicode[STIX]{x1D703}}=8200$ , located upstream of the end of the computational domain in order to avoid outflow effects.

This study will focus on the data at the streamwise position corresponding to $Re_{\unicode[STIX]{x1D703}}=4430$ ( $Re_{\unicode[STIX]{x1D70F}}=1324$ ), however the cases corresponding to $Re_{\unicode[STIX]{x1D703}}=2240$ ( $Re_{\unicode[STIX]{x1D70F}}=704$ ) and $Re_{\unicode[STIX]{x1D703}}=8200$ ( $Re_{\unicode[STIX]{x1D70F}}=2370$ ) will also be analysed in order to evaluate the effect of the Reynolds number in the predictions and to perform the off-design studies.

The number of collocation points employed in the simulation is $13\,824\times 513\times 1152$ in the streamwise, wall-normal and spanwise directions. As shown by Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014) this resolution corresponds to a very finely resolved LES, in excellent agreement with DNS studies of ZPG TBL. The spectral code SIMSON (Chevalier, Lundbladh & Henningson Reference Chevalier, Lundbladh and Henningson2007) is used for this simulation, and was previously used in a number of studies by this group (Schlatter et al. Reference Schlatter, Örlü, Li, Brethouwer, Fransson, Johansson, Alfredsson and Henningson2009; Schlatter & Örlü Reference Schlatter and Örlü2010).

A total of 19 410 time samples is used, with a spacing in time of $\unicode[STIX]{x0394}t^{+}=0.5$ , which is appropriate for convergence of the statistics and is also adequate to correctly resolve the near-wall (which corresponds to a period $\unicode[STIX]{x1D706}_{t}^{+}\approx 100$ ) and outer ( $\unicode[STIX]{x1D706}_{t}\approx 10\unicode[STIX]{x1D6FF}_{99}$ ) peaks. The first half of the time samples is used to design the modelling methodologies and the other half to test them. The spanwise discretization results in 768 points separated by $\unicode[STIX]{x0394}z^{+}=12$ , which is appropriate to capture the near-wall (with a typical spanwise wavelength $\unicode[STIX]{x1D706}_{z}^{+}\approx 100$ ) and outer ( $\unicode[STIX]{x1D706}_{z}\approx \unicode[STIX]{x1D6FF}_{99}$ ) structures. For further details concerning the convergence of statistics of this simulation the reader is referred to Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014).

3 Pre-multiplied spectra and coherence

Prior to the construction of the transfer functions, an evaluation of the behaviour of the two-dimensional spectra and the coherence between sets of positions separated along the wall-normal direction was performed. The objective of such assessment was to evaluate the suitability of models for the prediction.

The focus will be on three locations in the wall-normal direction, $y^{+}=15$ , 50 and 100, which correspond to positions in the near-wall region, just outside the buffer layer and representative of the logarithmic region, respectively. This may be observed in figure 1 which shows the mean flow profile, scaled in inner units, for the four Reynolds numbers evaluated in this work. Regions I, II, III and IV represent the viscous sublayer, buffer and logarithmic layers and wake region, respectively. The linear and logarithmic profiles are shown to better highlight the behaviour of the mean velocity profile; the inner and outer regions (the latter shown for $Re_{\unicode[STIX]{x1D703}}=8200$ ) are also shown on the upper part of the plot.

Figure 1. Inner-scaled mean velocity profiles for the four Reynolds numbers evaluated in the current study. Regions I, II, III and IV represent the viscous sublayer ( $y^{+}\lessapprox 7$ ), buffer ( $7\lessapprox y^{+}\lessapprox 30$ ) and logarithmic layers ( $30\lessapprox y^{+}\lessapprox 220$ ) and wake region ( $220\lessapprox y^{+}$ until the edge of the boundary layer, for the $Re_{\unicode[STIX]{x1D703}}=8200$ case). The pink dashed lines correspond to the linear ( $U^{+}=y^{+}$ ) and logarithmic profiles ( $U^{+}=1/\unicode[STIX]{x1D705}log(y^{+})+B$ , where $\unicode[STIX]{x1D705}=0.41$ and $B=5.1$ ), respectively. The upper limit of the outer region corresponds to the $Re_{\unicode[STIX]{x1D703}}=8200$ case, extending until the edge of the boundary layer.

The inner-scaled, two-dimensional, pre-multiplied spectrum is defined as $E_{uu}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6FD}/u_{\unicode[STIX]{x1D70F}}^{2}$ , where $E_{uu}=\langle \hat{\hat{u} }\hat{\hat{u} }^{\ast }\rangle$ , with the double hat representing the double Fourier transform of $u$ in the spanwise direction and time, the brackets $\langle \rangle$ denote an ensemble averaging as defined in appendix A, $\ast$ defines the complex conjugate, $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D6FD}$ are the frequency and spanwise wavenumbers and $\unicode[STIX]{x1D706}_{t}^{+}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$ , $\unicode[STIX]{x1D706}_{z}^{+}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FD}$ , given in inner units. Spectra have been averaged using Welch’s method, using segments with 256 snapshots with 75 % overlap and a time discretization $\unicode[STIX]{x0394}t^{+}=0.5$ , with a triangular window leading to a frequency resolution of $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}^{+}=0.05$ . Figure 2 shows the time signal and pre-multiplied spectra of the streamwise velocity fluctuation at three locations, $y^{+}=15$ , 50 and 100, for the $Re_{\unicode[STIX]{x1D703}}=4430$ , (corresponding to $Re_{\unicode[STIX]{x1D70F}}=1324$ ), case. Velocity fluctuations are normalized using the friction velocity $u_{\unicode[STIX]{x1D70F}}$ in the results shown in figure 2, and also for all time series presented later in the paper.

Figure 2. Behaviour of the streamwise velocity fluctuations and inner-scaled pre-multiplied two-dimensional spectra at different wall-normal positions for $Re_{\unicode[STIX]{x1D703}}=4430$ . The crosses correspond to reference locations for the inner ( $(\unicode[STIX]{x1D706}_{t}^{+},\unicode[STIX]{x1D706}_{z}^{+})\approx (100,100)$ ) and outer peaks ( $(\unicode[STIX]{x1D706}_{t},\unicode[STIX]{x1D706}_{z})\approx (10\unicode[STIX]{x1D6FF}_{99}/U_{\infty },\unicode[STIX]{x1D6FF}_{99})$ , $(\unicode[STIX]{x1D706}_{t}^{+},\unicode[STIX]{x1D706}_{z}^{+})\approx (1400,1300)$ ). The limits of the colour bar were kept the same between different plots to highlight the different amplitudes. (a) Instantaneous streamwise velocity fluctuations at $y^{+}=15$ . (b) Pre-multiplied spectra at $y^{+}=15$ . (c) Instantaneous streamwise velocity fluctuations at $y^{+}=50$ . (d) Pre-multiplied spectra at $y^{+}=50$ . (e) Instantaneous streamwise velocity fluctuations at $y^{+}=100$ . (f) Pre-multiplied spectra at $y^{+}=100$ .

The dominant spanwise wavelength of the fluctuations grows as one moves further in the wall-normal direction, which leads to a wider scale separation between the large-scale streaks in the outer region and the near-wall streaks, a feature which may be observed in the two-dimensional pre-multiplied spectra of figure 2. For $y^{+}=15$ , representative of the buffer layer, the classical near-wall streaks with spacing $\unicode[STIX]{x1D706}_{z}^{+}\approx 100$ are observed, with a characteristic time $\unicode[STIX]{x1D706}_{t}^{+}\approx 100$ . In the outer region a secondary peak, which scales in outer units, is observed with characteristic scales $\unicode[STIX]{x1D706}_{z}\approx \unicode[STIX]{x1D6FF}_{99}$ and $\unicode[STIX]{x1D706}_{t}\approx 10\unicode[STIX]{x1D6FF}_{99}/U_{\infty }$ . In order to extract the relation between near-wall and outer positions, the coherence between two sets of synchronized signals separated in the wall-normal direction was computed using time series of the streamwise velocity fluctuation. For a statistically two-dimensional ZPG TBL, which is homogeneous in the spanwise direction, we define:

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}^{2}(y_{in},y_{out},\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD})=\frac{|\langle \hat{\hat{u} }(y_{in},\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD})\hat{\hat{u} }^{\ast }(y_{out},\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD})\rangle |^{2}}{\langle |\hat{\hat{u} }(y_{in},\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD})|^{2}\rangle \langle |\hat{\hat{u} }(y_{out},\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD})|^{2}\rangle },\end{eqnarray}$$

where $y_{in}$ , $y_{out}$ correspond to the wall-normal position of the input/output pair of time series that are considered in the analysis. The coherence is a quantity analogous to the correlation, but defined in the frequency domain (here also in wavenumber domain, taking advantage of homogeneity of the spanwise coordinate) and it represents a metric of linear behaviour between two signals, for a given $(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD})$ pair. Its value is normalized between zero and one, where the former indicates a complete random behaviour and the latter an exact linear relation between the two signals (namely an amplitude change and a phase shift). Linear, time and span invariant models with single input and single output lead to unit coherence between input and output (Bendat & Piersol Reference Bendat and Piersol2011); hence, coherence values close to one suggest the possibility of modelling the observed behaviour with a linear model.

The objective that we will pursue in the next sections is to predict the behaviour of the fluctuations in the near-wall region, such that the output position will be fixed at $y^{+}=15$ , where a clear peak is observed in the pre-multiplied spectra shown in figure 2. Sample inputs are chosen at $y^{+}=50$ and $100$ , where the larger scales become prominent and in a third case using the wall-shear stress. The objective of the latter is to capture the finer-scale motions close to the wall.

Figure 3 shows the coherence for three pairs of positions. It is noticeable that for the higher wall-normal positions the long spanwise wavelengths are coherent. This observation is related to the wall footprint of large-scale motions that are attached to the wall, as in Marusic et al. (Reference Marusic, Mathis and Hutchins2010). We see that larger-scale structures at outer positions have a significant footprint at the buffer layer, with coherences ranging between 0.6 and 0.9 for wavelengths and periods corresponding to the outer peak. However, the near-wall peak $(\unicode[STIX]{x1D706}_{z}^{+},\unicode[STIX]{x1D706}_{t}^{+}=100,100)$ is seen to be incoherent with outer-layer disturbances. In turn, the coherence using wall-shear stress as input exhibits values close to unity at the lower wavelengths, corresponding to the near-wall peak in the pre-multiplied spectra of figure 2; this indicates that shear-stress fluctuations may be an interesting candidate as an input measurement for prediction of buffer-layer fluctuations (Örlü & Schlatter Reference Örlü and Schlatter2011).

Figure 3. Coherence between two positions separated in the wall-normal direction, for $Re_{\unicode[STIX]{x1D703}}=4430$ . The crosses correspond to reference locations for the inner ( $(\unicode[STIX]{x1D706}_{t}^{+},\unicode[STIX]{x1D706}_{z}^{+})\approx (100,100)$ ) and outer peaks ( $(\unicode[STIX]{x1D706}_{t},\unicode[STIX]{x1D706}_{z})\approx (10\unicode[STIX]{x1D6FF}_{99}/U_{\infty },\unicode[STIX]{x1D6FF}_{99})$ , $(\unicode[STIX]{x1D706}_{t}^{+},\unicode[STIX]{x1D706}_{z}^{+})\approx (1400,1300)$ ). (a) Coherence between streamwise velocity fluctuations at $y^{+}=50$ and 15. (b) Coherence between streamwise velocity fluctuations $y^{+}=100$ and 15. (c) Coherence between wall-shear stress and streamwise velocity fluctuations at $y^{+}=15$ .

In order to better assess the behaviour of the coherence and its relation to the peak in the two-dimensional spectra, figure 4 presents the one-dimensional pre-multiplied temporal spectrum ( $E_{uu}^{1D}\unicode[STIX]{x1D714}/u_{\unicode[STIX]{x1D70F}}^{2}$ ) as a function of the wall-normal direction overlaid on the one-dimensional coherence,

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{1D}^{2}(y_{in},y_{out},\unicode[STIX]{x1D714})=\frac{|\langle \hat{\hat{u} }(y_{in},\unicode[STIX]{x1D714})\hat{\hat{u} }^{\ast }(y_{out},\unicode[STIX]{x1D714})\rangle |^{2}}{\langle |\hat{\hat{u} }(y_{in},\unicode[STIX]{x1D714})|^{2}\rangle \langle |\hat{\hat{u} }(y_{out},\unicode[STIX]{x1D714})|^{2}\rangle },\end{eqnarray}$$

which is calculated between a pair of positions at the same spanwise location, one is at a fixed wall-normal location and the other one is moved. The reference input position is in the viscous sublayer, i.e. at $y_{in}^{+}=5$ and two Reynolds numbers were considered, $Re_{\unicode[STIX]{x1D703}}=4430$ ( $Re_{\unicode[STIX]{x1D70F}}=1324$ ) and 8200 ( $Re_{\unicode[STIX]{x1D70F}}=2370$ ). For both cases it is noticeable that the larger scales exhibit significant values of coherence with the near-wall region, illustrating the fact that the largest scales are indeed attached to the wall. It is also observed that, as the Reynolds number is increased, the larger scales remain coherent with the near-wall region until much further away from the wall, including the outer spectral peak. Wall-normal positions up to $y^{+}\approx 500$ continue to exhibit coherence levels of 0.4 with the near wall indicating that the footprint of such large scales is clearly present for the $Re_{\unicode[STIX]{x1D703}}=8200$ case. This enables the use of data in the outer region to predict specific flow features close to the wall. As for the finer scales, they present high levels of coherence only close to the wall, which permits the use of local wall-shear-stress data or even universal signals to predict the inner-peak behaviour.

The high values of coherence reported here motivate the use of predictive techniques for application in the turbulent zero-pressure-gradient boundary layer. Through this paper, such models will be defined in terms of linear and nonlinear transfer functions, which will be introduced in the next section.

Other studies using different methodologies, and based on experimental databases (Marusic et al. Reference Marusic, Mathis and Hutchins2010; Baars et al. Reference Baars, Hutchins and Marusic2016), perform analyses up to $Re_{\unicode[STIX]{x1D70F}}=19\,000$ and observe that the footprint of the larger scales on the near-wall region remains. Particularly, the work of Baars et al. (Reference Baars, Hutchins and Marusic2016) shows quantitative results based on the one-dimensional (1-D) coherence spectrum up to $Re_{\unicode[STIX]{x1D70F}}=13\,300$ . These works supply evidence that the assumptions of the modelling methodologies herein developed would remain valid at considerably higher Reynolds numbers. Nevertheless, the footprint of the larger scales over the near-wall region and corresponding high coherence values are necessary ingredients for the use of the methods outlined in this work.

4 Linear and nonlinear transfer functions

Figure 4. One-dimensional coherence (solid lines) overlaid on the one-dimensional pre-multiplied temporal spectrum as a function of the wall-normal direction for $Re_{\unicode[STIX]{x1D703}}=4430$ and 8200. Coherence levels of 0.4 until 0.9 are highlighted.

4.1 Single-input linear transfer function

The basic supporting idea for the different identification methods that will be considered in this study is that there exists a function $f(I(z,t))$ that maps the relation between two wall-normal positions in the TBL, namely $I(z,t)$ and $O(z,t)$ . These two sets of data, taken along the spanwise direction and time, will be referred to as the input and output of the system, respectively. A qualitative description of such relation may be given in terms of the block diagram shown in figure 5. It should be noted that no restrictions have been made regarding the variables to be considered or the linearity of the function  $f$ .

The methods therefore assume the availability of spanwise- and time-resolved measurements, which are available in a numerical simulation but could be more difficult to obtain experimentally. Such turbulence measurements above the wall would require the use of time-resolved particle image velocimetry (PIV), such as in the work of Cuvier et al. (Reference Cuvier, Srinath, Stanislas, Foucaut, Laval, Kähler, Hain, Scharnowski, Schröder, Geisler, Agocs, Röse, Willert, Klinner, Amili, Atkinson and Soria2017). An alternative is to employ wall-shear stress sensors, as performed in closed-loop applications, such as in Lundell (Reference Lundell2007).

Figure 5. Block diagram describing the relationship between the input and output signals within the TBL.

The assumption in this section will be that there is a linear relationship between the two sets of measurements, such that the output may be predicted by means of a double convolution (along time and the spanwise direction) between the input measurement and a kernel,

(4.1) $$\begin{eqnarray}O(z,t)=\int _{-z_{max}/2}^{z_{max}/2}\int _{-\infty }^{\infty }g_{IO}(\unicode[STIX]{x1D701},\unicode[STIX]{x1D70F})I(z-\unicode[STIX]{x1D701},t-\unicode[STIX]{x1D70F})\,\text{d}\unicode[STIX]{x1D701}\,\text{d}\unicode[STIX]{x1D70F},\end{eqnarray}$$

which can be written in short as,

(4.2) $$\begin{eqnarray}O(z,t)=g_{IO}(z,t)\otimes I(z,t),\end{eqnarray}$$

where $\otimes$ represents the double convolution. Here $g_{IO}(\unicode[STIX]{x1D701},\unicode[STIX]{x1D70F})$ is the convolution kernel, which, in accordance with the common nomenclature of the field, will be referred to as a transfer function in the remaining of this paper. The dummy variables $(\unicode[STIX]{x1D701},\unicode[STIX]{x1D70F})$ , which are analogous to $(z,t)$ , were introduced for the calculation of the convolution. In order to obtain $g_{IO}(\unicode[STIX]{x1D701},\unicode[STIX]{x1D70F})$ , the problem is formulated in the frequency/spanwise wavenumber space, where the optimal frequency response, in the least squares sense, may be defined from the auto- and cross-spectra of the input and output signals (Bendat & Piersol Reference Bendat and Piersol2011; Sasaki et al. Reference Sasaki, Piantanida, Cavalieri and Jordan2017, Reference Sasaki, Morra, Fabbiane, Cavalieri, Hanifi and Henningson2018a ):

(4.3) $$\begin{eqnarray}G_{IO}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})=\frac{S_{IO}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})}{S_{II}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})},\end{eqnarray}$$

where $S_{II}$ and $S_{IO}$ are the auto- and cross-spectra of the input and output, respectively. These quantities are calculated from the expected values of $\hat{\hat{I} }(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\hat{\hat{I} }^{\ast }(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})$ and $\hat{\hat{I} }(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\hat{\hat{O} }^{\ast }(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})$ , obtained from an ensemble averaging. The details of the method are outlined in appendix A. It should be noted that since the spanwise direction is periodic in the simulation, the ensemble averaging along this direction corresponds directly to the Fourier transform from $z$ to $\unicode[STIX]{x1D6FD}$ . Equation (4.3) is referred to as the $H_{1}$ estimator of the system (Rocklin, Crowley & Vold Reference Rocklin, Crowley and Vold1985), and it minimizes the error due to measurement noise in the output. Other formulations, such as the $H_{2}$ or $H_{\unicode[STIX]{x1D708}}$ estimators, exhibit different performances in terms of sensor noise minimization. They are, however, expected to perform equally well for this type of estimation, which does not consider the presence of measurement uncertainties.

One interesting property of the $H_{1}$ estimator is that it leads to a prediction error which is linearly uncorrelated with the input (Rocklin et al. Reference Rocklin, Crowley and Vold1985; Bendat Reference Bendat1993). Consider the input/output relationship, given in the frequency–wavenumber domain,

(4.4) $$\begin{eqnarray}\hat{\hat{O} }(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})=\hat{{\hat{G}}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\hat{\hat{I} }(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})+\unicode[STIX]{x1D716}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}),\end{eqnarray}$$

where $\unicode[STIX]{x1D716}$ will represent the error of the linear estimator. Post-multiplying equation (4.4) by the conjugate transpose of the input and taking the expected values, we obtain

(4.5) $$\begin{eqnarray}S_{IO}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})=\hat{{\hat{G}}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})S_{II}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})+S_{\unicode[STIX]{x1D716},I}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}).\end{eqnarray}$$

Inserting the $H_{1}$ estimator of (4.3) leads to $S_{\unicode[STIX]{x1D716},x}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})=0$ . Remaining correlated errors are related either to measurement noise, which can be considered negligible in this problem, or spectral leakage (Schoukens, Rolain & Pintelon Reference Schoukens, Rolain and Pintelon1997), which is minimized by use of long time series.

Performing an inverse Fourier transform of $G_{IO}$ as:

(4.6) $$\begin{eqnarray}g_{IO}(z,t)=\int _{-\infty }^{+\infty }\int _{-\unicode[STIX]{x1D6FD}_{n}}^{+\unicode[STIX]{x1D6FD}_{n}}G_{IO}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\text{e}^{\text{i}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D701}}\text{e}^{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}}\,\text{d}\unicode[STIX]{x1D6FD}\,\text{d}\unicode[STIX]{x1D714},\end{eqnarray}$$

where $\pm \unicode[STIX]{x1D6FD}_{n}$ represent the spanwise Nyquist wavenumbers calculated based on the spanwise discretization, leads to $g_{IO}(\unicode[STIX]{x1D701},\unicode[STIX]{x1D70F})$ . This represents an empirical, linear, time-domain transfer function, which can be used to estimate the evolution of the fields separated in the wall-normal direction.

Figure 6. Block diagram describing the relationship between the multiple inputs and output of a system.

4.2 Multiple-input linear transfer function

In this section the construction of a model with multiple inputs will be introduced. This method applies to problems where multiple time series are available to be used in the estimation of the output. A qualitative schematic of the estimation is shown in the block diagram of figure 6. In this figure a system with three inputs is considered, however the extension to any number of inputs is straightforward. The linear hypothesis is inherent to this formulation since the output is assumed to be determined directly from the superposition of the contributions from each input, $I_{i}(z,t)$ , via the function $f_{i}(I_{i}(z,t))$ , which is estimated in terms of a convolution kernel, $g_{I_{i}O}$ . The estimation of the output is therefore obtained as follows:

(4.7) $$\begin{eqnarray}O(z,t)=\mathop{\sum }_{i=1}^{n}I_{i}(z,t)\otimes g_{I_{i}O}(z,t).\end{eqnarray}$$

As in the case of the single-input problem, the transfer functions are found in the $(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD})$ space and then inverse Fourier transformed to the physical domain. For a system with $n$ inputs, the transfer functions are obtained from the solution of the linear system:

(4.8) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}S_{I_{1}O}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\\ S_{I_{2}O}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\\ \vdots \\ S_{I_{n}O}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\end{array}\right)=\left(\begin{array}{@{}cccc@{}}S_{I_{1}I_{1}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}) & S_{I_{1}I_{2}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}) & \ldots & S_{I_{1}I_{n}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\\ S_{I_{2}I_{1}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}) & S_{I_{2}I_{2}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}) & \ldots & S_{I_{2}I_{n}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\\ \vdots & \vdots & \ddots & \vdots \\ S_{I_{n}I_{1}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}) & S_{I_{n}I_{2}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714}) & \ldots & S_{I_{n}I_{n}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\end{array}\right)\left(\begin{array}{@{}c@{}}G_{I_{1}O}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\\ G_{I_{2}O}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\\ \vdots \\ G_{I_{n}O}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\end{array}\right),\end{eqnarray}$$

which is performed for each frequency independently. The terms $S_{I_{i}O}$ define the cross-spectral density between the inputs and output, whereas $S_{I_{i}I_{j}}$ are the auto-spectral density of the inputs, when $i=j$ , or the cross-spectral density between different inputs, when $i\neq j$ . The auto- and cross-spectral densities are also calculated by means of ensemble averaging, as described in appendix A.

Note that it could also be advantageous to deal with incoherent inputs, such that $S_{I_{i}I_{j}}=0$ if $i\neq j$ , and only the diagonal terms in (4.8) remain. If the inputs are not correlated, each multiple-input–single-output (MISO) transfer function can be calculated directly as per (4.3), $G_{I_{i}O}=S_{I_{i}O}/S_{I_{i}I_{i}}$ , without the need to solve the matrix problem in (4.8). In order to decorrelate the inputs, the spectral conditioning technique (Bendat & Piersol Reference Bendat and Piersol2011) may be used. This has been used in all multiple-input cases studied here and the technique to do so is outlined in appendix B.

Each transfer function has to be inverse Fourier transformed, as per (4.6), in order to perform the estimation via the sum of convolutions in (4.7).

4.3 Nonlinear transfer function

Now we introduce the concept of nonlinear spectral analysis, via the use of nonlinear transfer functions, to the prediction of the interaction between the input and output signals. One of the few uses of nonlinear estimation methods in fluid mechanics may be found in the work of Naguib, Wark & Juckenhöfel (Reference Naguib, Wark and Juckenhöfel2001), who made use of stochastic estimation using nonlinear flow sources. On the other hand, nonlinear spectral analysis has been applied in previous works by means of a Volterra series to the derivation of nonlinear output frequency responses, with the objective of deriving properties of nonlinear systems (Guillaume, Pintelon & Schoukens Reference Guillaume, Pintelon and Schoukens1992; Lang & Billings Reference Lang and Billings1996; Peng, Lang & Billings Reference Peng, Lang and Billings2007). This technique assumes the expansion of the data into a power series, using an $N$ -integral kernel, leading to higher-order spectra. It presents drawbacks related both to the prohibiting computational demand, which limits the size of the system, and the large errors which can only be avoided with very large amounts of data.

Here we will follow a somewhat different strategy, by considering a priori the existence of a nonlinear model between the defined input and output of the system and working on determining such model. The exact shape of the nonlinearities is unknown and they are sought with the objective of reproducing the observed behaviour. This method is thoroughly discussed in the works of Rice & Fitzpatrick (Reference Rice and Fitzpatrick1988, Reference Rice and Fitzpatrick1991) and Bendat (Reference Bendat1993), for time-dependent problems; in this study we extend the same ideas to application in a boundary layer with periodicity in the spanwise direction.

The nonlinear transfer-function technique considered here consists in replacing a single-input nonlinear system by a multiple-input linear system, where the linear identification techniques of the previous section may be applied. Consider a finite-memory, single-input–single-output nonlinear function $g(\cdot )$ , which acts on $I(z,t)$ , mapping it into the output $O(z,t)$ :

(4.9) $$\begin{eqnarray}O(z,t)=g(I(z,t)),\end{eqnarray}$$

as shown schematically in the block diagram of figure 7. Such nonlinear system may be represented by a superposition of linear and nonlinear operations performed on the input, as given by:

(4.10) $$\begin{eqnarray}O(z,t)=\mathop{\sum }_{i=1}^{N_{TF}}g_{i}(I(z,t)),\end{eqnarray}$$

where $g_{i}$ represents an operation performed on the signal $I(z,t)$ and $N_{TF}$ is the number of transfer functions considered. Each one of the nonlinear operators $g_{i}$ is now replaced by a nonlinear operator $h_{i}$ and a linear operator $A_{i}$ , such that:

(4.11) $$\begin{eqnarray}g_{i}(I(z,t))=A_{i}(h_{i}(I(z,t))).\end{eqnarray}$$

We now define the signal $\tilde{I_{i}}(O,t)$ , as the result of the nonlinear operation $h_{i}$ applied to $I(O,t)$ . If all nonlinearities are encapsulated in the operators $h_{i}$ , the resulting system may now be treated as having inputs $\tilde{I} _{i}(O,t)$ and output $O(z,t)$ , such that linear frequency-domain techniques are applied to the determination of the convolution kernels $A_{i}$ , in the same manner as derived in the previous section, for a multiple-input–single-output (MISO) problem. Figure 8 illustrates the steps in obtaining this new system on a sample case where three functions were considered in the determination of the output. Note that the extension to any number of input functions is straightforward.

Figure 7. Block diagram corresponding to the nonlinear operation applied to $I(z,t)$ .

Figure 8. Summary of steps to transform a single-input nonlinear problem into a multiple-input linear system. (a) Nonlinear function is replaced by a sum of operations. (b) Each nonlinear function is broken into linear and nonlinear contributions. (c) The nonlinear operations are treated as known and the system is converted into the more usual MISO problem.

It should also be noted that when the exact shape of the nonlinearities is unknown, as will be the case here, the usual strategy is to deal with polynomial-like functions, with a degree sought such that a given performance is achieved; hence the method becomes more data driven than the previous approaches.

6 Off-design predictions

As outlined in the previous sections of this work, the on-line prediction of the designed transfer functions is made in two steps. The convolution kernels are firstly built using simultaneous time- and spanwise-resolved signals at the positions corresponding to inputs and output of the system, which are separated along the wall-normal direction. Once this system has been built, the prediction is made exclusively from measurements of the inputs.

This therefore raises questions accounting for the robustness of the methods for mild Reynolds-number variations, which would certainly be present in experimental applications. Furthermore, the fact that simultaneous measurements have to be performed as an initial step may lead to constraints in practical applications, where the correspondingly high Reynolds number would require a very fine resolution in the wall-normal direction, which may be difficult to obtain. Such constraints are shared by other methods available in the literature which require an initial set of data to build the model (also referred to as a training dataset in the context of machine learning methods).

One solution would therefore be to perform a large-eddy or direct numerical simulation of the geometry to be studied experimentally, and build the transfer functions a priori, which would only require the capability of measuring the inputs in the experiment (i.e. wall-shear stress and the fluctuations outside the buffer layer). Since high-fidelity simulations can typically only be performed at Reynolds numbers lower than those in certain experimental studies, it is of interest to design a model using a moderate Reynolds number, obtained from a high-fidelity simulation, and apply it to an experimental implementation at higher $Re$ . This possibility is explored in this section, where the transfer functions will be designed at $Re_{\unicode[STIX]{x1D703}}=880$ and tested at $Re_{\unicode[STIX]{x1D703}}=8200$ , spanning approximately one order of magnitude in Reynolds number. Three cases will be considered here: linear SISO and MISO using wall shear and data outside the buffer layer, and the quadratic TF using wall-shear stress, which exhibits mild improvements in comparison to the linear case. The same input/output positions for the two Reynolds numbers will be considered in terms of the corresponding inner variables, which showed better performance than if the outer non-dimensionalization was utilized. Figure 21 shows the resulting correlations and variances for these three methods, for the original and off-design cases. For all the cases, the correlations in the off-design configurations are degraded, where the SISO implementation is the most affected. It should also be noted, however, that the prediction of the near-wall peak is not much affected, where compelling correlations are seen for the SISO using wall-shear stress and MISO methods. The prediction of the MISO transfer function is also less affected by the Reynolds-number variation, where it should be noted that for the case with four inputs, the off-design prediction remains with compelling correlations (above 60 %) throughout the evaluated range of wall-normal positions. The fact that the near-wall peak prediction from wall-shear data was practically unaffected supplies evidence that the near-wall region is directly correlated to the wall. The increase in $Re$ has an effect on the large-scale modulation, modifying the footprint of the large-scale structures near the wall which in turn leads to the high degradation of the prediction of the near-wall peak when SISO transfer functions using data above the wall are used. We have therefore demonstrated the feasibility of designing SISO and MISO transfer functions at lower Reynolds numbers and implementing them at higher $Re$ (one order of magnitude higher in this particular case), along with a considerable robustness of the methodology.

Figure 21. (a,c,e) Correlations between prediction and LES data and (b,d,f) comparison of the resulting variances using the SISO, MISO and nonlinear methods. Solid lines correspond to the same Reynolds number being used for design and test, whereas dotted lines are used for off-design cases, i.e. the model is built at $Re_{\unicode[STIX]{x1D703}}=880$ and tested at $Re_{\unicode[STIX]{x1D703}}=8200$ . (a) Correlations using linear SISO. (b) Variances using linear SISO. (c) Correlations using the linear MISO. (d) Variances using the linear MISO. (e) Correlations using quadratic SISO. (f) Variances using quadratic SISO.

Although not shown for the sake of brevity, variations of up to 50 % in Reynolds number had no significant impact on the accuracy of the prediction. This feature is expected to be appealing for the application of the developed methods for experimental implementations, where uncertainties in the measured parameters are expected to be present.

7 Comparison with other methods

The methodologies developed in this work bear similarities with the recent studies by Baars et al. (Reference Baars, Hutchins and Marusic2016) and Illingworth et al. (Reference Illingworth, Monty and Marusic2018). The different focus of the first, which is on the characterization of the statistics, and application of the second, a channel at $Re_{\unicode[STIX]{x1D70F}}=1000$ , prevent quantitative comparisons of performance between these methods and the one proposed here. The objective of the current section is therefore to highlight differences in the approach, applicability, limitations and computational cost of the methods.

Illingworth et al. (Reference Illingworth, Monty and Marusic2018) considered a turbulent channel flow of $Re_{\unicode[STIX]{x1D70F}}=1000$ , obtained from DNS. Their approach is based on a Reynolds decomposition of the Navier–Stokes operator and gathering the nonlinear terms from the fluctuations into a forcing to an otherwise linear system. A Fourier transform is then applied to the homogeneous coordinates, corresponding to stream- and spanwise directions. The resulting inhomogeneous Orr–Sommerfeld squire (OSS) system is expressed in the form of an input/output state-space system, where the input corresponds to the forcing (nonlinear terms) and the output to the three velocity components. The idea is then to design a state observer, i.e. a linear operator that permits the estimation of the quantities based on a measurement of the three velocity components at a given height. This is made by means of an optimization problem, which enables the consideration of noise in the measurement sensors. The size of this problem is of the same order as the discretized Orr–Sommerfeld squire (OSS) operator, which could correspond to a square matrix of a few hundred lines. The computation of the model is made ‘off line’, without access to the simulated/experimental quantities. An eddy viscosity is added to the linear operator, and therefore a model is used for the determination of its behaviour along the wall-normal direction. Once the observer has been designed, the prediction is made by measuring the three velocity components at a given wall-normal plane and integrating the observer equations. The Fourier transform applied to the homogeneous directions implies that such measurements correspond to a whole plane in the span and streamwise directions, which is of course more difficult to accomplish experimentally than the proposal of the current work, in which we consider a spanwise line at a fixed streamwise position.

Baars et al. (Reference Baars, Hutchins and Marusic2016) proposed a modification of the inner/outer interaction model by Marusic et al. (Reference Marusic, Mathis and Hutchins2010), where the superposition term, which was originally calculated by a linear operation over the input signal, is replaced by a linear stochastic estimator (LSE). Their flow case is a ZPG TBL of $Re_{\unicode[STIX]{x1D70F}}$ up to 19 000. The LSE is an entity calculated in the frequency domain and it requires measurements of the input and output quantities which are Fourier transformed to build a frequency-domain transfer function. The model is built using time measurements only, which is a configuration simpler to implement experimentally than the one proposed in this study. However, note that while we aim at predicting instantaneous quantities, the work of Baars et al. (Reference Baars, Hutchins and Marusic2016) is only focused on turbulence statistics. Once the LSE is available, a prediction is performed from Fourier-transformed measurements of the input. Since the predictions are obtained in the frequency domain and then inverse Fourier transformed, certain characteristics such as the causality of the prediction cannot be evaluated. The second part of the inner/outer interaction model is unchanged and its calculation involves an iterative procedure to find a demodulated ‘universal signal’ which ensures that the statistics of the resulting predicted signal are recovered. This also means that the full model, which exhibits very good agreement with experimental data, cannot be used for unsteady (time-domain) predictions, such as those necessary for control.

Finally, the method proposed here could be considered as an alternative to these two methods, depending on the desired application. Our results in terms of the predicted statistics are expected to be inferior to those of Baars et al. (Reference Baars, Hutchins and Marusic2016), with a moderately higher computational cost, since the method presented in that work avoids the use of a two-dimensional convolution. However, for time-domain predictions, the method outlined here is expected to outperform the one of Baars et al. (Reference Baars, Hutchins and Marusic2016), particularly when it is enhanced with multiple inputs which are complementary in terms of their coherent content with the output, a feature which, to the best of the authors’ knowledge, is not available in any other model in the literature for this particular application. The current method produces correlation coefficients on the same level as those of Illingworth et al. (Reference Illingworth, Monty and Marusic2018), with an expected lower computational cost and a simpler scheme to implement in an experimental application. The characteristics of the methods by Baars et al. (Reference Baars, Hutchins and Marusic2016) and Illingworth et al. (Reference Illingworth, Monty and Marusic2018), together with those of the method proposed in this study, are summarized in table 2.