2 Description of the methodology

Figure 1 illustrates schematically the construction of the model employed to estimate the unsteady forces. The time mean flow $\boldsymbol{u}_{0}(\boldsymbol{x})$ and a pointwise measurement $\boldsymbol{u}(\boldsymbol{x}_{0},t)$ of the velocity history are the inputs, corresponding to the leftmost blocks. In principle, mean flow and probe information could be obtained independently from either experiments or simulations. A spectral analysis of the probe signal $\boldsymbol{u}(\boldsymbol{x}_{0},t)$ identifies the active frequencies $\unicode[STIX]{x1D714}_{i}$ to be explored in the resolvent analysis of the mean flow. The dominant resolvent modes $\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D714}_{i},1}$ arising from the resolvent analysis corresponding to the active frequencies are calibrated with the probe signal to obtain the amplitude coefficients $a_{\unicode[STIX]{x1D714}_{i},1}$ . A linear combination of the weighted resolvent modes then provides an approximation of the fluctuating velocity and associated pressure fields.

Figure 1.

Diagram of the construction of the model. The mean flow and pointwise measurement inputs are in the leftmost blocks. A fast Fourier transform (FFT) of the probe signal provides the active frequencies $\unicode[STIX]{x1D714}_{i}$ to be explored in the resolvent analysis of the mean flow. The dominant resolvent modes $\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D714}_{i},1}$ corresponding to the active frequencies are calibrated with the probe signal to obtain the amplitude coefficients $a_{\unicode[STIX]{x1D714}_{i},1}$ . A linear combination of the weighted resolvent modes provides an approximation of the fluctuating velocity and pressure. The fluctuating pressure is included in the fluctuating velocity vector $\boldsymbol{u}=(u,v,w,p)^{\text{T}}$ for convenience.

3 Obtaining resolvent modes using time stepping

The most computationally demanding step in the creation of the model, depicted in figure 1, is obtaining the resolvent modes. As such, efficient methods to evaluate the resolvent operator are essential for the feasibility of the present methodology for the estimation of aerodynamic forces. We note from (2.4) that the resolvent operator can be evaluated via a shift and inversion of the Jacobian. Hence, matrix-free or matrix-forming methods, typically employed for global stability analysis (Paredes et al. Reference Paredes, Hermanns, Le Clainche and Theofilis2013), can be used to study the resolvent operator.

Although the simplest way to obtain the resolvent modes numerically is to assemble the resolvent operator and perform an SVD, this is difficult in practice due to the massive computational requirements resulting from the large dimensionality of the operator associated with flows with three non-homogeneous spatial directions. Thus, matrix-free methods are preferred for the study of these flows. In what follows, we describe a matrix-free methodology based on time stepping to obtain resolvent modes, potentially able to cope with three-dimensional flows.

The main idea of matrix-free methods is that the singular values of the resolvent operator are the eigenvalues of $\mathscr{H}_{\unicode[STIX]{x1D714}}\mathscr{H}_{\unicode[STIX]{x1D714}}^{\ast }$ , and the response and forcing singular vectors correspond to the eigenvectors of $\mathscr{H}_{\unicode[STIX]{x1D714}}\mathscr{H}_{\unicode[STIX]{x1D714}}^{\ast }$ and $\mathscr{H}_{\unicode[STIX]{x1D714}}^{\ast }\mathscr{H}_{\unicode[STIX]{x1D714}}$ respectively. As such, the action of the resolvent operator and its conjugate transpose on a forcing vector could enable a matrix-free iterative power method to obtain the SVD of $\mathscr{H}_{\unicode[STIX]{x1D714}}$ .

Following the rank-1 hypothesis posed in § 2.1, only the first singular vectors at each active frequency are required for the construction of the reduced-order model in (2.8). In this context, Monokrousos et al. (Reference Monokrousos, Åkervik, Brandt and Henningson2010) and Lu & Papadakis (Reference Lu and Papadakis2014) showed that obtaining the dominant singular vectors of the resolvent is equivalent to finding the optimal harmonic forcing of the forced linearized Navier–Stokes equations

The long-time integration of (3.1) leads to the harmonic relation (2.3), if all transient effects vanish. Consequently, the optimal forcing and response vectors of $\mathscr{H}_{\unicode[STIX]{x1D714}}$ can be obtained from an optimization problem in a time-stepping context. A Lagrange function for the optimal forcing in (3.1) can be constructed as

where the objective function is the fluctuation energy represented by the energy norm $(\boldsymbol{u}_{\unicode[STIX]{x1D714}},\boldsymbol{u}_{\unicode[STIX]{x1D714}})$ . The first Lagrange multiplier $\boldsymbol{v}_{\unicode[STIX]{x1D714}}$ enforces that the response and forcing satisfy the forced linearized Navier–Stokes equations (3.1), and the second Lagrange multiplier $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}}^{2}$ enforces a unit energy norm to the optimal forcing vector. Variations of $\boldsymbol{L}$ with respect to the fluctuating velocity $\boldsymbol{u}_{\unicode[STIX]{x1D714}}$ and to the forcing $\boldsymbol{f}_{\unicode[STIX]{x1D714}}$ yield respectively

The combination of the outcome of the optimization problem (3.3) and (3.4) with (2.3) leads to the eigenvalue problem

where the fact that the optimal forcing solution corresponds to the first (most amplified) forcing mode $\unicode[STIX]{x1D74D}_{\unicode[STIX]{x1D714},1}$ has been exploited. An iterative matrix-free power method can be then applied to the eigenvalue problem (3.5) in order to obtain the dominant eigenvalue. However, as described by Monokrousos et al. (Reference Monokrousos, Åkervik, Brandt and Henningson2010), time stepping does not directly provide a solution of (2.3) or (3.3), and the forced linearized and adjoint equations need to be integrated for long enough such that the transient dynamics vanish and the response is harmonic. It should be noted that the long integration of the adjoint of (3.1),

provides a harmonic solution for (3.3).

Although the algorithm originally proposed by Monokrousos et al. (Reference Monokrousos, Åkervik, Brandt and Henningson2010) is successful in the sense that it is matrix-free, a number of difficulties arise in practice. The time integration required of (3.1) until a harmonic response is obtained can be very long; hence, the method may not be computationally affordable. The method is particularly slow at low frequencies. Moreover, the complex velocity response $\boldsymbol{u}_{\unicode[STIX]{x1D714}}$ is obtained via a Fourier transform of $\boldsymbol{u}(t)$ during one period; hence, it could have a different phase at each iteration step. Finally, the response of (3.1) can be susceptible to large transient growths depending on the initial condition. Even if $\boldsymbol{u}_{\unicode[STIX]{x1D714}}$ and $\boldsymbol{f}_{\unicode[STIX]{x1D714}}$ are close to convergence, the transient can be significant if the two vectors do not have their relative phase fixed by the resolvent.

A modification of the time-stepping approach of Monokrousos et al. (Reference Monokrousos, Åkervik, Brandt and Henningson2010) that is able to cope with the abovementioned limitations is proposed by allowing the velocity field to take on complex values in the forced direct and adjoint linearized Navier–Stokes equations. The use of a complex variable has the advantage that it is not necessary to perform a Fourier transform of the response during one period to obtain $\boldsymbol{u}_{\unicode[STIX]{x1D714}}$ . Once the response of (3.1) is harmonic, any snapshot of the complex vector $\boldsymbol{u}^{\prime }(\boldsymbol{x},t)$ represents its Fourier transform $\boldsymbol{u}_{\unicode[STIX]{x1D714}}$ at some particular phase. Furthermore, if the time integration is taken as an integer number $n$ of periods $T=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$ , the response and force vectors obtained remain with the same relative phase imposed by the resolvent during a complete iteration. This feature avoids possible sources of transient effects during iterations.

In addition, the integration time can be limited to a small number of periods $n$ such that the strongest transients have vanished but the flow is not yet exactly periodic. This permits an estimation of the amplification value and of the response and forcing vectors in relatively short integrals. The estimated $\unicode[STIX]{x1D70E}$ is then employed between each direct and adjoint iteration to successively generate a better initial condition at each iteration. As a result, the transient effects are damped and the harmonic response is achieved after a few iterations.

All of these features permit the most amplified modes to be obtained using a smaller computational effort than with the original algorithm of Monokrousos et al. (Reference Monokrousos, Åkervik, Brandt and Henningson2010). The modified algorithm is listed in table 1. In practice, the harmonic response is also assessed by monitoring the fluctuation energy. If the difference between the maximum and minimum fluctuation energies is less than a given tolerance within one period, the response is considered to be harmonic and the integration stops. Thus, each iteration may require a number of integration periods less than  $n$ .

Table 1.

Time-stepping algorithm for computing resolvent modes in non-trivial domains.

4 Application to the flow around an inclined square cylinder

The flow past a two-dimensional inclined square cylinder may serve as a model for the flow around a non-trivial bluff body, and it is a good compromise between computational affordability and complex flow features (Sohankar, Norberg & Davidson Reference Sohankar, Norberg and Davidson1998; Yoon, Yang & Choi Reference Yoon, Yang and Choi2010). Beyond the critical Reynolds number, the wake becomes unsteady, presenting asymmetric vortex shedding; thus, this flow is interesting for the investigation of unsteady lift and drag forces.

Although the present methodology seems to be more appealing for experimental works, the mean flow and pointwise measure inputs to the model have been obtained via direct numerical simulation (DNS) using a spectral-element solver (Blackburn & Sherwin Reference Blackburn and Sherwin2004). A rectangular computational domain defined in $[-16,20]\times [-14,14]$ has been discretized with 236 spectral elements. The square cylinder has a unit side length and its centroid is located at $(x,y)=(0,0)$ . Temporal and spatial convergence has been achieved with a polynomial expansion of order 11 in each element and using a second-order temporal scheme with $\unicode[STIX]{x0394}t=8.5\times 10^{-3}$ . A constant velocity $(u,v)=(\cos \unicode[STIX]{x1D6FC},\sin \unicode[STIX]{x1D6FC})$ is imposed at the inlet of the domain, a no-slip boundary condition at the cylinder wall and Neumann boundary conditions at the outlet. The Reynolds number based on the cylinder side length $D$ and the modulus of the inlet velocity is fixed to $Re=100$ . The angle of attack is set to $\unicode[STIX]{x1D6FC}=10^{\circ }$ .

The inputs to the model corresponding to the leftmost block in figure 1 obtained via DNS, mean flow and a single pointwise measurement of the velocity are shown in figure 3. The probe is located at a random point in the wake of the cylinder where the shear is non-zero, and it serves to identify a single dominant frequency $\unicode[STIX]{x1D714}_{1}=0.908$ ; hence, only one resolvent mode corresponding to that frequency needs to be computed. The dynamics of the self-interaction of this mode in this kind of flow could be relevant (Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003); thus, it may be convenient to also consider the first harmonic $\unicode[STIX]{x1D714}_{2}=2\unicode[STIX]{x1D714}_{1}$ , although in the present case the dynamics of the first harmonic does not show any significant influence on the aerodynamic forces, as will be shown.

The iterative time-stepping algorithm described in § 3 has been employed to obtain the two resolvent modes required to construct the reduced-order model. The forced linearized Navier–Stokes equations and their adjoint version have been solved with the same spectral-element solver as employed for DNS. The boundary conditions for each equation are described by Barkley, Blackburn & Sherwin (Reference Barkley, Blackburn and Sherwin2008); however, two sets of boundary conditions at the inlet have been tested for the forced adjoint equations: (i) an extended domain $[-40,20]\times [-14,14]$ with Dirichlet boundary conditions and (ii) the same domain as employed for the DNS with a non-physical forcing $-m(\boldsymbol{x})\boldsymbol{u}(\boldsymbol{x},t)$ applied at $-16<x<12$ to force a zero amplitude of the forcing mode at $x=-16$ . Both boundary conditions have provided similar results; thus, the former has been adopted on account of smaller computational requirements.

A comparison of the vorticity fields obtained from DNS and the resolvent-based model via calibration of the amplitude against the probe data in figure 3(b) is shown is figure 4. It is remarkable that the flow in the region around the cylinder, and hence the unsteady separation, is accurately predicted by the model, despite the probe being located approximately three side lengths from the cylinder. On the other hand, the structure of the wake far from the cylinder presents obvious discrepancies. Although the general feature of vortex shedding is also reproduced, the DNS presents additional features that the present resolvent-based model does not capture, like a consecutive and opposite vertical displacement of the vortex cores. An inspection of the spatial structure of the forcing mode that drives the resolvent modes, shown with blue isolines in figure 3(a), reveals that their maxima are located within the boundary layer of the cylinder; thus, the discrepancies between the model and the DNS can be attributed to this observation. As such, we presume that the additional dynamics of the far wake is governed by additional subdominant resolvent modes associated with the shear in the wake not considered in the present model.

Figure 3.

Inputs to the model corresponding to the leftmost block in figure 1. (a) Streamwise velocity contours of the mean flow, coloured from black to white. Blue contour lines indicate 40 % and 80 % of the maximum of the kinetic energy corresponding to $\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D714}_{1},1}$ . (b) Evolution of the fluctuating streamwise velocity in the wake at the location $\boldsymbol{x}_{0}=(3,1.5)$ , measured from the DNS and fitted to the model. The probe position is highlighted in (a) with a yellow circle.

Figure 4.

Comparison of the fluctuating vorticity fields obtained from (a) DNS and (b) the resolvent-based model via calibration of the amplitude against the probe data in figure 3(b). The coloured contours represent $\pm 1/3$ of the maximum and minimum vorticity.

Figure 5.

(a) Comparison of the temporal evolutions of the unsteady lift and drag coefficients calculated via DNS and predicted by the present methodology. (b) Representation of the fluctuating pressure distribution along the sides of the cylinder at a random instant calculated via DNS and predicted by the present methodology. Each cylinder side acts as an $x$ -axis, while its normal direction indicates the relative value of the unsteady pressure at that location.

The present approach is validated by measuring the unsteady lift and drag coefficients, defined as

where $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{c}$ denotes the boundary of the cylinder, $\boldsymbol{n}$ is the normal vector around the square cylinder, while $\boldsymbol{e}_{x}$ and $\boldsymbol{e}_{y}$ are the unit vectors in the streamwise and normal directions respectively. At the present value of $Re$ , the contributions of the viscous forces are negligible. Nevertheless, they could be taken into account using the present methodology. A good agreement between the temporal evolutions of the unsteady lift and drag coefficients calculated via DNS and predicted by the present methodology is shown in figure 5(a). This result is consistent with the accurate prediction by the model of the near field around the cylinder shown in figure 4. To provide further insight, the fluctuating pressure distribution along the sides of the cylinder at a random instant is represented in figure 5(b). The resemblance between the pressure distributions corresponding to the DNS and the resolvent-based model supports the good agreement between the temporal evolutions of the unsteady lift and drag forces. Finally, the same approach has been carried out using different pointwise measurements. As remarked by Gómez et al. (Reference Gómez, Blackburn, Rudman, Sharma and McKeon2016), similar results are obtained provided that the probe is always positioned at a location where the fluctuating velocity is significant. As such, care is needed in the selection of the probe locations. For example, a probe far from the wake would present a negligible fluctuating velocity or a probe in the wake centreline far from the cylinder would only show dynamics corresponding to the first harmonic frequency $\unicode[STIX]{x1D714}_{2}$ . However, this potential pitfall can be easily overcome by selecting more than one probe location, as shown by Gómez et al. (Reference Gómez, Blackburn, Rudman, Sharma and McKeon2016). Furthermore, the locations where the fluctuating velocity is significant can be inferred from the spatial structure of the resolvent modes.