key: cord-0504566-bj05babp
authors: Zhang, Kai; Hayostek, Shelby; Amitay, Michael; He, Wei; Theofilis, Vassilios; Taira, Kunihiko
title: On the formation of three-dimensional flows over finite-aspect-ratio wings under tip effects
date: 2019-09-10
journal: nan
DOI: nan
sha: 5b6401c15a7c1055f6d5566ca3aefc1a4c2aee53
doc_id: 504566
cord_uid: bj05babp

We perform DNS of flow over finite-aspect-ratio NACA 0015 wings to characterize the tip effects on the wake dynamics. This study focuses on the development of three-dimensional separated flow over the wing, and discuss flow structures formed on the wing surface as well as in the far field wake. Vorticity is introduced into the flow in a predominantly two-dimensional manner. The vortex sheet from the wing tip rolls up around the free end to form the tip vortex. At its inception, the tip vortex is weak and its effect is spatially confined. As the flow around the tip separates, the tip effects extend farther in the spanwise direction, generating noticeable three dimensionality in the wake. For low-aspect-ratio wings, the wake remains stable due to the strong downwash over the entire span. Increasing the aspect ratio allows unsteady vortical flow to emerge away from the tip at sufficiently high angles of attack. These unsteady vortices shed and form closed vortical loops. For higher-aspect-ratio wings, the tip effects retard the near-tip shedding process, which desynchronizes from the two-dimensional shedding over the mid-span region, yielding vortex dislocation. At high angles of attack, the tip vortex exhibits noticeable undulations due to the strong interaction with the unsteady shedding vortices. The spanwise distribution of force coefficients is related to the three-dimensional wake dynamics and the tip effects. Vortical elements in the wake that are responsible for the generation of lift and drag forces are identified through the force element analysis. We note that at high angles of attack, a stationary vortical structure forms at the leading edge near the tip, giving rise to locally high lift and drag forces. The analysis performed here reveals how the vortical flow around the tip influences the separation physics, the global wake dynamics, and the spanwise force distributions.

Separated flows over lifting surfaces have been studied extensively due to their critical importance in aerodynamics and hydrodynamics. Substantial researches have been dedicated to the understanding of separated flows over canonical airfoils to reveal the † Email address for correspondence: kzhang3@ucla.edu arXiv:1909.04729v1 [physics.flu-dyn] 10 Sep 2019 separation mechanism and their influences on the aerodynamic characteristics (Anderson 2010) . In many of the fundamental studies, the assumption of two-dimensional (or quasitwo-dimensional) flow was applied to assess the performance of airfoils, offering valuable insights (Abbott & Von Doenhoff 1959; Pauley et al. 1990; Huang et al. 2001; Yarusevych et al. 2009; He et al. 2017a; Rossi et al. 2018) .

Over the past couple of decades, there have been a large number of studies focusing on three-dimensional post-stall flows at angles of attack much higher than what were traditionally considered. These studies were performed to examine fundamental aspects of flows related to biological fliers and swimmers as well as small-scale aircraft that experience large-amplitude perturbations in flight (Ellington et al. 1996; Dickinson et al. 1999; Eldredge & Jones 2019) . The lifting surfaces for these flows have some uniqueness. The wings have low-aspect-ratio planforms and generate prominent tip vortices at high angles of attack (Taira & Colonius 2009; Lentink & Dickinson 2009 ). For these reasons, the resulting wakes are highly three-dimensional with complex nonlinear dynamics generated by the interactions among the wake vortices. Under the tip effects, the wake dynamics of finite wings can be significantly different from its two-dimensional analogue. The three-dimensional nature of the wakes requires global analysis of the separation dynamics.

Separated flows over finite-aspect-ratio wings in pure translation have been examined with surface oil visualizations (Gregory et al. 1971; Wang 1976; . From these visualizations, proposed a flow field model for a post-stall rectangular wing. This model consists of a pair of counter-rotating tip vortices at the two ends of the wing and the mushroom-like three-dimensional separation bubble. Later, by smoke visualization, Freymuth et al. (1987) revealed the intricate threedimensional wake structures behind finite wings. On the aerodynamic characterizations of low-aspect-ratio wings at low Reynolds numbers, Pelletier & Mueller (2000) , Torres & Mueller (2004) and Ananda et al. (2015) obtained a wealth of force and pitching moment data for a large collection of aspect ratios, angles of attack, and planform geometries.

Direct numerical simulations have provided detailed insights into separated flows over finite-aspect-ratio wings in translation. The wake structures behind an impulsively started plate were analyzed by Taira & Colonius (2009) using direct numerical simulation at Reynolds numbers of 300 and 500. In the initial stage of wake development, the wake vortices share the same structures for all aspect ratios. At large time, the tip vortices significantly influence the vortex dynamics and the corresponding forces on the wings. The long-term stability of the wake was found to be dependent on the angle of attack, aspect ratio and the Reynolds number. Lee et al. (2012) applied the force element theory (Chang 1992) to the impulsively started flows around low-aspect-ratio wings and identified the wake regions responsible for the force generation. They have shown that the tip vortices exert lift force during the start-up maneuver. More recently, Garmann & Visbal (2017) conducted high-fidelity large eddy simulations of flow over finite-aspectratio wings at Re = 2 × 10 5 . Intricate details of the vortical structures near the wing tip were revealed.

Stability analysis have also been carried out to study the instabilities in the wake of finite-aspect-ratio wings. He et al. (2017b) conducted linear stability analysis on the flows over finite elliptical wings. Two classes of linearly unstable perturbations were identified, namely the highly-amplified symmetric modes and the weakly-amplified antisymmetric disturbances. Navrose et al. (2019) performed transient growth analysis on the wake of finite NACA 0012 wings. The optimal linear perturbation was found to be located near the surface of the wing in the form of chord-wise periodic structures whose strength decreases from the root towards the tip. Edstrand et al. (2018a) carried out parabolized stability analysis of the tip vortex generated from a finite NACA 0012 wing at the Reynolds number of 1000. The insights obtained from their analysis have enabled effective design of active control techniques for the attenuation of the tip vortex (Edstrand et al. 2018b) .

We note that three-dimensional separated flows over finite-aspect-ratio wings have also been extensively studied for wings under unsteady maneuvers, including surging (Chen et al. 2010; Mancini et al. 2015) , rotation (Lentink & Dickinson 2009; Harbig et al. 2013; Jones et al. 2016) , pitching (Buchholz & Smits 2006; Green & Smits 2008; Yilmaz & Rockwell 2012; Jantzen et al. 2014) , heaving (Visbal et al. 2013) , plunging (Akkala & Buchholz 2017) , and flapping (Birch & Dickinson 2001; Dong et al. 2006; Shyy et al. 2010 ). These works have closely examined the formation of large-scale vortical structures, with particular emphasis on the leading-edge vortex and its influence on added lift (Eldredge & Jones 2019) .

Despite the vast literature, there are still fundamental questions left unanswered concerning three-dimensional separated flows over finite-aspect-ratio wings in pure translation. It is generally known that the three-dimensional wake dynamics described in the previous works are closely related with the tip effects. However, the detailed process of how the tip-imposed three dimensionality is introduced into the flow remains elusive. The formation of the wake vortical structures under the tip effects, as well as the dynamical relations between those structures requires further analysis.

To address the aforementioned questions, we carry out a large number of direct numerical simulations of flows over finite-aspect-ratio NACA 0015 wings. By focusing on the development of the three-dimensional wake, we characterize the tip effects on separation and wake dynamics. In what follows, we present the computational set-up and its validation in §2. The results are discussed in §3 starting with the flow physics on the wing surface in §3.1. A detailed look at the three-dimensional wake structures is provided in §3.2, followed by the spectral analysis of the wake in §3.3. A classification of the wake states for cases with different angles of attack and aspect ratios is provided §3.4

The aerodynamic forces and their relationship with the wake dynamics are presented in §3.5. We conclude this study by summarizing our findings in §4.

We numerically analyze incompressible flow over a NACA 0015 wing with a straight cut tip. The wing is situated in uniform flow with velocity U ∞ as shown in figure 1. The three-dimensional flow over the finite wing is studied by numerically solving the Navier-Stokes equations

where u = (u x , u y , u z ) is the velocity vector and p is the pressure. An incompressible flow solver Cliff (in CharLES software package, Cascade Technologies, Inc.) based on the finite-volume formulation with second-order accuracy in time and space (Ham & Iaccarino 2004; Ham et al. 2006) is used in the present study. For non-dimensionalization, we normalize the spatial variables by the chord length c, velocities by the freestream velocity U ∞ , and time by c/U ∞ . The Reynolds number in the momentum equation (2.1) is defined as Re ≡ U ∞ c/ν, where ν is the kinematic viscosity. In this work, we set the Reynolds number to Re = 400 so that the flow remains laminar. As we are interested in the three-dimensional flow physics imposed by the tip effects, we simulate the half-wing model by prescribing the symmetry boundary condition along the mid-span. Denoting the length of the half wing as b, the (half-wing) aspect ratio is defined as AR = b/c and is varied from 1 to 6. The angle of attack α is varied from from 0 • to 30 • .

The computational domain covers (x, y, z) ∈ [−13, 25] × [−15, 15] × [0, 16], where x, y and z are the streamwise, crossflow and spanwise coordinates. This results in a maximum blockage ratio of 0.8% for the case of (α, AR) = (30 • , 6). The origin of the Cartesian coordinate system is placed at the leading edge of the wing on the mid-span plane (z = 0) plane. The computational domain is discretized with a C-type grid with the mesh refined in the vicinity of the wing as well as its wake. The adequacy of the grid resolution is examined in §2.2.

The inlet and far-field boundaries are prescribed with freestream u = (U ∞ , 0, 0). The convective boundary condition (∂u/∂t + U ∞ ∂u/∂x = 0) is applied at the outlet to allow wake structures to leave the domain without disturbing the near-field flow. The wing surface is treated with the no-slip boundary condition. The simulations are started with uniform flow, and statistics are recorded only after the initial transients are flushed out of the computational domain (typically 25 convective time units). Flow statistics are collected with over 100 convective time units to ensure statistical convergence.

To verify the convergence of the numerical results with respect to the grid resolution, we examine the case of (α, AR) = (22 • , 6) with two sets of meshes described in table 1. For the refined mesh, the temporal resolution is also increased by restricting the maximum CFL number to be half of the medium mesh. We compare the forces on the wing evaluated from the two meshes to assess the mesh requirement. The time history of

Np Ns Nz NCV CFL medium 60 80 120 14.8 × 10 6 1.0 refined 80 100 150 20.6 × 10 6 0.5 Table 1 . Setups for the different meshes for the case of (α, AR) = (22 • , 6). Np and Ns are the numbers of grid points on the pressure and suction side of the airfoil. Nz is the number of grids along the wing span. NCV is the total number of control volumes, and CFL is the Courant number used in the simulations. the lift force obtained from these meshes are reported in figure 2. Throughout this study, the lift and drag forces are reported in their non-dimensional forms through

where S denotes the wing surface, n is the unit normal vector pointing outward from the wing surface. Unit vectors e y and e x are in the lift and drag directions. Moreover, τ = µω × n is the skin-friction vector (µ is the dynamic viscosity of fluid and ω = ∇ × u is the vorticity vector). Once the initial transient settles (from the use of uniform flow as the initial condition), the lift coefficients exhibit quasi-periodic oscillations with lowfrequency beating. This is resolved well with both mesh resolutions even at large times. Thus, the medium resolution grid is deemed adequate to yield accurate results and is used throughout this study. To validate our computational setup, we compare the numerical results against the stereoscopic particle image velocimetry (SPIV) measurements from water tunnel experiments conducted at the Rensselaer Polytechnic Institute. Detailed descriptions of the experimental setup are documented in Hayostek & Amitay (2018) . In contrast to the symmetry boundary condition used for the mid-span in the present simulations, a wall is present at the root plane in the water tunnel. For this reason, comparison between the numerical and experimental results is focused near the tip region for the case of (α, AR) = (22 • , 6), where the influence of the root boundary condition is small. The timeaveraged streamwise velocity fields are shown in figure 3 for three streamwise locations at x = 1, 2 and 3. While the water tunnel experiments are conducted at a larger Reynolds number of 600, the flow remains laminar and is in good agreement with that from the present numerical simulation, confirming the fidelity of the numerical results.

We further validate the present numerical results against those from another flow solver Nektar++ that is based on the spectral/hp element method (Cantwell et al. 2015) . The computation is carried out with an independent numerical setup (He et al. 2019) . We compare the wake vortical structures obtained by Nektar++ (top) and Cliff (bottom) for the same case of (α, AR) = (22 • , 4) at t = 9.25 in figure 4. Shown are the iso-surfaces of ω = 1 (vorticity magnitude) and Q = 1 (second invariant of velocity gradient tensor). The vortical structures obtained from the two solvers agree in an excellent manner with each other even for the smallest flow details. Based on the above verification and validation, we have ensured that the current computational setup is reliable. We now proceed to present our results in the next section. 

The wake behind a finite-aspect-ratio wing exhibits rich three-dimensional flow features due to the tip effects, as shown in figure 5. In this section, we begin the discussions by focusing on the wing surface to look for the source of vorticity. Next, we describe in detail the three-dimensional wake dynamics. The aerodynamic forces are also discussed with an emphasis on their relationship with the three-dimensional wake structures.

We begin our analysis by examining the introduction of vorticity to the flow from the wing surface. In incompressible flow, vorticity is only generated on the solid surface and diffuses into the flow (Morton 1984; Hornung 1989; Wu et al. 2007) . The rate at which vorticity is created at the surface is given by the wall-normal boundary vorticity flux

where Σ = (∇ω) 0 is the vorticity gradient tensor at the surface. Let us show the vorticity flux along the wing surface for a representative case of (α, AR) = (22 • , 4). The time-averaged wall-normal boundary vorticity flux components σ x , σ y and σ z are plotted in figure 6(a) to identify the local generation of vorticity. As expected, the spanwise vorticity flux σ z dominates over the surface of the wing. Negative σ z is contained within a narrow region along the leading edge of the wing. This region injects a large amount of negative ω z into the flow, forming the leading-edge vortex sheet. Compared to σ z , the other two components of the vorticity flux appear negligible for the majority of the span. It is only at the close vicinity of the wing tip that considerable level of σ x and σ y can be found. The creation of vorticity on the finite-aspect-ratio wing is not strongly affected by the tip effects and is predominantly two-dimensional along the span.

Once the vorticity is generated at the wing surface, it is diffused into the flow and is convected by the freestream. The developed flow imprints the skin-friction field τ on the wall, as shown in figure 6(b). Three regions can be recognized from the skinfriction field. The region close to the leading edge away from the tip is occupied by the nearly two-dimensional boundary layer, which separates from the wing in a uniform manner. Post separation, a region of reversed flow appears on the suction side of the wing where the skin-friction lines exhibit spanwise component of the flow. The dominant three dimensionality is observed over the wing-tip region with a three-dimensional swirl pattern, as reported in the experimental work of .

To assess the influence of the angle of attack on the near-surface flow pattern, we visualize the time-averaged skin-friction lines for the cases of AR = 4 with different α in figure 7. For higher angles of attack, the flow separation occurs more upstream, and the spanwise extent of the boundary layer that is affected by the tip effects becomes smaller. The skin-friction lines in the recirculation region becomes more three-dimensional at higher α.

The development of the three-dimensional flow can be envisioned from the boundary vorticity flux and skin-friction lines. Vorticity is introduced from the wing surface in a predominantly two-dimensional manner. The vortex sheets emanate from the leading edge and the wing tip. While the roll up of the leading-edge vortex sheet is twodimensional, the roll-up from the tip is fundamentally three-dimensional yielding a streamwise tip vortex. As the tip vortex develops, its influence grows from being locallyconfined to imposing spanwise variations, as evident in figure 5. We discuss the threedimensional wake dynamics in the following sections.

3.2. Three-dimensional wake structures

Let us continue to examine the three-dimensional wake structures for the representative case of (α, AR) = (22 • , 4). A top view of the wake vortical structures is shown in figure  8 . As also seen in figure 5 , the wake can be divided into the tip vortex region at the free end and the unsteady vortex shedding region near the mid-span. The tip vortex appears practically steady over time in the current case. The unsteadiness in the tip vortex becomes noticeable for cases with smaller AR and larger α, as discussed later in §3.2.2 and §3.2.3. The unsteady region of the wake can be further divided into two subregions. Adjacent to the mid-span (z 1.5), the wake is dominated by vortex cores that are aligned predominantly with the spanwise direction. These vortices are formed by the roll-up of the vortex sheets from the leading and trailing edges, and shed alternatively into the wake. On the other hand, the shedding vortices at 2 z 2.5 feature braidlike structures comprised of both the streamwise vorticity ω x and crossflow vorticity ω y . These vortical structures shed at the same frequency with the spanwise vortices.

The complex vortical structures in the unsteady region can be better understood by visualizing the vorticity lines, as depicted by the black solid lines in figure 8. Each vortex core in the braid-like region connects with a pair of counter-rotating spanwise vortical structures. These structures together form a closed vortex loop (considering the symmetry boundary condition at the mid-span). The unsteady wake region is constituted by a series of such vortex loops arranged in a zig-zag manner. It is thus clear that the braid-like structures play the role of closing the spanwise vortex system. This is in accordance with Helmholtz's second law which states that a vortex filament can not terminate in a fluid (Batchelor 1967) . Similar vortex line models have also been proposed to describe the wakes of finite-span circular cylinders (Taneda 1952; Levold 2012) .

Separated flows behind finite-aspect ratio wings are strongly influenced by the aspect ratio AR and the angle of attack α. Let us show some representative wakes for a range of aspect ratios and angles of attack in figure 9 . The short wings with AR = 1 generates different wake features compared with those discussed in §3.2.1. Stable three-dimensional flow is observed at α = 14 • , while the analogous two-dimensional flow at the same angle of attack exhibits periodic unsteady shedding. This indicates again that the tip effects stabilize the local flow, which in the current case of AR = 1 is the entire span. As the angle of attack is increased to 22 • and 30 • , the wake instability becomes sufficiently strong to overcome the stabilizing effect from the tip, resulting in periodic shedding of the wake vortices. Originated from the leading edge with predominately spanwise vorticity, the shedding vortices on the suction side gradually morph into the hairpin vortices with their legs extending far upstream. The tip vortex cannot be clearly observed due to its strong interaction with the unsteady shedding vortices. Similar vortical features have also been reported in the wake of a flat-plate wing at similar Reynolds numbers by Taira & Colonius (2009) . For AR = 2, the increased spanwise extent alleviates the strong interaction between the tip vortex and unsteady shedding vortices as the case in AR = 1 wing. This enables the tip vortex to remain relatively steady despite the unsteadiness presented in the rest of the flow. The wake at α = 14 • exhibits weak periodic shedding according to a close examination of the lift force, although it may not be clearly visible in figure 9. At α = 22 • , the unsteady shedding vortices mainly take the form of the braid-like structures, without clear spanwise vortex cores. As the angle of attack further increases to α = 30 • , the unsteady vortex shedding grows stronger. Under its effect, the tip vortex exhibits noticeable undulations along the streamwise direction.

Increasing the aspect ratio to AR = 4 reveals the spanwise vortex shedding region. The spanwise vortical structures are quite weak at α = 14 • and appear only in the direct vicinity of the mid-span. With the increase in angle of attack, the unsteady shedding vortices grow stronger and extend farther into the wake. The tip vortex elongates compared with that at lower angles of attack. At α = 30 • , small-scale streamwise vortical structures emerge in the spanwise shedding region, making the wake complex. The streamwise undulations in the tip vortex become less significant compared with that in the case of (α, AR) = (30 • , 2).

As we consider higher aspect ratio of 6, the spanwise shedding takes place over an expanded region for α = 14 • and 22 • . An additional feature that becomes clear at such high aspect ratio is the vortex dislocation, which divides the spanwise shedding region into two parts. This phenomenon is responsible for the low-frequency beating in the lift coefficients as shown in figure 2, and will be examined in detail in §3.3. At α = 30 • , 

We characterize the tip vortex in detail here for the representative case of (α, AR) = (22 • , 4). The time-averaged streamwise vorticity fields in the tip region are shown in figure 10 (a) at x = 2.5, 5 and 7.5. The contours of ω x = −1.0, −1.5 and −2.0 are shown to capture the shape of the vortex core. At x = 2.5, the three contour lines are close to each other, depicting a mango-like shape. At x = 5, the regions enclosed by the contours of ω x = −1.5 and −2.0 shrink significantly. The vorticity level inside the vortex core also decreases. Further downstream at x = 7.5, the tip vortex core becomes more diffused, and its shape relaxes to a circle.

The strength of the tip vortex can be quantitatively assessed through its streamwise circulation Γ x = C u · dl, where the integration contour C is taken as one of the three contours shown in figure 10(a) . Evaluating the circulations along the streamwise direction, we obtain three |Γ x |-x profiles in figure 10(b) . The strengths of the tip vortices, evaluated with the three contour levels, decay linearly in the x direction. The line associated with smaller |ω x | measures extended length and reduced decay rate.

To characterize the strengths of the tip vortices, let us use the contour of ω x = −1.5 and calculate the length L and the decay rate d|Γ x |/dx for different cases, as summarized in figures 10(c) and (d ). We show in Appendix A that the choice of the contour level does not affect the insights drawn here. For the majority of the cases (AR 2, α 26 • ), the length of the tip vortex increases almost linearly with the angle of attack. The decay rate remains mostly constant at d|Γ x |/dx ≈ 0.05 due to viscous diffusion. Significant increase of the decay rate is observed at α ≈ 26 • − 30 • with AR 2 and α 20 • with AR = 1. In these cases, the strong wake vortex shedding induces unsteady motion of the tip vortex, as observed in figure 9 . This unsteadiness serves as an additional mechanism (to diffusion) to reduce the time-averaged strength of the tip vortex. As a result, the L−α relationship as shown in figure 10(c) no longer follows the linear growth as at smaller α. At AR = 1 and 2, the time-averaged tip vortex length appears to even decrease at large angles of attack.

As discussed in §3.2.2, for large-aspect-ratio wings, vortex dislocation occurs, dividing the spanwise vortex shedding region. Here, we focus on the particular case of (α, AR) = (22 • , 6). The crossflow velocity u y along (x, y, z) = (2, −0.3, 0 − 7) over time is presented in figure 11 (a). For z 5, which corresponds to the tip vortex region, the velocity remains steady over time. For 0 z 5, the variations in the velocity indicate the occurrence of the periodic vortex shedding. Along with the vortical structures in figure 9 for (α, AR) = (22 • , 6), vortex dislocation is clearly observed from the discontinuities in the spanwise vortices. With the point of dislocation translating repeatedly from z ≈ 3.5 to z ≈ 1, the wake also changes over time, as observed from the iso-surfaces of Q = 1 in figure 11 (c). At t = 95, the spanwise vortical structures appear intact. Around t = 100, vortex dislocation starts to develop in the near wake. By t = 105, the formation of two shedding cells are clearly visible.

To gain additional insights, we apply Fourier transform to the crossflow velocity u y at each spanwise section, as shown in figure 11(b) . We useû y to denote the power spectral density (PSD) of the velocity fluctuation. The two spanwise shedding structures are found to be associated with different frequencies. The structures closer to the mid-span shed with a frequency of St 1 = 0.439 (St ≡ f c/U ∞ , where f is the dimensional shedding frequency). On the other hand, the structures closer to the wing tip shed with a lower frequency of St 2 = 0.372. To identify the vortical structures associated with these two frequencies, we perform dynamic mode decomposition (Schmid 2010; Rowley et al. 2009 ) on the Q criteria. As shown in figure 11(d ), the modal structures corresponding to St 1 is composed of the columnar vortex cores aligned with the spanwise direction. In contrast, the mode with frequency St 2 presents three-dimensional modal structures comprised of skewed vortex cores as well as the streamwise vortical structures. These two modes overlap over 1 z 3. As time evolves, the time-varying phase between these two modes results in the changes in the vortical structures, as observed in figure 11(c) .

The dominant frequencies identified in the wake for (α, AR) = (22 • , 6) are also detected in the corresponding lift spectrum, as shown in figure 12 . Also overlaid are the lift spectra for other aspect ratios with α = 22 • . The dominant frequency St 1 for AR = 6 is close to the frequency in the two-dimensional case. This suggests the two-dimensional nature of the corresponding modal structure (mode 1) as shown in figure 11(d ) . On the other hand, St 2 almost coincides with the primary frequency in the cases of AR = 2 and 4, hinting that the lower shedding frequency is caused by the slow down of shedding from the tip effects. The downward induced velocity from the tip vortex slows the shedding of the neighboring vortical structures. For the much lower AR = 1 case, the strong interaction between the tip vortex and the unsteady vortices further slows down the vortex shedding, resulting in a lower shedding frequency.

The lift spectra for different AR reveal how the tip effects influence the wake dynamics along the span. For AR 4, the tip effects dominate over the whole span of the wing, as discussed in §3.2.1 and §3.2.2. For a higher aspect ratio, the tip effects become weaker over the extended region of z 4. This allows nominally two-dimensional structures to develop near the mid-span and compete with the three-dimensional shedding near the tip, giving rise to vortex dislocation. Based on the insights obtained from spectral analysis, we are able to identify mainly three wake regions along the spanwise direction for largeaspect-ratio wings: (i) the tip vortex region as a direct consequence of the wing tip, (ii) the nominally two-dimensional shedding region near the mid-span (mode 1) where the influence of the tip vortex is weak, and (iii) the three-dimensional shedding region (mode 2) that are influenced by the interaction between regions (i) and (ii).

We summarize our wake characterization over aspect ratio and angle of attack in figure  13 . At Re = 400, the wake remains stable for α 12 • for the aspect ratios considered in this study. For AR = 1, the flow can remain steady for even a higher angle of attack of α = 16 • due to the stabilizing tip effects. We shade in blue the conditions for stable wake in figure 13 . The transition from steady to unsteady flow is attributed to the Hopf bifurcation associated with the periodically shedding vortical structures (Ahuja et al. 2007; Taira & Colonius 2009 ).

With the increase in angle of attack, the vortex sheet generated from the leading and trailing edges form shedding vortices, as highlighted by the yellow region in figure 13 . At this stage, the downwash provided by the tip vortex is no longer able to keep these wake vortices stationary. The wake here is characterized by the presence of two types of vortices: the unsteady shedding vortices and the tip vortex. The unsteady vortices take the form of spanwise vortex cores with their braid-like structures forming closed loops as described in §3.2.1.

For higher-aspect-ratio wings, the unsteady wake exhibits an additional feature of vortex dislocation. With the tip vortex imposing downwash on the shedding vortices near the tip, their shedding frequency becomes lower than that of the shedding structures over the mid-span region. When the discrepancy between these two frequencies become significant, the spanwise shedding processes become desynchronized, resulting in vortex dislocation as shown in figure 11(b) . This phenomenon is observed for AR 4 and is shaded in purple in figure 13 . This type of vortex dislocation is ubiquitous for bluff body flows with spanwise inhomogeneity (Eisenlohr & Eckelmann 1989; Williamson 1989; Noack et al. 1991; Techet et al. 1998; Lucor et al. 2001 ).

As we consider higher angles of attack, the tip vortex exhibits noticeable undulation in the streamwise direction due to its strong interaction with the unsteady shedding vortices. Such cases are observed over the red region in figure 13 . We notice that the transition curve separating the unsteady tip vortex region from the yellow and purple regions is similar to the periodic/aperiodic transition in Taira & Colonius (2009) . This suggests that the strong interaction between the tip vortices and the unsteady shedding vortices provides a route for the flow to transition to an asymmetric wake, if the symmetry plane is removed. Such transition is likely caused by the out-of-phase oscillations of the tip vortices at the two ends of the full wing. Unsteady tip vortex can also interact with the spanwise vortices experiencing vortex dislocation to create a complex wake as indicated by the top right subregion of figure 13. This complex wake generates finer-scale vortical structures from the rich interactions among the tip vortices and the large shedding vortices.

The three-dimensional wake imposes spanwise variations in the sectional distributions of the forces on the wing. The sectional lift and drag forces for AR = 4 are presented for various angles of attack in figure 14 . For the cases of α = 8 • and 14 • , C l gradually decreases from the mid-span to the wing tip. The drag coefficient remains almost constant over the majority of the span except for the slight increase at the wing tip.

For angles of attack above 20 • , considerable temporal fluctuations of C l appear in the unsteady shedding region, especially near z ≈ 2. The temporal fluctuations in the drag coefficients are significantly lower than the lift coefficients. As the flow becomes steady towards the wing tip region, the temporal fluctuation of C l gradually disappears. We note here that the time-averaged sectional C l exhibits a peak at z ≈ 3.5 before it drastically drops at the wing tip. This is in contrast to the cases with lower angles of attack, where C l decays monotonically towards the tip. This difference in the lift profile will be related to the difference in the near-tip vortical structures in the discussion below.

Let us now resort to the force element theory (Chang 1992; Lee et al. 2012) to identify the wake structures that exerts aerodynamic forces on the wing. In this theory, we utilize an auxiliary potential φ L (take the lift force for example) that satisfies the boundary condition −n · ∇φ L = n · e y . By taking the inner product of the Navier-Stokes equation (2.1b) with the potential velocity ∇φ L and integrating over the entire fluid domain V , the pressure force term S pn·e y dS in equation (2.2) can be expressed as a volume integration of wake variables and a surface integration. The lift force can then be expressed as

The integrands in the first and second terms on the right hand side are called the volume lift and the surface lift elements, respectively. The drag force elements can be obtained in a similar fashion by replacing the velocity potential φ L with φ D and enforce −n · ∇φ D = n · e x on the wing surface. The lift elements for the case of (α, AR) = (26 • , 4) at two time instances, at which the instantaneous lift coefficient is lowest and highest, are shown in figures 15(a) and (b), respectively. In the current study, the contribution of the viscous surface lift element is below 15% of the total force. For this reason, we only report below the dominant volume lift element ω × u · ∇φ L . Due to the spatially compact nature of the of the auxiliary potential gradient ∇φ L , the lift elements are concentrated in the vicinity of the wing. The two free shear layers emanated from both the suction and pressure sides contribute positively to lift, whereas the boundary layer on the pressure surface contributes adversely to the total lift. Similar observations have also been reported on an impulsively started finite-plate wing by Lee et al. (2012) . Comparing the snapshots between two time instances, the difference is observed in the unsteady shedding region. The difference in the lift elements between the two time instances at z = 1 is subtle. However, at z = 2, significant amount of negative lift elements appear in figure 15 (a) and disappear in figure 15(b) . This explains the larger lift fluctuations at z ≈ 2 than the other regions, as observed in figure 14(a) .

Near the tip at z = 3.5, the local flow is steady and the distributions of lift elements at two time instants are indistinguishable from each other. Noteworthy here is the large distribution of positive lift element near the wing tip giving rise to increase of sectional lift at around z = 3.5 in figure 14 . Inward from the tip vortex is a vortical structure that forms by rolling up the leading-edge vortex sheet, which is seen from figure 16 for higher angles of attack. This also explains the absence of such local lift increase at lower angles of attack, as the rolled-up vortical structure from the leading edge does not form next to the tip vortex.

Farther towards the tip at z = 3.9, negative lift elements emerge in the near wake and the positive lift elements diminish, giving rise to the drastic drop in the sectional lift force. It is noteworthy that the trailing part of the tip vortex barely contributes to the aerodynamic forces. This is not only because of the fast decay of the velocity potential ∇φ L , but also due to the fact that in the core of the tip vortex, the velocity vector and the vorticity vector are both aligned with the streamwise direction. The cross-product ω × u is small in magnitude, thus making the integrand in equation (3.2) negligible.

The volume drag elements for the same case are shown in figure 15 (c). In general, the drag elements share similar patterns with that of the lift elements, where the two separated free shear layers contribute positively to the total drag and the boundary layer on the pressure side contributes negatively. The iso-surfaces of drag elements, compared with the lift element at the same time instances, is devoid of the small-scale structures. Moreover, the intensity of the drag elements is lower compared to the lift elements. Thus, the integrated drag is smaller than lift in both temporal fluctuations as well as time-averaged value, as depicted in figure 14 .

We further examine the effect of aspect ratio on the sectional lift coefficients for α = 22 • in figure 17 . Here, the force coefficients are presented by aligning the distributions from the wing tip. With the increase in AR, the lift coefficients increase, although they remain below the two-dimensional lift value at the same angle of attack. Focusing in the tip region, the peak in sectional lift near the tip is observed for finite-aspect-ratio wings. This is caused by the presence of the leading-edge vortical structure near the tip as discussed above. Such structure gives rise to the increase in lift about 0.3 to 0.4 chords away from the tip. Similar to the lift coefficients, the drag coefficients increase with aspect ratio. However, the spanwise variations for the drag coefficients are relatively small compared to lift. The increase of the sectional C d near the tip is also observed for all aspect ratios at the angle of attack of α = 22 • . At last, let us summarize our aerodynamic force data for the finite-aspect-ratio wings. The time-averaged lift coefficients C L , drag coefficients C D , and the lift-to-drag ratios C L /C D are plotted in figure 18 . Overall, we notice that the lift and drag coefficients on the wing are smaller for lower-aspect-ratio wings and tend toward the two-dimensional limit for higher-aspect-ratio wings.

For the lift coefficients shown in figure 18 (a), we observe that the tip effects reduce lift even at low angles of attack. For post-stall flows, the influence of the tip vortices on lift becomes more noticeable with wider variations in the lift values with respect to AR. At low Reynolds number of 400, lift coefficients monotonically increase for the range of angles of attack considered. The enhanced lift at higher angles of attack is created by the vortical structures that form above the wing as seen in figure 15 (Taira & Colonius 2009; Lee et al. 2012) . We can observe that above α ≈ 12 • , there is a change in the lift slope which is caused by the emergence of unsteadiness in the wake, as shown in figure  13 .

Next, let us examine the variation of drag coefficient over the angle of attack. As in the case with lift, we observe that the drag coefficients become higher as the aspect ratio increases. This is due to the flow being able to pass around the tip to reduce the drag coefficient for lower-aspect-ratio wings. However, we notice a striking feature that the drag coefficient is almost independent of the aspect ratio up to α ≈ 12 • . This suggests that the drag force is primarily due to viscous effects even though the flow is separated at the moderate angles of attack. Once the wake becomes unsteady beyond α ≈ 12 • , pressure effects contribute significantly to drag, which then is influenced by the aspect ratio and the resulting unsteady wake profiles. These observations are also in agreement with the findings for finite-aspect-ratio flat-plate wings from Taira & Colonius (2009) .

At this low Reynolds number, the lift-to-drag ratio is limited to O(1), as evident from figure 18(c). Another unique feature for such low-Re aerodynamics is the high angle of attack at which the peak lift-to-drag ratio is achieved. These ratios are generally in agreement with those for finite-AR flat-plate wings from (Taira & Colonius 2009 ), but achieve their peaks at slightly higher angles of attack for the NACA 0015 airfoil. Such observations are reported in Torres & Mueller (2004) for higher Reynolds number of 10 5 . If the finite wings are to be operated in steady condition, it can still achieve a high lift-to-drag ratio for α 12 • . The ratio can be further increased but unsteadiness in the forces from the wake dynamics needs to be considered.

We performed direct numerical simulations of the flow over unswept NACA 0015 wings with straight cut tip at a Reynolds number of 400 with the objective of characterizing the tip effects on the three-dimensional separated flows. The present study considered a half-span wing model with symmetry boundary condition prescribed at the mid-span. An extensive parametric study was performed for aspect ratios of AR = 1 to 6 and angles of attack of α = 0 • to 30 • .

Based on these simulations, we have made several observations on the emergence of three-dimensional separation and wake patterns. The inception of three dimensionality is closely tied to the formation of the tip vortex. Vorticity is introduced to the flow from the wing surface in a predominantly two-dimensional manner. As the vortex sheet rolls up around the wing tip to form the tip vortex, three dimensionality is introduced to the flow but is confined to the near-tip region of the boundary layer. Once the tip vortex reaches its maturity, the stronger downwash influences the roll-up of the leading and trailing-edge vortex sheets, generating three-dimensional vortical structures in the wake and imposing a three-dimensional flow pattern in the vicinity of the wing surface.

The three-dimensional wake is dependent on the aspect ratio and angle of attack of the finite-aspect-ratio wing. For wings with low AR, the strong downwash from the tip vortex inhibits the roll-up of the leading and trailing-edge vortex sheets over the entire span. For this reason, the wake is able to remain steady for moderate angles of attack. For wings with higher AR, the relatively weakened tip effects near the mid-span are no longer able to keep the leading-and trailing-edge vortex sheets from rolling up. As a result, unsteady shedding takes place away from the tip. These shedding vortices are in the forms of spanwise vortex cores and braid-like structures. Together they close the vortex loops as they shed into the wake. Further increasing the AR brings out the locally two-dimensional vortical structures at the mid-span. Their competition with the neartip three-dimensional shedding vortices results in vortex dislocation. For high angles of attack, the wake becomes complex with interactions among the wake vortices. With streamwise oscillations appearing at higher angles of attack, the wake may transition to asymmetric flow about the mid-span in the absence of symmetry boundary condition. These wake features are mapped over a range of aspect ratio and angles of attack.

The three-dimensional wake dynamics has a significant impact on the sectional distribution of aerodynamics forces. By using the force element analysis, we identified the dominant vortical structures responsible for exerting lift and drag forces on the wing. Noteworthy here are the high sectional force coefficients near the tip for cases at high angles of attack. Such local increases in force coefficients are attributed to the stationary vortical structure formed inward from the tip vortex. From the simulations, the aerodynamic force coefficients and the lift-to-drag ratio were also summarized for a wide range of aspect ratios and angles of attack. Through a comprehensive analysis of the forces and wake dynamics, the relationship between the aspect ratio and the unsteadiness in the aerodynamic characteristics of the wing were revealed.

This study offered a detailed look into the influence of the tip effects on the formation of three-dimensional wake behind finite-aspect-ratio wings. The insights obtained from this study forms a foundation for further understanding the complex flow dynamics at higher Reynolds number and wake generation by unsteady wing maneuvers. Moreover, the knowledge gained here can serve as an important reference to active flow control strategies to stabilize or take advantage of vortical forces.

Theory of wing sections, including a summary of airfoil data

Lowdimensional models for control of leading-edge vortices: equilibria and linearized models

Vorticity transport mechanisms governing the development of leading-edge vortices

Measured aerodynamic characteristics of wings at low Reynolds numbers

Fundamentals of aerodynamics

Spanwise flow and the attachment of the leading-edge vortex on insect wings

On the evolution of the wake structure produced by a low-aspect-ratio pitching panel

Nektar++: An open-source spectral/hp element framework

Effect of ωx contour level on the length and decay rate of the tip vortex

Potential flow and forces for incompressible viscous flow

The leading-edge vortex and quasisteady vortex shedding on an accelerating plate

Floquet stability analysis in the wake of a NACA0015 airfoil at post-stall angles of attack

Wing rotation and the aerodynamic basis of insect flight

Wake topology and hydrodynamic performance of low-aspect-ratio flapping foils

A parallel stability analysis of a trailing vortex wake

Active attenuation of a trailing vortex inspired by a parabolized stability analysis

Vortex splitting and its consequences in the vortex street wake of cylinders at low Reynolds number

Leading-edge vortices: mechanics and modeling

Leadingedge vortices in insect flight

Further visualization of combined wing tip and starting vortex systems

Investigation of the unsteady tip vortex structure on a NACA 0012 wing at fixed incidence

Effects of three-dimensionality on thrust production by a pitching panel

Progress report on observations of three-dimensional flow patterns obtained during stall development on aerofoils, and on the problem of measuring two-dimensional characteristics. Her Majesty's Stationery Office

Energy conservation in collocated discretization schemes on unstructured meshes

Accurate and stable finite volume operators for unstructured flow solvers

Reynolds number and aspect ratio effects on the leading-edge vortex for rotating insect wing planforms

Three-dimensional separation on finite aspect ratio swept back wings

Wake dynamics of finite aspect ratio wings. Part III: Triglobal linear stability analysis

Linear instability of low Reynolds number massively separated flow around three NACA airfoils

Linear instability in the wake of an elliptic wing

Vorticity generation and transport

Surface flow and vortex shedding of an impulsively started wing

Vortex dynamics around pitching plates

Characterizing a burst leading-edge vortex on a rotating flat plate wing

Vorticity forces on an impulsively started finite plate

Rotational accelerations stabilize leading edge vortices on revolving fly wings

Viscous flow around finite length circular cylinder. Master's thesis

Vortex dislocations and force distribution of long flexible cylinders subjected to sheared flows

Unsteady aerodynamic characteristics of a translating rigid wing at low Reynolds number

The generation and decay of vorticity

Transient growth in the near wake region of the flow past a finite span wing

On cell formation in vortex streets

The structure of two-dimensional separation

Low Reynolds number aerodynamics of low-aspectratio, thin/flat/cambered-plate wings

Multiple bifurcations of the flow over stalled airfoils when changing the Reynolds number

Spectral analysis of nonlinear flows

Dynamic mode decomposition of numerical and experimental data

Recent progress in flapping wing aerodynamics and aeroelasticity

Three-dimensional flows around low-aspect-ratio flat-plate wings at low Reynolds numbers

An experimental study on the structure of the vortex street behind a circular cylinder of finite length

Vortical patterns behind a tapered cylinder oscillating transversely to a uniform flow

Low aspect ratio aerodynamics at low Reynolds numbers

Three-dimensional vortex formation on a heaving low-aspect-ratio wing: Computations and experiments

Separation of three-dimensional flow

Oblique and parallel modes of vortex shedding in the wake of a circular cylinder at low Reynolds numbers

The effects of leading edge modifications on the post-stall characteristics of wings

Flowfield model for a rectangular planform wing beyond stall

On vortex shedding from an airfoil in low-Reynolds-number flows

Flow structure on finite-span wings due to pitch-up motion

We acknowledge the US Air Force Office of Scientific Research (Program Managers: Dr. Douglas Smith and Dr. Gregg Abate, Grant number: FA9550-17-1-0222) for funding this project.

In §3.2.3, the evaluation of the length L and decay rate d|Γ x |/dx of the tip vortex involves the selection of ω x contour level. We have selected ω x = −1.5 as a representative contour level and compiled the characteristics of the tip vortex for various cases. Here, we examine two other contour levels to show that the overall findings are not influenced by the choice of contour level. The lengths and decay rates of the tip vortex based on ω x = −1.0 and −2.0 are shown in figure 19 . As qualitatively observed in figure 10(b) , the absolute values of L and d|Γ x |/dx are dependent on the choice of contour level. The larger the selected value of |ω x | is, the shorter the length and the larger the decay rate become. However, the variations of L and d|Γ x |/dx with respect to α, AR and Re within each contour level follow the same trends with those based on ω x = −1.5. This suggests that observations we have made in §3.2.3 is independent of the choice of the contour level considered in the present study.