Geometrically-exact unsteady model for airfoils undergoing large amplitude maneuvers

A semi-analytical, geometrically-exact, unsteady potential ﬂow model is developed for airfoils undergoing large amplitude maneuvers. For this objec-tive, the classical unsteady theory of Theodorsen is revisited relaxing some of the major assumptions such as (1) ﬂat wake, (2) small angle of attack, (3) small disturbances to the mean ﬂow components, and (4) time-invariant free-stream. The kinematics of the wake vortices is simulated numerically while the wake and bound circulation distribution and, consequently, the associated pressure distribution are determined analytically. The steady and unsteady behaviors of the developed model are validated against experimental and computational results. The model is then used to determine the lift frequency response at diﬀerent mean angles of attack. Both qualitative and quantitative discrepancies are found between the obtained frequency response and that of Theodorsen at high angles of attack.

Circulation pitching moment at the hinge point Γ j jth vortex at wake sgn(•) Signum function Re(•) Real part of complex value Im(•) imaginary part of complex value F X , F Z Horizontal and vertical components of the resultant forces, respectively Aspect ratio H Amplitude of oscillation k Reduced frequency ω Reduced frequency α ef f Effective angle of attack

Introduction
The unsteady nonlinear aerodynamic characteristics emerge in many engineering applications such as rapidly maneuvering aircraft, flow over wind turbines, and flapping-wing flight.While computational methods can be used to accurately estimate the flow quantities and the associated aerodynamic loads, there is always a need for low-to-medium fidelity aerodynamic model.Such models can be used in the preliminary design phase (where thousands of design alternatives are investigated), optimization and sensitivity analyses, and flight dynamics and control analyses.
Wagner [35] determined the time-response of the lift due to a step input (indicial response problem).Theodorsen [31] determined the frequency response of the lift (i.e., lift response due to harmonically oscillating input) and applied it to the flutter problem of fixed-wing aircraft.Von Karman and Sears [34] developed an integral equation governing the spatial distribution of the wake circulation.In fact, these three seminal works represent the basis for most of the unsteady aerodynamic models that were developed later.Particularly, motivated by evaluating the unsteady aerodynamic loads due to arbitrarily moving airfoils and assessing the associated flight dynamics and aeroelasticity, researchers developed finite-state unsteady aerodynamic models.Most of the finite-state models are obtained by approximating the curves of either the Wagner function in the time-domain, or the Theodorsen function in the frequency-domain (e.g., see the efforts of Garrick [9], Jones [12], Jones [13], and Vepa [33]).In contrast to these approximations, there are finite-state models that were developed from the basic principles (e.g., Peters and Karunamoorthy [19,20], Peters et al. [21], and Peters [17]).• Wagner [35] Challenges: Quasi-steady • Theodorsen [31] • Unsteadiness aerodynamics • Peters [17] • Nonlinearity and new phenomena Three-Dimensional: • Efficiency (computational burden) • Reissner [26] • Unsteady vortex lattice [14] Table 1 shows a taxonomy of the unsteady flow regimes and the applicability of the above referenced models.For a low reduced frequency (k < 0.1), quasi-steady aerodynamics may still be applicable.For a relatively high reduced frequency, the unsteady aerodynamic characteristics cannot be neglected.If small-amplitude maneuvers are to be investigated, then the attached flow assumption and, hence, linearity of the flow dynamics may be justifiable.There are many unsteady, linear aerodynamic theories that can be applied successfully and efficiently in this regime for both of the two-dimensional and three-dimensional applications.However, if large-amplitude maneuvers are concerned, the nonlinearity of the flow can no longer be neglected and non-conventional phenomena and/or lift mechanisms emerge (e.g., leading edge vortex in flapping flight).In this regime, there is a challenge of capturing both of the prominent unsteadiness and flow nonlinearity with a feasible computational burden so as to perform optimization, sensitivity analysis, dynamics, and control synthesis.
Because of the many associated applications, there has been a recent flurry of modeling the unsteady nonlinear aerodynamics associated with large amplitude maneuvers and non conventional lift mechanisms.Ansari et al. [1,2] extended the work of Von Karman and Sears [34] to include the leading edge vortex effect on flapping wings by shedding vortices from both leading and trailing edges.They derived two nonlinear integral equations for the shed wake and leading edge vortices.The non-compact representation of Ansari's model dictates quite complicated numerical implementation which was the main reason for very high frequency components of the flow quantities when used for a coupled aeroelastic-flight dynamic simulation of flapping wings by Su and Cesnik [28].Su and Cesnik reported that their coupled aeroelasticflight dynamic simulation could not complete one cycle because of such a problem.Ramesh et al. [23] developed an inviscid theoretical approach to account for non-planar wakes.However, the developed model still considers the induced velocities as small perturbations to the airfoil velocity, which considerably influence the aerodynamic forces at high angles of attack.In an earlier work, Taha et al. [29,30] extended the application of Duhamel's superposition principle, usually applied in linear unsteady aerodynamics, to account for non-conventional lift mechanisms.The main objective was to capture the leading edge vortex effect on flapping wings in an unsteady fashion.The developed model comprises two internal aerodynamic states for two-dimensional applications and four internal aerodynamic states for threedimensional applications using strip theory.The strength of the model of Taha et al. is that its asymptotic steady behavior can match any arbitrarily given steady lift variation with the angle of attack.However, its underpin-ning flow dynamics is that of the Wagner's response.That is, it assumes that the linearity of the flow dynamics still hold even for high angles of attack.Therefore, there is still a need to develop similar efficient models while capturing the nonlinearity of the lift evolution dynamics.
The objective of this work is to develop a hybrid analytical-numerical approach to determine the lift coefficient associated with unsteady aerodynamics that involve high angles of attack.For this purpose, we revisit the classical model of Theodorsen [31] relaxing the small angle of attack and the flat wake assumptions, aiming for a geometrically-exact unsteady aerodynamic theory for two-dimensional applications.Due to the challenges encountered when relaxing these assumptions, full analyticity of the model is not expected (e.g., the lift deficiency due to wake cannot be analytically represented via some explicit expression such as Theodorsen's function).In contrast, the developed model is semi-analytic; that is, the vortex kinematics is simulated numerically.But the wake and bound circulation distribution and, consequently, the associated pressure distribution are computed analytically.Hence, there will be no need to solve an algebraic system of equations to determine the bound circulation as usually done in the implementation of discrete vortex models such as the unsteady vortex lattice method.The developed model is applied to some canonical pitch motions such as the ones introduced by Eldredge et al. [7].The obtained results are then validated against the numerical and experimental results of Ramesh et al. [23].

Problem Formulation
We consider a two-dimensional flat plate undergoing horizontal, vertical, and rotational motions, as shown in Fig. 1.The horizontal velocity U , the vertical velocity ḣ, and the angular velocity α of the flat plate are taken positive to the left, downward, and clockwise directions, respectively.Two axis-systems are considered; the inertial frame X-Z and the plate-fixed frame x-z.While the origins of the two systems (O and O ′ ) are located at the mid-chord point, the plate is rotating about the hinge point which lies at a distance ab along the x-axis ahead of the mid-chord point.Here, b denotes the half chord length.
In local coordinates (x-z system), the relative wind speed is expressed as: where α is the angle of attack and u ′ and w ′ represent the components of the induced velocity (disturbance/perturbation velocity) in the x and z directions, respectively.Therefore, the velocity potential is written as: where φ ′ is the disturbance velocity potential that satisfies the following expression: ∂φ ′ ∂x = u ′ and ∂φ ′ ∂z = w ′ .Assuming rigid plate and applying the no-penetration boundary condition (i.e., the total flow velocity perpendicular to the plate is zero), we write Figure 2 shows the conventional Joukowski conformal mapping between the plate (ξ-plane) and the cylinder (η-plane), see [3,25] for example.The mapping is written as where, ξ = x + iz and η = x + iẑ = re iθ .Since the velocity vectors in the two planes are related via the transformation relation then it can be shown that u ′ (x, z = 0, t) and w ′ (x, z = 0, t) are related to q ′ r (r = b/2, θ, t) and q ′ θ (r = b/2, θ, t) via Also, taking the datum of the disturbance velocity potential at the leading edge (i.e., φ ′ (r = b/2, π, t) = 0), then the distribution of φ ′ on the cylinder surface is written as The aerodynamic loads are calculated by integrating the pressure distribution over the surface, which can be determined using the following unsteady Bernoulli's equation (see [3,6] for example): where p and ρ are the pressure and density of the air, and C(t) is a spatiallyconstant quantity, which can be calculated from the far field conditions.Thus, the pressure distribution is related to the disturbance velocities and velocity potential via where p ∞ is the free-stream pressure.
Integrating the pressure distribution to calculate the resultant aerodynamic loads may be deficient because it does not account for the leading edge suction force.We use Blasius theorem [14] to remedy this deficiency.According to Blasius theorem, the components (F x and F z ) of the aerodynamic loads acting on an immersed body are given by where u and w are the components of the fluid velocity vector in the body axes.Since the perturbations u ′ and w ′ near the leading edge outweigh by far the other contributions, one can write the components of the leading edge suction force as 1 sin θ dθ = π, we determine the components of the leading edge suction force as where F s is the suction force.Next, we follow Theodorsen's work [31] relaxing some of the simplifying assumptions such as (1) flat wake, (2) small angle of attack, (3) small disturbances to the mean flow components, and (4) time-invariant free-stream.Some of the applicable Theodorsen's results will be stated here without proof and emphasis will be given to the extended parts.Following Theodorsen's approach, we divide the problem into non-circulatory and circulatory contributions.

Non-circulatory Contributions
The non-circulatory contributions are included by considering a timevarying source and sink distributions on the cylinder's upper and lower surfaces of strengths H + (θ, t) and H − (θ, t), respectively.Thus, the component q ′ r of the induced velocity is written as Similar to Theodorsen's work, it can be shown that the component q ′ r of the induced velocity on the cylinder's surface is given by and that the component q ′ θ is given by[3] Using Eq. ( 7), the non-circulatory component of the disturbance velocity potential on the cylinder's upper surface is obtained as Substituting by the no-penetration boundary condition, Eq. ( 3), into Eq.( 16) and carrying out the integrations, the non-circulatory disturbance velocity potential is calculated as Let the superscripts (u) and (l) denote the upper and lower surfaces of the plate, respectively.Since w ′(u) = w ′(l) and φ ∂t ), then the non-circulatory distribution of the pressure difference across the plate can be derived from Eq. ( 9) as Substituting by φ N from Eq. ( 17) into Eq.( 18) and integrating the outcome over the cylinder's surface, we obtain the non-circulatory normal force and pitching moment at the hinge point as In addition to allowing for a time-varying free stream, Eq. ( 19) represents a geometrically-exact extension of Theodorsen's result [31] N Equation ( 19) is consistent with the fact that the non-circulatory force is the force required to accelerate a mass of air cylinder of radius b with the acceleration of the mid-chord point.However, it should be noted that the acceleration term between the brackets is not the inertial acceleration of the mid-chord point (i.e., not the acceleration with respect to the still fluid), but rather the time rate of change of the velocity of the mid-chord point with respect to the body axes.Also, it is noteworthy to mention that the non-circulatory lift and drag forces are given by L N = N N cos α and D N = N N sin α, respectively.

Circulatory Contributions
The Kutta condition dictates finite velocity at the trailing edge (θ = 0).To satisfy this condition, Eq. ( 6) yields q ′ θ (r = b/2, 0, t) = 0.However, Eq. ( 15) indicates the inability to achieve this requirement for arbitrary motion of the plate which is represented by w ′ via the no-penetration boundary condition.Hence, an additional contribution (circulatory one) to the tangential velocity is invoked such that both contributions could satisfy the Kutta condition together; that is Thus, a circulation distribution is added to represent the wake of the flow past the flat plate (and the associated circular cylinder), as shown in Fig. 3.
For each vortex −Γ j in the cylinder's wake, there has to exist an image vortex Γ j to maintain zero total circulation, as shown in Fig. 4, see Von Karman and Sears [34].Since one of the main objectives of this work is to account for wake deformation, the position of the wake vortex −Γ j is considered arbitrarily free.In order to keep the no-penetration boundary condition satisfied, the radial components of the velocities induced by Γ j and −Γ j must cancel each other.Thus, if the wake vortex is arbitrarily  located at η j = r j e iϕ j , then its image has to lie inside the cylinder at b 2 4r j e iϕ j , as shown in Fig. 4.
Also, as shown in Fig. 4, the j th wake vortex −Γ j and its image Γ j induce velocity vectors q − and q + , respectively, at an arbitrary point p on the cylinder.The tangential components of the total velocity vector q + + q − on the upper and lower surfaces are given by q ′(u) Recalling the non-circulatory contribution of q ′ θ from Eq. ( 15) and discretizing the wake sheet of vorticity into a finite number of point vortices, we write the unsteady Kutta condition as where N is the total number of discrete wake vortices.Substituting by w ′ from the no-penetration boundary condition, Eq. ( 3), and carrying out the integration in the left hand side, we write where w 3/4 is the component of the plate velocity relative to the air normal to the plate surface at the three-quarter chord point.Equation ( 23) will be used to determine the strength Γ j of the newly shed wake vortex at each time step in the numerical implementation, as will be shown later.
After solving for the wake vortex distribution at each time step, the circulatory loads are determined by using the unsteady Bernoulli's equation.Using Eq. ( 7), we determine the disturbance velocity potential on the cylinder due to the j th vortex −Γ j and its image Γ j as follows: Using the above result to determine u ′(u) and u ′(l) , knowing that w ′(u) = w ′(l) = 0, and substituting into Eq.( 9), we determine the circulatory pressure distribution on the lower and upper surfaces.Integrating the pressure difference over the plate's surface, we determine the normal force and hinge moment due to the jth vortex −Γ j and its image Γ j as where x j is the x-coordinate of the wake vortex −Γ j in the plate domain, r j and ϕ j are its polar coordinates in the cylinder's domain, sgn(•) is the signum function, Re(•) and Im(•) denote the real and imaginary parts of their arguments, respectively, and the terms p 1 , p 2 , p 3 and p 4 are written as As for the circulatory contribution of the pitching moment about the hinge point, it is written as The details of the derivation of Eqs.(25,26) are given in Appendix A.
The terms in Eq. (25,26) that are proportional to Γ 2 j stem from the nonlinear term u ′2 in Bernoulli's equation (9), which has been usually neglected in the previous studies.If we neglect this term, considered a small angle of attack, constant free stream, and a flat wake (ϕ j = 0), then the shedding velocity will be given by ẋj = U and, hence, Theodorsen's result [31] is recovered Equations (25,26) are more than a geometrically exact extension of Theodorsen's result, Eq. ( 27).In addition to considering the geometric nonlinearities, Eqs.(25,26) account for a time-varying (non-uniform) free stream, non-flat wake deformation, shedding by the local velocity ẋj rather than the free stream one, and finally for the higher order perturbation u ′2 -effect.As for the total suction force on the flat plate, we substitute by the noncirculatory and circulatory contributions of q ′ θ into Eq.( 12) to obtain

Numerical Implementation
In this section, we show the implementation of the unsteady model developed in the previous sections.The developed model is semi-analytical in the sense that closed form algebraic equations are written for the aerodynamic loads, though the wake convection needs to be performed numerically.As stated before, we discretize the wake sheet of vorticity into a finite number of point vortices.As shown in Eqs.(25,26), the coordinates of the discrete wake vortices and their velocities with respect to O ′ are required to determine the unsteady aerodynamic loads.The velocity of each wake vortex −Γ j comprises three components: (1) a component due to the undisturbed flow, (2) a circulatory component due to the other separated vortices and their images, and (3) a non-circulatory component due to the source/sink distribution over the cylinder/plate.All the components are transformed to the inertial axes (X-Z system) where they are added together and marched forward to predict the location of each point vortex in the next time step.
The first contribution to the components of the convection velocity of a wake vortex −Γ j that is due to the undisturbed flow, as shown in Fig. 1, is written as Figure 5(a) shows the induced velocity on a wake vortex −Γ j due to another wake vortex −Γ k and its image.Then, the components of the total induced velocity on a wake vortex −Γ j due to the other wake vortices and their images are given by q (2) Similarly, Fig. 5(b) shows the induced velocity on a wake vortex −Γ j due to an element of the source/sink distribution.Hence, the non-circulatory components of the total induced velocity on a wake vortex −Γ j are given by q (3) According to the transformation between the plate domain (ξ-plane) and the cylinder domain (η-plane), as presented in Eq. ( 5), the components of the above induced velocity vectors q (2) and q (3) in the ξ-plane are written as where p = 2, 3.These velocity components are then transformed to the inertial axes (X − Z) as (a) Induced velocity on a wake vortex −Γ j due to another wake vortex −Γ k and its image.
(b) Induced velocity on a wake vortex −Γ j due to the source/sink distribution.
Figure 5: The circulatory and non-circulatory contributions to the induced velocity on a wake vortex −Γ j .
Thus, the total velocity components of the wake vortex −Γ j with respect to the plate's center O ′ in the inertial frame are determined as We use a simple first-order Euler scheme to predict the new position of the wake vortex as Finally, the velocities ( ẋ and ż) and positions (x and z) in the body frame, as required in Eqs.(25,26), can be calculated from their inertial counterparts according to the inverse of Eq. (33).
As for the aerodynamic loads, the non-circulatory contributions are given directly in Eq. (19).After determining the positions and the velocities of the wake vortices, the circulatory loads induced by each wake vortex can be determined from Eqs. (25,26).Finally, the total suction force is determined from Eq. ( 28).Thus, the horizontal and vertical components of the resultant force are given by

Validation
Although the model is mainly developed for unsteady applications, it is of interest to validate its asymptotic steady behavior first in the following subsection, while the unsteady behavior will be validated in the next subsection.

Steady Behavior
The asymptotic steady behavior of the developed model is validated against the wind tunnel experimental results of Brandon [4] on the F-18 wing experiencing high angles of attack (5 • − 75 • ).We divide the wing into 100 segments, apply the developed aerodynamic model on each segment, and integrate over the span to obtain the total lift force L.Then, we use to account for the induced downwash due to the three dimensional effects in a simple way through the Extended Lifting Line Theory (see [27]).In this equation, A R denotes the aspect ratio of the wing.Figure 6 shows a comparison between the lift coefficient using Eq.(37), with the lift L determined from the proposed model, and the experimental data given in Ref. [4].Two model predictions are presented in Fig. 6; with and without considering the leading edge (LE) suction force.The plots show that the variation of the static lift coefficient with the angle of attack (C L -α curve) using the proposed model accounting for the LE suction force exactly matches that of the classical wing theory.This matching indicates that both have the same nonlinearity (sin α), that is induced by the no-penetration boundary condition.Moreover, since both of these curves closely match the experimental data up to α = 20 • , we can conclude that the nonlinearity induced by the no-penetration boundary condition is the dominant one in this range of the angle of attack (as well as the flow remains attached).At large angles of attack, the separation effects are more pronounced.Consequently, the contribution of the LE suction force is expected to diminish as shown by Dickinson and Gotz [5] and Usherwood and Ellington [32].Therefore, within the framework of potential flow, we account for the flow separation by neglecting the LE suction force.The fact that the predicted C L -α curve, ignoring the LE suction force, closely matches the experimental data shows that the effects of the geometric nonlinearities may dominate other separation effects.Although the F-18 wing section (NACA 65A005) is a typical cambered, thick airfoil with a rounded LE, it is envisaged that the high sweep angle of the wing helps stabilizing a LE vortex, consequently, eliminates abrupt stall from taking place and leads to a smooth lift-curve.This may explain the good matching between the F-18 wing experimental data and the potential flow results when ignoring the LE suction (clearly, the LE suction force and the LE vortex cannot coexist [22]).
Since the developed model is two-dimensional, it is still of interest to validate its steady behavior versus two dimensional experimental results.We choose the experimental results of Dickinson and Gotz [5] for validation.They built an experimental setup for a two dimensional flat plate undergoing an impulsive start from rest to some specific forward speed through almost a constant acceleration at a constant angle of attack.They spanned the whole positive range of the angle of attack and a range for low Reynolds numbers (between 79-236).Dickinson and Gotz found that for large angles of attack, a strong leading edge vortex (LEV) forms with a large diameter (of the order of a chord length).After about four chord lengths (4c), the vortex becomes extended rearwards and eventually shed.As this vortex is shed, a vortex of opposite sign is formed at the trailing edge.This alternating pattern continues, forming Von Karman street, which is similar to the numerical results of Wang [36].Dickinson and Gotz provided two polar plots for the measured lift and drag coefficients at 2c and 7c referring to them as "early" and "late" measurements, respectively.The late data correspond to the inception of the Von Karman street and the load oscillations while the early results correspond to the steady loads reached after the transient response.Fig. 7 shows a comparison between the lift and drag coefficients as determined by the proposed model and the early results of Dickinon and Gotz [5].
Figure shows that the developed model ignoring the LE suction gives a satisfactory trend for the lift and drag variation in comparison to the experimental data of Dickinson and Gotz.Is should be noted that, in contrast to the previous validation case, the experiment of Dickinson and Gotz is performed at a very low Reynolds number.This explains the deviation of the experimental data for the lift coefficient from the current results even at small angles of attack.As expected, considering the LE suction force in a potential flow framework predicts no drag, as shown in Fig. 7(b).The difference between the potential flow results for drag and the experimental data at zero angle of attack is attributed to the skin friction drag that is not captured within the potential flow framework.

Unsteady Behavior
The unsteady predictions of the proposed model are validated against the computational and experimental results of Ramesh et al. [23].They solved the two-dimensional Navier-Stokes equations on a flat plate undergoing large amplitude canonical pitch maneuvers that are presented in Ref. [7] and shown here in Fig. 8.The plate is pitching around the mid point uniformly from 0 • up to a specific angle of attack ( 25• and 45 • ), retains this angle for a while, and then returns back uniformly to 0 • with a total maneuver non-dimensional time of six.
Figures 9 and 10 show comparisons of the time variations of the resulted lift coefficients among the proposed model (with and without considering the LE suction force), the computational results and experimental data of Ramesh at el. [23], and the classical unsteady model of Leishman and Nguyen [15] for the 25 • and 45 • maneuvers, respectively.Noting that the Reynolds number of the experiment and simulation of Ramesh et al. [23] is 10,000 where delayed stall is evidenced [5], considerable flow separation is not expected to occur during the relatively small angle of attack maneuver ( 25 • ).Hence, it is expected that the role of the LE suction force would be significant in this maneuver.This explains the agreement between the experimental data and the lift curve predicted by the proposed model (accounting for the LE suction) as well as that of the classical unsteady theory, as shown in Figure 9 for the 25 • maneuver.However, when the maneuver amplitude becomes larger, separation effects become more pronounced.Hence, the LE suction is not expected to play a role.As shown in Fig. 10, the lift curves predicted by the proposed model (accounting for the LE suction) as well as the classical unsteady theory deviate considerably from the experimental data, particularly at the times where the large angles of attack are encountered ( U t c = 3 − 4).On the other hand, the lift curve of the proposed model, neglecting the LE suction, closely matches the experimental data.
For each of the two considered maneuvers, we provide a quantitative comparison among the models discussed above.We consider the experimental data of Ramesh at el. [23] as a benchmark and calculate the deviation of each of the other results with respect to it, as shown in Table 2. Two error metrics are calculated.These are the root mean squared error in the lift coefficient in percent of the maximum lift coefficient e rms = √ max(Cexp) × 100, and the mean value of the deviation in the lift coefficient in percent of the maximum lift coefficient e = |C L −Cexp| max(Cexp) .We use the over bar to denote an average quantity.Because low Reynolds number (10,000) is encountered, the associated delayed stall indicates that the flow will remain attached during the relatively small angle of attack maneuver (25 • ).Thus, in such a maneuver, the two error metrics of the proposed model (considering the LE suction) and the classical unsteady theory are comparable to that of the Figure 9: A comparison for the unsteady lift evolution among the proposed models (with and without considering the leading edge suction force), the computational results and experimental data of Ramesh at el. [23], and the classical unsteady model of Leishman and Nguyen [15], for the 25 • -maneuver.Figure 10: A comparison for the unsteady lift evolution among the proposed models (with and without considering the leading edge suction force), the computational results and experimental data of Ramesh at el. [23], and the classical unsteady model of Leishman and Nguyen [15], for the 45 • -maneuver.computational results.But when considerable separation effects are encountered in the large amplitude maneuver, the two error metrics of the proposed model (considering the LE suction) and the classical unsteady theory are too large in comparison to the that of the proposed model, without considering Table 2: Error metrics quantifying the deviations from the experimental results of Ramesh et al. [23] of the results of the proposed model with and without including the LE suction force, the computational results of Ramesh et al. [23], and the results of the classical unsteady model of Leishman and Nguyen [15].
× 100 e = the LE suction, and that of the computational results.It is noteworthy to mention that in addition to the efficiency of the proposed model (in terms of computational cost and number of degrees of freedom), it resulted in the minimum error metrics in comparison to all of the other models for the high angle of attack case.That is, the proposed unsteady aerodynamic model is efficient enough to be used in multi-disciplinary analyses (e.g., dynamics, control, and optimization) and also rich enough to cover the gaps that the classical unsteady theory cannot cover.
The capability of the developed model to capture the effects of unsteady free stream is validated against the predictions of Isaac's theory [11], Greenberg's theory [10], and Peters finite state theory [18] for the case of a sinusoidal variation of the free stream.The unsteady free stream is written as where α is kept constant.Figure 11 shows a comparison for the resulted normalized lift L 2πρbU 2 0 α (unsteady to quasi-steady ratio), for the case µ = 0.8 and k = 0.2, between the current model and the models discussed above.The simulation is performed at a small angle of attack (in the linear range) in order to have a meaningful comparison with the linear theories.As such, both models (ignoring and including the LE suction) give satisfactory results in comparison to the linear theories for this test case.A comparison for the normalized lift (unsteady to quasi-steady ratio) among the proposed models (with and without considering the leading edge suction force), Isaac's theory [11], Greenberg's theory [10], and Peters finite state theory [18] for a sinusoidally varying free stream.

Frequency Response at Different Angles of Attack
Having validated the proposed model for a high angle of attack maneuver, it is interesting to investigate how the frequency response (Theodorsen's main result) changes as the angle of attack increases.Theodorsen's model is based on a linear approximation for the flow dynamics, which results in a frequency response that is independent on the operating condition and/or the amplitude of the aerodynamic input (airfoil motion).The developed model accounts for the geometric and nonplanar-wake nonlinearities.As such, it will not result in a single frequency response.Rather, a different frequency response (i.e., linearized flow dynamics) will be obtained at different operating conditions (angles of attack).In this section, we determine the frequency response of the circulatory normal force at different angles of attack, namely 5 • , 20 • , and 40 • .At each value of these angles of attack, the following plunging maneuver is simulated: where H is the amplitude of oscillation normalized by the half-chord length (b) and ω is the frequency of oscillation.That is, the effective angle of attack is given by where α 0 is the mean/steady angle of attack (5 o , 20 o and 40 o ) and k is the reduced frequency.For each value of α 0 , we simulate the above plunging maneuver at various frequencies to determine the frequency response at this value of α 0 .For each frequency, the value of H is adjusted such that the amplitude of the oscillating angle of attack (k H) is 5 • ; that is, we simulate small angle maneuvers around different mean values of α ef f to determine the frequency response of the linearized aerodynamic system around the considered α 0 .
For each combination of α 0 and k, the above plunging maneuver is simulated and the steady state time history of the non-circulatory normal force is obtained.We then simulate the same maneuver (i.e., same α ef f ) at almost steady conditions (at very small value of k).Thus, the magnitude of the frequency response at this combination of α 0 and k is simply the ratio of the unsteady amplitude to the almost steady one and its argument (phase) is the phase shift between the two signals far out in time.
Figure 12 shows the magnitude, phase, and polar lots of the calculated frequency response at the considered three values of α 0 along with that of Theodorsen.For α 0 = 5 • , the obtained frequency response closely matches that of the Theodorsen function.If α 0 is increased to 20 • , good matching is still being seen with a bit of discrepancy at large reduced frequencies.A noteworthy observation is the large deviation at α 0 = 40 • , not only quantitatively but also qualitatively.The magnitude of Theodorsen's frequency response does not go below 1  2 .In fact, it asymptotically reaches such a value as k goes to infinity.However, Fig. 12(a) shows that the magnitude of the frequency response decreases for large values of α 0 and approaches zero as k goes to infinity.Moreover, unlike the phase of the Theodorsen function, which asymptotically reaches zero (in phase) as k goes to infinity, the phase angle of the frequency response decreases considerably at large angles of attack so that it asymptotically reaches −180 • (out of phase) as k goes to infinity.This observation is worth discussing.Note that because the phase of the Theodorsen function approaches zero at high reduced frequencies (i.e., the input and output signals are exactly in phase), it invokes modeling the high-frequency, oscillatory flows with quasi-steady means with some magnitude drop [29].However, the obtained frequency response (which accounts for the nonlinear wake deformation) refutes this concept because the phase shift approaches −180 • rather (i.e., the input and output signals are completely out of phase).

Model Validity, Limitations and Extensions
Although the developed model provides an efficient, satisfactory means of analyzing unsteady aerodynamics associated with large amplitude maneu-vers, it is still lying within the framework of potential flow and, hence, cannot capture some physical aspects.These include viscous friction, LE separation, dynamic stall.One important extension is to study the appropriate conditions for switching between including and ignoring the LE suction.In fact, there have been relevant recent efforts in Refs.[16,24].Ramesh et al. [24] introduced a new criterion for LE separation based on LE suction.Note that according to Garrick [8], the flow velocity at the LE, V LE , can be written as √ x is a measure of the LE suction force and γ is the bound vorticity distribution that is written in Fourier series as where x = c 2 (1 − cos θ).All the fourier series terms vanish at the LE except the A 0 -term, which is infinite at the LE (having 1 √ x singularity).Thus, S is finite and equal to S = √ cU A 0 (t), where c is the chord length and U is the free stream velocity.As such, the A 0 parameter is a measure of the flow velocity at the LE and the LE suction force.Since the LE separation is very related to the flow conditions at the LE, Ramesh et al. [24] proposed the A 0 parameter to be a measure for LE separation and called it "Leading Edge Suction Parameter(LESP)".If A 0 exceeds a certain limit (LESP cr ) (that depends on the airfoil type and the Reynolds number), a LE vortex is shed.They determine the strength of the newly shed LE vortex such that A 0 = LESP cr .Morris and Rusak [16] also studied the onset of steady LE stall on thin airfoils using the method of matched asymptotic expansion.In their formulation, the problem outer region is the flow around most of the airfoil chord and is modeled using thin airfoil theory.The inner region is the flow in the vicinity of the LE and is determined by solving the incompressible viscous Navier Stokes equation on a semi-infinite parabola.The two solutions are then made to asymptotically match each other.They showed the existence of a critical angle of attack (that depends on Reynolds number) that is responsible for the inception of LE stall.It should be noted that, for steady conditions, A 0 = 1 2 sin α.Thus, the LESP criterion of Ramesh et al. [24] is equivalent to the angle of attack criterion of Morris and Rusak [16].However, for unsteady cases, the LESP criterion is more representative than the angle of attack criterion because it accounts for the effective angle of attack over the whole airfoil surface.
For the high Reynolds number of Brandon's experiment on the F-18 wing [4], the critical steady angle of attack α cr ∼ 15 • according to Morris and Rusak [16].This is consistent with our results of the F-18 wing (shown in Fig. 6); including the suction force closely matches the experimental results up to α cr ∼ 15 • at which separation occurs.Above 15 • , since separation occurs, ignoring the LE suction closely matches the experimental results.On the other hand, considering the unsteady validation cases (shown in Figs. 9 and 10), LESP cr ∼ 0.18 [24], which corresponds to an angle of attack of α cr ∼ 21.1 • .This is also consistent with our results.Since most of the 25 • -maneuver lies below the critical condition, including the LE suction gives satisfactory results.On the other hand, during the 45 • -maneuver, considerable separation occurs (and consequently LE suction diminish).As such, ignoring the LE suction yields very comparable results to the computational and experimental data.
Another straightforward extension to the developed model is to allow emanation of a vortex sheet from the LE.The strength of the LE vorticity can be determined to satisfy the Kutta condition at the LE in the same way it is satisfied here at the TE.That is, the following two equations will be solved simultaneously to determine the strengths of the shed vorticity at both edges: q ′ θN (r = b/2, θ = 0, t) + q ′ θT E (r = b/2, θ = 0, t) + q ′ θLE (r = b/2, θ = 0, t) = 0 q ′ θN (r = b/2, θ = π, t) + q ′ θT E (r = b/2, θ = π, t) + q ′ θLE (r = b/2, θ = π, t) = 0 where the subscripts N , T E, and LE refers to the non-circulatory, trailing edge wake, and leading edge wake contributions, respectively.This extension will be more suitable to capture the LE separation effects in comparison to the mere manipulation of the LE suction force inclusion.

Conclusions
We developed a hybrid analytical-numerical approach to determine the lift coefficient associated with unsteady aerodynamics that involve high angles of attack.For this purpose, we revisited the classical Theodorsen's frequency response model and relaxed the major simplifying assumptions that led to limited region of applicability of Theodorsen's model such as (1) flat wake, (2) small angle of attack, (3) small disturbances to the mean flow components, and (4) time-invariant free-stream.By relaxing these assumptions, we managed to develop a geometrically-exact potential flow model.In the developed model, the vortex kinematics were determined numerically.However, unlike the discrete vortex models, the circulation distribution and the associated aerodynamic loads were determined analytically after solving for the vortex kinematics.
The asymptotic steady behavior of the developed model was validated against two-dimensional experimental data and on the F-18 wing showing a good matching up to 75 • angle of attack.The unsteady behavior of the developed model was validated against some experimental and computational results of canonical large-amplitude pitch maneuvers.The model also showed a good agreement with the experimental results in comparison to the classical unsteady theory without requiring high computational burden.The developed model was then used to determine the lift frequency response at various angles of attack.For small angles of attack, the obtained frequency response closely matches that of Theodorsen function.However, for high angles of attack (40 • ), both qualitative and quantitative discrepancies are observed between the obtained frequency response and that of Theodorsen.The obtained frequency response at the high angle of attack approaches zero magnitude and −180 • phase at large values for the reduced frequency, which is in contrast to the 1  2 magnitude and 0 • phase approached by Theodorsen frequency response.The developed model is efficient enough to be used in multi-disciplinary applications (e.g., dynamics and control) and also rich enough to cover some gaps that the classical theory of Theodorsen cannot cover.

Figure 1 :
Figure 1: A schematic diagram showing the flat plate motion and axis-systems.Note that the translational motion of the plate is considered in the wind speed

Figure 3 :
Figure 3: Wake of the flat plate and the associated circular cylinder.

Figure 4 :
Figure 4: A schematic diagram showing a vortex −Γ j in the cylinder's wake and its image inside the cylinder along with their induced velocities on the cylinder.

Figure 6
also shows the static lift coefficient using the conventional extended lifting line theory; that is C

Figure 6 :
Figure 6: Comparison between the static lift coefficient using the proposed model (with and without suction force), the classical wing theory, and the experimental data[4].
considering suction force Proposed model considering suction force Experimental data (a) Lift Coefficient.
suction force Proposed model considering suction force Experimental data (b) Drag Coefficient.

Figure 7 :
Figure 7: Comparison between the two-dimensional static lift and drag coefficient using the proposed model (with and without suction force) and the experimental results of Dickinson and Gotz [5].

Figure 11 :
Figure11:A comparison for the normalized lift (unsteady to quasi-steady ratio) among the proposed models (with and without considering the leading edge suction force), Isaac's theory[11], Greenberg's theory[10], and Peters finite state theory[18] for a sinusoidally varying free stream.
Magnitude of the frequency response.

oo
Proposed model around: 40 o (b) Phase angle of the frequency response.Proposed model around: 40 o (c) Polar plot of the frequency response.

Figure 12 :
Figure 12: Magnitude and phase plots along with the polar plot for the obtained frequency response at different angles of attack along with that of Theodorsen.
Velocities in the x and z directions, respectively u ′ , w ′ Perturbation velocity in the x and z directions, respectively b