Triglobal resolvent analysis of swept-wing wakes

Abstract

Through triglobal resolvent analysis, we reveal the effects of wing tip and sweep angle on laminar separated wakes over swept wings. For the present study, we consider wings with semi-aspect ratios from $1$ to $4$, sweep angles from $0^\circ$ to $45^\circ$ and angles of attack of $20^\circ$ and $30^\circ$ at a chord-based Reynolds number of $400$ and a Mach number of $0.1$. Using direct numerical simulations, we observe that unswept wings develop vortex shedding near the wing root with a quasi-steady tip vortex. For swept wings, vortex shedding is seen near the wing tip for low sweep angles, while the wakes are steady for wings with high sweep angles. To gain further insights into the mechanisms of flow unsteadiness, triglobal resolvent analysis is used to identify the optimal spatial input–output mode pairs and the associated gains over a range of frequencies. The three-dimensional forcing and response modes reveal that harmonic fluctuations are directed towards the root for unswept wings and towards the wing tip for swept wings. The overlapping region of the forcing–response mode pairs uncovers triglobal resolvent wavemakers associated with self-sustained unsteady wakes of swept wings. Furthermore, we show that for low-aspect-ratio wings optimal perturbations develop globally over the entire wingspan. The present study uncovers physical insights on the effects of tip and sweep on the growth of optimal harmonic perturbations and the wake dynamics of separated flows over swept wings.

1. Introduction

Understanding flow separation over finite swept wings is essential to the study of aircraft and biological flight (Videler, Stamhuis & Povel 2004; Lentink et al. 2007; Anderson 2010). The aspect ratio, angle of attack and sweep play important roles in influencing stall and wake characteristics (Zhang et al. 2020a,a). Although a number of studies have deepened our knowledge of laminar separated wakes around swept wings, coherent flow structures associated with three-dimensional ($3$-D) flow separation have not been characterized in a comprehensive manner. Such findings would be crucial to explain the role played by the perturbations in characterizing the wakes and support efforts to control flow separation around finite wings.

Previous studies have shown the effect of sweep on post-stall wake characteristics with focus on the role of spanwise flow over wings (Harper & Maki 1964). The spanwise flow induced by sweep delays the emergence of stall (Yen & Hsu 2007; Yen & Huang 2009) and reduces wake oscillations, as shown for high-Reynolds-number flows over transonic buffets in biglobal (Crouch, Garbaruk & Strelets 2019; Paladini et al. 2019; Plante et al. 2021) and triglobal (Timme 2020; He & Timme 2021) linear stability analysis. Similar observations have been made for flows around aircraft models in experiments (Masini, Timme & Peace 2020) and computations (Houtman, Timme & Sharma 2022).

At a low Reynolds number, direct numerical simulations (DNS) from Zhang et al. (2020a) showed that sweep angle can significantly alter the wake patterns. For wings with low sweep angles, vortex shedding develops near the wing tip, while unsteadiness is suppressed for flows over highly swept wings. Similar attenuation of flow unsteadiness was further studied for forward-swept wings (Zhang & Taira 2022) revealing that wing sweep has a strong effect on attenuating wake oscillations in laminar flows. Furthermore, linear instabilities around swept wings were examined for a variety of swept and unswept wings, showing that the sweep angle suppresses the emergence of wake modes (Burtsev et al. 2022; Ribeiro et al. 2022a).

The aspect ratio of the wing also affects the wake dynamics on separated flows due to the wing tip vortex in steady (Devenport et al. 1996; Torres & Mueller 2004; Taira & Colonius 2009) and unsteady wing motion (Buchholz & Smits 2006; Yilmaz & Rockwell 2012). For low-aspect-ratio wings, the tip vortex may suppress leading-edge vortex formation, reducing the wake unsteadiness (Taira & Colonius 2009). Tip vortices can also produce adverse effects on the wing, with induced drag and a reduced lift.

To alter the wake dynamics with a proper actuation input, we need to identify the optimal forcing structures that can be amplified in the flow field (Edstrand et al. 2018a,b). For this task, we may use modal analysis techniques (Taira et al. 2017,  2020) to study the dynamics of flow oscillations. Resolvent analysis is an attractive tool for the present study because it identifies the optimal input perturbations in the flow field, their energy amplification and the characteristics of their unsteady response (Trefethen et al. 1993; Jovanović & Bamieh 2005). Furthermore, with the diverse steady and unsteady wakes observed around swept wings, resolvent analysis can provide a comprehensive study of the input–output dynamics around wings with different aspect ratios, angles of attack and sweep.

Resolvent analysis has been used to study a broad range of fluid flows (Moarref et al. 2013; Schmidt et al. 2018; Thomareis & Papadakis 2018; Skene & Schmid 2019; Yeh et al. 2020; Ricciardi, Wolf & Taira 2022). This approach was initially formulated for steady base flows, to identify modal structures that can be amplified in stable flow regimes (Trefethen et al. 1993). This perspective on fluid dynamics was later extended to unstable systems by Jovanović & Bamieh (2005) and to unsteady and turbulent flows by McKeon & Sharma (2010). In these formulations, a time-averaged flow is used as a base state and nonlinear terms act as sustained forcing in the flow field. In both steady and unsteady flows, resolvent analysis identifies harmonic forcings that produce an amplified response in the flow.

In this study, we identify the optimal spatial input–output modes around the wing through a $3$-D global (triglobal) resolvent analysis, which assumes no spatial homogeneity. Moreover, we gain insights into the self-sustained fluctuations that support unsteadiness on laminar separated flows using resolvent wavemakers, which are similar in spirit to eigenvector-based wavemakers (Giannetti & Luchini 2007; Giannetti, Camarri & Luchini 2010). The resolvent wavemakers, also named as structural sensitivity, are obtained from the overlap of forcing and response modes (Qadri & Schmid 2017; Skene et al. 2022b). These findings provide a comprehensive analysis of the energy amplification mechanisms in flows around swept wings through an input–output process, identifying the optimal locations where perturbations can be introduced to alter the wake behaviour. Therefore, these findings are crucial for the development of efficient flow control strategies (Yeh & Taira 2019; Liu et al. 2021) that aim to improve the aerodynamic performance of swept wings experiencing massive flow separation.

The present paper on triglobal resolvent analysis is organized as follows. In § 2, we describe the problem set-up for the current work. In § 3, we discuss our main findings from triglobal resolvent analysis. We identify the emergence of wake unsteadiness caused by the overlap of optimal forcing and response modes in the near wake. Perturbations are directed towards the region where vortex shedding takes place. The locations of the optimal forcing and response modes over the wingspan also suggest that wakes of highly swept wings are more resilient to external perturbations. Furthermore, we find that low-aspect-ratio wings limit the growth of perturbations to global modes extending over the entire wingspan. Finally, our conclusions are presented in § 4.

2. Problem set-up

We consider laminar flows over untapered swept wings with NACA 0015 cross-sectional profile, as shown in figure 1. The spatial coordinates are defined with $(x,y,z)$ being the streamwise, transverse and spanwise directions, respectively, with the origin placed at the leading edge of the wing root. The NACA 0015 airfoil geometry is defined on the $(x,y)$ plane. The wingspan is formed by extruding the airfoil profile in the spanwise direction. The semi-aspect ratio is defined through the half-span length $b$ and the chord length $c$ as $sAR = b/c$, with values set between $1 \le sAR \le 4$. For swept wings, the $3$-D computational set-up is sheared in the $x$ direction and the sweep angle is defined between the $z$ direction and the leading edge. In the present work, we consider sweep angles $0^\circ \le \varLambda \le 45^\circ$. The angle of attack, $\alpha = 20^\circ$ and $30^\circ$, is defined between the streamwise direction and the airfoil chord line. To focus on the effects of wing tip and sweep in the wake dynamics, we analyse a half-span model with symmetry boundary conditions imposed at the root plane. The wings have a straight-cut tip and sharp trailing edge. For all flows analysed herein, we define the chord-based Reynolds number $Re_{c} = U_\infty c/\nu = 400$, where $U_\infty$ is the free-stream velocity and $\nu$ is the kinematic viscosity. The free-stream Mach number is set to $M_\infty = U_\infty / a_\infty = 0.1$, where $a_\infty$ is the free-stream speed of sound.

Figure 1. Set-up for finite swept wing simulation. In the grey boxes, the instantaneous flow field for $\alpha = 20^\circ$, $\varLambda = 15^\circ$ and $sAR = b/c = 4$, with $Q = 2$ isosurfaces coloured by instantaneous $u_x$. Mesh coloured by time-averaged $\bar {u}_x$. In the blue box, isosurfaces of the primary response mode with mesh in light grey.

2.1. Direct numerical simulation

We perform DNS with the compressible flow solver CharLES (Khalighi et al. 2011; Brès et al. 2017), which uses a second-order-accurate finite-volume method in space with a third-order-accurate scheme in time. With the origin at the leading edge of the airfoil $(x/c,\ y/c,\ z/c) = (0,0,0)$, the computational domain extends over $(x/c, y/c, z/c) \in [-20,25] \times [-20,20] \times [0,20]$. We build a C-type grid for each angle of attack with $\min ({\rm \Delta} x, {\rm \Delta} y, {\rm \Delta} z)/c = (0.005, 0.005, 0.04)$ applying mesh refinement near the airfoil and in the wake, as shown in figure 1.

Inlet and far-field boundaries are prescribed with Dirichlet boundary conditions $(\rho, u_x, u_y, u_z, p) = (\rho _\infty, U_\infty, 0, 0, p_\infty )$, where $\rho$ is density, $p$ is pressure and $u_x$, $u_y$ and $u_z$ are velocity components in the $x$, $y$ and $z$ directions, respectively. Variables with subscript $\infty$ denote free-stream values. For all set-ups considered herein, the velocity boundary conditions applied at the inlet and far field are aligned with the $x$ direction, which enforces the same streamwise flow over wings with different sweep angles. The airfoil surface is provided with adiabatic no-slip boundary condition. To simulate a half-wing model, we prescribe the symmetry boundary condition along the root plane. A sponge layer is applied at the outlet over $x/c \in [15,25]$ with the target state being the running-averaged flow over five convective time units $t \equiv c / U_\infty$ (Freund 1997). Simulations start with uniform flow and time integration is performed with a constant acoustic Courant–Friedrichs–Lewy number of $1$. After transients are flushed out of the computational domain, the time-averaged base flow $\bar {\boldsymbol {q}}$ is determined over $50$ convective time units. The present results were carefully verified and validated. Close agreement for instantaneous and time-averaged velocity components was achieved with those from Zhang et al. (2020a). We have further validated our computations for time-averaged drag and lift coefficients,

(2.1a,b)\begin{equation} C_D = \frac{F_x}{\dfrac{1}{2} \rho U_\infty^2 b c} \quad \text{and} \quad C_L = \frac{F_y}{\dfrac{1}{2} \rho U_\infty^2 b c}, \end{equation}

respectively, where $F_x$ is the drag and $F_y$ is the lift over the wing, as reported in table 1.

Table 1. Time-averaged lift and drag coefficients ($\overline {C_L}$ and $\overline {C_D}$) compared with those of Zhang et al. (2020a) for laminar separated flow over NACA 0015 wings with $sAR = 4$, $\alpha = 20^\circ$ and $\varLambda = 0^\circ$, $15^\circ$, $30^\circ$ and $45^\circ$.

A variety of wake patterns can be observed for different $\alpha$, $\varLambda$ and $sAR$, as summarized in figure 2. In figure 2(b), the flow over the wing with $(sAR,\alpha,\varLambda ) = (2,30^\circ,0^\circ )$ exhibits a quasi-steady streamwise-oriented tip vortex. This structure is characteristic of flows over unswept wings and also appears around wings with different $\alpha$ and $sAR$. For such wings, unsteady spanwise vortices develop at the root plane. Between the root and the wing tip, there is an intermediate zone with braid-like vortices.

Figure 2. Instantaneous isosurfaces of $Q = 2$ coloured by $u_x$ for (a) $\alpha = 20^\circ$ and (b) $\alpha = 30^\circ$. Unsteady shedding near wing root ($\boldsymbol {\diamond }$, red), unsteady shedding near wing tip ($\boldsymbol {\square }$, yellow), steady flow with root structures ($\boldsymbol {\blacktriangleright }$, green), steady flow with streamwise vortices ($\boldsymbol {\blacktriangledown }$, blue).

Wing sweep affects the wake dynamics and structures. At low sweep angles, a spanwise flow develops over the wing and advects unsteady vortices towards the wing tip. For instance, for $(sAR,\alpha,\varLambda ) = (4,20^\circ,15^\circ )$, spanwise vortices still appear. These structures are similar to the ones observed over unswept wings, although they form closer to the wing tip, and break into helical structures in the wake. For the flow around this wing, streamwise-oriented tip vortices are absent.

For higher $\varLambda$, wing sweep can attenuate wake oscillations. For $(sAR,\alpha,\varLambda ) = (4,30^\circ,30^\circ )$, near-wake unsteadiness is reduced and unsteady vortices appear further downstream in the wake. We notice that these structures are absent when the angle of attack is lowered to $20^\circ$. Similarly, increasing the sweep angle to $\varLambda = 45^\circ$ suppresses unsteady vortices on both angles of attack and the wake becomes steady. On such highly swept wings, ram-horn-shaped streamwise-oriented vortices develop from the root plane and extend into the wake.

For each $(\alpha,\varLambda )$ pair, the wake exhibits similar characteristics for wings with $sAR = 4$ and $2$. Reducing the semi-aspect ratio to $sAR = 1$ has a strong influence on the wake dynamics, as shown in figure 2. For such wings, tip effects can suppress the formation of leading-edge vortices at lower angles of attack. For instance, at $\alpha = 20^\circ$, wake unsteadiness is reduced and swept wings exhibit steady flows with root structures.

Unsteady vortices are observed in flows over $sAR = 1$ wings at $\alpha = 30^\circ$ for all considered sweep angles. The unsteadiness appears near the root for lower $\varLambda$. Further downstream, unsteady vortices appear over the entire wingspan. For higher $\varLambda$, vortices are generated near the wing tip and helical structures are observed in the wake. These observations agree with the characterizations by Zhang et al. (2020a). To deepen our insights into swept-wing wake dynamics, we now call for triglobal resolvent analysis.

2.2. Resolvent analysis

Let us consider the Reynolds decomposition of state variable $\boldsymbol {q} = \bar {\boldsymbol {q}} + \boldsymbol {q}^\prime$, where $\bar {\boldsymbol {q}}$ is the time-averaged flow and $\boldsymbol {q}'$ is the statistically stationary fluctuation component (McKeon & Sharma 2010). This decomposition along with spatial discretization is used to linearize the compressible Navier–Stokes equations about $\bar {\boldsymbol {q}}$ to yield

(2.2)\begin{equation} \frac{\partial \boldsymbol{q}^\prime}{\partial t} = \boldsymbol{\mathsf{L}}_{\bar{\boldsymbol{q}}}\boldsymbol{q}^\prime + \boldsymbol{f}^\prime, \end{equation}

where $\boldsymbol{\mathsf{L}}_{\bar {\boldsymbol {q}}}$ is the discrete linearized Navier–Stokes operator (Sun et al. 2017) and $\boldsymbol {f}^\prime$ accounts for the external forcing and nonlinear terms. With the Fourier representation

(2.3)\begin{equation} \left[\boldsymbol{q}^\prime(\boldsymbol{x},t), \boldsymbol{f}^\prime(\boldsymbol{x},t) \right] = \int_{-\infty}^{\infty} \left[ \hat{\boldsymbol{q}}_{\omega}(\boldsymbol{x}), \hat{\boldsymbol{f}}_{\omega}(\boldsymbol{x}) \right] \exp({-{\rm i} \omega t})\,{\rm d} \omega, \end{equation}

we obtain

(2.4)\begin{equation} -{\rm i}\omega \hat{\boldsymbol{q}}_{\omega} = \boldsymbol{\mathsf{L}}_{\bar{\boldsymbol{q}}} \hat{\boldsymbol{q}}_{\omega} + \hat{\boldsymbol{f}}_{\omega}, \end{equation}

where $\boldsymbol {x} = (x,y,z)$ and the triglobal response and forcing modes are $\hat {\boldsymbol {q}}_{\omega }$ and $\hat {\boldsymbol {f}}_{\omega }$, respectively, for a temporal frequency $\omega$. This expression leads to

(2.5)\begin{equation} \hat{\boldsymbol{q}}_{\omega} = \boldsymbol{\mathsf{H}}_{\bar{\boldsymbol{q}},\omega} \hat{\boldsymbol{f}}_{\omega}, \end{equation}

in which the resolvent operator $\boldsymbol{\mathsf{H}}_{\bar {\boldsymbol {q}},\omega } \in \mathbb {C}^{m \times m}$, with $m$ defined by the product of the number of state variables and the number of spatial grid points. For the present triglobal base flows, the linear operators have size $m$ between $3$ and $5 \times 10^6$. We analyse the resolvent operator with the singular value decomposition (SVD)

(2.6)\begin{equation} \boldsymbol{\mathsf{H}}_{\bar{\boldsymbol{q}}} = \left[ -{\rm i} \omega \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{L}}_{\bar{\boldsymbol{q}}} \right]^{{-}1} = \boldsymbol{\mathsf{Q}} \boldsymbol{\varSigma} \boldsymbol{\mathsf{F}}^*, \end{equation}

where $\boldsymbol{\mathsf{F}} = [\hat {\boldsymbol {f}}_1,\hat {\boldsymbol {f}}_2,\ldots,\hat {\boldsymbol {f}}_m]$ is an orthonormal matrix holding the forcing modes, $\boldsymbol {\varSigma } = {\rm diag}[\sigma _1,\sigma _2,\ldots,\sigma _m]$ is the diagonal matrix with singular values (gain) in descending order and $\boldsymbol{\mathsf{Q}} = [\hat {\boldsymbol {q}}_1,\hat {\boldsymbol {q}}_2,\ldots,\hat {\boldsymbol {q}}_m]$ is the orthonormal matrix comprised of the response modes (Trefethen et al. 1993; Jovanović & Bamieh 2005). For visualization purposes, we show only the real part of the complex-valued resolvent modes. Here, we employ the Chu norm (Chu 1965) incorporating it within the resolvent operator through a similarity transform $\boldsymbol{\mathsf{H}}_{\bar {\boldsymbol {q}}} \rightarrow \boldsymbol{\mathsf{W}}^{{1}/{2}} \boldsymbol{\mathsf{H}}_{\bar {\boldsymbol {q}}} \boldsymbol{\mathsf{W}}^{-({1}/{2})}$, where $\boldsymbol{\mathsf{W}}$ is the weight matrix that accounts for numerical quadrature and energy weights.

Resolvent analysis requires careful consideration of the eigenvalues of $\boldsymbol{\mathsf{L}}_{\bar {\boldsymbol {q}}}$. In the presence of unstable modes in the linear operator eigenspectrum, the asymptotic input–output relationship is buried under the unstable dynamics behaviour. The present resolvent analysis utilizes a time-averaged flow as the base state. Since such flow is not the equilibrium state, stability characterization cannot be performed in a strict sense. However, it is important to check the location of the eigenvalues in the complex plane to capture the growth rate of the most unstable modes of $\boldsymbol{\mathsf{L}}_{\bar {\boldsymbol {q}}}$.

To use resolvent analysis to study the wake dynamics of unstable base flows, we examine the dynamics through the lens of temporal discounting (Jovanović 2004). Discounting applies a temporal damping on forcing and response modes as $[\hat {\boldsymbol {q}}_{\omega },\hat {\boldsymbol {f}}_{\omega }]\,{\rm e}^{-\beta t}$, where $\beta >0$ is a time-discounting parameter defined within the discounted resolvent operator (Jovanović 2004). With the discounted resolvent analysis, we can examine amplification dynamics that takes place on a time scale shorter than that of the most unstable mode. Detailed discussions of our choice of $\beta$ are provided in Appendix A. Through the discounted resolvent analysis, valuable insights have been provided in past studies for the dynamics and control of flows over airfoils (Yeh & Taira 2019; Yeh et al. 2020; Ribeiro et al. 2022b; Ricciardi et al. 2022).

The $\boldsymbol{\mathsf{H}}_{\bar {\boldsymbol {q}},\omega }$ operators were discretized over $3$-D structured grids with the leading edge at the root positioned at $(x/c,\ y/c,\ z/c) = (0,0,0)$, extending over $(x/c, y/c, z/c) \in [-10,15] \times [-10,10] \times [0,10]$ with near-wake grids shown at the bottom left of figure 1. The computational grids used for resolvent analysis have a smaller domain size than those used for DNS. For the base flow on the mesh for resolvent analysis, we perform a linear interpolation of the flow field from DNS mesh to resolvent mesh. We prescribe homogeneous Neumann boundary conditions for $T'$ and homogeneous Dirichlet boundary conditions for the fluctuating variables $\rho '$ and $u'$ along the far field, airfoil surface and outlet. Sponges are applied far from the airfoil and in conjunction with the boundary conditions (Freund 1997).

For the large linear operators in the present work, efficient numerical tools are needed for SVD (Halko, Martinsson & Tropp 2011). We use the randomized resolvent analysis algorithm from Ribeiro, Yeh & Taira (2020), sketching $\boldsymbol{\mathsf{H}}_{\bar {\boldsymbol {q}},\omega }$ with $10$ random test vectors. Each entry of the test vectors is associated with a particular grid point and the five state variables, scaled by $[\| \boldsymbol {\nabla } \rho \|, \|\boldsymbol {\nabla } u_x\|, \|\boldsymbol {\nabla } u_y\|, \|\boldsymbol {\nabla } u_z\|, \|\boldsymbol {\nabla } T\|]$ at each spatial location for each state variable (Ribeiro et al. 2020; House et al. 2022). A convergence analysis of the randomized resolvent algorithm is provided in Appendix B.

The computation of resolvent modes for large linear operators can be challenging for the resolvent analysis of high-Reynolds-number flows that require a large grid. The bottleneck is related to the time and memory requirements of the linear system solvers within the SVD. Building an optimal basis to avoid linear system solvers is possible (Barthel, Gomez & McKeon 2022), although a generalization for complex geometries is still challenging. It is possible, however, to obtain accurate resolvent modes with time-stepping instead of direct solvers. Those methods tend to penalize the computational time costs, although a considerable reduction in memory requirements can be achieved (Barkley, Blackburn & Sherwin 2008; Monokrousos et al. 2010; Gómez et al. 2016). The computational time required by time-steppers can also be reduced by incorporating streaming discrete Fourier transforms (Farghadan et al. 2021; Martini et al. 2021). The use of iterative solvers has shown promising results for computing resolvent modes around a commercial aircraft model (Houtman et al. 2022).

In the present work, the direct and adjoint linear systems were directly solved using the MUMPS (multifrontal massively parallel sparse direct solver) package (Amestoy et al. 2001). Moreover, we incorporate adjoint-based sensitivity analysis to interpolate the resolvent norm over frequencies $\omega$ (Schmid & Brandt 2014; Fosas de Pando & Schmid 2017). This approach is used to calculate the gradient of $\sigma$ with respect to $\omega$, allowing an accurate interpolation among frequencies (Skene & Schmid 2019). The codes used to compute the resolvent modes are part of the ‘linear analysis package’ made available by Skene, Ribeiro & Taira (2022a).

3. Triglobal resolvent analysis

3.1. Forcing and response mode structures

Let us first examine the dominant gains, forcing and response modes for $(sAR,\alpha,\varLambda ) = (4,20^\circ,0^\circ )$, as shown in figure 3. The dominant resolvent modes are observed at $St = 0.14$, where

(3.1)\begin{equation} St = \frac{\omega}{2{\rm \pi}} \frac{c \sin \alpha}{U_\infty \cos \varLambda} \end{equation}

is the Fage–Johansen Strouhal number (Fage & Johansen 1927) with a $1/\cos \varLambda$ scaling that incorporates the influence of the sweep angle. This frequency scaling is inspired by the independence principle (Wygnanski et al. 2011) and collapses the spectral behaviour of the resolvent modes over different sweep angles (Ribeiro et al. 2022b). This frequency matches the peak frequency for the lift coefficient shown in the bottom panel of figure 3(a). The dominant frequency for $\sigma _1$ and $\hat {C}_L$ agrees for all unsteady flows presented herein.

Figure 3. (a) Resolvent gains and (b) forcing–response mode pairs for $(sAR,\alpha,\varLambda ) = (4,20^\circ,0^\circ )$. For each mode, forcing ($\hat {\boldsymbol {f}}$) is the top half while response ($\hat {\boldsymbol {q}}$) it the bottom half with isosurfaces of velocity $\hat {u}_y \in [-0.2,0.2]$, with free stream directed to the right. Bottom panel of (a) shows power spectrum density of lift coefficient $\hat {C}_L$.

The spatial structures of forcing–response mode pairs are shown in figure 3(b) for representative frequencies. For $St = 0.14$, primary modes exhibit modal structures near the root plane. The forcing mode appears near and upstream of the wing, while the response mode develops downstream in the wake. The modal structures for the primary forcing and response modes are aligned with the wingspan, with the response mode similar to the unsteady vortices revealed from DNS. At this frequency, the secondary modes are comprised of spanwise-aligned vortices near the root plane, similar to the primary modes.

As we increase the frequency, the resolvent gains decay in magnitude and $\sigma _1$ decays faster than $\sigma _2$. Their magnitudes become approximately the same at $St = 0.16$. At this frequency, the spatial characteristics of the primary and secondary forcing–response mode pairs exhibit distinct behaviour. The primary forcing and response modes are aligned with the wingspan and near the root plane, similar to those at lower $St$. The secondary modes, however, exhibit modal structures near the wing tip, in contrast to the secondary modes at lower frequencies which reside near the wing root.

For $St = 0.18$, the primary forcing–response mode pair appears near the wing tip, while the secondary mode pair develops at the root plane. Such behaviour persists as we increase the frequency to $St = 0.20$. For $St \ge 0.18$, primary modes are tip-dominated while secondary modes are root-dominated around this wing. This means that root and wing tip modes switch their order of amplification at $St \approx 0.18$, i.e. mode switching.

The mode switching phenomenon is also observed for swept wings with $\varLambda = 15^\circ$. For such wings, root-supported structures appear as the primary forcing–response pairs at $St = 0.14$, as shown in figure 4(a). A distinct mode switching is observed over this wing, as the forcing–response pairs gradually transition towards the root at $z/c \approx 0$ with an increase in $St$. This type of concentrated resolvent mode at the wing root also appears for the unswept wings at $St = 0.20$ as a secondary mode, as shown at the bottom of figure 3(b).

Figure 4. Resolvent gain distribution for the top three mode pairs and forcing–response mode pairs for selected frequencies for $(sAR,\alpha ) = (4,20^\circ )$ and $\varLambda =$ (a) $15^\circ$, (b) $30^\circ$ and (c) $45^\circ$. Primary forcing ($\hat {f}$) and response ($\hat {q}$) modes shown with isosurfaces of velocity $\hat {u}_y \in [-0.2,0.2]$.

For $\varLambda = 30^\circ$, mode switching also occurs towards the root with an increase in $St$, in an opposite trend to that of the unswept wings. The dominant response modes at lower frequencies appear at the wing tip, as shown in figure 4(b). There is a gradual transition to root-supported modes as the frequency increases. At a higher sweep angle, $\varLambda = 45^\circ$, no mode switching occurs. The region of dominance of the forcing and response modes is slightly invariant for the frequencies shown herein.

In contrast with the lower-sweep-angle wings, for $\varLambda = 45^\circ$, forcing and response modes are dominant at distinct wingspan locations, as shown in figure 4(c). Response modes are tip-dominated while forcing structures appear upstream near the root plane, extending over the wingspan aligned with the sweep angle. For all $(sAR,\alpha ) = (4,20^\circ )$ wings, the highest amplification is found for $\varLambda = 30^\circ$, at $St \approx 0.12$. At $\varLambda = 45^\circ$, the dominant gain is an order of magnitude lower. This finding suggests that it is challenging to perturb flows over $\varLambda = 45^\circ$ wings. These wings are steady because self-sustained flow disturbances cannot introduce sufficient energy into the wake to generate vortex shedding.

3.2. Resolvent wavemakers

To characterize the self-sustained unsteadiness in the flows over swept wings, we study the spatial overlap between the forcing and response modes that supports the continuous formation of vortical structures. Since the forcing modes show regions receptive to external perturbations and the response modes reveal the structures being excited due to the forcing, the region over which forcing and response modes overlap can be interpreted as a mechanism for self-sustained oscillations in the flow. This idea is similar to the wavemaker concept deduced from direct and adjoint eigenmodes presented in Giannetti & Luchini (2007).

Through the wavemaker analysis, previous studies identified critical points responsible for sustaining wake shedding on laminar wakes around cylinders (Strykowski & Sreenivasan 1990; Hill 1992) and regions associated with their primary and secondary instability modes (Giannetti & Luchini 2007; Giannetti et al. 2010). Moreover, wavemakers revealed the physical mechanisms responsible for tonal noise generation in high-Reynolds-number flows over airfoils (Fosas de Pando, Schmid & Sipp 2017) and self-sustained flow instabilities in transonic buffet regimes (Paladini et al. 2019).

In the aforementioned studies, wavemakers were derived from direct and adjoint global stability eigenmodes. Our formulation derives wavemakers from global resolvent modes and it is closely related to the structural sensitivity devised by Qadri & Schmid (2017) and to the resolvent wavemaker studied by Skene et al. (2022b). The present resolvent wavemaker is not identical to the eigenvector-based wavemaker. Using the time-averaged base flow, the present forcing terms encapsulate nonlinear effects as an internal feedback mechanism within the flow. Hence, the spatial overlap between forcing and response identifies regions responsible for self-sustained wake oscillations.

Herein, the resolvent wavemaker modes are directly obtained from the resolvent modes, as the Hadamard product of forcing and response modes:

(3.2)\begin{equation} \hat{\boldsymbol{w}} = \hat{\boldsymbol{f}} \circ \hat{\boldsymbol{q}}, \end{equation}

where $\hat {\boldsymbol {w}}$ is the resolvent wavemaker mode. The resolvent modes presented herein are defined with the five state variables, $\hat {\boldsymbol {f}} = [\hat {\boldsymbol {f}}_\rho,\hat {\boldsymbol {f}}_{u_x},\hat {\boldsymbol {f}}_{u_y},\hat {\boldsymbol {f}}_{u_z},\hat {\boldsymbol {f}}_T]$ and $\hat {\boldsymbol {q}} = [\hat {\boldsymbol {q}}_\rho,\hat {\boldsymbol {q}}_{u_x},\hat {\boldsymbol {q}}_{u_y},\hat {\boldsymbol {q}}_{u_z},\hat {\boldsymbol {q}}_T]$. We define our resolvent wavemaker gain $\xi$ as

(3.3)\begin{equation} \xi = \sigma^2 \int_{S} |\hat{\boldsymbol{w}}(\boldsymbol{x})|\,{\rm d}S, \end{equation}

which follows $\xi = \sigma ^2 | \langle \hat {\boldsymbol {f}},\hat {\boldsymbol {q}} \rangle |$, derived by Skene et al. (2022b). This is similar to that presented in Ribeiro et al. (2022b). Qualitatively, both definitions of $\xi$ result in similar discussions and interpretations. The present expression for the resolvent wavemaker gain provides a proper quantitative definition (Skene et al. 2022b). The resolvent wavemaker gain $\xi$ can also be computed for each spanwise slice and each frequency. To this end, we consider $S = S(x,y)$, as $z$-normal planes at different spanwise locations, to build the $\xi$ contours shown in figure 5. Through this analysis, we highlight the spatial support of the resolvent wavemaker over the wingspan.

Figure 5. Wingspan location of primary resolvent wavemakers with isocontours of $\xi$ for $0.05 \le St \le 0.25$ for $(sAR,\alpha ) = (4,20^\circ )$ wings with (a–d) $0^\circ \le \varLambda \le 45^\circ$. Resolvent wavemaker modes at $St = 0.14$ are shown with isosurfaces of $\hat {u}_y / \| \hat {u}_y \|_\infty = \pm 0.1$ and instantaneous $Q=1$ are grey-coloured isosurfaces.

Let us focus our resolvent wavemaker analysis on the flow over the unswept wing with $(sAR,\alpha,\varLambda ) = (4,20^\circ,0^\circ )$, as shown in figure 5(a). At $St=0.14$, triglobal resolvent wavemakers with high $\xi$ appear between $2 \lesssim z/c \lesssim 3$ in the near wake. The resolvent wavemakers at this region support the formation of unsteady root vortices that propagate downstream in the wake. This resolvent wavemaker region is also characterized by the formation of braid-like structures that connect to the root shedding as vortex loops (Zhang et al. 2020a). Resolvent wavemakers for $(sAR,\alpha,\varLambda ) = (4,20^\circ,15^\circ )$ also show similar shedding behaviour, as seen in figure 5(b).

Resolvent wavemakers are also revealed for steady flows. The overlap of forcing and response modes for flows over wings with $(sAR,\alpha,\varLambda ) = (4,20^\circ,30^\circ )$, as shown in figure 5(c), develops over the wing and extends into the wake aligned with the wing tip. These resolvent wavemakers extend over the entire wingspan, being stronger and larger than those exhibited around wings with lower sweep angles. The $\xi$ peak appears at the wing tip at $St = 0.12$, indicating that the tip region is more susceptible to develop unsteadiness around this wing.

For $\varLambda = 45^\circ$, shown in figure 5(d), we reveal that resolvent wavemakers emerge from the leading edge near the root plane towards the wing tip and downstream at the wake, overlapping the region where steady ram-horn-shaped vortices appear in the DNS. These resolvent wavemakers exhibit a region of the flow field with high receptiveness to amplify forcing structures and disturb the steady ram-horn vortex. Because the dominant resolvent wavemakers around the $\varLambda = 45^\circ$ wing have a low $\xi$, in spite of occupying a large region of the wake, the energy they introduce to the flow field is insufficient to disturb the wake.

The resolvent wavemakers further exhibit the root- and tip-dominated modal characteristics and the mode switching phenomenon in figure 5, in agreement with the forcing–response modal behaviour shown in figures 3 and 4. For instance, the resolvent wavemaker modes at the peak $\xi$ values for the unswept wing appear near $z/c \approx 2$, with a gradual transition from root-supported to tip-dominated modes as $St$ increases. Moreover, for the $\varLambda = 15^\circ$ wing, there is a transition in the dominant region of resolvent wavemaker support from $z/c \approx 1.5$ at lower frequencies to $z/c \approx 0$ at higher frequencies, as shown in figure 5(b). Lastly, for the $\varLambda = 30^\circ$ wing, there is a tip-to-root transition with an increase in $St$ while the peak resolvent wavemakers for $\varLambda = 45^\circ$ are invariant over the frequencies, appearing near $z/c \approx 2$.

3.3. Forcing-to-response dynamics

Let us further explain how perturbations emerge around swept wings, by analysing the overlap between the forcing and response modes in the spanwise direction. To this end, we integrate the norm of $\hat {\boldsymbol {f}}$ and $\hat {\boldsymbol {q}}$ over $z$-normal planes, analogous to the resolvent wavemaker mode analysis, as

(3.4a,b)\begin{equation} \boldsymbol{\varOmega}_{\hat{\boldsymbol{f}}}(z) = \int_{S(x,y)} \|\hat{\boldsymbol{f}}\|_2 \,{\rm d}S \quad \text{and} \quad \boldsymbol{\varOmega}_{\hat{\boldsymbol{q}}}(z) = \int_{S(x,y)} \|\hat{\boldsymbol{q}}\|_2 \,{\rm d} S, \end{equation}

where $\|\hat {\boldsymbol {f}}\|_2$ and $\|\hat {\boldsymbol {q}}\|_2$ are the $2$-norm of $[\hat {\boldsymbol {f}}_{\rho },\hat {\boldsymbol {f}}_{u_x},\hat {\boldsymbol {f}}_{u_y},\hat {\boldsymbol {f}}_{u_z},\hat {\boldsymbol {f}}_{T}]$ and $[\hat {\boldsymbol {q}}_{\rho },\hat {\boldsymbol {q}}_{u_x},\hat {\boldsymbol {q}}_{u_y},\hat {\boldsymbol {q}}_{u_z},\hat {\boldsymbol {q}}_{T}]$, respectively, at each grid point of the computational domain. By performing the integral over $S(x,y)$, we obtain $\boldsymbol {\varOmega }_{\hat {\boldsymbol {f}}}$ and $\boldsymbol {\varOmega }_{\hat {\boldsymbol {q}}}$ computed for each spanwise slice and for each frequency. Here, we plot their contours normalized by the maximum $\boldsymbol {\varOmega }_{\hat {\boldsymbol {f}}}$ and $\boldsymbol {\varOmega }_{\hat {\boldsymbol {q}}}$ at each $St$, to emphasize the spatial support of forcing and response over the wingspan, as shown in figure 6 for wings at $\alpha = 30^\circ$ with $0^\circ \le \varLambda \le 45^\circ$. The locations of the maximum strength of forcing and response modes are shown by the dot-dashed lines. Black arrows indicate the direction from the maximum forcing to the maximum response at $St = 0.15$. This analysis depicts the preferential direction in which optimal forcing is transferred to optimal response over the wingspan at each frequency.

Figure 6. Wingspan locations of dominant forcing (red) and response (blue) with contours of $\boldsymbol {\varOmega }_{\hat {\boldsymbol {f}}}$ and $\boldsymbol {\varOmega }_{\hat {\boldsymbol {q}}} \in [0.4,1.0]$ for $0.05 \le St \le 0.35$ for $sAR = 4$ wings with $\alpha = 30^\circ$ with (a–d) $0^\circ \le \varLambda \le 45^\circ$. Dot-dashed lines are polynomial fit of maximum $z/c$ of forcing and response at each $St$. Arrows show direction of optimal forcing to response at $St = 0.15$.

For unswept wings, shown in figure 6(a), the optimal forcing structures appear closer to the wing tip than the response modes, which are slightly shifted towards the root, suggesting that fluctuations are directed towards the root. Indeed, as seen in the DNS, unsteadiness is concentrated towards the root, as evident from figure 2, also in agreement with the results reported by Zhang et al. (2020b). In addition, the flow around the wing tip for unswept wings is characterized by an almost steady tip vortex, suggesting that it is likely hard to amplify flow oscillations near the tip.

For swept wings, fluctuations are directed towards the wing tip. For $\varLambda = 15^\circ$, both forcing and response modes appear near the wing root. At the vortex shedding frequency for this wing, $St \approx 0.15$, we observe forcing and response modes to be dominant at $z/c \approx 1$, with the forcing mode supported closer to the wing root than the response mode. This concurs with the flow field we observe in the DNS, as vortices are formed near the wing root and evolve towards the wing tip where spanwise vortices appear and propagate in the wake. For the $\varLambda = 30^\circ$ wing, the dominant forcing–response mode pair emerges near the wing tip at low $St$, as seen in figure 6(c). For this wing, low-frequency vortical structures emerge downstream in the wake aligned at the tip, as shown in figure 2.

For the $\varLambda = 45^\circ$ wing, the distance between the maximum forcing and response mode locations significantly increases. For this sweep angle, the region of forcing is centred at $z/c \approx 1$, while the response is supported mostly at $z/c \approx 3$. As the peak $\sigma _1$ is lower for this wing compared with lower-sweep-angle planforms, we argue that a significant amount of energy is required for an external forcing to perturb the wakes of highly swept wings. For all wings with $sAR=4$, this distance between the dominant forcing–response mode pairs is strongly associated with the sweep angle, while having a minor dependency on the angle of attack and presenting a gradual decrease with the frequency.

The direction from forcing to response revealed by the optimal triglobal resolvent modes suggests a spanwise advection of flow structures associated with the sweep angle. As shown previously, we can relate the forcing-to-response characteristics to the vortical fluctuations observed in the DNS. We can further relate these findings to the modal convective speed from biglobal stability analysis over swept wings (Crouch et al. 2019; Paladini et al. 2019; Plante et al. 2021). Triglobal resolvent modes also reveal the advection of perturbations over the wingspan related to the sweep angle, the attenuation of flow unsteadiness and the resilience to amplify perturbations at high sweep angles. Even for unswept wings, the triglobal analysis uncovers a preferential root direction for advection of oscillations.

3.4. Influence of the aspect ratio

High sweep angle and low aspect ratio restrict the emergence of fluctuations in flows over finite wings. As shown in figure 6, tip- and root-dominated modes may extend over one or two chord lengths over the wingspan. For this reason, for flows over wings with $sAR < 2$, the dominance of the global modes may not be associated with root or tip regions, as they extend over the entire wingspan.

For flows over $sAR = 2$ wings, we observe a gradual transition between root-dominated and tip-dominated forcing and response modes, as shown in figure 7. For $\varLambda = 15^\circ$, the optimal forcing–response mode pair appears near the root for lower frequencies and at the wing tip for higher frequencies, characterizing a root-to-tip mode switching. For $\varLambda = 30^\circ$, the trend is opposite, with wing tip modes at lower frequencies and root modes at higher frequencies, characterizing a tip-to-root mode switching. These features are similar to the mode switching observed for these sweep angles with $sAR = 4$, as shown in figure 6.

Figure 7. (a–d) Wingspan locations of dominant forcing (red) and response (blue) with contours of $\boldsymbol {\varOmega }_{\hat {\boldsymbol {f}}}$ and $\boldsymbol {\varOmega }_{\hat {\boldsymbol {q}}} \in [0.4,1.0]$ for $0.05 \le St \le 0.35$ for wings at $\alpha = 30^\circ$, $sAR = 2$ and $1$ and $\varLambda = 15^\circ$ and $30^\circ$.

For wings with a low aspect ratio, the growth of root-dominated and tip-dominated perturbations is constrained and mode switching does not occur for $sAR = 1$, as shown in figure 8. Distinguishing between root-dominated and tip-dominated modes may be challenging for flows over $sAR = 1$ wings as forcing and response mode pairs appear globally, extending over the entire wingspan, independent of the sweep angle. Therefore, flows around wings with $sAR = 1$ tend to exhibit similar wake characteristics over different sweep angles. Indeed, the wake patterns for flows over $sAR = 1$ wings at a particular angle of attack and sweep exhibit characteristics different from those of the flows over higher-aspect-ratio wings, e.g. $sAR = 2$ and $4$.

Figure 8. Resolvent gain distribution and forcing–response mode pairs over frequency for $(\alpha,\varLambda ) = (30^\circ,15^\circ )$ and (a–c) $1 \le sAR \le 4$. Forcing ($\hat {f}$) and response ($\hat {q}$) modes shown with isosurfaces of velocity $\hat {u}_y \in [-0.2,0.2]$. Mode switching is absent for $sAR = 1$ due to merging of root and wing tip perturbations on the wake.

For high-aspect-ratio wings, for instance, the flow around $(sAR,\alpha,\varLambda ) = (4,30^\circ,15^\circ )$ wings, we observe in the DNS that the wake shedding structures appear over the entire wingspan. The resolvent modes depict these structures in three different flow mechanisms. As shown in figure 8(c) for $sAR = 4$, there are two types of root-dominated modes, which were also previously identified for this wing at $\alpha = 20^\circ$, shown in figure 4(a). The first one is characterized by root-dominated structures and appears at $St = 0.15$, while the second type, with a high $\sigma _1$, develops at $St = 0.25$ with compact root-concentrated modes. The third type is comprised of tip-dominated modes that become primary as the frequency increases to $St = 0.28$. These modes were primary at $\alpha = 20^\circ$ and $St \approx 0.40$, as shown in figure 4, although for $\alpha = 30^\circ$ they present a higher amplification gain.

For $sAR = 2$, root-dominated modes are primary for $St < 0.20$. Root-concentrated modes are absent and tip-dominated modes are the primary perturbations for $St \ge 0.20$, as shown in figure 8(b), characterizing a root-to-tip mode switching. The overall mode switching for $sAR = 2$ is the same, with root-to-tip transition. For $sAR = 1$, shown in figure 8(a), mode switching is absent. Both primary and secondary modes develop over the entire wingspan for all frequencies, as shown for the primary modes at $St = 0.14$. Although mode switching is absent we can still reveal two distinct root- and tip-dominated mechanisms on a single mode over low-aspect-ratio wings. For instance, at $St = 0.24$, modes emerge from the leading edge at the root and from the trailing edge near the tip. Combined, these two types of flow unsteadiness yield a global mode that appears over the entire wingspan.

4. Conclusions

We presented the triglobal resolvent analysis of laminar separated flows over swept wings and characterized the effects of wing tip and sweep angle on the wake dynamics. We revealed the forcing and response structures that can be amplified from harmonic oscillations or external actuation over finite wings. In the present triglobal analysis, we have identified the wingspan locations where forcing structures can be amplified near the wing and the regions where the unsteady response develops. We have further characterized the region of dominance of modal structures over the wingspan, with forcing–response mode pairs appearing near the wing root or tip as a function of their characteristic frequency.

Through resolvent wavemakers, we studied the steady to unsteady flow characteristics over swept wings. We showed the regions where self-sustained unsteadiness appears over swept wings and related those to the vortex shedding structures observed in the DNS. We also revealed the most sensitive regions for perturbation amplification in steady wakes over highly swept wings. The forcing–response mode pairs further revealed the mechanisms of spanwise advection of flow structures, which is further related to the spanwise convective speed found for two-dimensional and 3-D resolvent analysis over swept wings and also associated with the nonlinear flow characteristics observed in DNS.

Lastly, we showed for low-aspect-ratio wings that localized perturbations with root- or tip-dominant characteristics are limited as modes evolve globally over the entire wingspan. In fact, we have shown that root- and tip-dominant structures can appear over the wing in a single mode for low-aspect-ratio wings. This behaviour explains the characteristics of the laminar flows around these wings, as observed in DNS, to be different from those of the flows over $sAR = 2$ and $4$ wings. These findings provide fundamental insights into future studies on flow separation over swept wings at higher Reynolds numbers, in which a wider spectrum of fluctuations is present.