Evolutionary understanding of airfoil lift

This review attempts to elucidate the physical origin of aerodynamic lift of an airfoil using simple formulations and notations, particularly focusing on the critical effect of the fluid viscosity. The evolutionary development of the lift problem of a flat-plate airfoil is reviewed as a canonical case from the classical inviscid circulation theory to the viscous-flow model. In particular, the physical aspects of the analytical expressions for the lift coefficient of the plate-plate airfoil are discussed, including Newton’s sine-squared law, Rayleigh’s lift formula, thin-airfoil theory and viscous-flow lift formula. The vortex-force theory is described to provide a solid foundation for consistent treatment of lift, form drag, Kutta condition, and downwash. The formation of the circulation and generation of lift are discussed based on numerical simulations of a viscous starting flow over an airfoil, and the evolution of the flow topology near the trailing edge is well correlated with the realization of the Kutta condition. The presented contents are valuable for the pedagogical purposes in aerodynamics and fluid mechanics.


Introduction
The recurring questions on how aerodynamic lift is generated might have arisen when people wonder how birds and bats could fly effortlessly. Since modern aviation has become relatively mature, people might think that how lift is generated seems such a trivial question that they could find a standard answer by just searching on the Internet. Unfortunately, people who are interested in flight are still misled by some misconceptions and even wrong "theories" in non-technical literature. This is evidenced by the widely popularized myth that the laws of aerodynamics prove that the bumblebee cannot fly [1]. There are various popular explanations for lift generation [2]. These explanations capture certain physical aspects of lift generation at different levels of fidelity, but they fail to reconstruct a complete and consistent picture with all the main physical processes. Therefore, it is still required to elucidate the physical mechanisms of lift generation.
The aerodynamic force (lift and drag) of an airfoil is generated as a result of interaction between the incoming flow and airfoil. To understand the physical mechanisms of lift generation, the phenomenological aspects of the flow over an airfoil should be described based on flow visualizations and computational fluid dynamics (CFD) simulations. Figure 1 illustrates the flow over an airfoil with streamlines, marker lines and aerodynamic force vectors. Streamlines are significantly influenced by the airfoil, and the cross sections of streamline tubes are pinched. The cross sectional area of streamline tubes on the upper surface is reduced near the leading edge and enlarged downstream. A direct evidence on the faster motion of fluid particles can be obtained in flow visualizations. As illustrated in Fig. 1, the flow field around the airfoil is visualized by vertical maker lines (oil-smoke lines or dye lines in experiments) generated upstream at sequential times. The two pieces of a marker line separated at the leading edge never meet at the trailing edge, and the marker line on the upper surface moves much faster. It is also observed that the incoming flow is deflected downward in the wake of the airfoil, while the incoming flow near the leading edge is induced upward. Of course, the well-known empirical fact is that surface pressure on the upper surface is lower than that on the lower surface, and the pressure difference between the upper and lower surfaces is the main source of lift. Since the above observed flow phenomena are related to lift, a good explanation of lift should naturally incorporate all these interrelated physical aspects into a unified framework. There are two classes of popular explanations of lift: one based on Bernoulli's equation relating pressure to velocity and another based on Newton's second and third laws applied to the flow momentum deflected downward [2]. Here, a brief review of the popular explanations is given to show why they are incomplete.
Bernoulli's equation along a streamline gives a relation between static pressure p and fluid velocity u, i.e., p + ρ |u| 2 /2 = const.. When velocity at the upper surface is larger than that at the lower surface, the pressure difference between the upper and lower surfaces generates positive lift in the upward direction. However, a critical question is why air on the upper surface moves faster. Bernoulli's equation itself cannot provide an answer. It is tempting to find a more intuitive answer to this question. A naive assumption is that two fluid particles that separate at the leading edge of an airfoil have to meet again at the trailing edge, inferring higher velocity on the upper surface due to the longer traveling distance. However, as illustrated in Fig. 1, this so-called "equal transit time" assumption is not supported by the experimental facts. Therefore, this interpretation may be called the "equal-transit-time" fallacy [2]. An alternative is that higher flow velocity on the upper surface could be achieved by compressing streamline tubes due to the so-called "obstacle effect" or the Venturi effect [3,4]. When the cross section of a streamline tube is pinched, local velocity is increased just like an incompressible flow in a flexible pipe. An interesting historical footnote is that this pipe-flow analogy was once adopted by the physicist Albert Einstein to design an airfoil called "cat's-back" [5]. However, how the compression process of streamline tubes is related to the airfoil geometry is unknown. The Bernoulli-based explanations imply that a velocity field is the cause of the corresponding pressure field. It is conceptually questionable to assume that the velocity field comes first and the pressure field is established as a result. In fact, there is a reciprocal (or circular) cause-and-effect relationship between pressure and velocity [2]. Therefore, the Bernoulli-based explanation is not complete and satisfactory. Another Bernoulli-based explanation uses Bernoulli's equation in the normal direction of a streamline, i.e., ∂p/∂n = ρ |u| 2 /R, where ∂p/∂n is the gradient normal to a streamline and R is the radius of streamline curvature [6,7]. On the upper surface with the larger curvature, the normal pressure gradient is larger, and thus pressure decreases more rapidly when approaching the surface. Since the upper surface is more curved than the lower surface, if streamlines are attached on the surface, it is inferred that pressure on the upper surface is lower than that on the lower surface and as a result lift is generated. Although this conjecture seems conceptually correct, several questions remain unclear. First, it is not justified why the flows are attached on the upper surface with a large curvature (someone may consider the Coanda effect as an answer to this question). Further, to calculate the surface pressure distribution of an airfoil, the integration along a path normal to streamlines depends on the whole velocity field. Therefore, the familiar reciprocal cause-and-effect question remains: why the velocity field comes first and the pressure field is established as a result.
In the momentum-based explanation, it is argued that an airfoil acts like a turning vane deflecting the incoming flow downward, and thus lift is generated due to the deflection of the flow momentum flux according to Newton's second and third laws. However, a question is what the physical mechanism is for the deflection of the fluid stream. On the lower surface at the positive angle of attack (AoA), the fluid flow could be pushed down by the inclined surface. To explain the flow turning on the upper surface, the Coanda effect was proposed as a physical mechanism, originally describing the tendency of a powered jet flow to attach to the adjacent solid surface [7,8]. The application of the Coanda effect to this case where no jet exists could be problematic. In particular, the two different unrelated mechanisms are proposed for the downward deflection of the fluid stream on the upper and lower surfaces, which fails to provide a unified consistent physical explanation [2]. In fact, the deflection of the fluid stream in the momentum-based explanation is intrinsically associated with a vortex system generated by an airfoil/wing [9,10].
On the other hand, the mathematical theory of lift must be quantitative and predictive, unlike the popular explanations of lift. In classical aerodynamics textbooks, the two-dimensional inviscid potential flow theory of airfoils is developed, in which lift is calculated by using the Kutta-Joukowski theorem (the K-J theorem) and the Kutta condition is applied to the airfoil trailing edge to determine the airfoil circulation [11][12][13][14][15][16][17][18]. The potential-flow theory has generated some neat and insightful results in some simple lift problems. The potential-flow methods have been widely used to design various airfoils before CFD methods become available in the era of computers. Interestingly, the effects of the fluid viscosity have been discussed mainly in the drag problem. However, the critical role of the fluid viscosity in lift generation is not sufficiently addressed, which could be attributed to the very subtle nature of lift itself.
In this paper, to clarify the physical mechanisms of lift generation as a viscous-flow phenomenon, we consider a flat-plate airfoil as an example to review the evolution of our understanding of lift generation. Although a flat-plate airfoil is the simplest from a perspective of geometry, the flow over it is non-trivial, exhibiting a singular behavior with very large velocity and low pressure at the sharp leading edge. Thus, the popular explanations have some difficulties for a flat-plate airfoil. For example, a flat-plate airfoil is a clear counterexample for the "equal transit time theory". The Venturi effect of pinched streamline tubes for higher velocity is not solidly grounded on most upper surface of the plate (except near the leading edge) since the cross section of streamline tubes is expanded there. The application of the Coanda effect is questionable at the sharp leading edge where the flow exhibits a singular behavior. Furthermore, the explanation based on Bernoulli's equation along the normal direction to a streamline has a difficulty since both the upper and lower surfaces have zero curvature and a singularity with |∂p/∂n| → ∞ as R → 0 occurs at the leading edge.
The objective of this paper is to elucidate the physical origin of aerodynamic lift of an airfoil using simple formulations and notations, particularly focusing on the critical effect of the fluid viscosity. The content will be valuable as a supplementary material for teaching in aerodynamics and fluid mechanics. From a theoretical standpoint, the analytical expressions for the lift coefficient of a flat-plate airfoil are available for comparison, and the underlying assumptions of these lift models are examined. Newton's sinesquared law is discussed as the direct application of the Newtonian mechanics of solid particles, and then Rayleigh's lift formula is described based on the potential-flow solution for a prescribed flow pattern over a flat-plate airfoil. Thin-airfoil theory is discussed as a classical model of the inviscid circulation theory of lift. Further, an expression for the aerodynamic force of a flat-plate airfoil in an incompressible viscous flow is given, which explicitly elucidates the critical role of the fluid viscosity in generating aerodynamic force (lift particularly). To explore the origin of lift from a general perspective, the viscous vortex-force theory is briefly described, in which the relevant concepts such as the circulation, downwash, form drag and Kutta condition are naturally incorporated. The vorticity generation at the surface is discussed based on Lighthill's viscous coupling equations. The formation of the circulation and generation of lift in a starting viscous flow over an airfoil are described, elucidating the relationship between the trailing-edge flow topology and the Kutta condition.

Newton's sine-squared law
In the direct application of classical Newtonian mechanics of solid particles, fluid is modeled as a collection of individual non-interacting solid particles that impinge on a surface and then slide frictionlessly along the surface after collisions. Thus, the momentum transferred from the impinging particles to the surface is in the normal direction to the surface. As shown in Fig. 2a, we consider a flat-plate airfoil with the area S (a product of the chord length and the unit span) and angle of attack (AoA) α at the freestream velocity U ∞ . The momentum flux impinging onto the plate is U ∞ (ρU ∞ S sin α). According to Newton's second law, the normal force per unit span to the plate is F N ¼ ρU 2 ∞ S sin 2 α, and the normal force coefficient is where q ∞ ¼ ρU 2 ∞ =2 is the dynamic pressure, and ρ is the fluid density. The sectional lift and drag coefficients are expressed as where L' and D' are the sectional lift and drag, respectively. The lift-to-drag ratio (L/D) is C l /C d = cot α. It is noted that Newton first found that the force of a segment of a curved surface was proportional to the squared sine function of the incidence angle [19]. Figure 3a shows C l as a function of α for Newton's sine-squared law. Newton's sinesquared law indicates that the lift coefficient at small AoA is so small that lift is not sufficient enough to air vehicle in flight. In order to increase lift, the wing area S must increase proportionally such that the wing would be too large for practical flight. Furthermore, when the AoA is large for increasing lift, as shown in Fig. 3b, the drag coefficient increases significantly such that L/D becomes too small for practical flight. Therefore, Newton's sine-squared law gave a pessimistic prediction for aircraft flight before modern aerodynamic was formulated.
Essentially, Newton's sine-squared law is one of quantitative forms of the momentum-based explanation. The contribution to lift and drag is solely from the momentum exchange between impinging particles and the lower surface of the plate. On the upper surface, there is no particle interacting with the surface. Therefore, the contribution to lift and drag from the upper surface is totally neglected, which is now known as the major part of lift. This explains why Newton's sine-squared law significantly under-estimate both lift and drag, as shown in Fig. 3. From a physical viewpoint, the main problem in Newton's sine-squared law is that fluid cannot be simply treated as a collection of individual non-interacting solid particles. Newton's second law should be applied to fluid as a continuum, which leads to the momentum equation in fluid mechanics, i.e., the Navier-Stokes equations (the NS equations). In principle, lift and drag of an airfoil in an incompressible viscous flow can be calculated by solving the NS equations coupled with the continuity equation. However, due to the tremendous mathematical difficulty in solving the NS equations for the lift problem, various approximations and assumptions have to be made at different levels of fidelity to develop predictive analytical aerodynamics models. Some physical mechanisms, particularly the viscous effect on lift generation lost in simplified aerodynamics models, are not sufficiently elucidated, which could lead to misunderstandings on lift generation.

Rayleign's lift formula
It has been recognized that the working lift models have to be developed based on a cascade of suitable approximations and assumptions since the NS equations are difficult to be directly solved analytically. First, the inviscid-flow approximation is made such that the NS equations are reduced to the Euler equations that have an integral: Bernoulli's equation providing an explicit relation between velocity and pressure along a streamline. Bernoulli's equation can be used as a useful tool to calculate a pressure distribution on an airfoil and further airfoil lift when a velocity field around an airfoil is given. To reconstruct a velocity field around an airfoil, a further approximation is that the flow is potential and irrotational. Therefore, the flow is governed by the Laplace equation of the velocity potential. Some elemental solutions of the Laplace equation can be used as building blocks due to its linear nature to reconstruct a velocity field around an airfoil. In particular, for a two-dimensional (2D) potential flow, a conformal transformation can be applied to the airfoil theory.
Rayleigh [21] applied a conforming transformation (Kirchhoff's transformation) to a flat-plate airfoil at a given AoA to reconstruct a potential flow with two discontinuous lines (free streamlines) originating from the leading and trailing edges. As shown in Fig.  2b, this flow pattern is consistent with the observation of a massively separated flow around a flat plate with a large AoA, where two discontinuous lines represent two shear layers shedding from the leading and trailing edges. The region confined by the two free streamlines behind the plate is the so-called "dead air" or "dead water" region that simulates the wake in a viscous flow. The pressure difference across the plate can be calculated using Bernoulli's equation. Rayleigh gave a formula for the normal force coefficient, i.e., The sectional lift and drag coefficients are expressed as The lift-to-drag ratio is C l /C d = cot α. Figure 3a shows C l as a function of α for Rayleigh's lift formula. Similar to Newton's sine-squared law, Rayleigh's lift formula still predicts a smaller value of C l than experimental and CFD data of the lift coefficient of a flat-plate airfoil. Figure 3b shows the predicted drag coefficient is more consistent with the CFD data. This underestimated lift results from the fact that Rayleigh's flow pattern is suitable for the drag problem rather than the lift problem. Particularly, Rayleigh's analysis did not incorporate the vortex-lift mechanism that is a key in airfoil lift generation. Nevertheless, Rayleigh's work had a significant impact on the application of conforming transformations to various inviscid flows [22]. Along this line, the Joukowski transformation and its extended forms were later developed to calculate airfoil lift and design airfoils [12][13][14].

Circulation and thin-airfoil theory
The study of the vorticity, vortex and circulation has a long history in theoretical hydrodynamics [22,23]. However, the intrinsic connection between these apparently abstract mathematical concepts and lift generation for flight had not been explored until Lanchester recognized the importance of vortices in lift generation [24,25]. Then, the advent of the Kutta-Joukowski theorem (the K-J theorem) had formally laid a rational foundation for the circulatory lift theory [26,27]. For a 2D flow over an airfoil, the K-J theorem gives lift per unit span, i.e., In Eq. (5), Γ is the circulation defined as where u is the fluid velocity around the airfoil, C is an arbitrary closed contour enclosing the airfoil, the contour integral is clockwise as shown in Fig. 1, and the vorticity ω is clockwise. Physically, the circulation is the measure of fluid rotation around an airfoil, which is positive when lift is an upward force (the positive lift). In the second equality of Eq. (6), according to the Stokes theorem, the line integral of the velocity u along C equals to the surface integral of the spanwise vorticity ω in a 2D domain V C enclosed by C, where ω is one component of the vorticity ω = ∇ × u. In fact, the K-J theorem is a compact expression of the total momentum flux across a far-field control surface enclosing an airfoil. In this perspective, as illustrated in Fig. 1, the circulation is directly proportional to the downwash momentum flux through the boundary of a domain and therefore lift. For a 2D flow, a direct relation between the circulation and downwash could be found [27]. The Kutta-Joukowski theory is referred to as the circulation theory of lift, bringing the vorticity and vortex dynamics to the center of classical aerodynamics.
Thin-airfoil theory developed by Munk [15] is a classical application of the circulation theory of lift. From a standpoint of aerodynamics, a thin airfoil is modeled as a straight vortex sheet with the strength distribution γ(x) on the chordline from the leading edge to trailing edge. The vortex-sheet strength is given as γðxÞ ≈ ½ u e þ ¼ u þ e −u − e , representing the velocity difference between the upper and lower surfaces denoted by the superscripts + and −, respectively. From this perspective, the velocity difference between the upper and lower surfaces of an airfoil is a natural result of the vortex sheet that is an idealized model of the boundary layers (see Section 6). Figure 2c shows the flow over a vortex sheet that simulates a flat-plate airfoil. Generally, this flow pattern represents a more reasonable model for a thin airfoil.
To determine γ(x), the thin-airfoil integral equation can be derived based on the nonpenetration condition on the vortex sheet [14,15]. A Fourier-series solution of the thin-airfoil equation for γ(x) can be obtained, where the Kutta condition γ = 0 at the trailing edge is imposed (the velocity or pressure difference at the trailing edge is zero). The circulation Γ is calculated by integrating γ(x) from the leading edge to trailing edge, and the sectional lift is calculated using the K-J theorem. Thin-airfoil theory gives the sectional lift coefficient where α L = 0 is the AoA at zero lift that represents the effect of the airfoil camber. Thinairfoil theory predicts the lift slope of dC l /dα = 2π, as shown in Fig. 3a. Compared to Newton's sine-squared law and Rayleigh's lift formula, thin-airfoil theory gives the lift coefficient that is more consistent with the CFD and experimental data. Although both Rayleigh's theory and thin-airfoil theory are based on the potential-flow theory, the success of thin-airfoil theory in predicting airfoil lift is a result of introducing vortex (or vortices) leading to vortex lift as a correct physical mechanism of lift generation. As shown in Fig. 4, at large Reynolds numbers, thin-airfoil theory gives an improved prediction of the lift coefficient for thin airfoils with round leading edges that avoid leading-edge separation. However, when the Reynolds number based on the chord length c decreases to Re c~1 0 2 − 10 3 , the measured lift curve usually exhibits a nonlinear behavior and the lift slope becomes smaller, as shown in Fig. 3a.
The circulation in the K-J theorem cannot be determined in the potential-flow framework alone since the vorticity cannot be physically generated in an inviscid flow. Due to D′Alembert's paradox, the integrated pressure force of a body is zero in a steady inviscid irrotational incompressible flow [12,14]. The important implication of D′Alembert's paradox is that lift and drag of a moving body must be a result of a viscous flow no matter how small the fluid viscosity is. To determine the circulation, the Kutta condition is imposed at the trailing edge, providing a phenomenological model of the viscous effect on lift generation [13,14]. The viscous origin of lift has been recognized, which has to be essentially found in the viscous-flow framework (see Section 6). It is Fig. 4 The sectional lift coefficient as a function of AoA for several NACA airfoils noted that thin-airfoil theory gives zero drag as a consequence of D′Alembert's paradox although it correctly predicts lift.

Viscous-flow lift formula
Lift of a flat-plate airfoil can be directly evaluated as a viscous-flow result without explicitly using the circulation and K-J theorem in classical aerodynamics [20]. This nontraditional attempt is based on the intrinsic relationship between the skin friction topology and surface pressure structure in a viscous flow. For an incompressible viscous flow over a stationary curved surface, an exact relation between the surface pressure gradient ∇p and the skin-friction vector τ can be derived from the NS equations [20,28,29,30,31]. In particular, for a flat surface where the curvature effect vanishes, this on-wall relation is written as where f Ω = μ ∂Ω/∂n is the boundary enstrophy flux (BEF), μ is the dynamic viscosity of fluid, Ω = |ω| 2 /2 is the enstrophy, ω is the vorticity, and ∂/∂n is the derivative along the unit normal outward vector n of the surface. The BEF is an intriguing quantity that is particularly related to the topological features in complex flows such as isolated critical points, separation lines, and attachment lines in a skin-friction field [28][29][30]. Along a skin-friction line, Eq. (8) is re-written as dp/ds = s ⋅ ∇ p = μ f Ω |τ| −1 , where s = τ/|τ| is the unit vector along the skin-friction line. Therefore, the surface pressure is given by the path integral along a skin-friction line, i.e., where P(x) and P 0 (x 0 ) denote a point and a starting point on a skin-friction line, respectively. In principle, since skin-friction lines are densely distributed, a set of the starting points P 0 (x 0 ) from which skin-friction lines originate could be selected such that the points P(x) can cover densely the whole surface. According to Eq. (9), surface pressure is related to the historical effect of the viscous diffusion flux of the boundary vorticity magnitude (μ ∂|ω|/∂n) along a skin-friction line. At the limit of μ → 0, the boundary layer thickness approaches to zero, i.e., δ → 0 while |ω| remains finite, leading to ∂|ω|/∂n → ∞. Thus, μ ∂|ω|/∂n could be finite. It is noted that the viscous flow with μ → 0 is not equivalent to the purely inviscid flow with μ = 0. The aerodynamic force is formally expressed as the following surface integral Eq. (10) explicitly describes the critical role of the fluid viscosity in generating the aerodynamic force. For a steady inviscid incompressible flow with μ = 0, we have F = 0 according to Eq. (10), and thus D'Alembert's paradox is naturally recovered. The first term in Eq. (10) is the integrated pressure force expressed by the integral of the historical effect of the viscous diffusion flux of the boundary vorticity magnitude that is explicitly related to the fluid viscosity. The second term is the integrated skin-friction force. Lift is mainly contributed by the pressure term (the first term). The main consequence of Eq. (10) is that lift and drag (including the pressure and skin-friction drags) must coexist as a result of the viscous flow over an airfoil. Physically, lift cannot be generated without the cost of generating the viscous drag at the same time.
To elucidate the connection between boundary layer and lift generation through Eq. (10), a flat-plate airfoil at a small AoA (α) is considered. As shown in Fig. 2d, the boundary layers with opposite vorticity are developed on the upper and lower surfaces of the plate. To model the boundary layers of the flat plate, the Falkner-Skan flow is used, in which the BEF and skin friction can be expressed in the analytical forms [32]. Therefore, by using Eq. (9), the pressure coefficient difference across the plate (the non-dimensional pressure loading of the plate) is given by where Δ p is the pressure difference between the lower and upper surfaces of the plate, x ¼ x=c is the chordwise coordinate normalized by the chord length c, is the ratio between the trailing-edge and freestream dynamic pressures, U ∞ is the freestream velocity, U ref ¼ a 1 ð ν=U ∞ Þ m 1 is a reference velocity, Re c = U ∞ c/ν is the Reynolds number, Δ C p, TE is the value of Δ C p at the trailing edge, and . The Δ C p -distribution is weakly dependent of Re c since |m 1 | ≈ α/ π < < 1. At the leading edge, as x→0, there is a singularity of ΔC p →x 2m 1 with m 1 = − α/(π − α), and the singularity is weakened as α → 0, which seems physically reasonable in a viscous flow. In contrast, in the classical thin-airfoil theory, the stronger singularity Δ C p →x −1=2 at the leading edge remains unchanged even as α → 0. Figure 5 shows the Δ C p -distributions on the flat plate predicted by Eq. (11) for different AoAs at Re c = 200. Eq. (11) gives the consistent profiles with the CFD data [20]. Interestingly, the profiles of Δ C p /α approximately collapse, indicating that the effect of α on Δ C p predicted by the viscous-flow model is essentially consistent with the distri- given by thin-airfoil theory. The value of Δ C p, TE decreases with α, and the classical version of the Kutta condition ΔC p = 0 is not exactly Fig. 5 The chordwise distributions of the surface pressure coefficient difference across the flat plate at different AoAs. From Liu et al. [20] satisfied at the trailing edge. Nevertheless, the Taylor-Sears condition that the net vorticity flux is zero at the trailing edge holds in this viscous flow, which is a generalized form of the Kutta condition (see Section 6). In this sense, the Kutta condition is a reduced and more restricted case of the Taylor-Sears condition.
The normal force of the flat plate is calculated by integrating Δ p from x = 0 and x = c. Then, the sectional lift and drag coefficients are expressed as where the non-linear factor is defined as and C d, 0 is the parasite drag coefficient (the value of C d at zero AoA). The lift formula Eq. (12a) is different from the linear relation C l = 2π α given by thin-airfoil theory due to the non-linear factor F(α) depending on the parameters R q and Δ C p, TE . The parameter R q is determined by the velocity ratio U ref /U ∞ and the Reynolds number Re c . In the limiting case as α → 0, the asymptotic behavior of F(α) → 2(U ref /U ∞ ) 2 /π 2 = 1 is inferred when Eq. (12a) is consistent with thin-airfoil theory with Δ C p, TE = 0 in this limiting case. Therefore, the parameter Figure 3a shows the sectional lift coefficient as a function of α. The C l -curve predicted by Eq. (12a) exhibits the non-linear behavior that is consistent with the CFD results at Re c = 200 [20] and measurement data of an aluminum (Al) foil wing with the aspect ratio of 6 at Re c = 7500 [33]. In contrast, thin-airfoil theory gives the linear C lcurve with the larger lift slope. This low-Reynolds-number flow example indicates that the nonlinear effect of the viscosity on lift generation could not be sufficiently simulated by the Kutta condition imposed in a simple linear potential-flow model. Furthermore, as shown in Fig. 3b, the C d -curve given by Eq. (12b) is in reasonable agreement with the data given by CFD and measurements [20,33], and C d − C d, 0 increases with α mainly due to the pressure drag.

Viscous Vortex force
To provide a feasible physical explanation of lift generation, the circulation or vortexforce theory has to be reconstructed in the viscous-flow framework, and the critical role of the fluid viscosity in lift generation must be sufficiently addressed by incorporating all the relevant physical features into the theory. The general force expressions have been derived from the Navier-Stokes equations [34][35][36][37]. When the Reynolds number is large (Re > > 1), the force of a body in the fluid control volume V f bounded externally by a control surface Σ is given by Wu et al. [38], i.e., where P = p + ρ |u| 2 /2 is the total pressure with P ∞ being its constant value at upstream infinity. The volume integral of the Lamb vector (l = u × ω) in Eq. (14) represents the vortex force, which was first found by Prandtl in developing the finite wing theory [39].
In an inviscid flow with P ∞ − P = 0 (Bernoulli's equation), the K-J theorem can be exactly reduced from the vortex force in 2D [40,41]. As illustrated in Fig. 6, the boundary layer on an airfoil is a key element for lift generation since the Lamb vector as the main force constituent is concentrated in a very thin layer in the boundary layer on the surface (within about 10% boundary-layer thickness). In Eq. (14), since the total pressure loss measured by P ∞ − P > 0 appears in the viscous wake, the integral of P ∞ − P over the wake plane represents the viscous form drag denoted by D form . On the other hand, the component of the vortex force in the freestream direction is the induced drag denoted by D in that coexists with the lift of a finite wing. Therefore, all the force constituents, including lift, skin-friction drag, form drag, and induced drag, are related to the vorticity generated at the surface in a 3D viscous flow.

Kutta-Joukowski lift and form drag
In a 2D steady viscous flow, the vortical wake must extend downstream unboundedly, and any contour surrounding an airfoil must cut through the wake. In this case, since Bernoulli's equation no longer holds across the viscous wake, the original derivation of Eq. (5) by Joukowski in an inviscid flow is not strictly applicable. The applicability of the K-J theorem in a viscous flow was first studied experimentally by Bryant and Williams [42]. To explain the experimental findings, Taylor [43] provided a theoretical account. The main results of Taylor are stated as follows. If the downstream face of the outer contour Σ is a wake plane denoted by of W, for Re > > 1, the sectional lift and form drag are given by where Γ is the circulation along the outer contour Σ. The condition imposed on Eq. (15) is that at the wake plane, the net vorticity flux must vanish, i.e., Z where u is the velocity projected on the freestream direction. According to Eq. (15), the K-J theorem naturally coexists with the form drag formula in the viscous-flow framework, which resolves the dilemma in the classical inviscid circulation theory associated with D'Alembert's paradox.
As illustrated in Fig. 6, Eq. (16) indicates that the positive and negative advective vorticity fluxes (uω) from the boundary layers on the upper and lower surfaces are cancelled out in the wake. Therefore, Eq. (16) is necessary for the circulation Γ to be independent of the position of the wake plane. Sears [44,45] further proved that Eq. (16) would be equivalent to the requirement that the pressures at the outer edges of the boundary layers on the upper and lower surfaces must be the same at the trailing edge. Eq. (16) is referred to as the Taylor-Sears condition providing a viscous-flowtheoretical foundation for the empirical Kutta condition.

Thin airfoil as a reduced case
To elucidate the effectively viscous origin of lift, as illustrated in Fig. 7, thin-airfoil theory is a simplified model, where the boundary layer on a thin airfoil can be idealized as a vortex sheet on the camber line [46]. As indicated in Eq. (14), the Lamb vector (l = u × ω) is the sole contributor to lift. The Lamb vector is concentrated in the boundarylayer domain V bl that forms a folded band with the width of δ and its outer contour C AB wrapped around the flat plate from the point A to the point B in the wake as Fig. 7 Idealization from a boundary layer to a vortex sheet to a bound vortex shown in Fig. 8. The vortex lift of a body in a viscous flow is expressed as a volume integral of the Lamb vector, i.e., where e l is the unit vector normal to the freestream velocity. Formally, the vortex lift of an airfoil is expressed in the K-J theorem, i.e., L ' = ρ U ∞ Γ, where Γ is the circulation based on the Lamb vector integral. Further, since the Lamb vector in a 2D boundary layer is u × ω = n uω, the sectional lift of a flat-plate airfoil calculated based on pressure at the boundary-layer edge (the contour C AB ) is related to the vortex lift by where χ ¼ ½uω δ 0 is the advective vorticity flux across the boundary layer and ½χ A B is the jump of χ across the points A and B in the wake.
Eq. (18) indicates that in the viscous flow the pressure lift equals the vortex lift given by the K-J theorem when the condition ½χ A B ¼ ½uω A B ¼ 0 is satisfied in the wake. In other words, the positive and negative advective vorticity fluxes from the boundary layers on the upper and lower surfaces should be cancelled out in the wake for the K-J theorem to be applicable to a viscous flow. This condition is just the Taylor-Sears condition. For a very thin boundary layer, as the wake plane approaches to the trailing edge, the Taylor-Sears condition would be reduced to the requirement that pressures at the outer edges of boundary layers of the upper and lower surfaces must be the same at the trailing edge, which is the Kutta condition ΔC p = 0 at the trailing edge.
To examine the Taylor-Sears condition in the flow over a flat-plate airfoil, the difference of the Lamb vector integrals across the boundary layers on the upper and lower surfaces, i.e., Δl ¼ ½uω þ − , is evaluated based on numerical simulations [20], where "+" and "−" denote the upper and lower surfaces, respectively. The Lamb vector difference Δl can be interpreted as the local loading on the flat plate, and at the same time Δl ¼ ½uω þ − represents the net advective vorticity flux across the boundary layers on the flat plate. Figure 9 shows the chordwise distributions of 2Δl Ã =α ¼ 2Δl=U 2 ∞ α on the flat plate at different AoAs in comparison with Δ C p /α given by thin-airfoil theory. It is found that Δl ¼ ½uω þ − ¼ 0 at the trailing edge, and therefore the Taylor-Sears condition holds in this viscous flow even though the Kutta condition ΔC p = 0 is not satisfied at the trailing edge (see Section 5).  [28] In general, when a thin airfoil is divided into N segments with the length of Δx k (k = 1, 2, ⋯N), the sectional lift of the airfoil at a small AoA can be approximately expressed as where Δp is the surface pressure difference between the lower and upper surfaces, Δl is the difference of the Lamb vector (l = uω) between the upper and lower surfaces, γ = Δp/ρU ∞ is the strength of the vortex sheet, and δ is the characteristic thickness of the Lamb-vector layer in a boundary layer. For the homogeneous segments, a local relation is Δp k ≈ ρ U ∞ γ k ≈ ρ Δl k δ k , indicating that the pressure loading of a segment of a thin airfoil approximately equals the local K-J lift or the local vortex lift of the segment. This directly elucidates the relationship between the pressure difference across the airfoil and the boundary-layer properties (vorticity and Lamb vector). In other words, lift is generated by the surface pressure difference that is intrinsically coupled with the boundary-layer properties (vorticity and Lamb vector) on the upper and lower surfaces and strained by the Taylor-Sears condition at the trailing edge. This provides the viscous-flow foundation of thin-airfoil theory. In a similar approach, unsteady thinairfoil theory can be developed in the viscous-flow framework [46], in which the classical von Kárman-Sears theory of unsteady inviscid flow over a moving thin airfoil is recovered as a reduced case [47].

Lift, circulation and downwash
The circulation theory and the momentum-based models are related. According to Newton's third law, the upward lift must be balanced by the downwash of a massive fluid body pushed by the wing. The momentum-based explanation claims that the origin of lift should be the downward turning of the flow or downwash. As pointed out by McLean [2], however, lift and downwash are in a reciprocal causal relationship. It is Fig. 9 The normalized chordwise distributions of the Lamb vector integral across the boundary layers on the flat-plate airfoil. From Liu et al. [20] therefore of interest to examine how lift and downwash are related. The force exerted by a wing/airfoil is proportional to the rate of change of the total momentum through a sufficiently large control surface Σ enclosing the wing/airfoil. When the viscous effect is neglected in the far flow field, the vertical projection of the force gives a liftdownwash relation, where n z and v are the components of the unit vector n and velocity u = (u, v) in the zdirection normal to the incoming flow, respectively. The first and second terms in the right-hand side (RHS) of Eq. (20) are the contributions by pressure and the downwash momentum flux projected in the vertical direction through Σ, respectively. For a 2D flow, a direct relation between the circulation and downwash was found by Wu et al. [38]. A sufficiently large rectangular control surface Σ enclosing an airfoil is considered. The front face of Σ is at far upstream with the freestream velocity U ∞ , and the rear face (a wake plane W) is at x W downstream of the trailing edge. The origin of the coordinate x is set at the leading edge. In the non-dimensional variables u ¼ u=U ∞ , x = x/c and z = z/c with c being the chord length, this relation is expressed as The first term in the RHS of Eq. (21) is the downwash momentum flux on the wake plane. The second term is the difference of u 2 −v 2 between the upper and lower faces of Σ that vanishes as these faces recede to infinity. It is clear that both the circulation theory and the lift-downwash argument are on the same basis.

Generation of vorticity and circulation
The standard argument on the origin of the circulation was given based on flow visualizations in a starting flow [11,15,48]. A suddenly accelerating airfoil generates a starting vortex observed in experiments. According to the total circulation conservation theorem, since the total circulation in the region enclosed by a sufficiently large control contour is zero initially, it is inferred that a bound vortex with the opposite sign and the same magnitude of the starting vortex should be generated in the airfoil at the same time. However, this argument is more like logical inference than direct experimental validation [10]. The true cause of the circulation is the boundary-layer vorticity on an airfoil. The boundary-layer theory actually provides the physical foundations for the circulation theory. The origin of the vorticity was studied by Lighthill [49,50]. To recapitulate his result, we use the right-hand orthonormal frame (n, t) on a streamline, as shown in Fig. 1, where t and n are the unit tangent and normal vectors of any streamline. Lighthill gave a pair of equations describing the viscous coupling between the surface pressure gradient and boundary vorticity flux, i.e., Eq. (22a) reveals that the tangential pressure gradient can generate new vorticity at the rate measured by the boundary vorticity flux (BVF) denoted by σ = μ ∂ω/∂n. Figure 10 illustrates a rolling-fluid-ball analogy for the vorticity-generating mechanism revealed by Eq. (22a). Consider a small fluid ball on a solid wall. It is assumed that at t = 0 the tangential pressure gradient (∂p/∂s < 0) suddenly occurs to push the ball to move toward in the right direction. Due to the no-slip condition, the ball cannot slide but only rolls along the wall and hence gains an angular velocity. Thus, the vorticity in the rolling ball is newly created by the tangential pressure gradient, which first occurs in the fluid layer adjacent to the wall and then diffuses into the fluid. In contrast, the second equation in Eq. (22b) has a very minor effect since its terms in both the sides are in the order of Re −1/2 when the Reynolds number Re is large. It is possible to identify the causality mechanisms embedded in Eq. (22a). When an airfoil moves in fluid, pressure p is generated as the first causal mechanism since it is a result of molecules collisions in an equilibrium state. This is an inviscid process with the timescale t p~c /a, where c is the chord length and a is the speed of sound. In an incompressible flow, t p → 0. Since it is realized that pressure alone can by no means fully explain lift generation, the appearance of the viscous shear stress is a necessary mechanism, which comes after a short relaxation time [49,50]. Strictly, according to Eq. (22a), the tangential pressure gradient causes the BVF, but they occur almost simultaneously since the molecule relaxation time can be neglected in the local equilibrium assumption used throughout the continuous fluid dynamics. Then, the vorticity generated by the BVF at the wall diffuses into the fluid and advects downstream by the tangential pressure gradient. The advection time scale is t a~c /U ∞ so that t p /t a~M , where M is the Mach number. The diffusion time scale is t ω~δ 2 /ν, where δ is the effective viscous diffusion distance. If the estimated boundary-layer thickness is δ $ c Re −n c , there is an estimate t p =t ω $ M Re 2n−1 c , where Re c = U ∞ c/ν is the Reynolds number and n is a positive exponent that depends on the boundary-layer characteristics. For an incompressible flow with M → 0, there must be t p /t ω < < 1. Therefore, the generation of the vorticity in a boundary layer is an accumulated effect of the BVF in both space and time.
To study the generation of the circulation and the realization of the Kutta condition, numerical simulations were conducted by Zhu et al. [51] on a laminar accelerating uniform incoming incompressible flow over a NACA-0012 airfoil at AoA of 6 o . A scaling Fig. 10 Schematic illustration of vorticity generation at a wall. Adapted from Lighthill [49,50] analysis gives a generic transient process of aircraft takeoff. In the aircraft coordinate system, takeoff can be considered as the transient flow described by the normalized velocity U(t)/U ref as a function of the normalized time t/T ref , where the reference quantities U ref and T ref are the takeoff velocity and time, respectively, at which an aircraft lifts off from the ground. Figure 11 shows the normalized velocity data in takeoff for a Cessna 182 [52], a hybrid flight vehicle [53] and a Gulfstream American AA-5B Tiger [54]. Interestingly, these data collapse well to the model function Zhu et al. [51] in their numerical simulations. The solution of the dynamical equation for takeoff of a Cessna 182 is also plotted in Fig. 11 for comparison.
The major events in the formation process of the circulation and the flow patterns near the trailing edge are shown in Fig. 12. Streamlines reveal the evolution of the flow topology particularly in the zoomed-in views near the trailing edge. Isolated critical points are marked in Fig. 12, where node and saddle are denoted by N and S, and seminode and semi-saddle by N' and S', respectively. The front stagnation point at the leading edge is a semi-saddle. The number of critical points in a 2D airfoil flow must obey the following topological rule [55,56].
where # denotes the total number of nodes, saddles, semi-nodes or semi-saddles. The number of critical points at different stages in Fig. 12 satisfies Eq. (23), confirming the topological consistency of the computed flow patterns. Figure 13 shows the generation of the circulation in the transient process at the corresponding moments.
To describe the flow patterns near the trailing edge, several versions of the Kutta condition are given. From a phenomenological standpoint, Glauert [15] stated "the flow must leave the trailing edge smoothly"; von Kármán and Burgers [40] stated "in the final steady flow, the rear stagnation point shall coincide with the trailing edge of the airfoil". The physical foundation of the Kutta condition can only be found in the viscous-flow framework. Sears [45] gave a quantitative condition that the total vorticity flux would be zero at the trailing edge (the Taylor-Sears condition), which can be reduced to the requirement that the surface pressure difference at the trailing edge vanishes. The statements by Glauert, Kármán and Burgers, and Sears are referred to as the Glauert, Kármán-Burgers, and Sears versions of the Kutta condition, respectively.
The following stages (a)-(e) in the process of the generation of the circulation and lift are summarized below.
(a) In a very short time period immediately after the flow starts up, as shown in Fig. 12a at t/T ref = 0.002, the pressure field and the tangential pressure gradient on the airfoil surface are instantly established, and thus the boundary vorticity flux (BVF) is generated almost at the same time through the pressure-vorticity coupling at the surface according to Lighthill's local dynamic equations. However, in such a short time period, the generated boundary vorticity is not yet diffused into the fluid since the diffusion timescale is much larger than the molecule relaxation time. Therefore, the flow field over the airfoil exhibits the typical topology of the inviscid irrotational flow, in which a semi-saddle (the rear stagnation point) occurs on the upper surface. At this moment, the circulation and lift are zero. (b) As the generated boundary vorticity diffuses into the vicinity of the wall, a near-wall separated bubble starts to form and grows rapidly in its size on the upper surface. After the bubble grows to a critical size, as shown in Fig. 12b at t/T ref = 0.029, two semi-saddles on the upper surface and upstream the trailing edge merge into a single saddle that starts to lift off from the wall, so the bubble breaks and the upper-surface transient boundary layer becomes fully attached, leaving a single semi-saddle at the trailing edge. The small segment of free vortex sheet starts to roll rapidly, which signifies the initial formation of a starting vortex of counterclockwise vorticity. Thus, the clockwise circulation around the airfoil and lift starts to appear. The von Kármán-Burgers version of the Kutta condition is realized since the rear stagnation point coincides with the trailing edge.   Fig. 12 are marked. Adapted from Zhu et al. [51] (e) As shown in Fig. 12e, after a sufficient time of evolution, the flow field approaches a steady state at t/T ref = 1.0. Since the net advective vorticity flux in the wake is zero near the trailing edge, the Sears version of the Kutta condition is realized. Therefore, the legitimacy of the application the K-J theorem to viscous attached flows is fully established. Eventually, steady lift is generated along with the airfoil circulation.
The evolution of the near-wall flow topological structures near the trailing edge in the above process is possible only in a viscous flow, which is directly correlated with the realization of the Kutta condition. In the history of airfoil's flow evolution, the pressure field is first established around the airfoil in a very short time before the viscous effect becomes significant, and thus the inviscid flow pattern appears. Then, according to Lighthill's coupling equations, the vorticity is generated by the tangential pressure gradient near the wall and the viscous shear layer is developed. Eventually, interaction between the pressure field and the near-wall viscous flow leads to a steady flow field over the airfoil where the flow leaves the trailing edge smoothly and the corresponding pressure field. Therefore, the steady circulation and lift are generated.

Conclusions
The generation of lift happens only in a viscous flow over an airfoil. The incompressible aerodynamic force is solely related to the vorticity field and its associated velocity field. In this sense, the boundary layer on the airfoil surface is a direct physical origin of lift. This point is elucidated by reviewing the evolutionary development of the lift problem of a flatplate airfoil, including Newton's sine-squared law, Rayleigh's lift formula, thin-airfoil theory, and viscous-flow lift formula. The Kutta-Joukowski inviscid circulation theory for lift, including the Kutta condition for determining the circulation, is rationalized in the viscousflow framework. The viscous vortex-force theory provides more consistent interpretations for the relevant concepts such as the circulation, form drag, Kutta condition, and downwash. To illustrate the generation of the circulation and lift as a viscous-flow phenomenon, the numerical simulations reveal the dynamic evolution of the trailing-edge flow topology in a starting flow over an airfoil. The flow topological structures near the trailing edge are well correlated with the realization of the different versions of the Kutta condition.
The relevant aspects of lift are summarized below.
(1) Lift is a result of the surface pressure difference between the upper and lower surfaces, which is intrinsically coupled with the boundary-layer development characterized by the relevant physical properties particularly vorticity, Lamb vector, and boundary vorticity magnitude flux.