The Kutta condition has been extensively used to determine aerodynamic loads of steady and unsteady flows over airfoils. Nevertheless, the application of this condition to unsteady flows has been controversial for decades. A viscous correction to the Kutta condition was recently developed by matching the potential flow solution with a special boundary layer theory that resolves the flow field in the immediate vicinity of the trailing edge: the triple-deck boundary layer theory. In this work, we utilize this viscous condition to extend two common numerical methods for unsteady aerodynamics to capture viscous effects on the dynamics of unsteady lift and pitching moment—we develop viscous versions of the traditional discrete vortex method and unsteady vortex lattice/panel method. The resulting aerodynamic loads obtained from the proposed numerical models are compared against higher-fidelity simulations of unsteady Reynolds-averaged Navier–Stokes equations for an airfoil undergoing step, harmonic, and complex maneuvers. The obtained results are consistently in better agreement with the unsteady Reynolds-averaged Navier–Stokes simulations in comparison to their potential flow counterpart. In conclusion, the developed numerical methods are capable of capturing (i) unsteady effects; (ii) viscous effects (e.g., viscosity-induced lag) on the dynamics of lift and moment at high frequencies and low Reynolds numbers; and (iii) wake deformation, for arbitrary time-varying motion.

## I. INTRODUCTION

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 models^{6}) 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 Falkner^{9}—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 flight^{22,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 body^{26}; hence, an auxiliary condition is needed. The Kutta condition^{27} 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 Rezaei^{34} 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 model^{34} 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 Rezaei^{34} 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.

## II. VISCOUS CORRECTION OF THE KUTTA CONDITION USING THE TRIPLE-DECK BOUNDARY LAYER THEORY

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 Blasius^{36} until the recent efforts of Cantwell^{37} and Majdalani and Xuan.^{38} Prandtl's boundary layer concept^{2} 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 2*b* in a free stream *U* at a steady angle of attack *α _{s}* is written as

^{40}follows:

where *x* is the plate coordinate (i.e., $\u2212b\u2264x\u2264b$). 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 $\Gamma s=2\pi bUBs$ at the center of the cylinder. In other words, the Kutta condition dictates *B _{s}* = 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

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

To find the strength *B _{s}* of the trailing edge singularity, Brown and Stewartson

^{40}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 $\u211c\u22121/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 layer

^{36}(satisfying the no-slip boundary condition on the plate) and the flow downstream of the edge with a Goldstein shear layer

^{41}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 Stewartson^{40} constructed an asymptotic solution ($\u03f5\u21920$ where *ϵ* is related to the Reynolds number $\u211c$ as $\u03f5=\u211c\u22121/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 Melnik^{42} solved this partial differential equation numerically, and provided the strength of the singularity

in terms of the normalized angle of attack

where $\lambda =0.334$ is the Blasius skin friction coefficient. The term $Be(\alpha 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 $\alpha e=0.47$ where trailing edge stall occurs according to the triple-deck results.^{40} This value of *α _{e}* corresponds to $\alpha s=3.1\xb0\u22124.2\xb0$ for $\u211c=104\u2212106$ beyond which the flow separates upstream of the trailing edge.

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

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 Rezaei^{34} as follows:

and the normalized angle of attack *α _{e}* is determined from

*α*and Reynolds number using Eq. (5). Here,

_{s}*α*represents an effective angle of attack: an equivalent steady angle of attack defined as

_{s}^{34}

Therefore, to determine the viscous correction *B _{v}*, the coefficients

*a*($n=0,\u20091,\u20092,\u2009\u2026$), representing the potential flow pressure distribution, must be evaluated. We write the plate normal velocity in the Fourier cosine series as follows:

_{n}The no-penetration boundary condition provides the coefficients *a _{n}* ($n=1,\u20092,\u2009\u2026$) as follows:

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 Laurmann^{43} 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\alpha \u0307(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=\omega b/U$ is the reduced frequency of oscillation; $\alpha (t)$ is the instantaneous angle of attack; and $v3/4(t)=V3/4ei\omega 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, *a*_{0} can be obtained from the inviscid pressure distribution (coming from a traditional UVLM/DVM) as

where $\Delta P=Pl\u2212Pu$, *P _{u}* and

*P*denote the pressure on the upper and lower airfoil surfaces, respectively. Alternatively, it can be determined from the potential flow lift coefficient as

_{l}## III. VISCOUS EXTENSION OF THE DVM

Based on the viscous correction discussed above, Taha and Rezaei^{34,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 Sears^{45} 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 $\Gamma 1$ at $\zeta =\zeta 1$, where the *ζ*-plane ($\zeta =\chi +i\phi $) 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 2*b* aligned with the *x* axis of the *z*-plane ($z=x+iy$). Using the method of images, we place an image vortex of intensity $\u2212\Gamma 1$ inside the cylinder at $\zeta =r2/\zeta 1\xaf$, where $\zeta 1\xaf$ 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 $\zeta =r$ (corresponding to the trailing edge in the plate domain), thereby violating the Kutta condition. At this point, von Kármán and Sears^{45} 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 *F*_{1} associated with this system of vortices is written as follows:

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

The addition of Γ_{k} violates Kelvin's conservation of circulation: the sum of all vortices in the system is equal to $\Gamma k\u22600$. In addition, the strength of the wake vortex $\Gamma 1$ remains unknown. Both problems are solved when the quasi-steady circulation around the airfoil, $\gamma 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 $\Gamma 0=\u222b\u2212bb\gamma 0(x)dx$. Hence, Kelvin's conservation of circulation requires $\Gamma k=\u2212\Gamma 0$, which results in the following equation for the strength $\Gamma 1$ of the nascent wake vortex:

where the nascent position of the new vortex at each time step is $\zeta =\chi 0$, with $\chi 0\u2208\mathbb{R}>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 Rezaei^{34} and presented in Sec. II is equivalent to the addition of a vortex $\Gamma v=2\pi bBv(t)/U$ at the center of the cylinder. Consequently, the conservation of circulation must be modified as follows:

which yields the strength $\Gamma 1$ of the wake vortex as

Equation (17) provides an extension of the theoretical model of von Kármán and Sears^{45} [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., $\Gamma v=0$), one can observe after some mathematical manipulation that Eq. (17) recovers the result obtained by Tchieu and Leonard^{14} in their Eq. (18).

Having obtained the strength $\Gamma 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 Sears^{45} classifies the aerodynamic loading into three parts: a contribution of the apparent mass, which is analogous to the non-circulatory lift introduced by Theodorsen^{6}; 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 Theodorsen^{6}—which corresponds to the sum of quasi-steady and wake components of von Kármán and Sears^{45}—is given by

where $\u2111$ denotes the imaginary part of a complex number.

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

Finally, all the wake vortices move in the physical domain according to the local induced velocity^{13} as follows:

## IV. VISCOUS EXTENSION OF THE UVLM

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

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 (*x _{v}*,

*y*) onto a point (

_{v}*x*,

_{p}*y*) is

_{p}where Γ is the vortex strength and $x\u0302$ and $y\u0302$ represent vectors along the *x* and *y* directions, respectively.

Similar to the DVM, one vortex $\Gamma 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 $\Gamma 1$. The additional equation is given through Kelvin's conservation of circulation as follows:

where $\Gamma b=\u2211i=1m\Gamma bi$ denotes the total bound circulation.

The replacement of the Kutta condition by the viscous correction proposed by Taha and Rezaei^{34} 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\theta =2Bv(t)/U$. Consequently, this change affects the tangential velocity over the plate as

which corresponds to an additional distribution of circulation as follows:

where $\Delta si$ is the length of the *i*th 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 $\Delta \Gamma bi$, we write

where $Ai,j$ are the influence coefficients depending solely on the airfoil geometry. The values of $\Delta \Gamma 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 $\Gamma bi$ and the nascent wake vortex, $\Gamma 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 Plotkin^{10} for the standard UVLM setup.

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

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

## V. COMPARISON OF THE METHODS

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 $\Delta 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 $\alpha s=3\xb0$ in the angle of attack at $\u211c=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.

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 $\u211c=105$. Figure 7(a) shows aerodynamic loads computed from the viscous UVLM. Figure 7(b) also presents a visualization of the airfoil wake.

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 $\u211c=105$. URANS simulations employed the $k\u2212\omega $ SST turbulence model. The reader is referred to Rezaei and Taha^{51} or Taha and Rezaei^{34} 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).

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

with *k *=* *1 at $\u211c=105$ and $\mathcal{A}$ designed to ensure the maximum angle of attack equal to $1\xb0$ 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 $\Delta \Gamma bi$ in Eq. (26).

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.

## VI. CONCLUDING REMARKS

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.

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Author Contributions

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

## DATA AVAILABILITY

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

## References

*June 28–July 2, 1976*(