The ability and rich legacy of potential flow theory have paved the path for aerodynamicists to exploit it in a wide range of aeronautical engineering problems. Although its history goes back to D'Alambert and Euler,1 the theory did not ripen until Prandtl.2 Since his seminal efforts, scholars have sought various analytical and numerical techniques to extend its applicability to capture the essential physics of many complex problems. These potential flow-based approaches (e.g., thin airfoil theory,3 lifting line theory,4 Wagner,5 and Theodorsen models6) typically result in integral equations with singular kernels. As discussed by Belotserkovskii,7 one of the most convenient techniques to solve singular integral equations is via discretization, leading to a linear algebraic system of equations.8 This approach constructs the basis of the so-called vortex lattice method (VLM) introduced by Falkner9—one of the most common efficient solvers in steady-state scenarios. In this approach, the wing is divided into panels (lattice), and a vortex is placed at each panel. The strengths of these bound vortices are determined to satisfy the no-penetration boundary condition, resulting in a linear system of equations. In the case of transient problems, such as those associated with unsteady pitching, plunging, or surging, the unsteady vortex lattice method (UVLM)10 is one of the most convenient tools, which allows for wake deformation and arbitrary time-varying motion. The main difference between the steady (VLM) and unsteady (UVLM) versions is that the latter account for wake effects: vortices are continuously/discretely shed at the sharp edge, and the solver has to track their free motion in the fluid to account for their effect on the no-penetration boundary condition on the airfoil.

Another vortex-based approach for numerical modeling of unsteady aerodynamics is the discrete vortex method (DVM).11–15 Unlike the UVLM, the bound circulation distribution over the airfoil is determined directly without dividing the airfoil into panels, hence avoiding solving a linear system at each time step. The method employs a conformal mapping between the airfoil and a circular cylinder whose potential flow solution is readily known. The main unknown in the DVM at each time instant is the strength of the vortex being shed at the sharp edge, which is typically determined by applying the Kutta condition at the edge. As such, this technique, while avoids the inversion of a linear system, does not allow for arbitrary body shapes (e.g., time-varying camber); a conformal map between a cylinder and the body must be known a priori.

Because of the ability of the potential flow-based methods, such as UVLM or DVM, to capture the essential macroscopic physics of the flow, they are adopted for variety of engineering and research problems. There have been several previous efforts to develop extensions of UVLM and DVM to regions where the primary potential flow assumptions fade out, hence broadening the applicability of the aforementioned methods. For instance, Wang and Eldredge,16 Ford and Babinsky,17 Hemati, Eldredge, and Speyer,18 Ramesh et al.,15 Darakananda and Eldredge,19 SureshBabu et al.,20 Epps et al.21 among others have developed extensions of the UVLM/DVM to high angle of attack maneuvers. Such regimes are commonly reached in modern applications of bio-inspired flight22,23 and dynamic stall.24,25

One of the main concerns about the potential flow framework (either the analytical or numerical version) in two dimensions is its inability to determine a unique value of the circulation around the body26; hence, an auxiliary condition is needed. The Kutta condition27 is widely used as the closure condition, which dictates smooth flow off the sharp trailing edge where separation happens and the wake begins. The Kutta condition was first employed to determine the lift of a shell of infinite span and cross section of a circular arc at a small angle of incidence with respect to the flow.28 Using conformal mapping to transform the flow around the shell into the flow around a cylinder, Kutta observed that the potential velocity at the extremities of the shell was infinitely large. To avoid such a “disturbing fact”—in his words “störende Umstand”—he set the circulation to remove the singularity at one of the edges (the rear edge) and assumed that the geometry in reality would be rounded enough at the other edge (the leading edge) to remove the other singularity. He observed good agreement between the lift calculated using this consideration and Lilienthal's experimental values.27

The Kutta condition is generally accepted as a proper closure condition for steady attached flows at high Reynolds numbers, which may also be applicable to unsteady problems at a relatively small oscillation frequency. However, the discontent with the Kutta condition at highly unsteady flows (high reduced frequency) or more viscous ones (low Reynolds numbers) has been reported in numerous studies.19,28–33 The reader is referred to the recent work by Taha and Rezaei34 for more references reporting flaws with the application of Kutta condition in unsteady flows. In the same paper, the authors developed a viscous extension of the classical potential flow unsteady aerodynamics (with emphasis on the lift frequency response problem: Theodorsen's problem). They relaxed/corrected the Kutta condition by matching the potential flow solution with a special boundary layer theory that pays close attention to flow details in the immediate vicinity of the trailing edge: the triple-deck boundary layer theory.34 Their model introduces a singularity at the trailing edge, which must vanish according to the Kutta condition. In contrast, its strength is obtained from the triple-deck boundary layer theory, which provides dependence on the Reynolds number and a viscous unsteady effective angle of attack. Consequently, unlike the Theodorsen inviscid response, which is only a function of the motion frequency, the new model provides dependence on both the frequency and Reynolds number. Their model34 is called the viscous theory in this work. This viscous theory is semi-analytical; it relies on several simplifying assumptions commonly adopted in theoretical/analytical models of unsteady aerodynamics, namely flat plate and flat wake assumptions. Also, its focus was on the frequency response problem (i.e., considering harmonic motions only).

The main contribution of this paper revolves around the integration of the recently developed unsteady viscous theory by Taha and Rezaei34 and numerical vortex methods (UVLM and DVM) to allow for

• arbitrary body shapes (not necessarily flat plate),

• wake deformation, and

• arbitrary motion kinematics (not necessarily harmonic),

hence making it useful for many applications in aeronautical engineering. That is, the objective is to develop viscous UVLM and DVM methods. This goal is achieved by replacing the Kutta condition in the vortex methods with a viscous condition based on the triple-deck boundary layer theory.

In Sec. II, we present the viscous correction of the Kutta condition from the triple-deck boundary layer theory. We then develop a viscous DVM in Sec. III, and a viscous UVLM in Sec. IV. Finally, results and validations in terms of step response (i.e., Wagner function); frequency response (i.e., Theodorsen function); and lift history for an arbitrary time-varying (e.g., non-harmonic) input are presented.

The dynamics of fluid flows are governed by a set of nonlinear partial differential equations (Navier–Stokes) representing conservation laws. In the early 1900, Prandtl derived an approximation of the Navier–Stokes equations in the vicinity of a flat plate: the so-called boundary layer equations.35 This formulation has formed the basis for important advances in fluid mechanics from the early efforts of Blasius36 until the recent efforts of Cantwell37 and Majdalani and Xuan.38 Prandtl's boundary layer concept2 assumes that at high Reynolds numbers and small disturbances, vorticity effects are confined to a thin layer around the surface, and the flow outside can be assumed irrotational. An immediate consequence of such assumption simplifies the nonlinear Navier–Stokes equations to the linear Laplace equation, thereby allowing superposition of solutions. For example: the flow over a circular cylinder can be constructed from the superposition of a free stream and a doublet.10,39

The flow over a cylinder can be conformally mapped to the flow over an arbitrarily given shape (thanks to Riemann mapping theorem), which allows solving the flow field over more complex shapes, such as airfoils. However, no lift force can be computed from this framework, unless a vortex is placed at the center of the cylinder, which retains the cylinder as a streamline yet generates circulation or lift force through the Kutta–Joukowsky theorem. The potential flow framework requires an auxiliary condition to determine the strength of the added vortex, typically taken as the Kutta condition if the body has a sharp trailing edge. Formulating the problem in the cylinder domain, the Kutta condition dictates the value of the circulation by forcing the point on the cylinder corresponding to the body trailing edge to be a stagnation point. As such, the Kutta condition removes the velocity (and pressure) singularity at the trailing edge.26 It should be noted that there is also a singularity at the leading edge, which remains intact. Kutta left it to be removed via a smooth geometry around the edge, which is typically the case.

If the Kutta condition is relaxed, the velocity u on a flat plate of length 2b in a free stream U at a steady angle of attack αs is written as40 follows:

$uU=1±αsb−xb+x∓Bs1−(x/b)2,$
(1)

where x is the plate coordinate (i.e., $−b≤x≤b$). The first row of signs corresponds to the suction side and the second row corresponds to the pressure side of the plate. The first two terms on the right-hand side of Eq. (1) represent the potential flow solution, satisfying the Kutta condition. However, the third term represents a correction to the Kutta condition. This term is equivalent to the velocity induced by a vortex $Γs=2πbUBs$ at the center of the cylinder. In other words, the Kutta condition dictates Bs = 0 so that the trailing edge singularity is removed.26 It should be noted that similar to the common leading edge singularity, the new trailing edge singularity is an integrable singularity: although the velocity (or pressure) is singular, the integrated lift and moment are finite. Thus, by applying Bernoulli's equation to Eq. (1) and integrating over the plate, the aerodynamic coefficients for lift and moment at the leading edge can be written as

$cl=2π(αs−Bs) ,$
(2)
$cm=0.5π(αs−2Bs) .$
(3)

It can easily be shown that after relaxing the Kutta condition, the new stagnation point sits at $θst=arcsinBs$, where θ is defined as $x=b cos θ$.

To find the strength Bs of the trailing edge singularity, Brown and Stewartson40 adopted the triple-deck boundary layer theory to the case of a flat plate with a non-zero angle of attack on the order of $ℜ−1/16$. This theory investigates flow events in the vicinity of the plate trailing edge. In this region, two boundary layers interact with each other, as illustrated in Fig. 1: the flow upstream of the edge with a Blasius boundary layer36 (satisfying the no-slip boundary condition on the plate) and the flow downstream of the edge with a Goldstein shear layer41 in the wake region (which satisfies a zero-stress condition on the wake centerline). To resolve the rapid variations of flow quantities in the tangential direction in the singular region around the trailing edge, the triple-deck theory applies scaling and asymptotic analysis also in the tangential direction (differently from Prandtl's boundary layer theory,35 which performs asymptotic analysis only in the normal direction). By doing so, the triple-deck theory accounts for the discontinuity of the viscous boundary condition: from a zero slip on the plate to a zero stress on the wake centerline. In this theory, three decks divide the fluid domain: (i) the upper deck, representing irrotational flow outside the boundary layer; (ii) the main layer, which constitutes an inviscid layer (though rotational); and (iii) the lower deck—a sub-layer inside the boundary layer where Prandtl's nonlinear boundary layer equations apply subject to a discontinuous viscous boundary condition. The interaction between the decks is such that asymptotic solutions from one deck approaching a neighboring one must match.

Brown and Stewartson40 constructed an asymptotic solution ($ϵ→0$ where ϵ is related to the Reynolds number $ℜ$ as $ϵ=ℜ−1/8$) to model the local interactions between the Blasius boundary layer and the shear layer in the wake in the immediate vicinity of the trailing edge. Their approach resulted in a nonlinear boundary value partial differential equation; and they provided a solution for its linear approximations. Six years later, Chow and Melnik42 solved this partial differential equation numerically, and provided the strength of the singularity

$Bs=2αsε3λ−5/4Be(αe),$
(4)

in terms of the normalized angle of attack

$αe=αsϵ−1/2λ−9/8,$
(5)

where $λ=0.334$ is the Blasius skin friction coefficient. The term $Be(αe)$ is a nonlinear function that is obtained from the numerical solution of Chow and Melnik,42 as shown in Fig. 2; the function possesses a vertical asymptote at $αe=0.47$ where trailing edge stall occurs according to the triple-deck results.40 This value of αe corresponds to $αs=3.1°−4.2°$ for $ℜ=104−106$ beyond which the flow separates upstream of the trailing edge.

As for the unsteady cases, Taha and Rezaei34 introduced an additional (beyond Kutta's) vortex at the center of the cylinder, which modifies the pressure distribution over the airfoil as follows:

$P(θ,t)−P∞=ρ[12a0(t) tan θ2+∑n=1∞an(t) sin nθ+12Bv(t)(cotθ2+a0v(t) tan θ2)] ,$
(6)

where ρ is the air density; the first two terms represent the inviscid pressure distribution, which satisfies the Kutta condition; $Bv(t)$ represents the viscous correction, which is given by Taha and Rezaei34 as follows:

$Bv(t)=−2ϵ3λ−5/4(12a0(t)+2∑n=1∞nan(t))Be(αe) ,$
(7)

and the normalized angle of attack αe is determined from αs and Reynolds number using Eq. (5). Here, αs represents an effective angle of attack: an equivalent steady angle of attack defined as34

$αs(t)=1U2|12a0(t)+2∑n=1∞nan(t)| .$
(8)

Therefore, to determine the viscous correction Bv, the coefficients an ($n=0, 1, 2, …$), representing the potential flow pressure distribution, must be evaluated. We write the plate normal velocity in the Fourier cosine series as follows:

$vp(θ,t)=12b0(t)+∑n=1∞bn(t) cos nθ .$
(9)

The no-penetration boundary condition provides the coefficients an ($n=1, 2, …$) as follows:

$an(t)=b2nḃn−1(t)+Ubn(t)−b2nḃn+1(t).$
(10)

The above expression can be numerically determined for any arbitrary motion including a flexible one-dimensional structure (e.g., time-varying camber). The only unknown term at this point is $a0(t)$. Robinson and Laurmann43 presented an integral equation to determine $a0(t)$ for any arbitrary motion. In the case of harmonic motions of a flat plate, ignoring wake deformation $a0(t)=U[2v3/4(t)C(k)+bα̇(t)]$, where C(k) is the Theodorsen's lift frequency response function, $C(k)=H1(2)(k)H1(2)(k)+iH0(2)(k)$, written in terms of Hankel functions of the second kind, $k=ωb/U$ is the reduced frequency of oscillation; $α(t)$ is the instantaneous angle of attack; and $v3/4(t)=V3/4eiωt$ is the normal velocity at the three-quarter-chord point (positive upward). However, to allow for arbitrary motion and to account for wake deformation, a0 can be obtained from the inviscid pressure distribution (coming from a traditional UVLM/DVM) as

$a0(t)=∫−bbΔPπρbdx−a1(t) ,$
(11)

where $ΔP=Pl−Pu$, Pu and Pl denote the pressure on the upper and lower airfoil surfaces, respectively. Alternatively, it can be determined from the potential flow lift coefficient as

$a0(t)=−U2cl(t)π−a1(t) .$
(12)

Based on the viscous correction discussed above, Taha and Rezaei34,44 provided a viscous version of the celebrated Theodorsen's inviscid lift frequency response. Building up on this effort, the objective of this work is to develop a computational tool for unsteady aerodynamic simulation of airfoils undergoing arbitrary motion, not necessarily harmonic. We achieve this following the seminal theoretical effort of von Kármán and Sears,45 replacing the classical Kutta condition in their effort with the viscous correction provided by Taha and Rezaei.34 The model proposed by von Kármán and Sears45 assumes wake singularities convecting along the x axis only (i.e., rigid flat wake)—a convenient assumption for small amplitudes. However, for relatively large-amplitude maneuvers, wake deformation becomes relevant (allowing wake vortices to move with the local fluid velocity) to capture some of the essential flow physics; for the importance of wake deformation effects, the reader is referred to Refs. 46 and 47.

Figure 3 illustrates a schematic of the modified von Kárman and Sears approach. Consider a wake vortex $Γ1$ at $ζ=ζ1$, where the ζ-plane ($ζ=χ+iφ$) is a plane containing a circular cylinder of radius $r=b/2$ centered at the origin. A standard conformal mapping then transforms the cylinder in the ζ-domain to a plate of chord 2b aligned with the x axis of the z-plane ($z=x+iy$). Using the method of images, we place an image vortex of intensity $−Γ1$ inside the cylinder at $ζ=r2/ζ1¯$, where $ζ1¯$ is the complex conjugate of ζ1, so the resulting flow field (due to the two vortices) is tangential to the cylinder everywhere (i.e., it does not disturb the no-penetration boundary condition). However, this pair of vortices induces a non-zero tangential velocity at the point $ζ=r$ (corresponding to the trailing edge in the plate domain), thereby violating the Kutta condition. At this point, von Kármán and Sears45 satisfied the Kutta condition by introducing a vortex Γk at the center of the cylinder. Note that this vortex does not disturb the no-penetration boundary condition on the cylinder either. Therefore, the complex potential F1 associated with this system of vortices is written as follows:

$F1=iΓ12π[ln (ζ−ζ1)−ln (ζ−r2ζ1¯)]+iΓk2πln ζ .$
(13)

Enforcing the Kutta condition (the total tangential velocity at the point $ζ=r$ must vanish), we obtain

$Γk=−Γ1r(1r−ζ1−1r−r2ζ1¯) .$
(14)

The addition of Γk violates Kelvin's conservation of circulation: the sum of all vortices in the system is equal to $Γk≠0$. In addition, the strength of the wake vortex $Γ1$ remains unknown. Both problems are solved when the quasi-steady circulation around the airfoil, $γ0(x)$, is taken into account. This distribution, by construction, already satisfies the no-penetration boundary condition (accounting for the airfoil motion), and also the Kutta condition. The resultant quasi-steady bound circulation is defined as $Γ0=∫−bbγ0(x)dx$. Hence, Kelvin's conservation of circulation requires $Γk=−Γ0$, which results in the following equation for the strength $Γ1$ of the nascent wake vortex:

$Γ1=−rχ0−rχ0+r[Γ0r−∑n=2NΓn(1r−ζn−1r−r2ζn¯)] ,$
(15)

where the nascent position of the new vortex at each time step is $ζ=χ0$, with $χ0∈ℝ>0$ and the N − 1 previously shed vortices remain constant while moving with the local fluid velocity (the Kirchhoff velocity) according to Helmholtz laws of vortex dynamics.48

On the other hand, the approach formulated by Taha and Rezaei34 and presented in Sec. II is equivalent to the addition of a vortex $Γv=2πbBv(t)/U$ at the center of the cylinder. Consequently, the conservation of circulation must be modified as follows:

$Γk=−Γ0−Γv ,$
(16)

which yields the strength $Γ1$ of the wake vortex as

$Γ1=−rχ0−rχ0+r[Γ0+Γvr−∑n=2NΓn(1r−ζn−1r−r2ζn¯)] .$
(17)

Equation (17) provides an extension of the theoretical model of von Kármán and Sears45 [Eq. (11) in their work] to account for (i) wake deformation, considering discrete vortices; and more importantly to the objective of this paper (ii) viscous effects through the Γv-term. If the vortices are only allowed to move along the x axis and the Kutta condition is applied (i.e., $Γv=0$), one can observe after some mathematical manipulation that Eq. (17) recovers the result obtained by Tchieu and Leonard14 in their Eq. (18).

Having obtained the strength $Γ1$ of the nascent vortex, the aerodynamic lift force can be computed from conservation of momentum.13,14,16,45 The original work by von Kármán and Sears45 classifies the aerodynamic loading into three parts: a contribution of the apparent mass, which is analogous to the non-circulatory lift introduced by Theodorsen6; a contribution of quasi-steady aerodynamics, which is directly related to (solely determined by) the airfoil motion; and finally a contribution from the wake. The non-circulatory lift is assumed to remain the same in the present work. However, following the derivation by Michelin and Smith,13 one can note that the circulatory lift coefficient in the language of Theodorsen6—which corresponds to the sum of quasi-steady and wake components of von Kármán and Sears45—is given by

$clc=1U2bℑ[eiα(t)∑n=1NΓniddt(ζn−r2ζn¯)] ,$
(18)

where $ℑ$ denotes the imaginary part of a complex number.

Similarly, following Tchieu and Leonard,14 the circulatory moment coefficient at the mid chord is

$cmc=−clc4−Γ0+Γv4Ub .$
(19)

Finally, all the wake vortices move in the physical domain according to the local induced velocity13 as follows:

$dxdt−idydt=1z′(ζ){Ue−iα(t)(1−r2ζ2e2iα(t))−limζ→ζj[∑n=1NiΓn2π(1ζ−ζn−1ζ−r2ζn¯)−iΓj2π(1ζ−ζj−1ζ−r2ζj¯)]}.$
(20)

Unlike the DVM, which relies on conformal mapping between a cylinder and the airfoil shape, the UVLM works directly in the airfoil domain. Hence, it allows more freedom in handling arbitrary shapes and flexible (time-varying) cambers. Consider the standard UVLM, which divides the airfoil camber in panels, with a vortex placed at 1/4 of the panel length and a control point at 3/4 of the panel length, as illustrated in Fig. 4.

The Neumann boundary condition (no-penetrability condition) is applied at the control points in contrast to the DVM where the no-penetration boundary condition is satisfied everywhere over the airfoil. As such, we have

$(U+vb+vw)·n=0 ,$
(21)

where U is the local freestream velocity, $vb$ is the total velocity induced by bound vortices, $vw$ is the velocity induced by the wake vortices, and n is the normal vector to each panel. The velocity induced by a vortex at (xv, yv) onto a point (xp, yp) is

$v=Γ2π(yp−yv)x̂−(xp−xv)ŷ(xp−xv)2+(yp−yv)2 ,$
(22)

where Γ is the vortex strength and $x̂$ and $ŷ$ represent vectors along the x and y directions, respectively.

Similar to the DVM, one vortex $Γ1$ is shed in the wake at each time step, typically from the trailing edge. If the airfoil is divided into m panels, the no-penetration boundary condition offers m equations. However, one more equation is required because we have m +1 unknowns: m bound circulations and the newly shed vortex $Γ1$. The additional equation is given through Kelvin's conservation of circulation as follows:

$Γb(t)−Γb(t−Δt)+Γ1=0 ,$
(23)

where $Γb=∑i=1mΓbi$ denotes the total bound circulation.

The replacement of the Kutta condition by the viscous correction proposed by Taha and Rezaei34 must consider the effects of the additional vortex Γv at the center of the cylinder that is conformally mapped to the airfoil (as presented in Sec. III). To account for its effects in the UVLM framework without relying on conformal mapping, we follow von Karman and Sears ideas as follows: The additional vortex, Γv, induces a tangential velocity distribution over the cylinder by the spatial constant $vθ=2Bv(t)/U$. Consequently, this change affects the tangential velocity over the plate as

$u=−vθ21−(x/b)2 ,$
(24)

which corresponds to an additional distribution of circulation as follows:

$ΔΓbi=2uiΔsi ,$
(25)

where $Δsi$ is the length of the ith panel.

As such, a simple modification in the UVLM system of equations provides a viscous extension of the Kutta condition in the problem. Recalling the original system of equations, m no-penetration boundary conditions applied at the control points and Kelvin's condition, and adding the effects of $ΔΓbi$, we write

$(A1,1A1,2…A1,mA1,wA2,1A2,2…A2,mA1,w⋮⋮⋱⋮⋮Am,1Am,2…Am,mA1,w11…11)({Γb1Γb2⋮ΓbmΓ1}+{ΔΓb1ΔΓb2⋮ΔΓbm0})=−{(U+vw)1·n1(U+vw)2·n2⋮(U+vw)m·nm∑j=2NΓj−∑i=1mΔΓbi} ,$
(26)

where $Ai,j$ are the influence coefficients depending solely on the airfoil geometry. The values of $ΔΓbi$ are previously computed from the calculation of $Bv(t)$ at each time step. As a consequence, the only unknowns remaining in the problem are the values of the circulations $Γbi$ and the nascent wake vortex, $Γ1$. The other (previously shed) wake vortices are assumed to move freely with the local fluid velocity (Kirckhhoff velocity) while maintaining their strength (dictated at the shedding time) to satisfy Helmholtz laws of vortex dynamics.49,50 The reader may refer to Katz and Plotkin10 for the standard UVLM setup.

The pressure difference on each panel is given by the unsteady Bernoulli's equation as follows:

$ΔPi=ρUΓbidi+ρ∂∂t∑j=1iΓbj,$
(27)

and the integration of this pressure difference provides the lift force and pitching moment.

A flow chart is presented in Fig. 5 to illustrate how to proceed when applying viscous corrections in UVLM and DVM (according to the model developed by Taha and Rezaei34).

To verify and validate the proposed methods (viscous DVM and viscous UVLM), three different cases are considered: the step response; the frequency response (harmonic motion) to verify the proposed numerical methods with reference to the theoretical results of Taha and Rezaei;34 and the response due to a multi-frequency non-harmonic signal. In the viscous DVM simulations, the nascent vortices are shed at $x0=1.01b$ to avoid a trivial result in Eq. (17). The viscous UVLM simulations are performed with m =40 and $Δt=0.02b/U$ was selected for both UVLM and DVM.

Figure 6 shows the time-history of lift and moment due to a step change $αs=3°$ in the angle of attack at $ℜ=104$, also known as the Wagner response.5 The inviscid and viscous solutions approach different asymptotes; expectedly, the traditional UVLM and DVM approach Kutta's potential flow loading whereas the modified versions approach the viscous values given by Eqs. (2) and (3). In fact, these results verify that the coupling between the traditional UVLM/DVM and the viscous theory is correctly achieved.

FIG. 6.

Wagner (step) response at $ℜ=104$ and $α=3°$: (a) lift coefficient and (b) moment coefficient at the leading edge.

FIG. 6.

Wagner (step) response at $ℜ=104$ and $α=3°$: (a) lift coefficient and (b) moment coefficient at the leading edge.

Close modal

An important advantage of the viscous UVLM, in contrast to DVM, is its ability to handle arbitrary airfoil shapes (i.e., not confined to airfoils that are conformally mapped to a cylinder). To demonstrate this point, we simulate a step change in the camber line from zero (i.e., a flat plate) to that of a National Advisory Committee for Aeronautics (NACA) 2412 at a zero angle of attack and $ℜ=105$. Figure 7(a) shows aerodynamic loads computed from the viscous UVLM. Figure 7(b) also presents a visualization of the airfoil wake.

FIG. 7.

Step change of the camber from zero (flat plate) to that of a NACA 2412 at $ℜ=105$ and $α=0°$: (a) lift coefficient and moment coefficient at the leading edge and (b) wake visualization from the viscous UVLM.

FIG. 7.

Step change of the camber from zero (flat plate) to that of a NACA 2412 at $ℜ=105$ and $α=0°$: (a) lift coefficient and moment coefficient at the leading edge and (b) wake visualization from the viscous UVLM.

Close modal

Figure 8(a) provides a comparison in terms of the circulatory lift frequency response from Theodorsen, the viscous theory of Taha and Rezaei,34 the proposed viscous DVM and viscous UVLM, in reference to URANS simulations of a harmonically pitching NACA 0012 at $ℜ=105$. URANS simulations employed the $k−ω$ SST turbulence model. The reader is referred to Rezaei and Taha51 or Taha and Rezaei34 for more details on the setup in the computational fluid dynamics (CFD) solver. The viscous UVLM and DVM capture the additional phase lag between the airfoil motion (represented by the quasi-steady lift) and the resulting unsteady circulatory lift in contrast to Theodorsen's inviscid response, as discussed by Taha and Rezaei.34 Moreover, the amplitude of the circulatory lift slightly decreases with the relaxation of the Kutta condition due to viscous effects. The amplitude and phase of the viscous approaches are in consonance with CFD data. The wake computed by the viscous UVLM considering k =1 is also presented in Fig. 8(b).

FIG. 8.

Analyses of harmonic oscillations: (a) frequency response between the circulatory and quasi-steady lift ($ℜ=105$) and (b) wake visualization from the viscous UVLM (k = 1).

FIG. 8.

Analyses of harmonic oscillations: (a) frequency response between the circulatory and quasi-steady lift ($ℜ=105$) and (b) wake visualization from the viscous UVLM (k = 1).

Close modal

The proposed viscous DVM and UVLM also have the ability to simulate arbitrary time-varying motion, not necessarily harmonic oscillations. Figure 9 provides a comparison of the total lift and pitching moment predicted by the potential flow theory, URANS simulations and the two viscous methods presented in this work due to the non-harmonic maneuver

$α(t)=𝒜[e sin (ωt)−1] ,$
(28)

with k =1 at $ℜ=105$ and $𝒜$ designed to ensure the maximum angle of attack equal to $1°$ during the maneuver. The wake simulated by the viscous UVLM for this case is presented in Fig. 9(c). The lift coefficient predicted by the viscous UVLM and DVM are basically the same and they are closer to the CFD results than the results from the potential theory. The latter clearly over-predicts the unsteady loading. The resulting moment coefficient from the viscous approaches (DVM and UVLM) are even closer to the URANS computations. Therefore, the developed numerical techniques are capable of handling (i) unsteady effects; (ii) viscous effects (e.g., viscosity-induced lag) on the lift dynamics at high reduced frequencies and small Reynolds numbers; (iii) wake deformation; and (iv) arbitrary time-varying motion. In addition, the viscous UVLM allows handling arbitrary shapes and flexible time-varying cambers. It should be noted that the proposed two methods provide these features without requiring significant changes of existing potential flow codes. Given a traditional DVM code, one only needs to replace the Kutta condition in Eq. (15) with the viscous correction in Eq. (17); and given a potential flow UVLM code, one only needs to add the additional column of $ΔΓbi$ in Eq. (26).

FIG. 9.

Response to $α(t)=𝒜[e sin (ωt)−1]$ at $ℜ=105$: (a) lift coefficient, (b) moment coefficient at the quarter-chord, and (c) wake visualization from the viscous UVLM.

FIG. 9.

Response to $α(t)=𝒜[e sin (ωt)−1]$ at $ℜ=105$: (a) lift coefficient, (b) moment coefficient at the quarter-chord, and (c) wake visualization from the viscous UVLM.

Close modal

Finally, it is worth mentioning that the proposed viscous DVM has a strong dependence on the position of the nascent wake vortices. Such a limitation is not observed in the viscous UVLM. Therefore, the latter may represent a better (though more expensive) numerical method for viscous simulation of unsteady aerodynamics.

This article proposes two different numerical viscous extensions of commonly used computational methods for unsteady aerodynamics: the discrete vortex method (DVM) and the unsteady vortex lattice/panel method (UVLM). The extensions replace the largely used Kutta condition by a condition derived from a special boundary layer theory that resolves the flow field in the immediate vicinity of the trailing edge: the triple-deck boundary layer theory. The traditional DVM, which relies on conformal mapping between a cylinder and the body, is modified by introducing an additional vortex beyond Kutta's at the center of the cylinder whose strength is determined from the triple-deck boundary layer theory based on the airfoil motion and Reynolds number. Similarly, the classical potential flow UVLM is modified by introducing additional circulation over each panel such that the corresponding tangential velocity distribution, when mapped to a cylinder, matches the one dictated by the triple-deck additional vortex. As such, the proposed numerical techniques keep the same layout of the potential flow solvers (DVM and UVLM) while applying viscous corrections for any generic maneuver. The obtained results are validated against URANS simulations for harmonic, step, and complex maneuvers. The resulting lift and pitching moment are consistently in better agreement with the URANS simulations in comparison to their potential flow counterparts.

C. R. dos Santos acknowledges the sponsorship of the Brazilian agency (Grant Nos. 2017/09468‐6 and 2018/05247‐8), São Paulo Research Foundation (FAPESP). H. E. Taha would like to acknowledge the support of the National Science Foundation, USA (Grant No. CBET-2005541).

C.R.d.S. and A.S.R. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
R.
Dugas
,
A History of Mechanics
(
Courier Corporation
,
2012
).
2.
L.
Prandtl
, “
Über flussigkeitsbewegung bei sehr kleiner reibung
,”
Verhandl. III, Internat. Math.-Kong., Heidelberg, Teubner, Leipzig
1904
,
484
491
.
3.
W.
Birnbaum
, “
Die tragende wirbelfläche als hilfsmittel zur behandlung des ebenen problems der tragflügeltheorie
,”
ZAMM-J. Appl. Mathematics Mech./Z. Angew. Mathematik und Mechanik
3
,
290
297
(
1923
).
4.
L.
Prandtl
,
Gesammelte Abhandlungen Zur Angewandten Mechanik, Hydro-und Aerodynamik
(
Springer
,
1918
), Vol.
1
.
5.
H.
Wagner
, “
Über die entstehung des dynamischen auftriebes von tragflügeln
,”
ZAMM-J. Appl. Math. Mech./Z. Angew. Math. Mech.
5
,
17
35
(
1925
).
6.
T.
Theodorsen
, “
General theory of aerodynamic instability and the mechanism of flutter
,” Technical Report No. 496 (
NACA
,
1935
).
7.
S. M.
Belotserkovskii
, “
Study of the unsteady aerodynamics of lifting surfaces using the computer
,”
Annu. Rev. Fluid Mech.
9
,
469
494
(
1977
).
8.
R. M.
James
, “
On the remarkable accuracy of the vortex lattice method
,”
Comput. Methods Appl. Mech. Eng.
1
,
59
79
(
1972
).
9.
V. M.
Falkner
, “
,” Technical Report No. 1910 (
Aeronautical Research Council, London
, UK,
1943
).
10.
J.
Katz
and
A.
Plotkin
,
Low-Speed Aerodynamics
(
Cambridge University Press
,
New York
,
2001
), Vol.
13
.
11.
T.
Sarpkaya
, “
Computational methods with vortices–the 1988 freeman scholar lecture
,”
J. Fluids Eng.
111
,
5
(
1989
).
12.
L.
Cortelezzi
, “
On the unsteady separated flow past a semi-infinite plate: Exact solution of the Brown and Michael model, scaling, and universality
,”
Phys. Fluids
7
,
526
529
(
1995
).
13.
S.
Michelin
and
S. G. L.
Smith
, “
An unsteady point vortex method for coupled fluid–solid problems
,”
Theor. Comput. Fluid Dyn.
23
,
127
153
(
2009
).
14.
A. A.
Tchieu
and
A.
Leonard
, “
A discrete-vortex model for the arbitrary motion of a thin airfoil with fluidic control
,”
J. Fluids Struct.
27
,
680
693
(
2011
).
15.
K.
Ramesh
,
A.
Gopalarathnam
,
K.
Granlund
,
M. V.
Ol
, and
J. R.
Edwards
, “
Discrete-vortex method with novel shedding criterion for unsteady aerofoil flows with intermittent leading-edge vortex shedding
,”
J. Fluid Mech.
751
,
500
538
(
2014
).
16.
C.
Wang
and
J. D.
Eldredge
, “
Low-order phenomenological modeling of leading-edge vortex formation
,”
Theor. Comput. Fluid Dyn.
27
,
577
598
(
2013
).
17.
C. P.
Ford
and
H.
Babinsky
, “
,”
J. Fluid Mech.
720
,
280
313
(
2013
).
18.
M. S.
Hemati
,
J. D.
Eldredge
, and
J. L.
Speyer
, “
Improving vortex models via optimal control theory
,”
J. Fluids Struct.
49
,
91
111
(
2014
).
19.
D.
Darakananda
and
J. D.
Eldredge
, “
A versatile taxonomy of low-dimensional vortex models for unsteady aerodynamics
,”
J. Fluid Mech.
858
,
917
948
(
2019
).
20.
A.
SureshBabu
,
K.
Ramesh
, and
A.
Gopalarathnam
, “
Model reduction in discrete-vortex methods for unsteady airfoil flows
,”
AIAA J.
57
,
1409
1414
(
2019
).
21.
B. P.
Epps
,
B. T.
Roesler
,
R. B.
Medvitz
,
Y.
Choo
, and
J.
McEntee
, “
A viscous vortex lattice method for analysis of cross-flow propellers and turbines
,”
Renewable Energy
143
,
1035
1052
(
2019
).
22.
A. A.
Hussein
,
A. E.
Seleit
,
H. E.
Taha
, and
M. R.
Hajj
, “
Optimal transition of flapping wing micro-air vehicles from hovering to forward flight
,”
Aerosp. Sci. Technol.
90
,
246
263
(
2019
).
23.
M. Y.
Zakaria
,
H. E.
Taha
, and
M. R.
Hajj
, “
Design optimization of flapping ornithopters: The pterosaur replica in forward flight
,”
J. Aircr.
53
,
48
59
(
2016
).
24.
M. Y.
Zakaria
,
H. E.
Taha
, and
M. R.
Hajj
, “
Measurement and modeling of lift enhancement on plunging airfoils: A frequency response approach
,”
J. Fluids Struct.
69
,
187
208
(
2017
).
25.
C. R.
dos Santos
,
D. R. Q.
Pacheco
,
H. E.
Taha
, and
M. Y.
Zakaria
, “
Nonlinear modeling of electro-aeroelastic dynamics of composite beams with piezoelectric coupling
,”
Compos. Struct.
255
,
112968
(
2021
).
26.
D. G.
Crighton
, “
The Kutta condition in unsteady flow
,”
Annu. Rev. Fluid Mech.
17
,
411
445
(
1985
).
27.
M. W.
Kutta
, “
Über eine mit den Grundlagen des Flugproblems in Beziehung stehende zweidimensionale Strömung
,”
Sitzungsberichte der Königlich Bayerischen Akademie der Wissenschaften
2
,
1
58
(
1910
).
28.
F. S.
Archibald
, “
Unsteady Kutta condition at high values of the reduced frequency parameter
,”
J. Aircraft
12
,
545
550
(
1975
).
29.
C.
Ho
and
S.
Chen
, “
Unsteady Kutta condition of a plunging airfoil
,” in
(
Springer
,
1981
), pp.
197
206
.
30.
D. R.
Poling
and
D. P.
Telionis
, “
The response of airfoils to periodic disturbances-the unsteady Kutta condition
,”
AIAA J.
24
,
193
199
(
1986
).
31.
S. A.
Ansari
,
R.
Żbikowski
, and
K.
Knowles
, “
Non-linear unsteady aerodynamic model for insect-like flapping wings in the hover. Part 2: Implementation and validation
,”
Proc. Inst. Mech. Eng. Part G: J. Aerosp. Eng.
220
,
169
186
(
2006
).
32.
M.
La Mantia
and
P.
Dabnichki
, “
Unsteady panel method for flapping foil
,”
Eng. Anal. Boundary Elem.
33
,
572
580
(
2009
).
33.
X.
Xia
and
K.
Mohseni
, “
Lift evaluation of a two-dimensional pitching flat plate
,”
Phys. Fluids
25
,
091901
(
2013
).
34.
H.
Taha
and
A. S.
Rezaei
, “
Viscous extension of potential-flow unsteady aerodynamics: The lift frequency response problem
,”
J. Fluid Mech.
868
,
141
175
(
2019
).
35.
H.
Schlichting
and
K.
Gersten
,
Boundary-Layer Theory
(
Springer
,
2016
).
36.
H.
Blasius
,
Grenzschichten in Flüssigkeiten Mit Kleiner Reibung
(
Druck von BG Teubner
,
1907
).
37.
B. J.
Cantwell
, “
Integral measures of the zero pressure gradient boundary layer over the Reynolds number range
$0≤ℜτ<∞$,”
Phys. Fluids
33
,
085108
(
2021
).
38.
J.
Majdalani
and
L.-J.
Xuan
, “
On the Kármán momentum-integral approach and the Pohlhausen paradox
,”
Phys. Fluids
32
,
123605
(
2020
).
39.
F. M.
White
,
Fluid Mechanics
(McGraw-Hill,
2016
).
40.
S. N.
Brown
and
K.
Stewartson
, “
Trailing-edge stall
,”
J. Fluid Mech.
42
,
561
584
(
1970
).
41.
S.
Goldstein
, “
Concerning some solutions of the boundary layer equations in hydrodynamics
,”
Math. Proc. Cambridge Philos. Soc.
26
,
1
30
(
1930
).
42.
R.
Chow
and
R. E.
Melnik
, “
Numerical solutions of the triple-deck equations for laminar trailing-edge stall
,” in
Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics, Twente University
,
Enschede, the Netherlands
, June 28–July 2, 1976 (
Springer
,
1976
), pp.
135
144
.
43.
A.
Robinson
and
J. A.
Laurmann
,
Wing Theory
(
Cambridge University Press
,
2013
).
44.
H. E.
Taha
and
A. S.
Rezaei
, “
,” in
AIAA Scitech 2021 Forum
(
2021
), p.
1830
.
45.
T.
von Kármán
and
W. R.
Sears
, “
Airfoil theory for non-uniform motion
,”
J. Aeronaut. Sci.
5
,
379
390
(
1938
).
46.
R. A.
and
S. S.
Davis
, “
Visualization of quasiperiodic flows
,”
AIAA J.
17
,
1164
1169
(
1979
).
47.
J. M. R.
Graham
, “
The lift on an aerofoil in starting flow
,”
J. Fluid Mech.
133
,
413
425
(
1983
).
48.
P. G.
Saffman
,
Vortex Dynamics
(
Cambridge University Press
,
1995
).
49.
H.
Helmholtz
, “
On Integrals of the hydrodynamical equations, which express vortex-motion
,”
London, Edinburgh, and Dublin Philos. Mag. J. Sci.
43
(
226
),
485
512
(
1867
).
50.
P. G.
Saffman
,
Vortex Dynamics
(
Cambridge University Press
,
1992
).
51.
A. S.
Rezaei
and
H. E.
Taha
, “
Computational study of lift frequency responses of pitching airfoils at low Reynolds numbers
,” in
55th AIAA Aerospace Sciences Meeting
(
AIAA
,
2017
), p.
0716
.