2.1. Direct numerical simulations

To study the flows over swept NACA 0015 airfoils, we perform DNS with CharLES, a finite-volume-based compressible flow solver with second- and third-order accuracies in space and time, respectively (Khalighi et al. Reference Khalighi, Ham, Nichols, Lele and Moin2011; Brès et al. Reference Brès, Ham, Nichols and Lele2017). We position the leading edge of the airfoil at $(x^\prime /L_c,\ y/L_c) = (0,0)$. The C-shaped computational mesh extends over $(x^\prime /L_c, y/L_c, z^\prime /L_c) \in [-20,25] \times [-20,20] \times [0,4]$. We build the grids with $(\min {\rm \Delta} x, \min {\rm \Delta} y, \min {\rm \Delta} z)/L_c = (0.005, 0.005, 0.0625)$, with mesh refinement near the airfoil and wake. This yields a mesh with $100\,000$ to $200\,000$ cells on the root plane, extruded in the spanwise direction with equally spaced cell elements. We have verified our computational set-up and validated the results with Zhang et al. (Reference Zhang, Hayostek, Amitay, Burstev, Theofilis and Taira2020a). Our simulations obtained close agreement for the instantaneous and time-averaged velocity components, skin friction lines and pressure coefficients over the wing surface.

We prescribe Dirichlet boundary conditions at the inlet and far-field boundaries as $(\rho, u_{x'}, u_y, u_{z'}, p) = (\rho _\infty, U_\infty \cos \varLambda, 0, U_\infty \sin \varLambda, p_\infty )$, where $\rho$ is density, $p$ is pressure, $u_{x^\prime }$, $u_y$ and $u_{z^\prime }$ are velocity components in the ($x^\prime,y,z^\prime$) directions, respectively, $\rho _\infty$ is the free stream density and $p_\infty$ is the free stream pressure. On the airfoil surface, we prescribe the adiabatic no-slip boundary condition. For the outflow, a sponge layer (Freund Reference Freund1997) is applied over $x^\prime /L_c \in [15,25]$ with the target state being the running-averaged flow over five acoustic time units. For time integration, a fixed acoustic Courant–Friedrichs–Lewy number of 1 is utilized. We start the simulations with uniform flow. The initial transients are flushed out of the computational domain for $80$ convective time units, after which statistics are recorded over at least $100$ convective time units.

The current results are carefully validated by examining the lift, drag and pressure coefficients, defined as

(2.1a–c)\begin{equation} C_L = \frac{F_y}{\frac{1}{2}\rho U_\infty^2 L_c}, \quad C_D = \frac{F_x}{\frac{1}{2}\rho U_\infty^2 L_c} \quad \text{and} \quad C_p = \frac{p - p_\infty}{\frac{1}{2}\rho U_\infty^2 }, \end{equation}

respectively; where $F_x$ and $F_y$ are the $x$ and $y$ force components, respectively. The forces for $\alpha = 20^\circ$ are in close agreement with those reported by Zhang et al. (Reference Zhang, Hayostek, Amitay, Burstev, Theofilis and Taira2020a), as shown in table 1. The flow fields were also compared and exhibit agreement validating the current set-up.

Table 1. Time-averaged lift and drag coefficients ($\overline {C_L}$ and $\overline {C_D}$) compared with Zhang et al. (Reference Zhang, Hayostek, Amitay, Burstev, Theofilis and Taira2020a) for laminar separated flow over NACA 0015 airfoils with $\alpha = 20^\circ$, $\varLambda = 0^\circ$, $15^\circ$, $30^\circ$ and $45^\circ$.

2.2. Resolvent analysis

To analyse the perturbation dynamics over swept wings, let us decompose the flow state $\boldsymbol {q}$ as

(2.2)\begin{equation} \boldsymbol{q}= \bar{\boldsymbol{q}} + \boldsymbol{q}' ,\end{equation}

where $\bar {\boldsymbol {q}}$ is the time- and spanwise-averaged flow and $\boldsymbol {q}'$ is the fluctuating component. With this Reynolds decomposition, we can express the compressible Navier–Stokes equations as

(2.3)\begin{equation} \frac{\partial \boldsymbol{q}'}{\partial t} = \mathcal{L}_{\bar{\boldsymbol{q}}}\boldsymbol{q}' + \boldsymbol{f}' , \end{equation}

where $\boldsymbol {f}'$ gathers all nonlinear terms. The spanwise periodicity and statistical stationarity allows for $\boldsymbol {q}'$ and $\boldsymbol {f}'$ to be represented with temporal and spanwise Fourier modes

(2.4)\begin{equation} [ \boldsymbol{q}'(x',y,z',t), \boldsymbol{f}'(x',y,z',t) ] = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} [ \hat{\boldsymbol{q}}_{k_{z^\prime},\omega}(x',y), \hat{\boldsymbol{f}}_{k_{z^\prime},\omega}(x',y) ] \,{\rm e}^{{\rm i} (k_{z^\prime} z' - \omega t )} {\rm d} k_{z^\prime} \,{\rm d} \omega , \end{equation}

where $k_{z^\prime }$ is the spanwise wavenumber, and $\hat {\boldsymbol {q}}_{k_{z^\prime },\omega }$ and $\hat {\boldsymbol {f}}_{k_{z^\prime },\omega }$ are the biglobal modes for spanwise wavenumber $k_{z^\prime }$ and temporal frequency $\omega$. By substituting these modal expressions into (2.3), we obtain

(2.5)\begin{equation} -{\rm i} \omega \hat{\boldsymbol{q}}_{k_{z^\prime},\omega} = \mathcal{L}_{\bar{\boldsymbol{q}}}\hat{\boldsymbol{q}}_{k_{z^\prime}, \omega} +\hat{\boldsymbol{f}}_{k_{z^\prime},\omega} , \end{equation}

which is an inhomogeneous linear differential equation describing the perturbation evolution with input $\hat {\boldsymbol {f}}_{k_{z^\prime },\omega }$. This formulation leads us to

(2.6)\begin{equation} \hat{\boldsymbol{q}}_{k_{z^\prime},\omega} = \mathcal{H}_{\bar{\boldsymbol{q}}(k_{z^\prime},\omega)} \hat{\boldsymbol{f}}_{k_{z^\prime},\omega},\end{equation}

where

(2.7)\begin{equation} \mathcal{H}_{\bar{\boldsymbol{q}}(k_{z^\prime},\omega)} = [ -{\rm i} \omega I -\mathcal{L}_{\bar{\boldsymbol{q}}} (k_{z^\prime}) ]^{{-}1}\end{equation}

is known as the resolvent (Reddy & Henningson Reference Reddy and Henningson1993; Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). This operator acts as a state-space transfer function between the input forcing $\hat {\boldsymbol {f}}_{k_{z^\prime },\omega }$ and the output response $\hat {\boldsymbol {q}}_{k_{z^\prime },\omega }$ with respect to the base flow $\bar {\boldsymbol {q}}$ (Jovanović & Bamieh Reference Jovanović and Bamieh2005). To construct a discrete linear operator $\boldsymbol{\mathsf{L}}_{\bar {\boldsymbol {q}}}$ and a discrete resolvent $\boldsymbol{\mathsf{H}}_{\bar {\boldsymbol {q}}} \in \mathbb {C}^{m \times m}$, we use the time-averaged flow $\bar {\boldsymbol {q}}$ (McKeon & Sharma Reference McKeon and Sharma2010). Here, $m$ is the resolvent operator size defined by the product of the number of spatial grid points used to discretize the domain and the number of state variables. Appropriate boundary conditions are embedded in $\boldsymbol {H}_{\bar {\boldsymbol {q}}}$.

The discrete resolvent is examined through singular value decomposition,

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

where $\boldsymbol{\mathsf{F}} =[\hat {\boldsymbol {f}}_1,\hat {\boldsymbol {f}}_2, \ldots,\hat {\boldsymbol {f}}_m]$ is the orthonormal matrix of right singular vectors representing the forcing modes, $\boldsymbol {\varSigma } = {\rm diag}[\sigma _1,\sigma _2,\ldots,\sigma _m]$ is the diagonal matrix of singular values ranked in descending order, and $\boldsymbol{\mathsf{Q}} = [\hat {\boldsymbol {q}}_1,\hat {\boldsymbol {q}}_2,\ldots,\hat {\boldsymbol {q}}_m]$ is the orthonormal matrix of left singular vectors representing the response modes (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Jovanović & Bamieh Reference Jovanović and Bamieh2005; McKeon & Sharma Reference McKeon and Sharma2010).

We can also incorporate temporal damping into forcing and response modes as $[\hat {\boldsymbol {q}}_{k_{z^\prime },\omega }, \hat {\boldsymbol {f}}_{k_{z^\prime }, \omega }]\, {\rm e}^{-\beta t}$ through a discounted resolvent operator, where $\beta$ is a time-discounting parameter (Jovanović Reference Jovanović2004; Schmid & Brandt Reference Schmid and Brandt2014). Moreover, the pseudospectral analysis is dependent on the norm (Trefethen & Embree Reference Trefethen and Embree2005). In this work, we use the Chu norm (Chu Reference Chu1965) which is incorporated into the resolvent 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.

The weighted resolvent is dependent on the temporal frequency $\omega$, spanwise wavenumber $k_{z^\prime }$ and the base flow $\bar {\boldsymbol {q}} = \bar {\boldsymbol {q}}(\alpha,\varLambda )$. These parameters define a large parameter space to characterize the effect of sweep angle on the wake dynamics. To facilitate this characterization, we employ an adjoint-based parametric sensitivity method for $\omega$ and $k_{z^\prime }$ (Schmid & Brandt Reference Schmid and Brandt2014; Fosas de Pando & Schmid Reference Fosas de Pando and Schmid2017). This approach is helpful when the parameter space is large, as well as capturing the sensitivity of the resolvent norm to specific geometrical and flow parameters (Skene & Schmid Reference Skene and Schmid2019).

To perform the resolvent analysis, we construct a discrete linear operator $\boldsymbol{\mathsf{L}}_{\bar {\boldsymbol {q}}}$ (Sun et al. Reference Sun, Taira, Cattafesta and Ukeiley2017). This operator is discretized over a 2-D unstructured grid, as shown in figure 1, with a reduced-size spatial domain $x/L_c \in [-10,15]$ and $y/L_c \in [-10,10]$ to alleviate the computational costs of performing resolvent analysis without affecting the accuracy of the resolvent modes. The structured DNS base flow solution is transferred to the unstructured grid via cubic interpolation. We enforce homogeneous Dirichlet boundary conditions for the fluctuating variables $\rho '$ and $u'$ and homogeneous Neumann boundary conditions for $T'$ along the far field and airfoil boundaries, as well as to all variables at the computational outlet. In addition, we apply sponges far from airfoil in conjunction with the boundary conditions.

In the present work, we construct $\boldsymbol{\mathsf{L}}_{\bar {\boldsymbol {q}}}$ with $m$ ranging between $150\,000$ and $200\,000$. The resolvent modes were computed using the randomized resolvent algorithm (Ribeiro, Yeh & Taira Reference Ribeiro, Yeh and Taira2020), sketching the operator with $10$ random test vectors weighted by the gradients of the baseflow $(\| \boldsymbol {\nabla } \rho \|, \|\boldsymbol {\nabla } u_x\|, \|\boldsymbol {\nabla } u_y\|, \|\boldsymbol {\nabla } u_z\|, \|\boldsymbol {\nabla } p\|)$. The resolvent norm converges to at least seven significant digits. For the spectral analysis of $\boldsymbol{\mathsf{L}}_{\bar {\boldsymbol {q}}}$, eigenmodes were computed using the Krylov–Schur method (Stewart Reference Stewart2002) with $128$ vectors for the Krylov subspace and tolerance residual of $10^{-10}$. The direct and adjoint linear systems were solved using the MUMPS (multifrontal massively parallel sparse direct solver) package. The codes used to compute the resolvent modes and eigenvalues are part of the ‘linear analysis package’ made available by Skene, Ribeiro & Taira (Reference Skene, Ribeiro and Taira2022).

3.1. Wake characterization

The wake structures are affected by sweep and incidence angles, as shown by the isosurfaces of $Q$ coloured by the vorticity $\omega _z$ in figure 2. For angles of attack $\alpha \le 20^\circ$, the flow is 2-D and the wake vortices are aligned with the sweep angle. For $\alpha \ge 26^\circ$, the vortical structures exhibit spanwise oscillations with a transition from 2-D to 3-D vortex shedding. At $\alpha = 26^\circ$, a sinusoidal pattern of oscillations appears in the spanwise direction and, at $\alpha = 30^\circ$, streamwise vortical structures emerge.

Figure 2. Instantaneous flow field visualization with isosurfaces of $Q$ criterion coloured by the vorticity component $\omega _z$. For $\alpha \ge 26^\circ$, the wake is 3-D. As $\varLambda$ increases, the vortices become slanted with the sweep angle and wake three-dimensionality is reduced for $\alpha \ge 26^\circ$.

Highly swept wings induce spanwise flows and attenuate wake three-dimensionality as evident from the flow visualizations. When the sweep angle is $\varLambda \le 15^\circ$ the wake is similar to the flow over unswept wings for all angles of attack. The wake is significantly altered for sweep angles $\varLambda \ge 30^\circ$, especially at high angles of attack, when spanwise oscillations are advected by the spanwise flow. For instance, at $\varLambda = 45^\circ$ and $\alpha = 26^\circ$, the spanwise oscillations are almost suppressed. The attenuation of spanwise fluctuations also occurs for $\alpha = 30^\circ$, as we observe a similar effect for $\varLambda \ge 30^\circ$.

Even though sweep attenuates three-dimensionality, the sustained unsteadiness of the wake suggests the existence of self-supported mechanisms that yield distinct vortex shedding patterns in swept wings at high incidence angles. We gain further insights into the separated flows by analysing the time-averaged flow field contours of streamwise and spanwise velocities, as seen in figure 3. In general, as we increase the sweep angle, the $\bar {u}_x = 0$ contour line approaches the wing surface. For $\alpha = 30^\circ$, we also notice a circular $\bar {u}_z$ profile appearing in the wake region where the spanwise flow is stronger.

Figure 3. Time-averaged flow field visualization with $z$-velocity component, $\bar {u}_z \in [0,0.5]$, in greyscale, for $0^\circ \le \varLambda \le 45^\circ$ and $\alpha = 16^\circ$ and $30^\circ$. Red solid contours mark the laminar separation bubble, with six equally distributed isolines of $x$-velocity component, $\bar {u}_x \in [0,0.5]$. Cross-flow component $\bar {u}_z$ strengthens with $\alpha$ and $\varLambda$.

The aerodynamic loads exerted on the wing are also affected by the sweep angle (Zhang et al. Reference Zhang, Hayostek, Amitay, Burstev, Theofilis and Taira2020a; Zhang & Taira Reference Zhang and Taira2022). However, it is possible to observe similarities between force characteristics with different sweep angles through the independence principle. This leads to the application of proper scaling factors to collapse aerodynamic properties for a variety of swept wings where adverse pressure gradients and spanwise fluctuations are negligible (Wygnanski et al. Reference Wygnanski, Tewes, Kurz, Taubert and Chen2011).

For the present flows, however, the adverse pressure gradients cannot be neglected due to the massive separation. In figure 4, we show that $\overline {C_L}$ and $\overline {C_D}$ differ for the same $\alpha$ if we analyse the flow variables in $(x,y,z)$. The coefficients collapse if we consider scaling the same force coefficients in ($x^\prime,y,z^\prime$). Here, the vector-valued variables in ($x,y,z$) aligned with $x$ are scaled with $\cos \varLambda$, and the effective chord length

(3.1)\begin{equation} L_{c}^\prime = L_c ( \cos^2 \alpha \cos^2 \varLambda + \sin^2 \alpha )^{1/2} \le L_c\end{equation}

is used to form the scaled $C_L^\prime$ and $C_D^\prime$ coefficients as

(3.2a–c)\begin{equation} C_L^\prime = \frac{F_y}{\frac{1}{2}\rho (U_\infty \cos \varLambda)^2 L_c^\prime}, \quad C_D^\prime = \frac{F_x \cos\varLambda}{\frac{1}{2}\rho (U_\infty \cos \varLambda)^2 L_c^\prime},\end{equation}

where $F_x \cos \varLambda = F_{x^\prime }$ is the $x^\prime$ component of the pressure and viscous forces integrated over the airfoil surface per unit depth. As shown in figure 4(d–f), these scaled coefficients collapse over the angles of attack. Deviations are noticed only for sweep angles $\varLambda \ge 30^\circ$ at high angles of attack $\alpha \ge 26^\circ$.

Figure 4. Time-averaged lift, $\overline {C_L}$, drag, $\overline {C_D}$ and coefficient ratio $\overline {C_L / C_D}$, for all $\alpha, \varLambda$ pair of the present study compared with 2-D incompressible results shown by Zhang et al. (Reference Zhang, Hayostek, Amitay, He, Theofilis and Taira2020b). Panels (d–f) show the scaled time-averaged coefficients where the flow is analysed in ($x',y,z'$), and the results collapse for each $\alpha$, for all sweep angles.

The shown scaling can also be applied to the Fage–Johansen Strouhal number in the $(x^\prime,y)$ plane as

(3.3)\begin{equation} St^\prime = \frac{\omega}{2{\rm \pi}} \frac{L_c^\prime \sin \alpha_{eq}}{U_\infty \cos \varLambda}. \end{equation}

As presented in figure 5, the PSD profiles of $C_L^\prime$ for swept wings exhibit peaks at their characteristic Strouhal number of the vortex shedding and its harmonics. When the adapted Fage–Johansen Strouhal number $St^\prime$ is considered, the spectral peaks collapse at the same frequencies as observed for the unswept wings. For $\alpha = 16^\circ$, the flow is characterized by a single 2-D vortex shedding and the $C_L^\prime$ spectra is smooth with distinct peak values. For $\alpha = 30^\circ$, the spectra exhibits secondary peaks for $\varLambda \le 30^\circ$ and is smooth for $\varLambda = 45^\circ$, when three-dimensionality is attenuated.

Figure 5. Scaled lift power spectrum density (PSD) for (a,c) $\alpha = 16^\circ$ and (b,d) $\alpha = 30^\circ$ at sweep angles $\varLambda = 0^\circ$, $30^\circ$ and $45^\circ$. In (c,d) the Fage–Johansen Strouhal number is analysed in the $(x',y)$ plane and the dominant and harmonic frequencies for swept wings collapse with the unswept wings.

We can also similarly normalize the pressure coefficients $C_p$ by considering the $U_\infty \cos \varLambda$ in place of $U_\infty$ in (2.1a–c). Indeed, large differences in $C_p$ distribution over the wing are shown in figure 6(a,b), however, if we consider the scaled-$C_p$, we reveal that the pressure distributions collapse for moderate angles of attack, even though these flows exhibit massive separation, as shown in figure 6(c,d). Although we can bring pressure coefficients closer using $C_p / \cos ^2 \varLambda$, we notice some deviations for higher angles of sweep and attack, as observed for $\varLambda = 45^\circ$ and $\alpha = 30^\circ$. This motivates us to further understand how massively separated streamwise flow and the strong spanwise flow may impose additional loads over the wing. These deviations indeed suggest that even when the three-dimensionality is attenuated, the wake over laminar swept wings can exert additional aerodynamic forces over the wing for $\varLambda \ge 30^\circ$ and $\alpha \ge 26^\circ$.

Figure 6. Time-averaged pressure coefficients (a,b) $C_p$ and (c,d) $C_p/\cos ^2 \varLambda$ over the airfoil surface at (a,c) $\alpha = 16^\circ$ and (b,d) $30^\circ$.

3.2. Force element analysis

To further understand the sources of deviations in the independence principle, we use the force element theory of Chang (Reference Chang1992) to relate the near-body vortical structures to aerodynamic forces. Through this method, we identify lift and drag force elements in the near-wake region of the current low-Reynolds number flows and analyse the distribution of force elements near the surface as the wing is swept. This analysis captures the emerging wake structures over swept wings at high angles of attack that exert nonlinear poststall forces onto the wing.

We observe that the emergence of these force elements is associated with a departure from the collapsed force and pressure coefficients. The sweep-angle dependent scaling in (3.2a–c) and (3.3) assumes independence of streamwise and spanwise flow components. The interaction between them can cause a departure from the collapsed force and pressure coefficients. By using the force element analysis, we uncover near-wake structures that are responsible for the extra forces at high angles of sweep and attack.

To perform this analysis, we start by defining an auxiliary potential with a specific boundary condition of $-\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\nabla } \phi _i = \boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {e}_i$ along the wing surface, where $\phi$ is the auxiliary potential, $\boldsymbol {n}$ is the unit wall-normal vector and $\boldsymbol {e}_i$ ar normal vectors in the $i$th direction. For a solenoidal velocity field, the force exerted on the wing in the $i$th direction can be written as

(3.4)\begin{equation} F_i ={-}\int_V \boldsymbol{u} \times \boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi_i \,\text{d}V + \frac{1}{Re} \int_S \boldsymbol{n} \times \boldsymbol{\omega} \boldsymbol{\cdot} (\boldsymbol{\nabla} \phi_i +\boldsymbol{e_i}) \, \text{d}S, \end{equation}

where the first integral term on the right-hand side is named the vortical elements force and the second integral term is called the surface element force. We can visualize the lift and drag force elements by the Hadamard product of the $\boldsymbol {\nabla } \phi _i$ and the Lamb vector ($\boldsymbol {\omega } \times \boldsymbol {u}$). The auxiliary potential field decays rapidly far from the wing surface, hence the force elements of $(\boldsymbol {u} \times \boldsymbol {\omega }) \circ \boldsymbol {\nabla } \phi _i$ are concentrated near the wing.

Sweep has a strong influence in limiting the validity of the independence principle at high angles of incidence and it favours the emergence of additional force elements near the wing surface. To show this, let us reveal the vortical structures that generate lift $(\boldsymbol {\omega } \times \boldsymbol {u}) \boldsymbol {\cdot } \boldsymbol {\nabla } \phi _y$ at the instance when the maximum lift is achieved, as shown in figure 7. Drag elements have a similar behaviour as the lift and are not shown for brevity. For $\alpha \le 20^\circ$, the force elements are located near the airfoil's leading and trailing edges, in the shear-dominated region of the flow, along the edge of the laminar separation bubble.

Figure 7. Characterization of force elements on swept wings. The symbols refer to $\boldsymbol {\circ }$ (yellow), force elements only near the leading and trailing edges; $\boldsymbol {\square }$ (blue), additional equally spaced small lift elements; and $\boldsymbol {\triangle }$ (red), large force structures observed on the suction side. On the right-hand side, we show isosurfaces of lift force elements, $(\boldsymbol {\omega } \times \boldsymbol {u}) \boldsymbol {\cdot } \boldsymbol {\nabla } \phi _L \in [-0.3,0.3]$ and isosurfaces of $Q$ values coloured by streamwise velocity component $\bar {u}_x$ with range $[0,1]$ for the time step with the highest lift coefficient $C_L^\prime$ for $\alpha = 30^\circ$ and $15^\circ \le \varLambda \le 45^\circ$.

Additional lift elements appear for higher angles of sweep and attack, as shown in figure 7. These lift elements are observed over the final quarter chord of the airfoil on the suction side, and as the sweep angle increases, they also increase in size. The force elements can be contrasted with the vortical structures in $Q$ criterion visualization in figure 7. As such coherent structures are present inside the laminar separation bubble, with size and shape similar to the force elements, they can be identified as the lift elements related to the larger deviations in figure 6(c,d).

We observe that sweep affects the coherent structures, time-averaged flow fields and aerodynamic loads and, although some similarities are perceived, sweep has a strong influence on the wake flow downstream at the higher angles of attack. This suggests that flow perturbations originated near the airfoil in the laminar flows over swept wings can be related to the features observed in the nonlinear simulations. To further understand how sweep alters the vortex dynamics we conduct global resolvent analysis to identify the sources of self-sustained mechanisms near the wing that affect the wake behaviour.

3.3. Resolvent analysis

To identify the existence of regions susceptible to the growth of perturbations in the flows over swept wings, we study the effect of sweep using resolvent analysis. Details on the stability of the linear operators and the usage of time-discounting are provided in the Appendix. We characterize the present flows through the dominant singular value $\sigma _1$ of the resolvent operator and its corresponding singular vectors $\hat {q}_{k_{z^\prime },\omega }$ and $\hat {f}_{k_{z^\prime },\omega }$. The influence of sweep on the vortex dynamics in the wake is analysed through the forcing and response modes in figure 8. The shown modal structures highlight the regions of the flow field which are more sensitive and responsive to the growth of perturbations.

Figure 8. (a) Forcing (red boxes) and (b) response (blue boxes) contours for the $|u_x| / \|u_x\|_\infty \in [0.1,1]$ in blue–green–red scale at the largest $\sigma _1$ with $k_{z^\prime } = 0$ at $\varLambda = 0^\circ$ and $45^\circ$. Forcing modes extend in the wake and response modes become closer to the airfoil in swept wings. Isosurfaces of dominant resolvent gain $\sigma _1$ in the $\varLambda$–$St^\prime$–$k_{z^\prime }$ space for $\alpha = 16^\circ$ to $30^\circ$.

Forcing modes are more localized than response modes, which are supported in the shear-dominated region of the flow, where perturbations can be highly amplified. On the other hand, the response modes are seen in the wake. For swept wings, the response modes are deformed spatially towards the airfoil surface and, for $\alpha = 30^\circ$, we notice the emergence of a characteristic responsive region near the airfoil leading edge. Such response structures appear over swept wings only at high angles of attack, as seen in figure 8(a). We observe the contours of the magnitude of modal streamwise velocity component $| \hat {u}_x |$ to reveal the regions of the flow that have more responsive to introduced perturbations. We visualize similar contours for the forcing counterpart and see that the modes extend slightly into the wake, over the laminar separation bubble. This behaviour reveals that the flow over swept wings can amplify optimal disturbances closer to the airfoil suction side, which can be used to alter the formation of the laminar separation bubble.

Furthermore, the present resolvent analysis predicts the formation of oblique vortex shedding, as observed in Mittal et al. (Reference Mittal, Pandi and Hore2021) and Zhang et al. (Reference Zhang, Hayostek, Amitay, Burstev, Theofilis and Taira2020a), even though the present study is performed on spanwise periodic wings. Previous works have shown that oblique coherent structures become spatially periodic for large aspect ratio wings, making spanwise periodic analysis valid to study such 3-D structures. Through resolvent analysis, we can explain how oblique coherent structures are advected by the flow stream using the spatiotemporal frequencies of the optimal response modes.

The frequency at maximum $\sigma _1$ for each spanwise wavenumber is a function of the sweep angle and is characterized by the convection speed of the optimal oblique modes. We compute this phase speed as $c = \text {d} \omega / \text {d} k_{z^\prime }$, the slope of the optimal response frequencies for each spanwise wavenumber. This value is a function of $\varLambda$ and $k_{z^\prime }$ (Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019; He & Timme Reference He and Timme2021; Plante et al. Reference Plante, Dandois, Beneddine, Laurendeau and Sipp2021). For flow regimes studied in the present work, $f \approx 0.25 U_\infty \tan \varLambda$ gives a reasonable prediction for the frequency of the maximum $\sigma _1$ for each spanwise wavenumber and sweep angle for all angles of incidence, as shown in figure 9. This function can be used to predict the optimal forcing and response modes for laminar separated flows over swept wings. Additionally, the mode shapes of the optimal forcing and response are similar for low $k_{z^\prime }$. The present results reveal the optimal actuation location and response as well as the spatiotemporal behaviour of the flow perturbations over laminar separated flows on swept wings.

Figure 9. The leading resolvent gain contours at $\alpha = 20^\circ$ and $30^\circ$ for $0^\circ \le \varLambda \le 45^\circ$. The dash–dotted slopes represent spanwise convection speeds. The green line exhibits the convection speed prediction with $0.25 U_\infty \tan \varLambda$.

In general, oblique modes are the most amplified optimal disturbances for all swept wings. The effect of sweep on $\sigma _1$, however, depends on the angle of attack, as shown in figure 8 and summarized in table 2. For $\alpha \le 20^\circ$, swept wings have higher amplification than unswept wings. This is a distinct behaviour when compared with $\alpha \ge 26^\circ$, in which swept wings have lower resolvent norm than unswept wings. This behaviour suggests that it is more challenging to perturb highly swept wings at high angles of attack.

Table 2. The maximum leading resolvent gain ${max}(\sigma _1)$ for each $\alpha,\varLambda$ pair. In the inset figure, we plot ${max}(\sigma _1)$ in $St^\prime$–$k_{z^\prime }$ space coloured in blue scale with respect to the minimum and maximum $\sigma _1$ for each $\alpha$.

The spatial and temporal frequency of the maximum resolvent gain ${max}(\sigma _1)$ in the $St^\prime$–$k_{z^\prime }$ space depends on the sweep angle, as shown in table 2. For unswept wings, the largest resolvent gain $\sigma _1$ is found for the 2-D setting of $k_{z^\prime } = 0$ associated with the temporal frequency of the characteristic vortex shedding. However, both $k_{z^\prime }$ and $St^\prime$ of the optimal disturbances ${max}(\sigma _1)$ increase with the sweep angle. Thus, the 3-D oblique modal structures are not only predicted by the present resolvent analysis but also found to be the most amplified flow mechanism in swept wings, as shown in figure 10.

Figure 10. (a,c) Forcing ($\,\hat {f}$, red boxes) and (b,d) response ($\hat {q}$, blue boxes) modes isosurfaces with $y$-velocity components $\hat {u}_y / \| \hat {u}_y \|_\infty \in \pm 0.2$ in red–blue colour scale for $\alpha = 20^\circ$ and $\varLambda = 45^\circ$. Wingspan length is $10 L_c$. Forcing and response modes are associated with the largest resolvent gain for each $k_{z^\prime }$ as shown in the $\sigma _1$ contours over $St^\prime$–$k_{z^\prime }$ plane.

Although all flows analysed in this work are spanwise periodic, and oblique shedding is absent in the DNS, the large $\sigma _1$ in $k_{z^\prime } > 0$ modes suggest that the spanwise flow over swept wings supports the formation and shedding of 3-D oblique vortices in infinite wings. To analyse the spatial behaviour of oblique modes, let us focus on the resolvent analysis at the angle of attack $\alpha = 20^\circ$ and sweep angle $\varLambda = 45^\circ$, as seen in figure 10. As noticed at $k_{z^\prime } = 0$, the 2-D forcing and response modes are aligned with the wingspan, however, the maximum resolvent gain $\sigma _1$ in the $St^\prime$–$k_{z^\prime }$ for this flow is found at $k_{z^\prime } = 0.3{\rm \pi}$ and $St = 0.19$, where modes are oblique with respect to the wingspan. For this reason, the flow over swept wings has a higher propensity to develop oblique shedding when compared with the flow over unswept wings and such characteristic is revealed through resolvent analysis to be associated with the sweep angle.

The flow mechanisms that are responsible for oblique shedding and the attenuation of the spanwise oscillations for swept wings were described as the growth of response modes towards the airfoil surface and the extension of forcing modes into the wake and over the laminar separation bubble. This phenomenon also creates an overlapping region of the flow where both forcing and response modes are supported. This overlap of forcing and response structures is more prominent at the higher sweep angles, although it is also present in unswept wings. The overlap of forcing and response modes and their associated resolvent gain can both be relevant to characterize how the flow over swept wings at high incidence gives rise to perturbations on the flow as we discuss such a phenomenon through the lens of wavemakers.

3.4. Wavemakers

Wavemakers have been described as regions of the flow field characterized by both high sensitivity and responsiveness to perturbation growth (Giannetti & Luchini Reference Giannetti and Luchini2007; Giannetti et al. Reference Giannetti, Camarri and Luchini2010; Fosas de Pando et al. Reference Fosas de Pando, Schmid and Sipp2017). Such regions are optimal for the introduction of self-sustained instabilities, acting as the source of global instabilities of the system, and motivate the analysis of structural sensitivity of the modal forcing and response structures.

In global stability analysis, wavemakers are generally derived with direct and adjoint modes. Here, we quantify the strength of the wavemakers with the inner product of the forcing and response modes $\langle \hat {q}_{k_{z^\prime },\omega }, \hat {f}_{k_{z^\prime },\omega } \rangle$ and visualize the corresponding wavemaker modes with their Hadamard product. Wavemaker modes exhibit a higher magnitude downstream the airfoil, as shown in figure 11, in the region where 3-D flow develops for $\alpha = 30^\circ$. Thus, the emergence of strong wavemakers near the leading edge and over the separation bubble highlights the presence of self-sustained oscillations in the flow field. As the sweep angle tends to empower forcing and response modes overlap, wavemakers tend to be spatially wider in swept wings.

Figure 11. Contours of $\xi = (\sigma _1 / {max}(\sigma _1)) \langle \hat {q}_i, \hat {f}_i \rangle$, the inner product between forcing and response modes, scaled by the ratio of $\sigma _1$ and the maximum $\sigma _1$ for each $\alpha$. Green line shows the convection speed $c = \text {d} \omega / \text {d} k_{z^\prime }$ for the optimal wavemakers. Spatial modes shown by the Hadamard product of $\hat {q}_i$ and $\hat {f}_i$, in magnitude, normalized by their maximum value, and coloured in purple scale.

We evaluate the strength of wavemakers in the $\varLambda$–$St^\prime$–$k_{z^\prime }$ space with the inner product of pseudomodes $\langle \hat {q}_{k_{z^\prime },\omega }, \hat {f}_{k_{z^\prime },\omega } \rangle$ and their associated resolvent gain $\sigma _1$. In this way, we avoid accounting for the wavemakers where $\sigma _1$ is too small to amplify perturbations. Hence, to understand which combination of temporal frequencies, spanwise wavenumber and sweep angle is most likely to generate wavemakers, we consider the coefficient $\xi = (\sigma / {max} (\sigma _1)) \langle \hat {q}_i, \hat {f}_i \rangle$, where ${max} (\sigma _1)$ is evaluated for each angle of attack over the $St^\prime$–$k_{z^\prime }$–$\varLambda$ space, with contours shown in figure 11 for $\alpha = 20^\circ$ and $30^\circ$.

The resolvent modes with higher values of $\xi$ are observed for swept wings at higher angles of attack. The wavemaker modes appear where vortex shedding develops in the wake. Hence a higher $\xi$ suggests perturbations are introduced with higher gain to be amplified in this region for flows over wings at high incidence. Furthermore, those disturbances feed the flow with self-generated disturbances that maintain three-dimensionality of the wake, as observed for instance at $\alpha = 30^\circ$. In unswept wings, the $St^\prime$–$k_{z^\prime }$ frequencies with strong wavemaker modes is found for 2-D wavemakers.

For swept wings, the modes with the highest $\xi$ coefficient for each angles of attack are located at $\varLambda = 45^\circ$ and non-zero $k_{z^\prime }$, hence being associated with oblique modes. This finding is in agreement with the previous observations on the overlap of forcing and response modes, in figure 8. Hence, even if the amplification gain $\sigma _1$ is reduced for swept wings at $\alpha = 30^\circ$, the overlap of forcing and response modes is stronger, which introduces wavemakers over swept wings that are stronger than wavemakers for unswept wings, which explains the three-dimensionality observed in these flow fields.

This finding, however, is in contrast with the DNS results that show an attenuation of spanwise oscillations with the sweep angle. To understand why such alleviation occurs, we must observe that optimal responses and optimal wavemakers also have an associated wave speed, characterized by their spatial and temporal frequencies, that is associated with the transport of disturbances over the periodic direction.

We note that an optimal wavemaker speed being faster than the optimal response wave speed leads to an attenuation on three-dimensionality. Wavemakers yield self-sustained instabilities in swept wings with an associated wavemaker phase speed $c_w = \text {d}\omega / \text {d} k_{z^\prime }$, which we characterize by the slope of the slash–dotted green lines in figure 11. As observed in table 3, when $c_w$ is large for high sweep angles and the optimal response $c$ is small, the reduction of spanwise oscillations is seen in the flow field. In such cases, a misalignment appears between optimal responses and wavemakers which cannot support spanwise oscillations. For this reason, wavemakers cannot sustain 3-D disturbances over swept wings.

Table 3. Convective speed $c = \text {d} \omega / \text {d} k_{z^\prime }$ for the optimal response and the optimal wavemakers for the 3-D flows at $\alpha = 26^\circ$ and $30^\circ$ and sweep angles $0^\circ \le \varLambda \le 45^\circ$.

Finally, even for $\alpha = 20^\circ$, in which the flow field is 2-D, resolvent analysis reveals the presence of wavemakers. Those are associated with the sustained unsteady 2-D vortex shedding. To sustain three-dimensionality, wavemakers must introduce sufficiently strong 3-D structures to the flow with high amplification gain. For swept wings at high incidence, even though optimal wavemakers have a high gain, they are advected faster than the optimal responses, which reduces the flow three-dimensionality.