Linear modal instabilities around post-stall swept finite-aspect ratio wings at low Reynolds numbers
Anton Burtsev, Wei He, Shelby Hayostek, Kai Zhang, Vassilios Theofilis, Kunihiko Taira, Michael Amitay
TThis draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Linear modal instabilities around post-stallswept finite-aspect ratio wings at lowReynolds numbers
Anton Burtsev , Wei He † , Shelby Hayostek , Kai Zhang , VassiliosTheofilis , Kunihiko Taira , and Michael Amitay Department of Mechanical, Materials and Aerospace Engineering, University of Liverpool,Brownlow Hill, England L69 3GH, United Kingdom Department of Mechanical, Aeronautical, and Nuclear Engineering, Rensselaer PolytechnicInstitute, Troy, NY 12180, USA Department of Mechanical and Aerospace Engineering, Rutgers University, Piscataway, NJ08854, USA Department of Mechanical and Aerospace Engineering, University of California, Los Angeles,CA 90095, USA Escola Politecnica, Universidade S˜ao Paulo, Avda. Prof. Mello Moraes 2231, CEP 5508-900,S˜ao Paulo-SP, Brasil(Received xx; revised xx; accepted xx)
Linear modal instabilities of flow over finite-span untapered wings have been investigatednumerically at Reynolds number 400, at a range of angles of attack and sweep ontwo wings having aspect ratios 4 and 8. Base flows have been generated by directnumerical simulation, marching the unsteady incompressible three-dimensional Navier-Stokes equations to a steady state, or using selective frequency damping to obtainstationary linearly unstable flows. Unstable three-dimensional linear global modes ofswept wings have been identified for the first time using spectral-element time-steppingsolvers. The effect of the wing geometry and flow parameters on these modes has beenexamined in detail. An increase of the angle of attack was found to destabilize the flow,while an increase of the sweep angle had the opposite effect. On unswept wings, TriGlobalanalysis revealed that the most unstable global mode peaks in the midspan region ofthe wake; the peak of the mode structure moves towards the tip as sweep is increased.Data-driven analysis was then employed to study the effects of wing geometry and flowconditions on the nonlinear wake. On unswept wings, the dominant mode at low angles ofattack is a Kelvin-Helmholtz-like instability, qualitatively analogous with global modes ofinfinite-span wings under same conditions. At higher angles of attack and moderate sweepangles, the dominant mode is a structure denominated the interaction mode . At highsweep angles, this mode evolves into elongated streamwise vortices on higher aspect ratiowings, while on shorter wings it becomes indistinguishable from tip-vortex instability.
Key words: † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b A. Burtsev et al.
1. Introduction
Our present concern is with linear global instability mechanisms associated with un-steadiness of laminar three-dimensional separated flows over finite aspect ratio, constant-chord unswept and swept wings at low Reynolds numbers. To-date the vast majorityof instability studies has focused on simplified models of laminar separation with nospanwise base flow component, as encountered in flows over two-dimensional profiles, orspanwise homogeneous flow over infinite-span three-dimensional wings, that have beenused as a proxy to understand fundamental mechanisms of separation in practical fixed-or rotary-wing applications. However, either of these approximations fails to addressthe essential three-dimensionality of the flow field (Wygnanski et al. et al. per se on the wing surface, or a deepunderstanding of the complex vortex dynamics arising from this instability on a finite-span wing, as a function of the angles of attack ( α ) and sweep ( Λ ). In fact, there is a void inthe literature that employs three-dimensional global (TriGlobal) linear instability analysisappropriate for the fully inhomogeneous three-dimensional flow field around a finiteaspect ratio wing. The present work aims to close this knowledge gap by documentingmodal instability mechanisms and their evolution on different wing geometries and flowconditions.A short review of existing literature on the subject sets the scene for the work performedherein. Laminar separation studies have extensively employed adverse pressure gradientsto generate nominally two-dimensional separation bubbles on a flat plate geometry.Although such bubbles were known to be structurally unstable (e.g. Dallmann 1988),Theofilis et al. (2000) showed that the physical mechanism leading to unsteadinessand three-dimensionalisation of a nominally two-dimensional laminar separation bubble,as well as to breakdown of the associated two-dimensional vortex, arises from self-excitation of a previously unknown stationary three-dimensional global mode. Soon afterthe latter work, global linear stability theory was applied to understand instabilitymechanisms in two-dimensional airfoils (Theofilis et al. et al. ∼ .
5% of backflow as a necessary condition for self-excitation of the nominally two-dimensional flow. Furthermore, linear superposition of the global mode discovered byTheofilis et al. (2000) upon the two-dimensional laminar separation bubble revealed thewell-known three-dimensional U-separation pattern (Hornung & Perry 1984; Perry &Chong 1987; D´elery 2013), while the surface streamlines topology induced by the globalmode resembled the characteristic cellular structures known as stall cells (SC), that areobserved experimentally to form on stalled wings.Stall cells consisting of pairs of counter rotating vortices on the suction side of a wingarise at angles of attack near maximum lift and have been observed experimentally for along time (Moss & Murdin 1968; Bippes & Turk 1980; Winkelman & Barlow 1980; Weihs& Katz 1983; Bippes & Turk 1984; Schewe 2001; Broeren & Bragg 2001). Application of inear modal instabilities around post-stall swept finite-aspect ratio wings et al. (2012).Steady massively separated spanwise homogeneous flow over stalled wings was revisitedby He et al. (2017 a ) using global linear modal and nonmodal stability tools. Flow overthree different NACA airfoils (symmetric 0009 and 0015 as well as cambered 4415) wasanalysed at 150 (cid:54) Re (cid:54)
300 and 10 ◦ (cid:54) α (cid:54) ◦ . Two different mechanisms wereidentified by linear modal instability analysis: a travelling Kelvin-Helmholtz (K-H) modedominating the flow at large spanwise periodicity length and a three-dimensional station-ary mode most active as the spanwise periodicity length becomes smaller. On all threeairfoils considered, nonmodal analysis showed that linear optimal perturbations evolveinto travelling K-H modes. Furthermore, the travelling K-H mode was the first to becomeunstable as Reynolds number or angle of attack were increased, with the flow becomingtwo-dimensionally unsteady prior to becoming three-dimensionally unstable. Secondaryinstability analysis of the time-periodic base flow ensuing linear amplification of theK-H mode, via temporal Floquet theory, revealed two secondary amplified modes withspanwise wavelengths of approximately 0.6 and 2 chords. These modes are reminiscentof the classic Mode A and B instabilities of the circular cylinder (Barkley & Henderson1996; Williamson 1996) although, unlike the cylinder, the short-wavelength perturbationwas the first to become linearly unstable over all three airfoils considered.The work of He et al. (2017 a ) showed that stall-cell like streamline patterns on thewing surface arise from linear amplification of this short-wavelength secondary instability.By contrast to the primary-instability based scenario proposed by Rodr´ıguez & Theofilis(2011), the mechanism discovered by He et al. (2017 a ) could explain the emergence ofstall cells at lower angles of attack. Zhang & Samtaney (2016) extended the analysisof He et al. (2017 a ) to study instability of unsteady flow over a NACA 0012 spanwiseperiodic wing at higher Reynolds numbers, 400 (cid:54) Re (cid:54) α = 16 ◦ . At Re = 800and 1000 these authors identified two oscillatory unstable modes corresponding to near-wake and far-wake instabilities, alongside a stationary unstable mode, while only oneunstable mode was found at the lower Re = 400 and 600. Ground-proximity effects onthe stability of separated flow over NACA 4415 at low Reynolds numbers were studiedusing BiGlobal theory with consideration of both flat (He et al. c ) and wavy groundsurfaces (He et al. b ). Finally, Rossi et al. (2018) considered incompressible flow overa NACA 0010 airfoil and a narrow ellipse of same thickness at a large angle of attack of30 ◦ (100 (cid:54) Re (cid:54) et al. (2019), who investigated cellular patterns on infiniteswept wings in subsonic stall (NACA 4412, Re = 3 . × , Mach = 0.2) and transonicbuffet (OALT25, Re = 3 × , Mach ≈ A. Burtsev et al. and numerical work of the latter configuration is discussed. Early experimental studieson finite aspect ratio ( AR ) wings are summarised in Boiko et al. (1996). More recentlyaerodynamic performance of small aspect ratio ( AR = 0 . −
2) wings that can be found insmall unmanned aerial vehicles has been studied experimentally (Torres & Mueller 2004)and computationally (Cosyn & Vierendeels 2006). Taira & Colonius (2009) used three-dimensional direct numerical simulation (DNS) to study impulsively translated flat-platewings ( AR = 1 −
4) of different planforms at a wide range of angles of attack and Reynoldsnumbers (300 (cid:54) Re (cid:54) et al. et al. (2017 b ) performed linear global instability analysis using spatial BiGlobaleigenvalue problem and linear PSE-3D disturbance equations in the wake of a low aspectratio three-dimensional wing of elliptic planform constructed using the Eppler E387airfoil at Re = 1750. Two types of disturbances were identified: symmetric perturbationscorresponding to instability of the vortex sheet connecting the trailing vortices andantisymmetric perturbations peaking at the vortex sheet and also in the neighbourhoodof the trailing vortex cores. A comparison with local stability theory has demonstratedthat local analysis has limited ability to accurately capture the physics of the instability.Edstrand et al. (2018 a ) carried out spatial and temporal stability analysis of a wakeand trailing vortex region behind a NACA 0012 finite wing at Re = 1000, α = 5 ◦ and AR = 1 .
25, documenting seven unstable modes with the wake instability dominating inboth temporal and spatial analyses. Unlike many stability analysis works focusing onlyon the vicinity of the tip vortex, the full half-span of the wing was considered. AlthoughBiGlobal stability analysis was employed, streamwise rather than spanwise homogeneitywas used exploiting the absence of large scale separation at the low angle of attack.This allowed capturing three-dimensional modes with structures in the tip and the wakeregions. Subsequent work of Edstrand et al. (2018 b ) on the same geometry employedparabolised stability analysis to guide the design on active flow control for tip vortex. Asubdominant fifth instability mode was found to counter-rotate with the tip vortex. Flowcontrol based on this fifth mode and forcing introduced at the trailing edge rather thanthe wing tip attenuated the tip vortex. Dynamic mode decomposition of the controlledflow was used to assess the effectiveness of such control.Navrose et al. (2019) conducted nonmodal stability analysis of a trailing vortex systemover a flat plate and NACA 0012 wing at a range of angles of attack, aspect ratios andReynolds numbers. Unlike in earlier studies their analysis included the tip vortex and flowover the wing. It was shown that the linear optimal perturbation is located near the wingsurface which advects into the tip vortex region during its evolution, which agrees with thefindings of Edstrand et al. (2018 b ). The displacement of the vortex core due to evolutionof the optimal perturbation was proposed as a possible mechanism behind trailing vortexmeandering. The results of this study point out the need to consider the entire complexthree-dimensional flow over the wing rather than focusing on a specific isolated regionin order to capture the full mechanics of the flow. These studies have demonstratedthat addressing the three-dimensionality of finite wing wake through stability analysisallows for enhanced understanding of the underlaying physical mechanisms. However, inear modal instabilities around post-stall swept finite-aspect ratio wings et al. et al. et al. et al. et al. et al. et al. et al. ram’s horn” vortex (Black 1956).As soon as local stall appears on a swept wing spanwise boundary layer flow alters thestall characteristics of sections of attached flow along the span (Harper & Maki 1964). Yen& Hsu (2007) have experimentally studied the surface flow patterns and wake structuresof a NACA 0012 with a sweep angle, Λ of 15 ◦ ( AR = 5 , Re = 2785) at a wide rangeof angles of attack: from fully attached flow all the way to bluff body like behaviourat α = 45 ◦ and above. Further experimental investigations at higher Reynolds numbersperformed over a range of sweep angles were used to classify the flows into seven boundarylayer flow regimes that were found to be closely related to aerodynamic performance(Yen & Huang 2009). However, the effect of the spanwise velocity component to laminarthree-dimensional flow separation and vortex dynamics, emphatically demonstrated byWygnanski et al. (2011, 2014) to be central to the prediction of turbulent boundary layerproperties on a swept wing, remained unexplored.In the framework of our present combined theoretical/numerical and experimentalefforts, Zhang et al. (2020 a ) employed direct numerical simulation to analyse the de-velopment of three-dimensional separated flow over unswept finite wings at a range ofangles of attack ( Re = 400 , < AR < AR > et al. (2020 b ) addressed swept wing flows at the same conditions.Several stabilisation mechanisms additional to those found in Zhang et al. (2020 a ) werereported for swept wings. At small aspect ratios and low sweep angles the tip vortexdownwash effects still stabilise the wake while the weakening of the downwash withincreasing span allows the formation of unsteady vortex shedding. For higher sweep anglesthe source of three-dimensionality was shown to transition from the tip of the wing tomidspan where a pair of symmetric vortical structures is formed. Their mutually induceddownward velocity stabilising the wake. Therefore, three-dimensional midspan effects A. Burtsev et al. leading to formation of the stationary vortical structures allow steady flow formationat higher aspect ratios which would not be feasible on unswept wings. At higher aspectratios the midspan effects weaken near the tip leading to unsteady vortex shedding in thewing tip region. Finally, for high aspect ratio highly swept wings steady flow featuringrepetitive formation of the streamwise aligned finger-like vortices along the span ensues.Hayostek et al. (2021) conducted an experimental study of unswept, cantilevered, lowaspect ratio NACA 0015 wings ( Re = 600 and 1000) under wall effects closely integratedwith numerical simulations. Secondary vortical structures along the span were observed,arising from the interaction between the tip and root regions and the separated flow overthe wing. A global stability analysis at AR = 2 indicated that the least stable mode,while having structures in the tip and horseshoe vortex regions is most pronounced inthe wake.Despite the substantial improvement of understanding of complex vortical structuresarising in separated flows over finite aspect ratio wings that our experimental and large-scale computational efforts have offered, several key questions remain open and motivatethe present work. The origin of the unsteadiness of the trailing edge separation line,observed in the simulations of Zhang et al. (2020 a ) and those performed herein, remainsunexplained and the conjecture that this unsteadiness arises on account of a presentlyunknown flow eigenmode needs to be examined. Further, the frequency content andspatial structure of this, and possibly other modes existing in the flow both during thelinear regime and at nonlinear saturation needs to be documented and classified, e.g. interms of the respective energy content. From an aerodynamic point of view, it is necessaryto understand the origin of the concentration of frequencies in the lift coefficient spectrum,reported in Zhang et al. (2020 a ), at particular nondimensional frequencies; again, itcan be conjectured that a flow eigenmode with that frequency is responsible for thisobservation. Last but not least, the stabilization effects observed by Zhang et al. (2020 b )on swept wings, their relationship with vortex-induced downwash and the resulting effectson the wake behaviour need to be further explored. Here, it may be conjectured thatthe magnitude of the spanwise velocity component along the wing leading edge, whichchanges with changing angle of sweep, is responsible for this observation; this is anotherconjecture put to test in our present work.We perform linear TriGlobal modal analysis of separated flow over finite aspect ratiothree-dimensional wings, followed by data-driven modal analysis (Taira et al. § §
3. Results are reported in § inear modal instabilities around post-stall swept finite-aspect ratio wings
2. Theory
TriGlobal linear stability theory
The flow under consideration is governed by the nondimensional incompressible Navier-Stokes equations and continuity equations: ∂ t u + u · ∇ u = −∇ p + Re − ∇ u , ∇ · u = 0 , (2.1)where the Reynolds number, Re ≡ U ∞ c/ν , is defined by reference to the free-streamvelocity, U ∞ , the chord, c , and the kinematic viscosity, ν . The flow field can be expressedon an orthogonal coordinate system as a function of the unsteady velocity componentsand pressure q ( x , t ) = ( u , p ) T , (2.2)which are decomposed into a base flow component ¯ q and a superimposed small pertur-bation ˜ q , such that q = ¯ q + ε ˜ q , ε (cid:28) . (2.3)The approach followed to obtain steady stable, or stationary unstable base flows will bediscussed in § O (1) andneglecting O ( ε ) terms leads to the linearised Navier-Stokes equations (LNSE) ∂ t ˜ u + ¯ u · ∇ ˜ u + ˜ u · ∇ ¯ u = −∇ p + Re − ∇ ˜ u , ∇ · ˜ u = 0 . (2.4)For the incompressible flow of interest the pressure perturbation can be related to thevelocity perturbation through ˜ p = −∇ − ( ∇ · ( ¯ u · ∇ ˜ u + ˜ u · ∇ ¯ u )). Now the LNSE can bewritten compactly as the evolution operator L forming an initial value problem (IVP) ∂ t ˜u = L ˜u . (2.5)For steady basic flows, the separability between time and space coordinates in (2.5)permits introducing a Fourier decomposition in time of the general form ˜ u = ˆ u ( x ) e − iωt .Depending on the number of inhomogeneous spatial directions in the base flow analysedand the related number of periodic directions assumed different forms of the ansatz for˜ u can be used (Theofilis 2003; Juniper et al. q and the perturbation ˜ u are inhomogeneous functions of all three spatial coordinates giving the following ansatz˜ u ( x, y, z, t ) = ˆ u ( x, y, z ) e − iωt + c.c.. (2.6)Here, ˆ u is the amplitude function, and c.c. is a complex conjugate to ensure real-valuedperturbations. Substituting (2.6) into (2.5) leads to the TriGlobal eigenvalue problem(EVP) A ˆ u = − iω ˆ u . (2.7)The matrix A results from spatial discretisation of the operator L and comprises ofthe basic state ¯ q ( x ) and its spatial derivatives, as well as the Reynolds number as aparameter.The IVP solution in (2.5) can be obtained for finite time horizons t → τ over a timeinterval τ or in asymptotically large time t → ∞ . The modal approach considers thelimit t → ∞ and solves the three-dimensional (TriGlobal) eigenvalue problem (2.7). Inthe nonmodal approach, a finite-time horizon is introduced to describe the disturbancebehaviour at arbitrary times and the solution of the IVP is sought. This allows transientgrowth analysis where the energy growth of perturbations over a finite time interval ismonitored. In view of the absence of modal analysis in the problem at hand, the present A. Burtsev et al. work focuses on modal stability analysis. The TriGlobal EVP (2.7) is solved numerically,using time-stepping tools implemented in nektar++ (Cantwell et al. ω and the corresponding eigenvectors ˆ u , which are referred to as theglobal modes. The real and imaginary components of the complex eigenvalue ω = ω r + iω i correspond to the frequency and the growth/decay rate of the global mode.2.2. Residuals algorithm
When the unsteady equations of motion are advanced in time and a steady state isreached, the eigenvalue problem pertaining to that state can be solved. Alternatively,linear global modes can be computed by employing the residuals algorithm (Theofilis2000) during linear decay of perturbations; this algorithm has also been used here, atconditions applicable to steady laminar three-dimensional flow. In its simplest version,when the DNS has reached the state of monotonic exponential decay of residuals, it isstraightforward to extract both the steady state to which the simulation marches andthe amplitude function of the least-damped global mode pertinent to the flow by simplealgebraic expressions. Given transient DNS flow field solutions q , at arbitrary times t ,and t during the exponential decay of the signal, the damping rate of a monotonically-decaying signal is obtained from the logarithmic derivative ω i = ln q − ln q ∆t , (2.8)where ∆t = t − t and q , represent field values of any flow quantity at arbitrary probesin the flow, extracted from the flow filed at times t and t . A system of equations forthe flow at the two instances in time can be written (cid:26) q = ¯ q + ε ˆ q e ωt , q = ¯ q + ε ˆ q e ωt , (2.9)where ¯ q is the steady base state toward which the simulation is marching, ˆ q is theamplitude function of the least-damped global mode and σ is the damping rate. Ex-pressions pertinent to oscillatory signals have been presented by Theofilis & Colonius(2003). Solving this system of equations at arbitrary times t and t allows predictingthe steady state (well ahead of convergence in time) and recover the spatial structureof the amplitude function of the global mode. Several examples of successful applicationof the algorithm can be found in the literature (e.g. G´omez et al. et al. § Proper orthogonal decomposition
Data-driven modal analysis is performed in the nonlinear saturation regime, with theaim of understanding the physics of massively separated flow past the primary flowbifurcation. Proper orthogonal decomposition (POD) provides a decomposition of theflow field data into a minimum number of modes that capture the most energy of theflow at any one time. POD, referred to in the literature also as principal componentanalysis or empirical eigenfunction decomposition, is applicable to signals taken eitherduring linear growth of perturbation, as well as during nonlinear signal saturation. Itwas first introduced in fluid dynamics by Lumley (1967) and its reviews can be found inHolmes et al. (1996); Berkooz et al. (1993) and Taira et al. (2017, 2020).A given flow field q ( x , t ) is stacked as column vectors, each representing a fluctuatingportion of a flow component (for example, u, v, w ) with the time-averaged value ¯ q ( x )removed. This is done at a number of snapshots chosen such that the important flow inear modal instabilities around post-stall swept finite-aspect ratio wings χ ( t ) = q ( x , t ) − ¯ q ( x ) ∈ R n , t = t , t , ..., t m . (2.10)These column vectors are combined into an n by m data matrix X , where n is the physicalsize of the problem and m is the number of snapshots. POD seeks optimal basis vectors φ j to represent the flow field data q ( x ) with the least number of modes. This can bedetermined by solving the eigenvalue problem: XX T φ j = λ j φ j φ j ∈ R n λ (cid:62) ... (cid:62) λ n (cid:62) . (2.11)known as the spatial or classical POD. The issue with this method is the fact that a n × n eigenvalue problem has to be solved, with n the total degree of freedoms used inthe spatial discretization. This makes the computational requirements of classic PODdirectly dependent on the size of the data field and would be impractical to analyse flowpast 3D wings, where a significant portion of the wake has to be included in the data.However the temporal correlation matrix yields the same dominant spatial modes, whilegiving rise to a much smaller m × m eigenvalue problem as pointed out by Sirovich (1987): X T X ψ j = λ j ψ j , ψ j ∈ R m , m (cid:28) n. (2.12)The nonzero eigenvalues of (2.12) are same as those of (2.11), while the eigenvectors of(2.12) can be related to POD modes through φ j = X ψ j (cid:112) λ j ∈ R n , j = 1 , , ..., m. (2.13)When POD is used to describe flow development during linear instability, it is straight-forward to identify in the decomposition˜ q ( x , t ) = (cid:88) j a j ( t ) φ j ( x ) . (2.14)Following the terminology used by Aubry et al. (1991) the coefficients a j are termed chronos , as they represent the time behaviour of the mode, as opposed to the spatialstructures φ j , which are referred to as topos . When the signal analysed is in the phase oflinear growth/decay, chronos and topos can be respectively identified with the exponentialtime factor and the amplitude functions shown in (2.6), and POD has been utilisedto cross-verify damped global modes computed by solution of the TriGlobal eigenvalueproblem.Numerically POD coefficients are obtained by projecting the snapshots onto the PODmodes a j = Φ T χ j , j = 1 , , ..., r, (2.15)where Φ = [ φ , φ , ..., φ r ] and r can be a subset of the total number of snapshots m .
3. Numerical work
Geometry and mesh
The geometry under consideration is a untapered wing based on the symmetric NACA0015 airfoil with a sharp trailing edge and a straight cut wing tip. Taking advantage ofthe symmetry of the problem, half of the wing is considered as shown in figure 1. Thechord-based Reynolds number Re = 400 is held constant, while the wing sweep ( Λ ),aspect ratio ( sAR ) and angle of attack ( α ) are varied. Here, we use the symmetry aspect0 A. Burtsev et al.
Figure 1: Problem setup showing wing and the computational domain. The symmetrycondition is applied at the BACK plane. The half wing model is shown in grey and isnot to scale. Light grey indicates the opposite side of the wing when mirrored in thesymmetry plane and is shown for visualisation purposes only.ratio defined as sAR = b/ c , where b is the wingspan defined from wing tip to wing tipand c is the wing chord.It is important to take into account the order of the operations performed to constructa swept wing at an angle of attack. Translating the wing tip along x first and then rotatingthe wing by α or first rotating and then translating the wing tip along x will give twodifferent geometries. The latter order of operations would also result in a wing that has adihedral angle. To obtain the wing geometry used here the first a two-dimensional meshwas generated and the airfoil was rotated by α . It was then extruded along a vector { x, y, z } = { b/ Λ cos α, − b/ Λ sin α, b/ } . This is equivalent to rotating the wingabout an axis parallel to the leading edge and achieves a swept back wing without dihedralangle.The computational extent is ( x, y, z ) = [ −
15 : 20] × [ −
15 : 15] × [0 : 15]. The half wingwas meshed using Gmsh (Geuzaine & Remacle 2009), with a structured C type mesharound the wing. Macroscopic elements for a typical sAR = 4 straight wing mesh areshown in figure 2( a ), the closeup in 2( b ) shows refinement near the wing. Within eachelement both spectral code (discussed in § c ) shows the mesh used for data-driven analysis that will be explained below. Severalcomputational meshes having different number of macroscopic elements were tested withdifferent polynomial order p to ensure spacial and temporal convergence. A combinationof 46735 hexahedra and prisms as macroscopic elements for an sAR = 4 wing andpolynomial order of 5 was selected.For analysing the effect angle of attack, the aspect ratio and sweep angle are keptconstant at sAR = 4 and Λ = 0 ◦ . While the effects of sweep are analysed at a constantangle of attack of 22 ◦ at which the flow is separated with sweep angle varied between 0 ◦ and 30 ◦ for wings of sAR = 4 and 2. Length and velocity is nondimensionalized by wing inear modal instabilities around post-stall swept finite-aspect ratio wings XYZ ( a ) XYZ ( b )( c ) Figure 2: Computational mesh. ( a ) DNS mesh showing full domain ( b ) Close up of theDNS mesh near the airfoil ( c ) section of the mesh used for modal analysis. For clarityonly the macroscopic elements are shown in ( a ) and ( b ), while the internal field and themesh resulting from a high-order polynomial fitting are not shown. Similarly, only everysecond grid line is shown in ( c ). WING NORTH SOUTH WEST EAST FRONT BACK¯ u D U U U O N S¯ v D U U U O N S¯ w D U U U O N Sˆ u D D D D N N Sˆ v D D D D N N Sˆ w D D D D N N S
Table 1: Boundary conditions for the base flow ¯ q and perturbation ˆ q components.chord c and U ∞ respectively. Time refers to nondimensional convective time normalisedby c/U ∞ and the Strouhal number is defined as St = f c sin( α ) /U ∞ . For modal stabilityresults shown in further section, each perturbation component is normalised by its ownabsolute maximum, for example u = u/ max( | u | ).3.2. Solvers and boundary conditions
Direct numerical simulation is used to solve equations of motion using either of the nek5000 (Fischer et al. nektar++ (Cantwell et al. nektar++ was used for computing artificially stationery base flows and triglobal stabilityanalysis via time stepping, while nek5000 was chosen to generate flow fields for data-driven modal analysis due to its high speed.In order to close the systems of equations solved, the boundary conditions defined intable 1 are used. In this table, D denotes homogeneous Dirichlet, N denotes homogeneousNeumann, S denotes symmetry boundary and O denotes robust outflow (Dong et al. nektar++ , while at the inlet uniform flow U =2 A. Burtsev et al.
Figure 3: Comparison of v velocity signal between nek5000 and nektar++ for( sAR, Λ, α ) = (2 , ◦ , ◦ ) at ( x, y, z ) = (4 , , Present results Zhang et al. (2020a)Case α C L POD DMD C L StsAR = 4 12 ◦ ◦ sAR = 2 12 ◦ ◦ ◦ Table 2: Comparison of mean lift coefficient over unswept NACA 0015 wings at Re = 400with literature.( U ∞ , , T is imposed. The base flow solutions obtained by both codes were compared toensure that identical results are achieved. Figure 3 shows good agreement in variation ofvertical velocity with time for a given wing geometry between the two codes. In general,for the configurations considered good agreement between the two codes is achieved whenusing time steps ∆t (cid:54) × − and polynomial orders p (cid:62) C L ), presented in table 2 are compared toresults of Zhang et al. (2020 a ) showing good agreement. Further comparisons betweenthe CharLES and nektar++ solvers have been presented in He et al. (2019 b ) and Zhang et al. (2020 a ). The frequencies of the dominant modes as will be discussed in detail in § et al. (2020 a ). Very goodagreement between the dominant frequencies of the C L signal and dominant modes canbe seen. inear modal instabilities around post-stall swept finite-aspect ratio wings Base flows and signal processing
As discussed in § § et al. et al. nektar++ , has been used to compute anartificially stationary, unstable base state that is used for subsequent modal analysis.This methodology was shown to recover amplified global modes of a sphere (He et al. a ). SFD uses filtering and control of unstable temporal frequencies in the flow, thetime continuous formulation in nektar++ can be expressed as (cid:26) ˙ q = N S ( q ) − γ ( q − ¯ q ) , ˙¯ q = ( q − ¯ q ) /∆ (3.1)where q represents the problem unknown(s), the dot represents the time derivative, N S represents the Navier-Stokes equations, γ ∈ R + is the control coefficient, ¯ q is a filteredversion of q , and ∆ ∈ R + ∗ is the filter width of a first-order low-pass time filter (Jordi et al. γ and ∆ affects the convergence to the steady-statesolution when q = ¯ q . If the dominant mode is known, for example through DMD, andspecified as input one can adjust the filter parameters to accelerate convergence.3.4. Modal analyses
TriGlobal instability analysis was performed using time-stepper algorithm imple-mented in nektar++ (Cantwell et al. − using thefully three-dimensional base flow. POD and DMD were carried out using an in houseFortran code that uses LAPACK (Anderson et al. nek5000 .Snapshots were recorded every 0.1 unit of characteristic time once the flow has reachedperiodic shedding, with 300 snapshots collected in the interval from 40 to 70 time-units. Each snapshot was interpolated using built in spectral interpolation of nek5000 onto the equidistant mesh shown in figure 2( c ) such that scaling due to changing cellvolume for each date point does not have to be taken into account. This simplifies theconstruction of the covariance matrix. The new mesh was generated by extruding anelliptic two-dimensional mesh around the airfoil, which resulted in areas of changing cellvolume and some mesh distortion directly above and below the wing, especially belowthe pressure side. Hyperbolic mesh could have offered better control over the cell shapeat the boundary but would not allow to define the far field boundary itself. It shouldbe emphasised that data-driven methods do not require computation of gradients, sothe effects of mesh non-orthogonality at the boundary is minimal, especially since thereare minimal changes in the base flow at these locations. Finally, the mode structuresassociated with massively separated flows such as the high- α wings considered hereare expected to lie in the wake where cell volume is constant. This mesh extended 7 c downstream of the training edge, 1 . c above and below the wing, and 1 c beyond the wing4 A. Burtsev et al.
Case Variables n x × n y × n z n Memory (GB) sAR = 2 POD u, v, w, p × ×
78 5 013 216 11.2DMD u, v, w sAR = 4 POD u, v, w, p × ×
130 8 355 360 18.7DMD u, v, w
Table 3: Mesh and data matrix sizes for modal analysis.
WakeInteractionTip α = 10 ◦ α = 14 ◦ α = 18 ◦ α = 22 ◦ Figure 4: Isocontours of Q -criterion Q = 1 coloured by streamwise vorticity showing theeffect of angle of attack on instantaneous DNS solution at sAR = 4, Λ = 0 ◦ .tip in z direction. With 26 points per characteristic length this resulted in 412 × × sAR = 4 and 412 × ×
78 mesh for sAR = 2 wings. Three velocity componentswere stacked to generate the data vector for DMD with the addition of pressure forPOD the resulting data matrix size is shown in table 3.
4. Results
Wake dynamics
The evolution of the flow over the unswept sAR = 4 wing with angle of attack isshown in figure 4. The vortical structures of the three-dimensional wake over unsweptwings is in agreement with the DNS results of Zhang et al. (2020 a ). For the separatedflows at high angles of attack, three regions can be identified behind the wing. These arethe wake region close to the symmetry plane containing mostly spanwise vortices, the tipregion containing the tip vortex and the interaction region in between where a pair ofcounter-rotating wake vortices are connected into a closed loop by braid-like structures(Zhang et al. a ).At α = 10 ◦ , shown in figure 4, the flow is steady with separation occurring at inear modal instabilities around post-stall swept finite-aspect ratio wings ( a ) ( b ) ( c ) Figure 5: Snapshots of DNS for (
AR, Λ, α ) = (4 , ◦ , ◦ ) at three different times ( a ) t = 50, ( b ) t = 60, ( c ) t = 71. Grey/white contour shows the surface-friction lines on thewing. Black lines show the streamlines in slices corresponding to spanwise location at z = 0, 1, 2 and 3. The separation zone is highighted by blue lines. Vortices are visualisedby isocontour of Q = 1 and its center location and rotational direction are highlightedby red lines and arrows, respectively.approximately two-thirds of the chord, independently of the angle of sweep. At α = 14 ◦ ,an unsteady wake is formed, the shed vortices being practically parallel to the trailingedge of the wing. The separation location moves approximately at half-chord and thespanwise region of the flow affected by the tip vortex is reduced, with the separationbubble extending closer to the tip. At the higher angles of attack of 18 ◦ and 22 ◦ , alsoshown in figure 4, the three distinct regions first identified by Zhang et al. (2020 a )develop: these regions are the wake, consisting of spanwise vortices near the symmetryplane, the essentially steady tip vortex, and the interaction region between the wake andtip characterised by the braid-like vortices.With increasing α the separation location moves closer to the leading edge and the tipvortex becomes stronger. At α = 18 ◦ and 22 ◦ the reattachment line close to the trailingedge of the airfoil starts to oscillate periodically, as can be seen in the surface streamlinesmovies available as online supplementary material. The origin of this motion is presentedin detail at a single set of parameters ( AR, Λ, α ) = (4 , ◦ , ◦ ), in figure 5. Vortexcenters and surface streamlines are plotted at three instances in time. When a counter-clockwise (on the xy plane, as viewed in the top row of figure 5) trailing edge vortex startsto separate from the wing near the symmetry plane the reattachment line (shown in blue)is pushed further away from the trailing edge of the wing as in figure 5( a ) . The sheddingof a clockwise leading edge vortex causes the reattachment location to return closer to thetrailing edge as in figure 5( b ). Alternating shedding of clockwise and counter-clockwisevortices that is non-uniform in the spanwise direction leads to the wavy-like motion of thereattachment line. Further downstream, each pair of counter-rotating spanwise vortices6 A. Burtsev et al. sAR = 4 sAR = 2 Λ = 0 ◦ Λ = 5 ◦ Λ = 10 ◦ Λ = 15 ◦ Λ = 20 ◦ Λ = 25 ◦ Λ = 30 ◦ Figure 6: Isocontours of Q -criterion Q = 1 coloured by streamwise vorticity showing theeffect of sweep on instantaneous DNS solution at Re = 400 and two wing aspect ratios.Top to bottom are Λ = 0 ◦ − ◦ with increment of 5 ◦ . For clarity Λ = 30 ◦ is shown with Q = 0 . (cid:54) z (cid:54) .
5, forming thevortex loops discussed by Zhang et al. (2020 a ).The frequency of these motions was measured by extracting a line parallel to thetrailing edge of the wing, located 0 . c above the trailing edge, and performing a fastFourier transform of the signal on this line. The dominant frequencies based on the timesignal of vertical velocity component were found to be St = 0 . .
139 and 0 .
143 for the inear modal instabilities around post-stall swept finite-aspect ratio wings (a) ( b ) Figure 7: Effects of spanwise and vortex induced flow ( a ) tip vortex region ( sAR, Λ, α ) =(4 , ◦ , ◦ ) one c downstream of wing trailing edge showing vortex center and core, ( b )comparison of components parallel to the wing leading edge of free-stream velocity U ∞ and vortex induced velocity U ind for low angles of sweep. sAR = 4 wing at sweep angles of Λ = 0 ◦ , 5 ◦ and 10 ◦ respectively. The spanwise locationsof these peaks of the power spectral density correspond well with the extremities of thereattachment line motion shown in figure 5.The effect of sweep on the flow over the longer sAR = 4 and the shorter sAR = 2 wingare shown in figure 6. As the longer wing is swept back, the braid-like region is movedcloser to the wing tip by the increased spanwise cross flow, which results in the tip vortexbecoming noticeably less steady. Interestingly, stronger periodic vortices are observed inthe near wake behind the wing at ( sAR, Λ ) = (4 , ◦ ), shown in figure 6, that meet atan oblique angle on the symmetry plane and form chevron-like patterns reminiscent ofthose reported in the wake of large aspect ratio cylinders by Williamson (1989). Thesevortices extend about four chord lengths downstream of the wing with four cores clearlyvisible before the wake changes to the two distinct regions seen behind the straight wing.In order to analyse the cause of this inherently two-dimensional shedding at Λ = 5 ◦ , thespanwise component of velocity (parallel to wing leading edge) caused by the sweep angleand the velocity induced by the tip vortex along the opposite direction were compared.The circulation Γ of the tip vortex was calculated from time-averaged flow at one chorddownstream of the wing. Figure 7(a) shows the location of the tip vortex relative to theprojections of wing LE and trailing edge (TE) as well as the vortex core size defined hereas the contour of maximum azimuthal velocity for the sweep angle of 10 ◦ . The vortex coreis located in the wake behind the wing with its spanwise position moving inboard furtherdownstream of the wing. For a point P outside the vortex core the induced velocity canbe approximated as U ind = Γ/ πr where r is the distance between vortex center C and P . The induced velocity by the tip vortex is compared to the spanwise component ofthe free-stream flow upstream of the wing ( U ∞ sin( Λ )) in figure 7( b ) as a function ofthe sweep angle Λ . With increasing sweep the strength of the tip vortex and hence U ind slightly decreases while the spanwise component of the free-stream that is opposite indirection increases significantly. Near Λ = 5 ◦ the magnitudes of the spanwise velocityof both components are similar, which supports the previous finding that at this sweepangle the wake evolves in a quasi-two-dimensional manner.The oscillations of the reattachment line seen on the unswept case are also observedwith sweep. However, at Λ = 10 ◦ a steady region of attached flow forms at the inboardtrailing edge of the wing. There is a qualitative change in the structure sAR = 4 wingwake as sweep angle reaches Λ = 15 ◦ . The periodic vortices passing though the symmetry8 A. Burtsev et al. plane are no longer visible, and the wake now consists of two series of braid-like vorticesforming behind the half-wing, that do not pass through the symmetry plane. The tipvortex is now less pronounced and clearly unsteady. As the sweep angle increases furtherto 20 ◦ , the tip vortex is no longer visible and the wake is dominated by the two braid-likeregions. The size of the attached flow region at the inboard root of the wing increaseswith sweep and at Λ = 20 ◦ the separation bubble splits into two halves on either side ofthe symmetric wing, with attached flow forming in the middle. Interestingly, the presenceof such region of attached flow at the root of a swept wing was also reported by Visbal& Garmann (2019) for turbulent flow at much higher Reynolds numbers. At Λ = 25 ◦ thebraid-like regions become narrower and vortices extending from the inboard section ofthe wing into the wake behind the tip are starting to form; these structures are sometimesreferred to as ”ram’s horn” vortices (Black 1956). The growth of these structures resultsin another fundamental change in the wake at the maximum sweep of Λ = 30 ◦ examinedin our work. A ”ram’s horn” vortex is generated on the suction side of the wing close tothe leading edge at the symmetry plane and a stronger counter-rotating vortex emanatesfrom the trailing edge. These two vortices form a closed structure further downstreamthat starts to shed with the shed vortices reminiscent of the hairpin-vortex. Overall higherangle of sweep has a stabilising effect on the flow. It was shown by Zhang et al. (2020 b )that as the sweep is further increased the flow wake will become steady at Λ ≈ ◦ .The effects of sweep are qualitatively analogous on the shorter wing. The breakdown ofvortices at the symmetry plane to two braid-like regions happens at Λ = 20 ◦ rather then15 ◦ . Horn-like vortices similar to the larger wing form at Λ = 30 ◦ . Unlike the sAR = 4case, these structures do not form hairpin vortices further in the wake, but rather mergeinto a spiral vortex, an effect attributed to the stronger tip effects on the shorter wing.4.2. Linear global modes
TriGlobal modal linear stability analysis was performed at conditions at which steadyflow naturally existed or could be computed using the selective frequency dampingmethod discussed in § Q -critetion.The most unstable mode for the unswept sAR = 4 wing takes the form of vorticalstructures localised at the symmetry plane. The components of the normalised pertur-bation velocity plotted in figure 9 show the presence of two branches in the ˆ u and ˆ w perturbation velocity components, originating from the leading and trailing edges. Theleading mode of the shorter sAR = 2 wing is more amplified but is qualitatively identicalto the sAR = 4 result. With increasing sweep angle, the peak of the mode structure movesalong the spanwise direction towards the tip. Qualitatively, the structures observed at Λ = 0 ◦ , 5 ◦ , and 10 ◦ represent the same wake mode as it is evolving with sweep angleand are highlighted in blue in figure 8. The difference in the non-dimensional frequencyof these modes remains below 10% in this range of sweep angles.There is a change in the structure of the leading mode at Λ = 15 ◦ . Although thestructures of vortices close to the wing still resemble the lower sweep results, as visible inear modal instabilities around post-stall swept finite-aspect ratio wings α = 22 ◦ .Most unstable or least stable modes are indicated with filled symbols with correspondingspatial structures visualised with contours of Q = ± St isdefined as St = ω r c sin α/ πU ∞ .in figure 8, closer examination of the velocity components reveals significant changes.While mode structures at the lower sweep angles were mostly streamwise periodic, withlittle change in the x direction, the leading mode at Λ = 15 ◦ transitions to a vorticalinstability further downstream of the wing as evident from the vorticity plots in figure11( a, c ) where positive ω x indicates clockwise rotation (counter-rotating with respectto the tip vortex). At the maximum considered sweep of Λ = 30 ◦ the leading modeis stable and takes the form of a vortical instability that grows with x downstream ofthe wing. Unlike the Λ = 15 ◦ mode that exhibits vortical structures centred aroundcounter rotating vortices of the SFD base flow, the Λ = 30 ◦ global mode clearly showsthe influence of the tip vortex, as can be seen in yz -plane at x = 20 shown in figure 11( d ).The change of amplification rate with sweep is significant, showing an overall trend oflower amplification and fewer unstable modes with increased sweep, until no amplifiedmodes at all are observed at Λ = 30 ◦ with one notable exception of Λ = 10 ◦ .In order to validate the above findings, full direct numerical simulations have beeninitialised using selected global modes as initial condition. As an example, the TriGlobalmode at ( sAR, Λ, α ) = (4 , ◦ , ◦ ) is superimposed at an amplitude O (10 − ) onto theSFD-obtained base flow. The time history of the ˆ v component of the perturbation at aprobe located at the symmetry plane is shown in figure 12( a ), where the linear growthregion leading to nonlinear saturation is clearly visible. The slope, measured from thelogarithmic plot in figure 12( b ) is equal to 0 . Comparison with local theory
The relative simplicity and low computational cost of local stability analysis make itan attractive tool to compare its results with those of the computationally challengingTriGlobal analysis. The obvious shortcoming of local theory is the assumption of flowhomogeneity in two spatial directions; however it might be argued that locally this0
A. Burtsev et al. sAR = 4 Λ : ˆ u ˆ v ˆ w ◦ ◦ ◦ sAR = 20 ◦ Figure 9: Most unstable TriGlobal modes for sAR = 4 and 2 wings at different angles ofsweep shown with contours of ˆ u , ˆ v , and ˆ w at ± . et al. b ; Paredes et al. u computed by the SFDmethod. The larger aspect ratio wing was chosen, and profiles were extracted at themidspan for the unswept wing, so as to stay away of the tip vortex and the three-dimensional vortical structures at quarter-span. For the low sweep angles of Λ = 5 ◦ and10 ◦ , spanwise locations of z/c = 1 and 2 were chosen based on the maximum velocitydeficit in the streamwise velocity profile. Due to high three-dimensionality of both thebase flow and the global mode at larger sweep angles, local analysis was not applied atthose conditions. The perturbation was assumed to be homogeneous in the streamwisedirection and the analysis was carried out at spawnwise wavenumber, β of zero.Spatial analysis results in figure 13( b ) show two unstable modes. For the unswept case inear modal instabilities around post-stall swept finite-aspect ratio wings Λ = 15 ◦ Λ = 30 ◦ ˆ u ˆ v ˆ v Figure 10: Most unstable TriGlobal mode for ( sAR, Λ, α ) = (4 , ◦ , ◦ ) and least stablemode for (4 , ◦ , ◦ ) ˆ u , ˆ v and ˆ w at ± . ( a ) ( b )( c ) ( d ) Figure 11: Streamwise vorticity of leading global modes for ( a, c ) ( sAR, Λ, α ) =(4 , ◦ , ◦ ) and ( b, d ) (4 , ◦ , ◦ ) plotted at x = 5 in ( a, b ) and x = 20 in ( c, d ). Thedashed line indicates the projection of the wing trailing edge onto the the slice.2 A. Burtsev et al. ( a ) ( b ) Figure 12: Growth of the perturbation for ( sAR, Λ, α ) = (4 , ◦ , ◦ ) ( a ) perturbationwith time showing linear growth and nonlinear saturation ( b ) close up plotted with logscale showing exponential growth. ( b ) Figure 13: Local spatial stability analysis results for small sweep Λ = 0 ◦ − ◦ showingvariation of the spatial amplification rate α i with frequency ω .mode 1 reaches maximum amplification at ω ≈ α r approaches zero while mode2 peaks at ω ≈ .
5. For Λ = 5 ◦ and 10 ◦ mode 1 reaches maximum amplification at ω ≈ ω = 1 . Λ = 5 ◦ and 10 ◦ and amplification rate of the global mode iswithin 20% and 24% respectively. On the other hand, mode 2 has an order-of-magnitudelarger amplification than that of the global mode. The frequency of the global modeslies approximately in between the frequencies corresponding to maximum amplificationof the local modes 1 and 2.In summary, and somewhat unsurprisingly due to the highly inhomogeneous nature ofthe base flows, this comparison shows that, despite some qualitative agreement between inear modal instabilities around post-stall swept finite-aspect ratio wings ( a ) ( b ) ( c ) Figure 14: Comparison of the spatial local analysis eigenfunctions with correspondingslices through the global modes at ω equal to the frequency of the global mode ( a ) Λ = 0 ◦ ω = 1 . b ) Λ = 5 ◦ , ω = 1 .
658 ( c ) Λ = 10 ◦ ω = 1 . ( a ) ( b ) ( c ) Figure 15: Modes of ( sAR, Λ, α ) = (4 , ◦ , ◦ ) wing ( a ) global mode obtained withresidual algorithm ( b ) most energetic POD mode based on sample 10 (cid:54) t (cid:54)
40 ( c )damped wake instability POD mode based on sample 40 (cid:54) t (cid:54)
70. Contours of spanwisevorticity plotted at a xz slice located at the wing trailing edge.global and local eigenfunction results, large differences in quantitative predictions of theamplification rates and frequencies exist. These differences underline the inability of localanalysis to adequately capture the physics of linear instability in this configuration evenfor locations close to the symmetry plane and away from the main sources of three-dimensionality. 4.4. Wake modes in the nonlinear saturation regime
Effect of angle of attack
Data-driven analysis was performed in the nonlinear saturation regime, in order toanalyse the effects of wing geometry on the wake. Both POD as described in § sAR, Λ, α ) = (4 , ◦ , ◦ ) is first considered. Two setsof snapshots were taken: one in the region of a exponential decay, at 10 (cid:54) t (cid:54) (cid:54) t (cid:54)
70. At this set of parameters theresidual algorithm is used and its results, shown in figure 15( a ), are compared to dominantmode delivered by POD, the latter applied in the two aforementioned time-windows. Itcan be seen that the leading POD mode during exponential decay corresponds to thedecaying wake instability/tip-vortex structure shown by the global mode. The dampedwake instability recovered in the α = 10 ◦ case on the sAR = 4 finite aspect ratio wingis a Kelvin–Helmholtz mode akin to the two-dimensional mode reported by He et al. (2017 a ) in their analysis of spanwise homogeneous base flow.The spatial structure of the dominant mode as the angle of attack is increased, is shownin figure 16. The results correspond to modes obtained by POD analysis, since DMD givespractically identical, in both frequencies, shown in table 4, and spatial structure of the4 A. Burtsev et al. α : ˆ u ˆ v ˆ w ◦ ◦ ◦ ◦ Figure 16: Evolution of the leading mode for sAR = 4 Λ = 0 ◦ wing with angle of attackfor α = 10 ◦ , 14 ◦ , 18 ◦ and 22 ◦ visualised with isosurfaces of ˆ u in ˆ v and ˆ w at ± . α = 10 ◦ these structures are parallel to the trailing edge of the wing, while as the angleof attack increases they are deformed, following the curvature of the vortex cores in theDNS results now meeting at an oblique angle to the symmetry plane. It can be seen thatat α = 14 ◦ the mode still resembles Kelvin–Helmholtz instability, while at higher anglesof attack the structure is increasingly more complex, as seen in the side views of ˆ u, ˆ v, ˆ w ,shown in figure 16. This increased three-dimensionality and the peaks in ˆ u, ˆ v developingat about 1 / b indicate the significant effect of the interaction region vortices on the mode.We refer to this structure as the interaction mode . This mode contains up to 47% of theflow kinetic energy as shown in table 4. It should be noted that POD modes are grouped inear modal instabilities around post-stall swept finite-aspect ratio wings ( a ) ( b ) Figure 17: Effects of angle of attack: ( a ) spectrum of POD and ( b ) variation of theinteraction mode frequency for unswept sAR = 4 wing.in pairs, with nearly identical frequencies and similar energies, the table shows only thefirst mode of the pair, similarly when classifying modes only the first member of the pairis shown. When combined, the two modes of the POD pair would actually contain twicethe E k quoted in this table.The dependence of the three most dominant modes at sAR = 4 on the angle of attackis visualised in figure 18, using contours of Q -criterion computed from the normalisedperturbation components. The interaction mode (IM) shown at α = 14 ◦ , ◦ and 22 ◦ infigure 18, clearly shows spanwise vortical structures in the wake. The pattern of theseIM structures near the wing from α (cid:62) ◦ shows a clear connection to the curvature ofvortex core lines shown in figure 5. It is also worth mentioning that the frequency of theIM is in close agreement with the dominant frequencies of the reattachment line motiondiscussed in § (cid:54) z (cid:54) . α inthe DNS results. Higher harmonics of the dominant IM also appear in the spectrum withthe most energetic one labelled IM2 shown in table 4. The wake mode (WM) appearsto be concentrated in the wake close to the symmetry plane, as also seen in 18. It has alower frequency compared to that of the IM and shows the least contribution from thebraid-like structures of the interaction region. Also, unlike the IM and its harmonic IM2,the structure of the WM mode does not follow the vortex line curvature. The frequenciesof IM and WM at α = 22 ◦ as shown in table 4 match very well with the two distinctfrequency peaks of the lift coefficient signal observed by Zhang et al. (2020 a ) at sameconditions on sAR = 4 wing.Finally, the dependence of POD eigenvalues on angle of attack is plotted in figure 17( a ).At α = 10 ◦ the results for analysis in both linear and nonlinear regimes show that thefirst few modes contain the majority of flow energy. The ratio λ j / (cid:80) m λ m between themode eigenvalue and the sum of eigenvalues decreases for a given mode with increasing α meaning that higher order modes become more important as their relative energycontribution increases. The Strouhal number of the dominant mode increases with angleof attack as shown in figure 17( b ) indicating an increase of dominant mode frequency.4.4.2. Effect of sweep
Turning our attention to the effect of angle of sweep, the effect of this parameter on thethree dominant modes, the interaction mode (IM), its second harmonic (IM2) and the6
A. Burtsev et al. α ◦ ◦ ◦ ◦ K-H (IM) St DMD g j -0.3962 0.0001 0.0000 -0.0003 St POD E k St DMD - 0.2434 0.2716 0.2789 g j - 0.0099 0.0001 -0.0011 St POD - 0.2433 0.2716 0.2794 E k - 2.0291 3.7739 3.4648WM St DMD - 0.1151 0.1489 0.1600 g j - -0.0194 -0.0043 -0.0093 St POD - 0.1159 0.1484 0.1608 E k - 0.8722 0.2957 2.5794 Table 4: Change of POD and DMD modes with angle of attack for sAR = 4 straightwing. Strouhal numbers according to POD and DMD, growth rate ( g j ), and percentageof kinetic energy ( E k ). IM WM St = 0 . St = 0 . α = 14 ◦ St = 0 . St = 0 . α = 18 ◦ St = 0 . St = 0 . α = 22 ◦ Figure 18: Modes of Λ = 0 ◦ , sAR = 4 wings visualised with contours of Q = 1 andcoloured by streamwise vorticity for α = 14 ◦ , 18 ◦ , and 22 ◦ , showing the interactionmode (IM) and the wake mode (WM). inear modal instabilities around post-stall swept finite-aspect ratio wings Λ ◦ ◦ ◦ ◦ ◦ ◦ ◦ IM St POD E k St POD E k St POD E k St POD - - - - 0.00884 0.0084 0.0085 E k - - - - 2.8337 3.2341 39.8751 Table 5: Change of POD modes with sweep angle for sAR = 4 wing showing Strouhalnumbers according to POD coefficients and percentage of kinetic energy.wake mode (WM) identified in the previous subsection, is monitored at a fixed angle ofattack α = 22 ◦ . Table 5 shows the Strouhal numbers obtained from POD coefficients andthe percentage of kinetic energy. The dominant interaction mode, having St = 0 .
122 at α = 14 ◦ , can be seen in figure 19, its structure reflecting the changes in the correspondingflow fields as the angle of sweep increases. At Λ = 10 ◦ the vortices of the IM no longerpass through the symmetry plane forming two structures on either side of the symmetricwing.The wake mode (WM), having St ≈ .
18, is shown in figure 19. In a manner analogousto the unswept wing, the WM mode exhibits periodic structures that are strongest inthe wake. At Λ = 5 ◦ it has the highest percentage of kinetic energy and overtakes IM2 asthe second most energetic mode. At Λ = 10 ◦ , WM structures seem to shrink and appearonly very close to the symmetry plane. This is the highest sweep angle at which thismode has been detected.At the sweep angle of Λ = 15 ◦ there is a change in the wake from shedding centredon the symmetry plane to shedding of the braid-like vortices closer to the tip vortex.This is reflected in the interaction mode shown in figure 19, as the mode structures areno longer passing through the symmetry plane. The periodic structures originating fromthe outboard side of the wing are terminated on either end by spiral vortical structuresthat correspond to the location of interaction region vortices in the DNS results. Furtherdownstream the periodic structures break down, leaving only the spiral vortices. Thefrequency spectrum at this angle of sweep is dominated by harmonics of the leadingmode, similar to the frequencies of sAR = 2 modes as will be seen in § Λ = 20 ◦ the spatialstructure of the dominant interaction mode (IM) is shifted closer to the wing tip, butis still qualitatively similar to the Λ = 15 ◦ case. In contrast, the mode at the highestsweep of Λ = 30 ◦ is fundamentally different and is dominated by tip effects. Note thatresults for the highest sweep angle are not to scale with the rest. This is because thestreamwise extent of the domain used for modal analysis for this case was doubled inorder to capture the mode structure resulting from the elongated vortices in the wake;these structures start to shed much further from the wing. A new low-frequency modewith St ≈ .
008 appears in the POD spectrum starting at Λ = 20 ◦ and rapidly gainsenergy from 2 .
83% at this angle of sweep to 39 .
86% at Λ = 30 ◦ , where it overtakes IMas the mode containing the highest energy. Since the structure of this mode consistsof streamwise vortices emanating from near the wing tip, as shown in 20, it is referred8 A. Burtsev et al.
IM WM St = 0 . St = 0 . Λ = 5 ◦ St = 0 . St = 0 . Λ = 10 ◦ St = 0 . Λ = 15 ◦ Figure 19: Modes of Λ = 0 ◦ , sAR = 4 wings visualised with contours of Q = 1 andcoloured by streamwise vorticity for α = 14 ◦ , 18 ◦ , and 22 ◦ , showing the interactionmode (IM) and the wake mode (WM).to as the streamwise mode (SM). A pair of elongated structures propagating from thewing downstream visible at Λ = 30 ◦ are associated with the counter rotating streamwisevortices observed in both the instantaneous and the time averaged DNS results. It shouldbe mentioned that it has not been possible to recover these low-frequency modes byDMD at Λ = 20 ◦ , ◦ ; DMD predicts the dominant mode at the highest sweep to be theinteraction mode.The dependence of the power spectral density of the time coefficients of the dominantPOD modes with wing sweep is plotted in figure 21. The interaction mode (IM), wakemode (WM) and the streamwise mode (SM) are shown, but the harmonics of IM are notincluded for clarity. Although POD modes can in principle contain multiple frequencies,in most of the cases examined modes have been found to correspond to a distinct singlefrequency, which explains the good overall agreement with DMD results. However, thewake mode on the sAR = 4 wing shows strong secondary peak in power spectral densityat the dominant frequency of the IM at sweep angle of 5 ◦ . Similarly, the IM shows asecondary peak at the dominant frequency of WM at these conditions as can be seen infigure 21( a ). The spatial structures of IM recovered by DMD are virtually identical tothat of POD while the WM structures are slightly different, showing less influence fromthe interaction region. This is to be expected as by definition DMD modes contain onlyone frequency and show the proper wake mode, whereas POD results were contaminatedby the interaction mode. Hence the structures WM shown in figures 18 and 19 are thatproduced by DMD. The power spectral density of SM based on the POD time coefficientsfor different sweeps is shown in figure 21( a ). This mode has predominantly low frequency inear modal instabilities around post-stall swept finite-aspect ratio wings IM SM St = 0 . St = 0 . Λ = 20 ◦ St = 0 . St = 0 . Λ = 25 ◦ St = 0 . St = 0 . Λ = 30 ◦ Figure 20: Modes of high Λ , sAR = 4 , α = 22 ◦ wings visualised with contours of Q = 1for Λ = 20 ◦ , 25 ◦ and 30 ◦ , showing interaction mode (IM) and wake mode WM. Λ = 30 ◦ results are not to scale with the rest of the figure. ( a ) ( b ) Figure 21: The frequency spectra of POD coefficients corresponding to dominant modeswith sweep ( a ) for sAR = 4 and ( b ) for sAR = 2. The spectrum for each sweep angle isnormalised by its peak and harmonics of the dominant IM are not plotted.content ( St ≈ .
01) with the exception of the sAR = 4 wing at sweep of 25 ◦ where apeak at St ≈ .
22 is visible.4.4.3.
Effect of aspect ratio
Finally, the effect of wing sweep on the leading modes is discussed at the sAR = 2.Frequencies of the three most dominant modes at low angles of attack are presented intable 6 and the corresponding spatial structures are shown in figure 23. At 0 ◦ (cid:54) Λ (cid:54) ◦ the overall structure of modes IM and its harmonic IM2 is analogous with that foundat sAR = 4. However, unlike the larger wing, no true wake modes have been found,0 A. Burtsev et al. ( b ) Figure 22: Lissajous curves of higer harmonics of the dominant POD mode for( sAR, Λ, α ) = (2 , ◦ , ◦ ).that would be independent of the interaction region vortices discussed in section 4.1.This can be attributed to the stronger downwash effects of the tip vortex on the shortwing. As a result, the spectrum of physically relevant modes (contributing to E k (cid:62) . sAR, Λ, α ) = (2 , ◦ , ◦ ). The first line, labelled IM,corresponds to the coefficient of the second member of the first POD pair plotted againstthat of the IM. Since POD mode pairs have the same amplitude and a phase differenceof 90 ◦ the resulting pattern is a circle. The ratio of frequencies determines the numberof ”lobes” in the curve with rational ratios resulting in closed curves. It can be seen thatlines corresponding to the second, third and fourth harmonics of IM have two, three andfour lobes respectively. The smaller relative size of the patterns corresponding to higherharmonics reflects the ratio of amplitudes, which are smaller for higher order modes.The spatial structures of sAR = 2 wing modes at high angles of sweep are visualised infigure 24, where higher harmonics of IM are not shown for clarity. In a manner analogousto the larger wing, the structures of the dominant mode show spiral vortices at spatiallocations where vortices are present in the flow; these vortices move towards the tip withincreased sweep angles. Also analogously to the sAR = 4 case, a low frequency SM modewith St ≈ .
008 appears at Λ (cid:62) ◦ and is shown in figure 20. At Λ = 30 ◦ the structureof SM reflects the presence of the ”ram’s horn” vortices in the flow. However, unlike the sAR = 4 wing this mode does not become the dominant at Λ = 30 ◦ . Instead, on the lowaspect ratio wing, the IM clearly shows strong tip effects and retains the largest portionof kinetic energy. This supports the point made in § sAR = 2 wingcan be compared with that on the larger wing in figure 21( b ). The dependence of thePOD eigenvalues with Λ at the two aspect ratios is shown in figure 25( a ). The spectrumat sAR = 2 is more ordered with each pair of modes representing significantly less energy inear modal instabilities around post-stall swept finite-aspect ratio wings IM IM2 St = 0 . St = 0 . Λ = 0 ◦ St = 0 . St = 0 . Λ = 5 ◦ St = 0 . St = 0 . Λ = 10 ◦ St = 0 . St = 0 . Λ = 15 ◦ Figure 23: Modes of low
Λ sAR = 2 wings visualised with contours of Q = 1 for Λ = 0 ◦ ,5 ◦ , 10 ◦ , and 15 ◦ showing interaction mode (IM) and its harmonic (IM2). IM SM St = 0 . St = 0 . Λ = 20 ◦ St = 0 . St = 0 . Λ = 25 ◦ St = 0 . St = 0 . Λ = 30 ◦ Figure 24: Modes of high Λ , sAR = 2 wings visualised with contours of Q = 1 for Λ = 20 ◦ , 25 ◦ and 30 ◦ showing interaction mode (IM) and streamwise mode (SM).2 A. Burtsev et al. Λ ◦ ◦ ◦ ◦ ◦ ◦ ◦ IM St POD E k St POD E k St POD E k St POD - - - - 0.0086 0.0082 0.0083 E k - - - - 13.5899 0.0983 0.8054 Table 6: Change of POD modes with sweep angle for sAR = 2 wing showing Strouhalnumbers according to POD coefficients and percentage of kinetic energy. ( a ) ( b ) Figure 25: Effects of sweep at α = 22 ◦ : ( a ) variation of the POD spectrum with increasing Λ for sAR = 4 and 2 ( b ) variation of the dominant interaction mode frequency with sweepangle and aspect ratio.compared to the larger wing. Interestingly, at sAR = 4 , Λ = 15 ◦ the POD eigenvaluesshow a pattern close to the sAR = 2 behaviour.The different behaviour of the dominant interaction mode with sweep for the twoaspect ratios is summarised in figure 25(b). For the larger sAR = 4 wing the frequencyof the dominant mode does not change significantly with sweep at first, staying within2% of the unswept wing value in the range 0 ◦ (cid:54) Λ < ◦ , but from Λ = 20 ◦ onwardit decreases systematically, dropping by approximately 25% at the extreme end of thesweep scale with the rate of decrease seemingly reducing between 25 ◦ and 30 ◦ . In thecase of the shorter wing, a similar decay of the dominant frequency of IM of about 22%is observed, however this starts at a lower sweep of Λ = 5 ◦ with a constant gradientmaintained up to Λ = 20 ◦ ; after this sweep angle there is an increase of frequency atto Λ = 25 ◦ . Both sAR cases show a small increase in the frequency of IM before thestart of the decrease stage, that appears at Λ = 5 ◦ and 10 ◦ for sAR = 4 and sAR = 2respectively. However, the increase is more pronounced in the sAR = 2 case. inear modal instabilities around post-stall swept finite-aspect ratio wings
5. Summary
Linear modal three-dimensional (TriGlobal) instability analysis of laminar three-dimensional separated flows over finite aspect ratio, constant-chord wings have beenperformed at Re = 400, two aspect ratios and a range of angles of attack and sweep.Monitoring the unsteady, three-dimensional base flows, the following observations weremade, as the angle of sweep increased. When 0 ◦ (cid:54) Λ < ◦ the three distinct regionsreported by Zhang et al. (2020 a ) were also observed, namely the tip vortex, wakeand the interaction regions with braid-like vortices. For 15 ◦ (cid:54) Λ < ◦ the braid-likevortices of the interaction region become dominant and absorb the tip vortex. Finally,at 25 ◦ (cid:54) Λ (cid:54) ◦ hair-pin vortices form in the wake, tip stall and ”ram’s horn” vorticesare present on the wing. The overall effect of an increasing sweep is flow stabilisation.Linear TriGlobal instability analysis, performed at conditions where steady or station-ary unstable base flows could be computed, revealed the leading eigenmodes of this classof flows for the first time. This analysis gave insight to the formation of the unstablewake, which is found to be caused by a unstable global wake mode at half-span for anunswept wing. Changes in the aspect ratio were seen to have little effect on the shape ofthe leading global flow eigenmode. As sweep increases from Λ = 0 ◦ to 10 ◦ the structuresof the unstable mode peak closer to the wing tip and at Λ = 15 ◦ the mode evolves intoa vortical instability. At the maximum considered sweep of Λ = 30 ◦ the leading mode isstable and takes the form of a vortical instability growing further away from the wing.Data-driven analysis employing POD, and to a lesser extent DMD, was subsequentlyconducted during the nonlinear saturation stage, with the objective of classifying thedominant structures of the separated flow and assessing the effects of wing geometry.On the larger aspect ratio wing, sAR = 4, at low angles of attack the dominant wakemode takes the form of Kelvin-Helmholtz instability. For higher α , this mode exhibitsmore complex three-dimensional features, changing into a structure denominated theinteraction mode. The interaction mode consists of both periodic structures, locatedat half-span of the wing associated with the wake, as well as vortical structures, thatcorrespond to the location of braid-like vortices in the wing wake, recently reportedby Zhang et al. (2020 a ). The interaction mode is also found to be associated withunsteady motion of the reattachment line, observed near the trailing edge in directnumerical simulation. The spatial characteristics of the interaction mode change withsweep, reflecting the changes of the underlying base flow. The mode eventually splitsinto separate structures on either side of the symmetric wing at Λ ≈ ◦ . At even higherangles of sweep the dominant structure unravelled in the data-driven analysis takes theform of elongated streamwise vortices on the larger sAR = 4 wing, which are associatedwith the recently discovered multiple streamwise-aligned vortices that form on the finiteaspect ratio wing at higher sweep angles (Zhang et al. b ). On the shorter, sAR = 2,wing the tip effects lead to the dominant mode becoming indistinguishable from a tip-vortex instability.The wake structures observed in Zhang et al. (2020 a ) have been analysed and classifiedinto two classes of modes, the interaction mode (IM) and the wake wode (WM), havingdistinct energy content at a range of angles of attack and sweep. The frequency contentand the spatial structures of the interaction mode could be linked to the oscillatory motionof reattachment line observed in large-scale simulations. The distinct frequency peaks ofthe lift coefficient signal observed by Zhang et al. (2020 a ) have also been identified inthe frequency spectra of the interaction and wake modes, leaving little doubt that theIM and WM are the major contributors to lift on finite aspect ratio wings. Finally, thestabilisation effects due to sweep and tip vortex discussed by Zhang et al. (2020 a , b )4 A. Burtsev et al. were quantified over a range of moderate sweep angles on the basis of the relativeimportance between the spanwise velocity component and the velocity induced by the tipvortex over at fixed values of the wing aspect ratio and angle of attack. The balance ofthese stabilisation mechanisms at Λ = 5 ◦ was found to lead to a quasi-two-dimensionalstructure of the wake close to the wing.Overall, the effect of sweep on the structure of the leading global flow eigenmodes, aswell as of the modes computed in the data-driven analysis is relatively mild, up to andincluding Λ = 10 ◦ , and mainly affects the spanwise location of the peak of the modestructures. On the other hand, higher sweep angle results, up to Λ = 30 ◦ examinedherein, show that flow over the sAR = 4 wing is dominated by the strength of thespanwise velocity on the suction side of the wing, which prevails over that of the tip-induced spanwise velocity. Conversely, the flow over the shorter sAR = 2 wing at thesame conditions is dominated by tip effects, with the tip-induced velocity countering thespanwise flow.Future work may address the effect of the last remaining unexplored parameter, namelywing taper, on the present findings. However, the low-Reynolds number results reportedhere establish a basis for understanding flow dynamics and instabilities on finite three-dimensional constant-chord wings at low Reynolds numbers, as a first step towardsunderstanding turbulent flow at higher Reynolds numbers and providing theoretically-founded, physics-driven flow control strategies. Acknowledgements.
Support of AFOSR Grant FA9550-17-1-0222 with Dr. Gregg Abate and Dr. DouglasSmith as Program Officers is gratefully acknowledged. The authors also acknowledgecomputational time made available on the UK supercomputing facility ARCHER viathe UKTC Grant EP/R029326/1 and on the DoD Copper supercomputer, via projectAFVAW10102F62 with Dr. Nicholas Bisek as Principal Investigator.
Declaration of Interests.
The authors report no conflict of interest.
Appendix A. Dynamic mode decompostion
Dynamic mode decomposition (DMD) can identify temporal and spatial coherentstructures purely from data. It is based on the best-fit linear operator and allows for thetime-resolved flow fileds to be decomposed into modes, each having a single characteristicfrequency and growth rate (Schmid 2010). The method is related to the Koopman (1931)operator, as discussed by Rowley et al. (2009) and was introduced into fluid dynamics bySchmid & Sesterhenn (2008) and and Schmid (2010). Several formulations and extensionsof the method have been created over the past years (see Kutz et al. X = [ χ ( t ) χ ( t ) ... χ ( t m )] ∈ R n × m , (A 1) X ∗ = [ χ ( t ) χ ( t ) ... χ ( t m +1 )] ∈ R n × m . (A 2)A reduced singular value decomposition (SVD) is performed on the data matrix X = U Σ V T where T denotes conjugate transpose. Optionally the SVD can be truncated byonly considering the first r columns of U and V and the first r rows and columns of Σ to obtain U r , Σ r , and V r . The relationship between the snapshots is approximated in a inear modal instabilities around post-stall swept finite-aspect ratio wings X ∗ = AX = AU r Σ r V Tr , (A 3)where A = X ∗ X + with X + denoting the pseudoinverse of X . Multiplying both side by U T and rearranging while making use of the fact that V Tr = V − r gives˜ A = U Tr X ∗ V r Σ − r ∈ R r × r , (A 4)where ˜ A = U Tr AU . The eigenvectors ˜ ν j and eigenvalues µ j can be found with:˜ A ˜ ν j = µ j ˜ ν j . (A 5)For a DMD eigenvalue (every nonzero µ j ) the corresponding DMD mode ν j is found by ν j = µ j − X ∗ V r Σ r − ˜ ν j . (A 6)Unlike with POD where the ordering of the modes already represents the energy con-tribution, it is more difficult to identify the physically relevant modes with DMD. Anumber of approaches to address this exist. In general, these include ranking using thenorm of each mode (Rowley et al. et al. et al. et al. et al. et al. b j = W + χ ( t ) , (A 7)where W is a matrix whose columns contain the DMD eigenvectors ν j , χ ( t ) is the firstsnapshot and + denotes the Moore-Penrose inverse. This approach requires the inverseof W which is a complex n × m matrix. Singular value decomposition is used for this W + = V Σ U T , (A 8)where Σ is a diagonal matrix containing the singular values and U and V are matricescontaining left and right singular vectors of W , respectively. Finally, the oscillationfrequency f j and the growth rate g j of a DMD mode can be determined from the DMDeigenvalue λ j using f j = ∠ λ j / (2 πδt ) , (A 9) g j = log | λ j | /δt, (A 10)where δt is the uniform sampling increment between snapshots. Appendix B. data-driven analysis validation
B.1.
Code validation
In order to verify that the modal results are independent of the number of snapshots m and the interval between them δt , a sensitivity analysis was performed on a singlecase ( sAR, Λ, α ) = (2 , ◦ , ◦ ). Three datasets were considered: 300 snapshots in theinterval of 40 (cid:54) t (cid:54)
70 with δt = 0 . δt = 0 .
2, and 150 snapshots inthe interval of 40 (cid:54) t (cid:54)
55 with δt = 0 . A. Burtsev et al. -4 -3 -2 -1 Figure 26: Sensitivity of POD to number of snapshots m and time interval δt . T1 S1BiGlobal EVP (Nektar++) ± . . i − . i POD, DMD (linear) ± . . i − . i POD, DMD (nonlinear) 0 . i − . i Table 7: Leading travelling and stationary modes for a reference case of two-dimensionalcylinder flow at Re D = 60.shows that the POD eigenvalues generated using these datasets are identical for the first15 modes, which collectively represent 99.995% of flow kinetic energy. It is thereforeconcluded that m = 300 δ = 0 . a ) that increasing δt to 0.2 (+ markers) leads to DMD loosing every other ofthe high initial amplitude modes. Note that the difference in the initial amplitudes ( b j )themselves is not a cause for concern since only the relative difference in b j is used todifferentiate the dominant modes. Reducing the time span of the data by half ( × markersdenoting m = 150, δt = 0 .
1) recovers dominant modes well but loses some accuracy atlower frequencies of
St < .
5. The two nonzero frequency modes at b j ≈ × − (markedwith arrows) that are not recovered by the shorter dataset are not physically significantas they lie far from the unit circle as shown in 27( b ). Finally, the effect of computinginitial amplitudes based on recovering less than m DMD modes is considered. Once theDMD eigenvalues are found, corresponding modes are computed using equation A 6,which is a computationally expensive step. Red square markers in 27( a ) show that b j computed based on 80 modes is accurate for the first 4 dominant modes. Note that thisdoes not affect the DMD eigenvalues in 27( b ). Therefore, the m = 300 δ = 0 . Relationship between global stability analysis and data-driven methods
In order to illustrate the relationship between operator based global analysis and data-driven POD and DMD consider a reference problem of a two-dimensional flow over acylinder at Re D = 60. The leading travelling and stationary modes obtained by BiGlobal inear modal instabilities around post-stall swept finite-aspect ratio wings ( a ) ( b ) Figure 27: Sensitivity of DMD to number of snapshots m and the interval between them δt ( a ) frequencies and normalised initial amplitudes showing the effect of recovering lessthan m modes, for clarity modes with b j < − are omitted ( b ) DMD eigenvalues onthe unit circle.stability analysis are compared to DMD performed on snapshots recorded in the lineargrowth regime and in the nonlinear saturated regime as shown in table 7. data-drivenanalysis shows good agreement with global stability in the linear regime but recoversvery different travelling mode when performed on the data taken in saturated regime.It is worth pointing out that POD gives identical results if the frequency and the decayrate of the temporal coefficients. REFERENCES˚Akervik, E., Brandt, L., Henningson, D. S., Hœpffner, J., Marxen, O. & Schlatter,P.
Physics of Fluids , 068102. Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz,J., Greenbaum, A., Hammarling, S., McKenney, A. & Sorensen, D.
LAPACKUsers’ Guide , 3rd edn. Philadelphia, PA: Society for Industrial and Applied Mathematics.
Aubry, N., Guyonnet, R. & Lima, R.
Journal of Statistical Physics (3-4), 683–739. Barkley, D. & Henderson, R. D.
Journal of Fluid Mechanics , 215–241.
Berkooz, G., Holmes, P. & Lumley, J. L.
Annual Review of Fluid Mechanics (1), 539–575. Bippes, H. & Turk, M.
Tech. Rep.
DFVLR-IB 251-80. DFVLR.
Bippes, H. & Turk, M.
Tech. Rep.
DFVLR-FB 84-20. DFVLR.
Black, J.
The Journal of the Royal Aeronautical Society (541), 51–60. Boiko, A. V., Dovgal, A. V., Zanin, B. Yu & Kozlov, V. V.
Thermophysics and Aeromechanics (1),1–13. Broeren, A. P. & Bragg, M. B.
AIAA Journal (9), 1641–1651. A. Burtsev et al.
Calderon, D. E., Cleaver, D. J., Gursul, I. & Wang, Z.
Physics of Fluids (7), 071907. Cantwell, C. D., Moxey, D., Comerford, A., Bolis, A., Rocco, G., Mengaldo, G.,Grazia, D. De, Yakovlev, S., Lombard, J. E., Ekelschot, D., Jordi, B., Xu, H.,Mohamied, Y., Eskilsson, C., Nelson, B., Vos, P., Biotto, C., Kirby, R. M.& Sherwin, S. J.
Computer Physics Communication , 205–219.
Citro, V., Luchini, P., Giannetti, F. & Auteri, F.
Journalof Computational Physics , 234 – 246.
Cosyn, P. & Vierendeels, J.
Journal of Aircraft (3), 713–722. Dallmann, U.
FluidDynamics Research (1), 183 – 189. D´elery, J.
Three-dimensional separated flow topology: critical points, separation lines andvortical structures . John Wiley & Sons.
Dong, H., Mittal, R. & Najjar, F. M.
Journal of Fluid Mechanics , 309–343.
Dong, S., Karniadakis, G. E. & Chryssostomidis, C.
Journal of Computational Physics , 83–105.
Edstrand, A. M., Davis, T. B., Schmid, P. J., Taira, K. & Cattafesta, L. N.
Journal of Fluid Mechanics , R1.
Edstrand, A. M., Schmid, P. J., Taira, K. & Cattafesta, L. N. a A parallel stabilityanalysis of a trailing vortex wake.
Journal of Fluid Mechanics , 858–895.
Edstrand, A. M., Sun, Y., Schmid, P. J., Taira, K. & Cattafesta, L. N. b Activeattenuation of a trailing vortex inspired by a parabolized stability analysis.
Journal ofFluid Mechanics , R2.
Eldredge, J. D. & Jones, A. R.
AnnualReview of Fluid Mechanics (1), 75–104. Elimelech, Y., Arieli, R. & Iosilevskii, G.
Physics ofFluids (2), 024104. Fischer, P. F., Lottes, J. W. & Kerkemeier, S. G.
Geuzaine, C. & Remacle, J. F.
International Journal for NumericalMethods in Engineering (11), 1309–1331. G´omez, F., Clainche, S. Le, Paredes, P., Hermanns, M. & Theofilis, V.
AIAA Journal (12), 2731–2743. Gursul, I., Cleaver, D. J. & Wang, Z.
Progress in Aerospace Sciences , 17–55. Gursul, I., Wang, Z. & Vardaki, E.
Progress in Aerospace Sciences (7), 246–270. Harper, C. W. & Maki, R. L.
Tech.Rep.
D-2373. NASA.
Hayostek, S., Zhang, K., Taira, K., Burtsev, A., He, W., Theofilis, V. & M., Amitay
Journal ofFluid Mechanics in review . He, W., Burtsev, A., Theofilis, V., Zhang, K., Taira, K., Hayostek, S. & Amitay,M. a Wake dynamics of finite aspect ratio wings. part III: Triglobal linear stabilityanalysis.
AIAA Paper 2019-1386 . He, W., Gioria, R., P´erez, J. M. & Theofilis, V. a Linear instability of low Reynoldsnumber massively separated flow around three NACA airfoils.
Journal of Fluid Mechanics , 701–741. inear modal instabilities around post-stall swept finite-aspect ratio wings He, W., Guan, Y., Theofilis, V. & Li, L. K. B. b Stability of low-Reynolds-numberseparated flow around an airfoil near a wavy ground.
AIAA Journal (1), 29–34. He, W., P´erez, J. M., Yu, P. & Li, L. K. B. c Non-modal stability analysis of low-re separated flow around a NACA 4415 airfoil in ground effect.
Aerospace Science andTechnology , 269 – 279. He, W., Tendero, J. A., Paredes, P. & Theofilis, V. b Linear instability in the wakeof an elliptic wing.
Theoretical and Computational Fluid Dynamics , 483–504. Holmes, P., Lumley, J. L. & Berkooz, G.
Turbulence, Coherent Structures, DynamicalSystems and Symmetry . Cambridge Monographs on Mechanics . Cambridge UniversityPress.
Hornung, H. & Perry, A. E.
Zeitschrift f¨ur Flugwissenschaft und Weltraumforschung ,77–87. Huang, Y., Venning, J., Thompson, M. C. & Sheridan, J.
Journal of Fluid Mechanics ,341–369.
Jantzen, R. T., Taira, K., Granlund, K. O. & Ol, M. V.
Physics of Fluids (5), 053606. Jones, A. R., Medina, A., Spooner, H. & Mulleners, K.
Experiments in Fluids (4), 52. Jordi, B. E., Cotter, C. J. & Sherwin, S. J.
Physics of Fluids (3), 034101. Juniper, M. P., Hanifi, A. & Theofils, V.
Appl. Mech. Rev. (2), 024804–1 – 024804–22. Keller, H. B.
Application of Bifurcation Theory (ed. P. Rabinowitz), pp. 359–384. New York: Academic.
Kim, D. & Gharib, M
Experiments in Fluids (1), 329–339. Kitsios, V., Rodr´ıguez, D., Theofilis, V., Ooi, A. & Soria, J.
Journal ofComputational Physics (19), 7181 – 7196.
Koopman, B. O.
Proceedingsof the National Academy of Sciences (5), 315–318. Kutz, J. N., Brunton, S. L., Brunton, B.W. & Proctor, J. L.
Dynamic ModeDecomposition: Data-Driven Modeling of Complex Systems . SIAM, Philadelphia, PA.
Lumley, J. L.
Atmospheric turbulenceand wave propagation pp. 66–178.
Mancini, P., Manar, F., Granlund, K. O., Ol, M. V. & Jones, A. R.
Physicsof Fluids (12), 123102. Manolesos, M. & Voutsinas, S. G.
Physics of Fluids (4), 045101. Medina, A., Eldredge, J. D., Kweon, J. & Choi, H.
AIAA Journal (9), 2607–2620. Moss, G. F. & Murdin, P. M.
Tech.Rep.
TR68104. Royal Aeronautical Establishment.
Navrose, Brion, V. & Jacquin, L.
Journal of Fluid Mechanics , 399–430.
Paredes, P., Gosse, R., Theofilis, V. & Kimmel, R.
Journal of Fluid Mechanics , 442–466.
Perry, A. E. & Chong, M. S.
Annual Review of Fluid Mechanics (1), 125–155. Plante, F., Dandois, J., Beneddine, S., Sipp, D. & Laurendeau, ´E.
AAAFAERO2019 . PARIS, France. A. Burtsev et al.
Rodr´ıguez, D. & Theofilis, V.
Journal of Fluid Mechanics , 280–305.
Rodr´ıguez, D. & Theofilis, V.
Theoretical andComputational Fluid Dynamics (1-4), 105–117. Rossi, E., Colagrossi, A., Oger, G. & Le Touz´e, D.
Journal of Fluid Mechanics , 356–391.
Rowley, C. W., Mezi´c, I., Bagheri, S., Schlatter, P. & Henningson, D. S.
Journal of Fluid Mechanics , 115–127.
Sayadi, T., Schmid, P. J., Nichols, J. W. & Moin, P.
Journal of Fluid Mechanics , 278–301.
Schewe, G.
Journalof Wind Engineering and Industrial Aerodynamics (14), 1267 – 1289, bluff BodyAerodynamics and Applications. Schmid, P. J.
Journalof Fluid Mechanics , 5–28.
Schmid, P. J. & Sesterhenn, J. . Schmid, P. J., Violato, D & Scarano, F.
Experiments in Fluids , 1567–1579. Sirovich, L.
Quarterlyof Applied Mathematics (3), 561–571. Smith, L. R. & Jones, A. R.
Journal of Fluid Mechanics , A22.
Son, O. & Cetiner, O.
Journal of Aerospace Engineering (5), 04017053. Taira, K., Brunton, S. L., Dawson, S. T. M., Rowley, C. W., Colonius, T., McKeon,B. J., Schmidt, Oliver T., Gordeyev, S., Theofilis, V. & Ukeiley, L. S.
AIAA Journal (12), 4013–4041. Taira, K. & Colonius, T.
Journal of Fluid Mechanics , 187–207.
Taira, K., Hemati, M. S., Brunton, S. L., Sun, Y., Duraisamy, K., Bagheri, S., M.,Dawson S. T. & Yeh, C. A.
AIAA Journal (3), 998–1022. Teixeira, R. de S. & Alves, L. S. de B.
Theoretical and Computational FluidDynamics (5-6), 607–621. Theofilis, V.
Tech. Rep.
F61775-99-WE090.
Theofilis, V.
Progress in aerospace sciences (4), 249–315. Theofilis, V., Barkley, D. & Sherwin, S.
The Aeronautical Journal (1968) (1065), 619–625.
Theofilis, V. & Colonius, T.
AIAA Paper 2003-4143 . Theofilis, V., Hein, S. & Dallmann, U.
Philosophical Transactions of the RoyalSociety of London. Series A: Mathematical, Physical and Engineering Sciences (1777),3229–3246.
Torres, G. E. & Mueller, T. J.
AIAA Journal (5), 865–873. Tu, J. H., Rowley, C. W., Luchtenburg, D. M., Brunton, S. L. & Kutz, J. N.
Journal of ComputationalDynamics (2), 391–421. Tumuklu, O., Levin, D. A. & Theofilis, V. inear modal instabilities around post-stall swept finite-aspect ratio wings laminar separated flows over a double cone geometry using a kinetic approach. Physics ofFluids (4), 046103. Tumuklu, O., P´erez, J.M., Theofilis, V. & Levin, D.A.
Visbal, M. R. & Garmann, D. J.
AIAA Journal (8), 3274–3289. Wan, Z-H, Zhou, L., Wang, B-F & Sun, D-J
European Journal of Mechanics - B/Fluids , 16 –26. Weihs, D. & Katz, J.
AIAAJournal (12), 1757–1759. Williamson, C. H. K
Journal of Fluid Mechanics , 579–627.
Williamson, C. H. K.
Annual Review of FluidMechanics (1), 477–539. Winkelman, A. E. & Barlow, J. B.
AIAA Journal (8), 1006–1008. Wygnanski, I., Tewes, P., Kurz, H., Taubert, L. & Chen, C.
J.Fluid Mech. , 336–346.
Wygnanski, I., Tewes, P. & Taubert, L.
Journal of Aircraft (1). Yen, S. C. & Hsu, C. M.
AIAAJournal (1), 228–236. Yen, S. C. & Huang, L. C.
Journal of Fluids Engineering (11), 111101.
Zhang, K., Hayostek, S., Amitay, M., He, W., Theofilis, V. & Taira, K. a On theformation of three-dimensional separated flows over wings under tip effects.
Journal ofFluid Mechanics , A9.
Zhang, K., Hayostek, S., M., Amitay, Burtsev, A., Theofilis, V. & Taira, K. b Laminar separated flows over finite-aspect-ratio swept wings.
Journal of Fluid Mechanics , R1.
Zhang, W. & Samtaney, R.
Physics of Fluids28