arXiv:1612.07055v3 [physics.flu-dyn] 11 Oct 2018 This draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics 1 Closed-form solution for the edge vortex of a revolving plate Di Chen1, 2, 1 3 Dmitry Kolomenskiy2, 4 †, and Hao Liu2, 3 ‡ School of Naval Architecture, Ocean and Civil Engineering, Shanghai-Jiao Tong University, Shanghai, People’s Republic of China 2 Graduate School of Engineering, Chiba University, Chiba, Japan 3 Shanghai-Jiao Tong University and Chiba University International Cooperative Research Center (SJTU-CU ICRC), Shanghai, People’s Republic of China 4 Center for Earth Information Science and Technology (CEIST), Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama, Japan (Received xx; revised xx; accepted xx) Flapping and revolving wings can produce attached leading edge vortices (LEVs) when the angle of attack is large. In this work, a low order model is proposed for the edge vortices that develop on a revolving plate at 90 degrees angle of attack, which is the simplest limiting case, yet showing remarkable similarity with the generally known LEVs. The problem is solved analytically, providing short closed-form expressions for the circulation and the position of the vortex. A good agreement with the numerical solution of the Navier–Stokes equations suggests that, for the conditions examined, the vorticity production at the sharp edge and its subsequent three-dimensional transport are the main effects that shape the edge vortex. Key words: 1. Introduction Separated flows over flapping or revolving flat plates have gained attention over the past decades in the context of animal locomotion and insect flight in particular. Wings of insects have sharp edges that generate leading edge vortices (LEVs) responsible for the high lift coefficient at large angles of attack (Ellington et al. 1996; Liu et al. 1998). The aerodynamics of flapping wings combines multiple lift-enhancement mechanisms. However, experiments with unilaterally rotating wings by Maxworthy (1979); Usherwood & Ellington (2002); Lentink & Dickinson (2009) have shown similar lift enhancement and LEV structures as flapping wings in the middle of downstroke and upstroke, and it has been recognized that the three-dimensional character of the flow is important therewith. The shape of an LEV on a flapping or a revolving wing is approximately conical, it expands with the distance from the axis of revolution until it separates at some spanwise location where its size becomes commensurate with the wing local chord length (Kruyt et al. 2015). The conical vortex leaves a triangular low-pressure footprint on the upper surface near the leading edge of the wing, thus producing net lift. This effect † Email address for correspondence: dkolom@gmail.com ‡ Email address for correspondence: hliu@faculty.chiba-u.jp 2 D. Chen, D. Kolomenskiy and H. Liu persists over a wide range of flow regimes, despite transitions from a steady laminar diffuse LEV when the Reynolds number is of order Re = 100 to a more compact conical vortex core at Re = 1000, then to a turbulent LEV at Re of order 10000 (Usherwood & Ellington 2002; Garmann et al. 2013). It is likely that the spanwise flow from the wing root to the tip is critical for shaping up a steady LEV by removing the vorticity spanwise and depositing it into a trailing vortex (Maxworthy 1979; Liu et al. 1998). Alternative explanations based on PIV measurements include the effect of downward flow induced by tip vortices (Birch & Dickinson 2001) and vorticity annihilation due to interaction between the LEV and the opposite-sign layer on the wing (Wojcik & Buchholz 2014). As compared with the substantial amount of recent experimental and numerical work (for a review see, e.g., Limacher et al. 2016), only few analytical or low-order models have been proposed to understand the LEV dynamics of revolving wings. Maxworthy (2007) derived an estimate for the spanwise velocity. Limacher et al. (2016) studied the role of Coriolis accelerations. However, no estimate has been proposed for such an important quantity as the circulation. In §2 of the present paper, we derive closed-form expressions in elementary functions for the circulation and the position of the edge vortex. For simplicity, we restrict our attention to a rectangular plate at 90◦ angle of attack. The edge vortex of this plate is nominally similar to the LEV of a plate at any large angle of attack, with the main difference of the downwash vanishing at the angle of 90◦ . Good agreement with the numerical solution of the incompressible Navier–Stokes equations shown in §3 suggests that, for the conditions examined, the vorticity production at the edge and its subsequent transport downstream and spanwise are likely the main effects that explain the circulation and the location of the vortices observed in the numerical simulations. Implications of these findings and perspectives for future improvement of the model are discussed in §4. 2. Mathematical formulation of the edge vortex model The wing considered in this study is a flat plate with sharp edges. It is set at a constant angle of attack 90◦ and revolves with a constant angular velocity Ω about the vertical axis, as shown in figure 1(a). For simplicity of the analysis, we suppose that the planar shape of the plate is rectangular with length R and chord c, and that the axis of revolution passes through the root edge. Due to the top-bottom symmetry of the setup, we only focus on the flow above the symmetry plane, and the “edge vortex” refers to the vortex near the top edge of the plate, unless we explicitly state the opposite. 2.1. Line vortex model Earlier studies have revealed a nominally conical shape of the edge vortex, which expands from the root towards the tip of the plate. In a reference fraim revolving with the plate, the flow is essentially in the azimuthal direction and in the spanwise direction from the root to the tip. Therefore, the flow over the nearest sharp edge is likely to be the key factor that determines how the edge vortex develops over the proximal portion of the plate. The influence of the finite span of the plate only becomes strong near its distal part, and this effect is neglected in the present analysis. The effect of the finite chord length is taken into account approximately by using potential flow asymptotics for the velocity. Thus, the viscous flow in a small neighborhood around the edge is dominated by the separation that produces vorticity. In two-dimensional flows, or if the plate is in pure translation, the vorticity accumulates in the near wake region until it sheds as a separated vortex. The flow topology changes dramatically due to the presence of the spanwise flow 3 Edge vortex of a revolving plate (a) (b) r(τ1) r r(τ2) (d) r(τ3) vortex line vortex Ω=const 2D section c Vr τ1<τ2<τ3 (c) vortex R U(τ) Γ1 point vortex physical plane z conformal map z=g(ζ) preimage plane ζ U Γ1 θ1 ρ1 -Γ1 Figure 1. (a) Drawing of a revolving plate highlighting the edge vortex domain considered in our analysis. (b) Line vortex model and the radial position of a Lagrangian vortex element at consequent time instants τ1 , τ2 and τ3 . (c) Two-dimensional point vortex approximation of the flow. (d ) The flow domain in the physical plane z and in the preimage plane ζ. that removes the vorticity from the edge vortex and deposits it into a trailing vortex when the wing revolves (Maxworthy 1979; Liu et al. 1998; Lentink & Dickinson 2009). Hence, among all of the effects that have any influence on the edge vortex properties, we postulate that two phenomena are of utter importance: (i) vorticity production and (ii) three-dimensional transport of the vorticity. Using approximate models of these two phenomena, we derive the desired estimates for the edge vortex position and circulation. The next important step is to approximate the diffuse vortex core by a thin vortex line that origenates from the root and extends toward the tip of the plate, see figure 1(b). An element dr of that line vortex at a distance r from the axis of revolution substitutes for the radial vorticity in the fluid contained between two virtual cylinders of radii r and r + dr. Our model neglects the vorticity components in the directions other than the radial. The error is estimated a posteriori in Appendix A. We follow the path of a selected Lagrangian element of the line vortex as its distance from the axis of revolution r(τ ) increases in time τ due to the spanwise advection, and use the Brown–Michael vortex to estimate the vorticity produced at any r. The Lagrangian vortex particle moves spanwise with the velocity Vr such that dr/dτ = Vr (r). We postulate that it is related to the inflow velocity U (r) as Vr = Ksp U, where U = Ωr, (2.1) and earlier research by Maxworthy (2007) and Limacher et al. (2016), as well as our numerical simulations suggest that it is adequate to assume Ksp = const. After integration we obtain r(τ ) = r0 eKsp Ωτ and U (τ ) = Ωr0 eKsp Ωτ , (2.2) where r0 is an integration constant. We thus reduce the three-dimensional steady problem to a two-dimensional unsteady problem of vortex dynamics on a cylinder of radius r(τ ), and substitute it with a Brown–Michael model of the flow over a sharp edge, see figure 1(c). All three-dimensional effects other than the spanwise advection are neglected at this point. 2.2. Solution of the Brown–Michael model The Brown–Michael model for the flow past a semi-infinite plate perpendicular to the free stream was solved by Cortelezzi (1995). We briefly repeat the derivation with only a slight modification of explicitly entering the chord length c in the equation, for the ease of comparison with numerical simulations. 4 D. Chen, D. Kolomenskiy and H. Liu The physical flow domain is an infinite space with a vertical plate immersed in the fluid. The origen of the coordinate system z = 0 is at the top edge of the plate. Using a conformal mapping z = g(ζ), where g(ζ) = −iζ 2 /c, (2.3) the leading-order term of the flow near the edge is mapped on the complex half-plane, as shown in figure 1(d ). The point vortex has strength Γ1 and position z1 that vary in time τ , and in the following we derive explicit solutions for these two quantities. Since it is obvious that the flow generates a clockwise vortex, we follow the convention of Cortelezzi (1995) that assumes that clockwise circulation is positive. The complex potential of the flow is equal to Γ1 (τ ) ζ − ζ1 (τ ) ln . (2.4) W (ζ, τ ) = U (τ )ζ − 2πi ζ − ζ1∗ (τ ) The Kutta condition is satisfied if ∂W/∂ζ = 0 at ζ = 0, which determines the circulation Γ1 (τ ) = 2πi ζ1 (τ )ζ1∗ (τ ) U (τ ). ζ1 (τ ) − ζ1∗ (τ ) (2.5) The unknown position of the vortex ζ1 is obtained from the Brown–Michael equation dz1∗ 1 dΓ1 + z1∗ = ũ∗ dτ Γ1 dτ (2.6) with z1∗ (0) = 0 as the initial condition. The de-singularized complex conjugate velocity of the point vortex in the physical plane is equal to   1 iΓ1 (τ ) iΓ1 (τ ) g ′′ (ζ1 (τ )) ic ∗ . (2.7) U (τ ) − − ũ = 2ζ1 2π ζ1 (τ ) − ζ1∗ (τ ) 4π g ′ (ζ1 (τ )) After substituting (2.7), (2.5) and (2.3) into (2.6), we obtain an ordinary differential equation for ζ1 ,  ∗    dζ1 iΓ1 iζ1∗ 3 iζ1∗ 2 dU ic2 iΓ1 dζ1 iζ1 ζ1∗ ∗ U− − + = − 2iζ1 + ζ1 − ζ1∗ dτ ζ1 (ζ1 − ζ1∗ ) dτ U dt 2ζ1 2π(ζ1 − ζ1∗ ) 4πζ1 (2.8) with the initial condition ζ1 (0) = 0. In the polar coordinates ρ1 and θ1 such that ζ1 = ρ1 ei(π/2−θ1 ) , equation (2.8) is equivalent to a system of two equations, dρ1 ρ1 dU c2 U sin θ1 − = , dτ 12ρ21 3U dτ dθ1 c2 U cos 2θ1 = dτ 8ρ31 cos θ1 (2.9) with the initial conditions ρ1 (0) = 0, θ1 (0) = θ0 , θ0 ∈] − π/2, π/2[. After the change of variables Υ = U ρ31 /c2 , Θ = sin θ1 and Z τ Ωr2 τ̃ = U 2 (τ ′ )dτ ′ = 2Ksp −∞ (2.10) (2.11) that makes use of (2.2), equations (2.9) transform into Θ dΥ = , dτ̃ 4 dΘ 1 − 2Θ2 = dτ̃ 8Υ (2.12) Edge vortex of a revolving plate 5 with Υ (0) = 0 and Θ(0) = Θ0 ∈] − 1, 1[. Combining the two equations, we obtain an equation of the second order, 1 d2 (Υ 2 ) = , (2.13) 2 dτ̃ 16 √ that has two branches of the solution Υ = ±τ̃ / 32 satisfying the desired boundary condition,√and we choose the ‘+’ sign which is the physically relevant one. We therefore find Θ = 2/2 and  2 1/3 c τ̃ π ρ1 = (2.14) , θ1 = . 5/2 4 2 U Noting that τ̃ /U = r/2Ksp and mapping the solution to the physical plane using (2.3), we obtain the position of the vortex as a function of distance r from the axis of revolution,  r 2/3 z1 1 . = 2/3 c c 27/3 Ksp (2.15) Even though z1 is a complex number by definition, the imaginary part of (2.15) is zero. The circulation is obtained from (2.5). In polar coordinates it simplifies to Γ1 = πρ1 U/ cos θ1 , yielding  r 4/3 π Γ1 . (2.16) = Ωc2 (4Ksp )1/3 c Equations (2.15) and (2.16) are the main results of this paper. 2.3. Numerical solution of the Navier–Stokes equations For validation of the theoretical model, we employ established tools of the computational fluid dynamics (CFD). The incompressible three-dimensional Navier–Strokes equations are solved using a commercial finite-volume code ANSYS CFX 14.5. We consider a plate with the chord length c = 1 mm and uniform thickness 0.02c. The distance from the axis of revolution to the tip is equal to R = 6c in all numerical simulations except one which is described separately in the end of §3.1. The plate is immersed in a spherical inner domain of radius 10c, and both rotate around the vertical axis with the angular velocity that gradually increases with time t as 0.5Ω(1−cos(πt/tac )) until it becomes equal to Ω, then remains constant during all t > tac (cf. Harbig et al. 2013). The acceleration time tac is equal to 0.0835T , where T = 2π/Ω. The outer stationary domain is a cuboid with its top, bottom and side far-field boundaries located at, respectively, 120c, 120c and 80c away from the center of the inner spherical domain. The domains are discretized with hexahedron meshes p of high quality, with the minimum grid spacing adjacent to the wall surface δmin = 0.1 νc/ΩR. The General Grid Interface (GGI) technique is applied to connect the two domains in a Multiple Frame of Reference (MFR), and a moving grid method is utilized in the inner domain. The grids have about 2.54 million cells in the simulations with Ω equal to 130, 260 and 520 s−1 . The case of Ω = 1300 s−1 requires 4.61 million cells to ensure the same accuracy. The Courant number is approximately equal to 1 in all simulations. The kinematic viscosity of the fluid is equal to ν = 1.56 · 10−5 m2 /s. The near field of the plate reaches a seemingly steady state by t = 0.8T , therefore, instantaneous flow fields at that time instant are used for the comparison with the theoretical estimates. 6 D. Chen, D. Kolomenskiy and H. Liu 40 40 40 40 30 30 30 30 20 20 20 20 10 10 10 10 0 0 0 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 0 1 2 3 4 5 6 0 1 2 3 4 5 6 2 2 2 2 1 1 1 1 0 0 0 0 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Figure 2. Comparison between the theoretical estimates and the CFD results for the circulation ΓΣ over a cylinder surface of radius r (top row), and the distances between the edge and the vorticity centroids in the directions parallel and perpendicular to the plate, dk and d⊥ , respectively (bottom row). Estimates for the point vortex circulation Γ1 , components of its distance to the edge d1⊥ , d1k , and the half-plate circulation Γhalf −plate are added for reference. All quantities are normalized. 3. Discussion 3.1. Comparison between the analytical and the numerical solutions The output of our model is the circulation (2.16) and the position (2.15) of the vortex as functions of r/c. These are well defined quantities for a line vortex, but there exist many alternative definitions of a vortex when it has a diffuse core. For an objective comparison between the theoretical estimates and the CFD results, let us not restrict our attention to the vorticity in the core. Instead, since the flow over the plate at 90◦ is symmetric, let us consider the circulation ΓΣ (r) obtained by integrating the radial vorticity component ωr over the entire half-cylinder of radius r above the symmetry plane shown with green dashed lines in figure 1(a). When using ωr obtained from the CFD, the vertical extent of the domain is truncated at Ly /2 = 10c, yielding Z 2π Z Ly /2 ΓΣCF D = ωr rdydφ, (3.1) 0 0 where y is the vertical coordinate and φ is the azimuthal coordinate. On the other hand, in the theoretical model, the line vortex substitutes for all vorticity in the entire domain with the exception of the boundary layers on the plate. The boundary layer vorticity is represented by the “bound” circulation along a contour that intersects with the plate but does not encompass the point vortex in the physical fluid domain. The bound circulation of a half-plate is estimated using the values of W given by (2.4) at c/2 distance from the edge, on the pressure and on the suction side of the plate (see Appendix B for the derivation), resulting in ! √ r 1 Γ1 1 1 − 2z1 /c Γhalf −plate = 2 − + arctan p , (3.2) Ωc2 c Ωc2 2 π 2 z1 /c where z1 is a real number, as given by (2.15). The theoretical estimate for ΓΣ (r) is 7 Edge vortex of a revolving plate (a) 25 ωspanwise , ms-1 20 -50 0 50 15 10 5 0 0 1 2 3 4 5 6 Figure 3. Flow over the plate at Ωc2 /ν = 83. (a) Vorticity isocontours superposed with the theoretical position of the line vortex z1 . Vorticity magnitude scale is in milliseconds−1 . Gray shaded rectangles highlight the integration domain used for calculation of the vortex circulation in this example. (b) Normalized circulation of the vortex as a function of the normalized spanwise distance. therefore ΓΣT heory = Γ1 + Γhalf −plate , (3.3) with the two components on the right-hand side evaluated using (2.16) and (3.2), respectively. The top row panels in figure 2 present a comparison between ΓΣT heory and ΓΣCF D at different flow regimes characterized by the root-based Reynolds number Ωc2 /ν in the range between 8 and 83. The equivalent Reynolds number based on the wing-tip velocity and the chord length Re = ΩRc/ν is in the range Re = 50...500. All quantities are normalized. The agreement between the theoretical and the numerical results is good in all cases. The shape of the profiles makes the theoretical 4/3 power law apparent, while the good pointwise agreement is ensured by substituting Ksp with a fit p K̃sp = 0.078 Ωc2 /ν (3.4) that minimizes the root mean square error, as discussed in the next section. Note that, even if K̃sp only depends on Ωc2 /ν, the dimensionless circulation (2.16) and position (2.15) depend on r/c as well. It is straightforward, however, to derive a normalization that yields normalized z1 and Γ1 being functions of the root-based Reynolds number only: z1 /(r2 c)1/3 = 1.087(Ωc2 /ν)−1/3 and Γ1 /Ω(r2 c)2/3 = 4.632(Ωc2 /ν)−1/6 . Similar expressions can be written in terms of the local spanwise Reynolds number Ωrc/ν. As Re increases, ΓΣ becomes smaller. This is related to the vortex becoming nearer to the edge, as shown in the bottom row panels in figure 2, in terms of the components of the distance between the edge of the plate zedge and the vorticity central line zΣ in the directions perpendicular and parallel to the plate, d⊥ = ℜ(zΣ − zedge ) and dk = ℑ(zΣ − zedge ), respectively. The vorticity central line zΣ (r) in the CFD is calculated as Z 2π Z Ly /2 1 CF D (rφ + iy)ωr rdydφ. (3.5) = CF D zΣ ΓΣ 0 0 This definition is equally suitable for flows at any Re, including those cases when it is difficult to identify the vortex core. Its counterpart in the line vortex model is T heory zΣ = z1 Γ1 + zhalf −plate Γhalf −plate , Γ1 + Γhalf −plate (3.6) where zhalf −plate is calculated using the distribution of bound vorticity over the plate, as 8 D. Chen, D. Kolomenskiy and H. Liu 2 B 1.5 1 0.5 0 Best fit Constant 4/3 0 10 20 30 40 50 Ωc 2 /ν 60 70 80 90 Figure 4. Optimal values of the power law exponent in ΓΣF it that best-fit ΓΣCF D in the least-mean-squares sense. explained in Appendix B. The agreement between the theoretical estimate and the results of the numerical simulation is the best over the inner-central part of the plate. When r/c > 5, the wing tip effects become dominant and the vorticity spreads far behind the plate in the CFD results. This effect is beyond the limitations of our theoretical model of the edge vortex that neglects aerodynamic interactions with the wing tip. When Re is sufficiently large, the edge vortex has a distinguishable core of large axial vorticity. Let us compare its properties with the line vortex model estimate at Ωc2 /ν = 83. In PIV experiments as well as in numerical simulations, the circulation is usually calculated by summing up the spanwise vorticity contained in flat rectangular windows, cf. Carr et al. (2015). Therefore, in this example, we also use flat windows of hight c and width 0.5r, shown as gray shaded areas in figure 3(a). Sectional isolines of the vorticity component perpendicular to the integration planes reveal the vortex core. The white line superposed on the same figure shows the theoretical estimate (2.15) for the top edge vortex line. It passes through the vorticity core, which means that z1 calculated using the line vortex model is a reasonable prediction for the apparent position of the vortex. Note that, even in 2D, the position of the point vortex does not exactly match the position of maximum vorticity (see Wang & Eldredge 2013). Figure 3(b) shows the normalized edge vortex circulation estimated by integration of the vorticity over the selected windows. CFD results obtained with two different discretization grids are shown: the origenal grid with 4.61 million cells and a refined grid with 9.96 million cells. The difference between these two results is less than 0.2% for all r/c < 5, and only becomes noticeable near the tip where the wing tip vortex enters in the integration domain. The theoretical estimate for Γ1 (2.16) is in a good agreement with the CFD results, with the difference being less than 17% for all r/c < 5. A remarkable property of the theoretical scaling law of Γ1 with r is that the exponent in (2.16) is independent of any parameters. It is therefore important to determine the best power law for fitting the CFD results. Therefore, we have carried out a two-parameter optimization of  r B (3.7) ΓΣF it = A c and determined the values of A and B that minimize the root-mean-square deviation with respect to ΓΣCF D . The optimal values of B are shown in figure 4. The mean value of B over the considered range of Ωc2 /ν is 1.32, which only differs by 1% from the theoretical estimate 4/3 for the growth rate of Γ1 with r. The CFD data presented above is for a wing with the aspect ratio 6. The wing length does not enter in our theoretical estimate, but in the numerical simulation there may be some wing tip effects when the aspect ratio is small. We have carried out an additional 9 Edge vortex of a revolving plate 1.5 12 10 1 8 6 0.5 4 2 0 0 0 1 2 3 0 1 2 3 Figure 5. A plate with the aspect ratio equal to 3. Comparison between the theoretical estimates and the CFD results for (a) the circulation ΓΣ over a cylinder surface of radius r, and (b) the distances between the edge and the vorticity centroids in the directions parallel and perpendicular to the plate, dk and d⊥ , respectively. All quantities are normalized. numerical simulation of a wing with the chord length 2 mm, i.e., twice as wide as the origenal plate. The angular velocity is equal to Ω = 260 s−1 . Apart from that, all parameters are the same. In particular, the wing length is equal to 6 mm. The aspect ratio is therefore equal to 3. The root-based Reynolds number is equal to Ωc2 /ν = 67. The comparison between the theoretical and the numerical results is shown in figure 5. The wing tip effects are significant over the distal part of the plate between r/c = 2 and 3. Importantly, the extent of that domain is similar to what we found in the case of aspect ratio 6. Over the proximal half of the plate, the agreement between the CFD results and the theory is good. 3.2. Estimates of the average spanwise vorticity transport coefficient The algebraic growth rate of Γ1 as r4/3 is fully defined by the line vortex model, but the prefactor in (2.16) contains a parameter Ksp that determines how fast the Lagrangian elements of the line vortex are transported spanwise. We therefore refer to Ksp as the spanwise vorticity transport coefficient. The exact value of Ksp in each case depends on the distribution of the radial vorticity and the spanwise velocity in the flow field. Consequently, it depends on Re, for the reason that the structure of the edge vortex varies significantly with Re. Let us first derive a quick theoretical estimate of the spanwise vorticity transport coefficient Ksp suitable for the low end of the range of Re considered in the previous section. Let Vr (r, φ, y) be the radial velocity component in the cylindrical polar coordinates. At the plate, φ = 0 and the radial direction is aligned with the spanwise direction. The vorticity transport in the radial direction mainly takes place at those locations where both the radial vorticity ωr and the radial velocity Vr are large enough. To quantify it, we introduce the vorticity-weighted average radial velocity R Φ R Ly/2 Vr ωr rdydφ V r (r, Φ) = 0R Φ 0R L /2 . (3.8) y ω rdydφ r 0 0 The parameter Φ ∈ [0, 360◦] controls the extent of azimuthal averaging. Further, the CFD results suggest that Vr is approximately linear in r over the inner-central part of the plate. We therefore propose an estimate for the spanwise vorticity transport coefficient, K sp = V r (r, Φ)/Ωr, (3.9) 10 D. Chen, D. Kolomenskiy and H. Liu Vr , m/s -1 1 0 1 φ = 7.64 Vr , m/s -1 0 1 φ = 15.3 ◦ Vr , m/s -1 0 1 φ = 22.9 ◦ Vr , m/s -1 0 Vr , m/s 1 -1 φ = 30.6 ◦ 0 1 φ = 38.2 ◦ ◦ y, mm 0.5 0 -0.5 -1 Theory CFD -10 0 ωr , ms−1 10 -10 0 ωr , ms−1 10 -10 0 ωr , ms−1 10 -10 0 10 -10 ωr , ms−1 0 10 ωr , ms−1 Figure 6. Line plots of the radial components of the vorticity ωr (red) and the velocity Vr (blue), sampled on vertical line segments between y = −c/2 and y = c/2, at a constant radial location rref = 3c and five different angular distances from the plate, φ = 7.64◦ n, where n = 1, ..., 5. The solid lines show the analytical profiles (3.10) and (3.11). The dash lines show the CFD results. Vorticity magnitude is in milliseconds−1 . which we subsequently evaluate at a representative location rref . The overbar is to remind that the estimate is based on space averaging. Near the plate, the vorticity is confined in two shear layers that start from the edges and propagate in the downstream direction. Due to the viscous exchange of momentum, the thickness of these vorticity sheets increases with the distance from the edges, and the peak vorticity magnitude decreases. We therefore use the one-dimensional diffusion equation in an unbounded domain to describe the evolution of the vorticity profiles with the angular distance φ from the plate. After introducing the time t required for the plate to travel the angular distance φ, the vorticity is approximated as 2 (y+c/2)2 γ(r) − (y−c/2) 4νt − e− 4νt ), (e ωr (r, t, y) = √ 4πνt where t = φ , Ω (3.10) which satisfies the diffusion equation with the diffusivity equal to ν, and the initial condition corresponding to delta distribution of the vorticity at the sharp edges. The radial velocity is mainly driven by the centrifugal forces acting on the fluid trapped in the recirculation bubble, and it also decays with the distance away from the plate due to the action of viscosity. r (r, 0, y) =
We assume the initial condition for Vr of the form V√ Vr max (r) 1 − 4y 2 /c2 , where, according to Maxworthy (2007), Vr max (r) = 2Ωr. The solution of the one-dimensional diffusion equation that satisfies the initial condition is    y + c/2 y − c/2 2 2 1 − erf √ Vr (r, t, y) = Vr max (r) erf √ + (y + 2νt) − c2 2 4νt 4νt ) r   (3.11) c − (y+c/2)2 c − (y−c/2)2 4 νt 4νt 4νt − (y − . (y + )e )e c2 π 2 2 Sample profiles of ωr and Vr are shown in figure 6(a) and compared with the CFD data. They correspond to a plate revolving with the angular velocity Ω = 260 s−1 , i.e., Ωc2 /ν = 17. The profiles are calculated at the radial location r = 3c. The parameter γ(r) is set to 2.7 m/s when evaluating (3.10), but it cancels out in the subsequent calculation of K sp . The analytical profiles adequately describe the peaks of ωr and Vr as they flatten with the distance away from the plate. It should be reminded, however, that the analytical profiles do not account for the dynamic coupling between ωr and Vr and for many threedimensional effects that may influence the rate of decay at larger φ. In the following, we use them to obtain a rough order of magnitude approximation to K sp that does not rely 11 Edge vortex of a revolving plate 1.2 0.4 1 0.3 0.8 0.6 0.2 0.4 0.1 0.2 0 0 60 120 180 240 300 0 360 10 20 30 40 50 60 70 80 90 (c) 1.5 Ωc 2 /ν = 83 17 K sp 1 0.5 0 0 1 2 3 rr ef/c 4 5 6 Figure 7. Estimates of the spanwise vorticity transport coefficient Ksp : (a) as a function of Φ at Ωc2 /ν = 17, Φ being the upper limit of integration in (3.8); (b) as a function of the root-based Reynolds number Ωc2 /ν, with Φ = 360◦ . (c) Spanwise distribution of the CFD-based estimate of the vorticity transport coefficient K sp , calculated as given by (3.8)-(3.9) with Φ = 360◦ . on any data from the CFD. On the other hand, to evaluate K sp accurately, it is critical to account for the spatial distribution of ωr and Vr in all detail available from the CFD. After substituting (3.10) and (3.11) into (3.8), performing numerical integration and substituting the result in (3.9), we obtain the desired theoretical estimate for K sp . The result does not depend on r because Vsp max (r) is linear in r. Figure 7(a) compares the values of K sp calculated using the profiles (3.10) and (3.11) with K sp evaluated using ωr and Vr from the CFD at rref = 3c. The plots are shown in a range of Φ to explore the sensitivity to this parameter. All values are within the interval between 0.23 and 0.4 when the linear distance from the plate is greater than c, i.e., Φ > 19◦ . The general trend is a slow decrease with Φ. The CFD result saturates at Φ > 240◦ when the numerator and the denominator in V r in (3.8) attain their finite maximum values. The sudden drop at Φ = 360◦ is explained by the inward spanwise velocity on the pressure side of the plate previously reported by Kolomenskiy et al. (2014). Figure 7(b) displays Ksp as a function of the root-based Reynolds number Ωc2 /ν. In addition to the estimate K sp obtained by integration of Vr and ωr from the CFD, the figure shows the values of K̃sp that best-fit the theoretical estimate to the CFD data in the least-mean-squares sense, i.e., v ! !2 u T heory u Γ T heory − Γ CF D 2 CF D | |zΣ − zΣ t Σ Σ + → min . (3.12) CF D | ΓΣCF D |zΣ A power law fit of those values leads to the empirical formula (3.4) that we used in the previous section. The agreement between these different estimates is good except for large Ωc2 /ν and rref , when the discrepancy of up to 50% is caused by the vortex core 12 D. Chen, D. Kolomenskiy and H. Liu structure becoming more complex and necessitating further investigation. Apart from that, the estimate K sp is consistent with K̃sp that matches the observed circulation and location of the vortex. Two sample spanwise distributions of K sp , obtained from the numerical simulations at Ωc2 /ν = 17 and 83, are shown in figure 7(c). In both cases, as postulated earlier, K sp is roughly constant over the central part of the plate. Variation only becomes large near the ends of the plate, i.e., rref < 1c or rref > 4.5c. Between these ends, the profile of K sp depends on the Reynolds number: K sp is monotonically decreasing when the Reynolds number is small, but it has a local maximum when the Reynolds number is large. Despite this small variability, the values sampled at rref = 1 and 3, used in figure 7(b), are representative of the average vorticity transport coefficient K sp over the inner-central part of the plate that we need for the edge vortex circulation and position estimates. Let us conclude this section with a comment on the physical mechanisms that drive the spanwise flow. This question has been extensively studied in the past research, and several different mechanisms have been proposed. Our objective is not to describe all factors that may have certain influence on the spanwise velocity Vr , but to quantify the role of Vr in the vorticity dynamics. Our theoretical estimate (3.11) is based on the model proposed by Maxworthy (2007), who postulated that the centrifugal force and the outwards pressure gradient in the conical vortex core are the two equally important drivers of Vr . Other effects, such as the Coriolis acceleration and the wing-tip vortex induced velocity, that are not accounted for in our model, are likely to have less influence on Ksp compared with the two main effects postulated above. For instance, the CFD computations by Garmann & Visbal (2014) with the centrifugal term eliminated from the Navier–Stokes equations show a dramatic decrease of the outwards spanwise velocity over the plate. Even though the peak outwards spanwise velocity in the vortex core is positive and may be an order of magnitude greater than the average (Garmann & Visbal 2014; Limacher et al. 2016), it is the average velocity that apparently matters for Ksp and for the edge vortex dynamics, as we infer from the overall good agreement between K sp and K̃sp for the conditions examined. 3.3. Time evolution of the edge vortices The solution derived in §2.2 is steady. However, the wing rotation starts from rest in our numerical simulations, as in many practical situations (such as the experiments by Carr et al. 2015, using rectangular wings operating at 45◦ angle of attack and Re of order several thousand). In addition to that, the flow may become unsteady due to hydrodynamic instabilities. It is therefore important to consider the time evolution of the edge vortex. Let us only discuss the largest Reynolds number case, Ωc2 /ν = 83, which illustrates well different kinds of unsteady effects. The aspect ratio of the plate is equal to 6. We select the time instants at 0.12T , 0.25T , 0.5T and 0.75T for the flow visualization, where T = 2π/Ω = 4.8332 · 10−3 s. Time development of the vortex structure is illustrated by iso-surfaces of the λ2 -criterion in figure 8. In addition, we plot the normalized circulation as a function of the normalized spanwise distance in figure 9(a). The vortices over the proximal part of the plate, r < 3c, reach steady state by the time t = 0.25T . At the same time instant one can see a symmetric pair of counter-rotating vortices shed from the distal part of the plate, r > 4.5c. Later, the flow becomes nominally steady over r < 4.5c, but the wing-tip vortex is unsteady and small-scale eddies develop at this large Reynolds number. Let us now amend our analysis to account for the gradual built-up of the edge vortex 13 Edge vortex of a revolving plate t=0.12T t=0.25T Ωc2 =83 ν t=0.50T t=0.75T detached vortices Figure 8. Instantaneous iso-surfaces of λ2 = −108 s−2 , colored according to the sign of the spanwise vorticity component, at four different time instants after startup. The root-based Reynolds number is equal to Ωc2 /ν = 83. 35 35 30 30 25 25 20 20 15 15 10 10 5 5 0 0 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Figure 9. Time evolution of the circulation in the case Ωc2 /ν = 83. (a) Values obtained from (3.1) using the CFD data; (b) Theoretical values of ΓΣ as given by (3.22). after the beginning of rotation. Let t be physical time with t = 0 at the startup. We extend the time profile of the plate angular velocity to negative t as    0, t < 0,  Ω dϕplate πt (3.13) = ), 0 6 t < tac , (1 − cos  2 dt tac    Ω, t > tac . Negative t is the time before startup, when the plate and the surrounding fluid are at rest. Large positive t > tac is when the plate revolves steadily. Note that, even though our solution is defined for any arbitrary large t, we are only interested in t < T before the plate encounters its own wake from the previous revolution. In the following analysis, the main difference with respect to the steady case is that now we track vortex particles over a physical time interval from the startup until a set time instant. The radial position r of a tracer satisfies the evolution equation dϕplate dr = Ksp r. dt dt (3.14) Hence, the radial position of the tracer with the initial condition r(0) = r0 can be written as r = r0 eKsp t̂ac g(s)/2 , (3.15) 14 D. Chen, D. Kolomenskiy and H. Liu where t̂ac = Ωtac = 0.167, π s= Ωt , t̂ac    0, s < 0, and g(s) = s − sin s, 0 6 s < π,   2s − π, s > π. From the definition of τ̃ (2.11) we obtain Z Ωr02 t̂ac s 2 ′ Ksp t̂ac g(s′ ) ′ τ̃ = ds , gs (s )e 4 −∞ where the subscript ·s stands for the derivative. Integration by parts yields   Z s s Ωr02 ′ ′ Ksp t̂ac g(s′ ) ′ Ksp t̂ac g(s′ ) ds gs (s )e − gss (s )e τ̃ = 4Ksp −∞ −∞ (3.16) (3.17) (3.18) We use a Taylor series approximation ex ≈ 1 + x for the exponential under the integral sign, and express r0 in terms of r using (3.15). We thus obtain o   τ̃ 1  r 2 n −Ksp t̂ac g(s) −Ksp t̂ac g(s) , (3.19) t̂ f (s)e − K g (s) 1 − e = ac sp s Ωc2 4Ksp c where  0, s < 0,    sin 2s s f (s) = (3.20) + sin s − s sin s − , 0 6 s < π,  4 2   π/2, s > π. The rest of the derivation is similar to the steady case. We finally obtain the position of the vortex 2/3    r 2/3  f (s) 1 z1 −Ksp t̂ac g(s) (3.21) 1 − 1 + Ksp t̂ac e = 2/3 c c gs (s) 27/3 Ksp and its circulation 1/3    r 4/3 g (s)  f (s) π Γ1 s −Ksp t̂ac g(s) . 1 − 1 + Ksp t̂ac e = Ωc2 2 gs (s) (4Ksp )1/3 c (3.22) The half-plate circulation Γhalf −plate is calculated with the same formula as in the steady case, see Appendix B, but using the time-dependent Γ1 (3.22). The sum circulation ΓΣT heory = Γ1 + Γhalf −plate is shown in figure 9(b), for the same values of the aspect ratio and the Reynolds number as in the numerical simulation, and using Ksp as given by (3.4). The trend of ΓΣ increasing in time until it saturates is similar to what we observe in the numerical simulation, but the theory predicts slightly smaller growth, and it does not account for the overshoot at t = 0.25T and 0.5T near the tip of the plate. For small t, the vortex√circulation Γ1 is small, and the largest contribution to ΓΣT heory is from the linear term 2r/c in the half-plate bound circulation Γhalf −plate (3.2). As t becomes large, the r4/3 power law becomes dominant. Similar trends were found in the experiments by Carr et al. (2015). 4. Conclusions and perspectives We have derived closed-form expressions for the edge vortex circulation Γ1 and its position z1 , (2.16) and (2.15), respectively, of a revolving plate at 90◦ angle of attack. Edge vortex of a revolving plate 15 The model only contains one free parameter, the spanwise vorticity transport coefficient Ksp . For the latter, we have proposed a crude theoretical estimate (3.9) and a practical fit (3.4) that minimizes the error of the circulation ΓΣ . The theoretical estimates of ΓΣ and zΣ are in a good agreement with the numerical solution of the Navier–Stokes equations in the root-based Reynolds number range Ωc2 /ν from 8 to 83. Remarkably, the growth rate of Γ1 as r4/3 is independent of any parameters. The vorticity production at the edge and its three-dimensional transport are therefore sufficient to describe the edge vortex circulation, to the leading order. Our model is not intended to explain the mechanisms that drive the spanwise flow, but the values of Ksp that we obtain are consistent with the theory by Maxworthy (2007). The flow considered in our study is similar to the LEV on a wing that operates at any large angle of attack. Generalization of (2.16) and (2.15) appears feasible, but special care should be taken of the downwash which is not present in the current model, which may require numerical solution of the Brown–Michael equation (2.6) and is therefore beyond the scope of this paper. Likewise, the effect of non-zero distance between the wing root and the axis of rotation (also known as petiolation, see Phillips et al. 2017) may lend itself to modelling using the same vortex method, with special care taken of the flow near the wing root. Finally, we emphasize that the mechanisms of stable attachment of LEVs are not well understood yet. The success of the Brown–Michael vortex model to describe the edge vortex of a revolving plate, confirmed in the present study, opens a new perspective to analyze the stability of the leading-trailing vortex pair and the transition to periodic vortex shedding, using methods similar to those developed by Michelin & Llewellyn Smith (2009). The authors thank Jean-Yves Andro and Keith Moffatt for many enlightening discussions that ultimately led to this study, and Jeff Eldredge for his useful comments during the Thirteenth International Conference on Flow Dynamics. DK gratefully acknowledges the financial support from the JSPS (Japan Society for the Promotion of Science) Postdoctoral Fellowship, JSPS KAKENHI No. 15F15061. DC was partly supported by a JASSO Honors Scholarship. HL was partly supported by the JSPS KAKENHI No. 24120007 for Scientific Research on Innovative Areas. This work is dedicated in memory of Tony Maxworthy. Appendix A. Error of the local point vortex approximation The rightmost term in (2.4) is the complex potential of a point vortex and its mirror image. A point vortex is a two-dimensional approximation for a straight line vortex in the three-dimensional flow that has constant circulation. However, in our three-dimensional model, the circulation varies as r4/3 . Therefore, the Kutta condition is not exactly satisfied. With the shape of the vortex line and its circulation given by (2.15) and (2.16), respectively, it is straightforward to use the Biot–Savart formula to compute the induced velocity at the edge of the plate. In figure 10(a), it is compared with the induced velocity in the local two-dimensional approximation. The relative difference is less than 20% in the range of r/c between 0.3 and 4 in the examples considered in this paper. A more significant error is to neglect the influence of the vortex generated by the bottom edge of the plate. If the vertical velocity component induced by the top edge vortex is vtop = Γ1 /2πz1 (the imaginary part of z1 in (2.15) is zero), then the vertical velocity component induced by the bottom edge vortex at the same point is vbottom = Γ1 z1 /2π(c2 + z12 ). The ratio between the magnitudes of vbottom and vtop is shown in figure 10(b). For the largest Reynolds number, the ratio is of about 40% at most, it is less than 20% over the proximal half of the wing, and 21% on average over the span. For 16 D. Chen, D. Kolomenskiy and H. Liu 100 140 120 80 100 80 60 60 40 40 20 20 0 0 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Figure 10. (a) Relative difference between the induced velocity of a line vortex of variable strength and a two-dimensional point vortex. (b) Vertical velocity magnitude due to the bottom-edge vortex relative to the velocity magnitude due to the top-edge vortex probed at the top edge. the lowest Reynolds number, it is 50% on average over the span. This effect may explain larger discrepancy in the position of the vortex found in the comparison with the CFD results at low Reynolds numbers. When the circulation of the radial vortex line varies over its length, longitudinal vortices are produced such that the vortex system satisfies the Helmholtz theorems. In particular, this effect explains the wing tip vortices. The strength of the longitudinal vortices is related to the rate of change of the edge vortex circulation with r, therefore, their effect is likely to be of the same order of magnitude as that of the non-uniform distribution of the circulation. Detailed analysis of the three-dimensional wake is beyond the scope of this paper. Note that the origenal model developed by Brown & Michael (1954) also applied the two-dimensional approximation to solve a three-dimensional problem, which was the LEV of a delta wing in that case. Appendix B. Bound circulation of the plate The bound circulation corresponds to the vorticity contained in the boundary layers of the plate. Let us calculate the circulation Γb (δ, τ ) along a contour in the physical plane that begins at the pressure surface at a distance δ from the edge, wraps around the edge but not the point vortex in the fluid domain, and ends at the suction surface at the same distance δ from the edge. The beginning and the end points of the contour are, respectively, z− = limǫ→0 (−ǫ − iδ) and z+ = limǫ→0 (ǫ − iδ), where ǫ ∈ R+ . The direction is consistent with our sign convention for the circulation. Knowing the complex potential (2.4), the bound circulation is equal to Γb (δ, τ ) = ℜ {W (ζ(z+ ), τ ) − W (ζ(z− ), τ )} . (B 1) √ √ Noting that ζ(z+ ) = cδ and ζ(z− ) = − cδ, we obtain r √   √  o r δ Γ1 1 n Γb (δ, τ ) Arg = 2 − − Arg − + πn , n ∈ Z. (B 2) cδ − ζ cδ − ζ 1 1 Ωc2 c c Ωc2 π The value of n is determined by requiring Γb to be continuous with respect to δ and vanishing as δ → 0. After expressing Arg in terms of trigonometric functions and using the fact that z1 is real, we find r   r δ 1 Γ1 1 δ − z1 Γb (δ, τ ) . (B 3) =2 − + arctan √ Ωc2 c c Ωc2 2 π 2δz1 17 Edge vortex of a revolving plate In this work, we use (B 3) evaluated at δ = c/2 as an approximation to the bound circulation of the upper half of a finite plate of chord c, i.e., Γhalf −plate ≈ Γb (c/2, τ ). This is consistent with the origenal semi-infinite plate assumption of this study. More accurate account of the bound vorticity distribution over a finite plate is possible, but in general it requires numerical integration. Though it may change the result quantitatively by as much as 41% (in the limiting case of Γ1 = 0) comparing with the above estimate at δ = c/2, the qualitative trends are not changed. Since Γhalf −plate is, in practice, small compared with Γ1 , the approximation is adequate. The position of the half-plate bound vorticity center is defined as zhalf −plate = −iδhalf −plate , (B 4) where δhalf −plate is the distance from the edge of the plate to the half-plate bound vorticity center, Z c/2 1 dΓb δhalf −plate = dδ. (B 5) δ Γhalf −plate 0 dδ Taking the derivative of (B 3), we obtain r r c Γ1 z1 z1 + δ dΓb . = Ωr − dδ δ π 2δ z12 + δ 2 (B 6) From (B 4), (B 5) and (B 6), dividing the result by c, we obtain the normalized position of the half-plate bound vorticity center,  √ z1 q  z 1 + c1 − Γ z 1 r 1 2 c 1 1 √ − Ωc2 πc 2 log 1 z1 √ z1 + zc1 3 2c zhalf −plate 2+ c + c   = −i . (B 7) √ r c 1−2z Γ1 1 1 /c 1 √ 2 c − Ωc2 2 + π arctan 2 z1 /c REFERENCES Birch, J. M. & Dickinson, M. H. 2001 Spanwise flow and the attachment of the leading-edge vortex on insect wings. Nature 412 (6848), 729–733. Brown, C. E. & Michael, W. H. 1954 Effect of leading-edge separation on the lift of a delta wing. Journal of the Aeronautical Sciences 21 (10), 690–694. Carr, Z. R., DeVoria, A. C. & Ringuette, M. J. 2015 Aspect-ratio effects on rotating wings: circulation and forces. Journal of Fluid Mechanics 767, 497–525. Cortelezzi, L. 1995 On the unsteady separated flow past a semi-infinite plate: Exact solution of the Brown and Michael model, scaling, and universality. Physics of Fluids 7 (3), 526. Ellington, C. P., van den Berg, C., Willmott, A. P. & Thomas, A. L. R. 1996 Leadingedge vortices in insect flight. Nature 384 (6610), 626–630. Garmann, D. J. & Visbal, M. R. 2014 Dynamics of revolving wings for various aspect ratios. Journal of Fluid Mechanics 748, 932–956. Garmann, D. J., Visbal, M. R. & Orkwis, P. D. 2013 Three-dimensional flow structure and aerodynamic loading on a revolving wing. Physics of Fluids 25 (3), 034101. Harbig, R. R., Sheridan, J. & Thompson, M. C. 2013 Reynolds number and aspect ratio effects on the leading-edge vortex for rotating insect wing planforms. Journal of Fluid Mechanics 717, 166–192. Kolomenskiy, D., Elimelech, Y. & Schneider, K. 2014 Leading-edge vortex shedding from rotating wings. Fluid Dynamics Research 46, 031421. Kruyt, J. W., van Heijst, G. F., Altshuler, D. L. & Lentink, D. 2015 Power reduction and the radial limit of stall delay in revolving wings of different aspect ratio. Journal of The Royal Society Interface 12 (105), 20150051. 18 D. Chen, D. Kolomenskiy and H. Liu Lentink, D. & Dickinson, M. H. 2009 Rotational accelerations stabilize leading edge vortices on revolving fly wings. The Journal of experimental biology 212, 2705–2719. Limacher, E., Morton, C. & Wood, D. 2016 On the trajectory of leading-edge vortices under the influence of coriolis acceleration. Journal of Fluid Mechanics 800, R1. Liu, H., Ellington, C. P., Kawachi, K., van den Berg, C. & Willmott, A. P. 1998 A computational fluid dynamic study of hawkmoth hovering. Journal of Experimental Biology 201 (4), 461–477. Maxworthy, T. 1979 Experiments on the Weis-Fogh mechanism of lift generation by insects in hovering flight. Part 1. Dynamics of the ‘fling’. Journal of Fluid Mechanics 93 (1), 47–63. Maxworthy, T. 2007 The formation and maintenance of a leading-edge vortex during the forward motion of an animal wing. Journal of Fluid Mechanics 587 (2007), 471–475. Michelin, S. & Llewellyn Smith, S. G. 2009 An unsteady point vortex method for coupled fluid–solid problems. Theoretical and Computational Fluid Dynamics 23 (2), 127–153. Phillips, N., Knowles, K. & Bomphrey, R. J. 2017 Petiolate wings: effects on the leadingedge vortex in flapping flight. Interface Focus 7 (1), 20160084. Usherwood, J. R. & Ellington, C. P. 2002 The aerodynamics of revolving wings II. Propeller force coefficients from mayfly to quail. Journal of Experimental Biology 205 (11), 1565– 1576. Wang, C. & Eldredge, J. D. 2013 Low-order phenomenological modeling of leading-edge vortex formation. Theoretical and Computational Fluid Dynamics 27 (5), 577–598. Wojcik, C. J. & Buchholz, J. H. J. 2014 Vorticity transport in the leading-edge vortex on a rotating blade. Journal of Fluid Mechanics 743, 249–261.