Airfoil Panel Method

Airfoil Panel Method

Airfoil Panel Method: Potential Flow About a Section

The tool computes the two-dimensional, incompressible, inviscid flow about an airfoil of any thickness and camber. The surface is divided into straight panels carrying vortex or source sheets whose strengths are found from two conditions: no flow passes through the surface, and the flow leaves the trailing edge smoothly (the Kutta condition). The result is the velocity and pressure everywhere on the surface, and from them the lift, the moment and the quantities derived from the pressure distribution. Thin-airfoil theory, which the Airfoil and Wing Lift tool uses, is the limit of zero thickness; the panel method keeps the thickness and so gives the pressure distribution itself.

Conventions

The chord lies on the x axis from the leading edge at x = 0 to the trailing edge at x = 1; the angle of attack α is measured from the chord line, and the free stream has unit speed. The nodes run clockwise: from the trailing edge along the lower surface to the leading edge and back along the upper surface. Moments are positive nose up. C_p = 1 − (v/V)² is the incompressible pressure coefficient. A plain flap is rotated about its hinge, trailing edge down positive, and the chord and angle of attack stay those of the unflapped section.

Geometry

NACA four-digit sections (NACA Report 460). The half thickness is ±y_t = 5 t (0.2969 √x − 0.1260 x − 0.3516 x² + 0.2843 x³ − 0.1015 x⁴) for a maximum thickness t; the mean line is two parabolas meeting at the point of maximum camber m at x = p, y_c = (m/p²)(2 p x − x²) ahead of it and (m/(1 − p)²)(1 − 2p + 2 p x − x²) behind it; the thickness is laid off perpendicular to the mean line, x_u = x − y_t sin θ, y_u = y_c + y_t cos θ with tan θ the slope of the mean line (the lower surface with the signs reversed). The first digit is m in percent, the second p in tenths, the last two t in percent.

NACA five-digit sections (NACA Report 537). The same thickness form on a mean line whose curvature falls linearly to zero at a juncture x = m slightly behind the maximum camber and stays zero from there to the trailing edge: y_c = (k₁/6)(x³ − 3 m x² + m²(3 − m) x) ahead of the juncture and (k₁ m³/6)(1 − x) behind it. Table I of the report gives m and k₁ for the camber positions 0.05 to 0.25 of the chord (mean lines 210 to 250) at a design lift coefficient of 0.3, and the reflexed lines 221 to 251, which add a second constant k₂ with k₂/k₁ = (3 (m − p)² − m³)/(1 − m)³ to bring the pitching moment to zero. The first digit scales the camber: the design lift coefficient is 0.15 times that digit, and k₁ is scaled in proportion. Mean line 211 is not in the report and is not offered. NACA 6-series sections are not generated from their names (their thickness forms are tabulated, not given by an equation): paste the coordinates.

NACA 6-series sections (NACA Report 824). The basic thickness forms of the 6-series were derived by conformal mapping from prescribed velocity distributions and have no closed equation; Report 824 gives them as tables of ordinates at 26 chord stations with the leading-edge radius. The tool carries those tables for the 63, 64, 65 and 66 families at 6, 8, 9, 10, 12, 15, 18 and 21 percent thickness (the 63 family has no 8 percent table). A thickness between two tables is interpolated, the ratio y/t taken linear in the thickness between the neighboring tables; predicting each printed table from its two neighbors in this way misses it by at most 0.013 percent of the chord across gaps of up to 6 percent, where scaling one table in proportion misses by 0.05 to 0.09 percent. Below 6 or above 21 percent the nearest table is scaled in proportion, which is no longer a member of the family; the result page says which of the three was done. Between the stations the half-thickness is a cubic spline in √x whose slope at the nose gives the printed leading-edge radius, y_t = √(2 r x). The trailing edge is closed as tabulated.

The mean line carries a uniform chordwise load from the leading edge to x = a and a load falling linearly to zero from there to the trailing edge; for the design lift coefficient c_li its ordinate is (Report 824, equation 6) y_c = c_li/(2π(a + 1)) { [½(a − x)² ln|a − x| − ½(1 − x)² ln(1 − x) + ¼(1 − x)² − ¼(a − x)²]/(1 − a) − x ln x + g − h x }, with g = −[a²(½ ln a − ¼) + ¼]/(1 − a) and h = [½(1 − a)² ln(1 − a) − ¼(1 − a)²]/(1 − a) + g, and its ideal angle of attack is α_i = −c_li h/(2π(a + 1)). For a = 1, the mean line assumed when none is named, y_c = −(c_li/4π)[(1 − x) ln(1 − x) + x ln x] and α_i = 0. The thickness is laid off perpendicular to the mean line as for the four-digit sections.

The designation is written with the hyphen: 64-212 is the 64 family, design lift coefficient 0.2, 12 percent thick, a = 1.0; 65(2)-415 a=0.5 names the low-drag range 2 (also accepted as 652-415 or 65_2-415) and the mean line a = 0.5; values that are not whole are written 64-(1.5)(13). The low-drag range is not used, because the thickness form depends only on the family and the thickness. Without the hyphen, 63012 is the five-digit section. Not built in: the 6A-series, the 67 family and the earlier individually derived sections that Report 824 writes with a comma (65,3-018); paste their coordinates.

Pasted coordinates. A list in Selig order (one loop from the trailing edge round the leading edge) or Lednicer order (upper surface, then lower surface, each from the leading edge) is recognized, shifted and scaled to a unit chord with the leading edge at the point farthest from the trailing edge, splined in arc length and re-paneled.

Paneling. NACA sections are paneled at the chord stations x = (1 − cos β)/2 with β equally spaced, half of the panels on each surface, which crowds them at both edges; pasted sections get the same spacing in arc length.

Trailing edge. The NACA thickness form ends with a small thickness (0.25 % of the chord for a 12 % section), and the panel methods need a closed body with a sharp edge. The default takes the gap out as a wedge: half of it is removed from each surface in proportion to x, which keeps the trailing-edge angle and thins the section slightly toward the rear (by 0.12 % of the chord at the trailing edge of a 12 % section). The second choice replaces the coefficient 0.1015 by 0.1036, which makes the thickness equation end at zero. Simply joining the two edge points was tried and is not offered: with cosine spacing the last panel of each surface then becomes a blunt base whose angle grows with the panel count, and both methods are sensitive to it.

Plain flap. The nodes behind the hinge line are rotated about a hinge point midway between the surfaces. A node is placed on the hinge line on each surface; on the side that opens a node is added on the arc swept by the hinge-line point, and on the side that closes the rotated nodes that would fall ahead of the hinge line are dropped. The gap is sealed.

The two panel methods

For a panel from z₁ to z₂ of length S and direction θ, write ζ = (z − z₁) e^−iθ for the field point in the panel's own axes and L = ln(ζ/(ζ − S)). Integrating the point source and the point vortex along the panel gives the complex velocity u − i v in the panel axes:

constant source q: (q/2π) L;    constant vortex γ: −(iγ/2π) L;    vortex varying linearly from γ₁ to γ₂: −(i/2π) [γ₁ (L − f) + γ₂ f],   f = ζ L/S − 1.

At the panel's own midpoint, approached from outside, L = −iπ. These closed forms are checked in the test against numerical quadrature of the point singularities.

Linear-strength vortex panels (the default). The vortex strength varies linearly along each panel and is continuous at the nodes, so N panels have N + 1 unknown node strengths. The N equations of no flow through the panel midpoints are closed by the Kutta condition γ₀ + γ_N = 0 at the trailing-edge node. This is the formulation of Kuethe and Chow. As it stands it leaves the strength at the trailing-edge node itself almost undetermined, because the two edge panels nearly coincide and carry equal and opposite strengths there: the same NACA 2412 paneled in two ways gave −0.66 and +2.6 for it, and on a cusped trailing edge the lift stayed 1.7 % low at any panel count. The tool therefore adds two weak equations (weight 0.01) asking the edge-node strength to continue linearly from the two neighboring nodes of each surface, and solves the N + 3 equations in least squares. The lift and the pressures away from the last panel are unchanged to the fifth digit on sections with a finite trailing-edge angle, and the cusped Joukowski section comes within 0.03 % at 200 panels. With the interior of the body at rest the surface velocity equals the local vortex strength, which is the value used here (the mean of the two node strengths at each midpoint).

Hess-Smith. Each panel carries a constant source strength and all panels share one vortex strength: N + 1 unknowns again. The N flow-tangency equations are closed by the Kutta condition that the tangential velocities at the midpoints of the two trailing-edge panels are equal and opposite. The surface velocity is the sum of the free stream and all the induced velocities at each midpoint.

Both systems are linear in cos α and sin α. They are solved once with the two right-hand sides, and the surface velocity at any angle of attack is v_i(α) = V_c,i cos α + V_s,i sin α. The sweep, the ideal angle and the scan for the least suction peak therefore cost nothing more.

Results

Lift. The Kutta-Joukowski theorem gives c_l = 2Γ/(V c) from the circulation Γ, the integral of the vortex strength round the section. Since Γ = Γ_c cos α + Γ_s sin α, the lift is exactly c_l = a₀ sin(α − α_L0) with a₀ = 2 √(Γ_c² + Γ_s²)/c the slope at zero lift and tan α_L0 = −Γ_c/Γ_s.

Moment and center of pressure. The pressure on each panel, taken constant at its midpoint value, gives the force and the moment about the leading edge; the moment about the quarter chord is c_m,c/4 = c_m,LE + c_n/4 and the center of pressure is at x_cp/c = −c_m,LE/c_n. The same integration gives a lift, which should agree with the circulation value, and a drag, which should be zero: both are printed as measures of the paneling error.

Stagnation point and suction peak. The stagnation point is where the surface velocity changes sign near the nose; the suction peak is the lowest C_p on the surface.

Ideal angle of attack. At the ideal (design) angle the stagnation point is on the leading edge, so the flow divides there without turning round the nose; the tool finds it in closed form as the angle at which the velocity at the leading-edge node vanishes, and reports the lift there as the design lift coefficient. It also scans for the angle at which the largest surface velocity is smallest: the angle of least peak suction, which is the angle of the highest critical Mach number. For a symmetric section both are zero; for a cambered section they differ, because the mid-chord velocity of the upper surface, not the nose, sets the peak near the design condition.

Thin-airfoil theory beside. For comparison the tool evaluates the thin-airfoil integrals of the mean line, taken as straight between closely spaced stations so that the integrals are exact sums: α_L0 = −(1/π) ∫ (dz/dx)(cos θ − 1) dθ, c_m,c/4 = (π/4)(A₂ − A₁), ideal angle (1/π) ∫ (dz/dx) dθ, a₀ = 2π. The thickness raises the panel method's slope above 2π, by a factor of about 1 + 0.77 t/c for the NACA four-digit thickness forms, and the camber's lift and moment with it. For the NACA 2412 the panel zero-lift angle is −2.15° against −2.08°, and the moment about the quarter chord at zero angle of attack is −0.055 against −0.053. In potential flow the aerodynamic center of a thick section lies a little behind the quarter chord (0.262 of the chord for the 2412), so the quarter-chord moment grows with lift, to −0.061 at 4°; measured aerodynamic centers lie at or slightly ahead of the quarter chord (0.239 to 0.247 for the 2412 in NACA Report 824), because the boundary layer thickens toward the trailing edge.

Compressibility

Below the critical Mach number the incompressible pressure coefficient C_p0 is corrected by the Karman-Tsien rule,

C_p = C_p0 / [√(1 − M²) + (M²/(1 + √(1 − M²))) C_p0/2],

or, more simply, by the Prandtl-Glauert rule C_p0/√(1 − M²), which also scales the lift and moment coefficients. The flow first becomes sonic where the pressure is lowest, when C_p there equals the critical value for isentropic flow from the free stream,

C_p* = (2/(γM²)) {[(2 + (γ − 1) M²)/(γ + 1)]^{γ/(γ − 1)} − 1},   γ = 1.4.

The critical Mach number is the root of the two expressions for the minimum C_p0 of the section at the angle of attack in question. For C_p0 = −0.43 the Prandtl-Glauert rule gives 0.737 and the Karman-Tsien rule 0.723. Above the critical Mach number the corrections do not apply.

Boundary-layer estimate

The option runs an integral boundary-layer calculation along each surface from the stagnation point, on the inviscid surface velocity u(s) (in units of the free stream, s in chords) at the chord Reynolds number Re. It is an estimate without viscous-inviscid interaction: the boundary layer does not change the pressure distribution, there is no separation bubble model and no wake, so the profile drag is indicative at low lift and too low near the stall, and the lift is not corrected.

Laminar part (Thwaites). θ² = (0.45/(Re u⁶)) ∫₀^s u⁵ ds, with the pressure-gradient parameter λ = Re θ² du/ds and the shape factor H(λ) and friction parameter from the usual curve fits to Thwaites' correlation. The laminar layer separates at λ = −0.09.

Transition. At the first station where Michel's criterion Re_θ > 1.174 (1 + 22,400/Re_x) Re_x^0.46 is met, at laminar separation if that comes first (treated as a short bubble: the layer continues turbulent with the same momentum thickness), or at an entered trip station.

Turbulent part (Head). The momentum equation dθ/ds = c_f/2 − (H + 2)(θ/u) du/ds and Head's entrainment equation d(u θ H₁)/ds = u F(H₁), F = 0.0306 (H₁ − 3)^−0.6169, with H₁ = 3.3 + 0.8234 (H − 1.1)^−1.287 for H ≤ 1.6 and 3.3 + 1.5501 (H − 0.6778)^−3.064 above, and the Ludwieg-Tillmann friction law c_f = 0.246 × 10^{−0.678 H} Re_θ^−0.268, integrated by a fourth-order Runge-Kutta scheme from H = 1.4 at transition. Turbulent separation is flagged when H reaches 2.4.

Profile drag (Squire-Young). c_d = Σ 2 θ_e (u_e/V)^{(H_e + 5)/2} over the two surfaces, from the momentum thickness, shape factor and velocity at the end of each. The end station is the last panel midpoint ahead of 99 % of the chord, because the inviscid velocity falls to the stagnation value at a trailing edge of finite angle; the result changes by about 1 % if 99.9 % is used instead.

On a flat plate the laminar routine gives θ √Re = 0.669 (Blasius 0.664) and the turbulent routine a drag within 7 % of Schlichting's 0.455/(log Re)^2.58 between Reynolds numbers of 10⁶ and 10⁷. For the smooth NACA 0012 and 2412 at a Reynolds number of 6 million the estimate is 0.0060 and 0.0061 at zero angle of attack, against measured minimum drag coefficients of about 0.006 in NACA Report 824 (read from its charts), and about 0.008 with transition forced at 2 % of the chord against about 0.010 measured with standard leading-edge roughness.

Accuracy

The test script compares both methods with the exact solution for Karman-Trefftz sections (a circle mapped to an airfoil with a chosen trailing-edge angle) and for the Joukowski section, whose trailing edge is a cusp. The linear-vortex lift is within 0.11 % of the exact value at 100 panels and 0.03 % at 200 on all of them, the cusp included, and its error falls fourfold when the panel count doubles. The Hess-Smith lift is within 1.3 % at 100 panels and 0.6 % at 200 for a finite trailing-edge angle; on the cusp it is 2.5 % low at 200 panels and 1.9 % low at 300, so it does not reach 1 % there. On a circle the Hess-Smith velocities are exact to rounding and the linear-vortex node strengths are within 0.0005 of the exact surface speed at 60 panels. The surface velocities of the NACA 0010, 0012 and 0015 at zero angle of attack agree with the tables of NACA Report 824 (computed in 1945 by Theodorsen's method) within 0.012, 0.012 and 0.014 of the free-stream velocity behind 1 % of the chord with the default wedge closure (the largest differences at 95 % of the chord, where the closure has thinned the section; with the two edge points joined instead and the Hess-Smith method the agreement is 0.008, 0.008 and 0.014).

For the 6-series, the tabulated ordinates of four cambered sections in Report 824 (the 63(2)-415, 64(1)-212, 65(2)-415 with a = 0.5 and 66-206; 400 stations and ordinates) are reproduced within 0.0025 percent of the chord from the thickness tables and the mean-line equation, and the mean-line ordinates and slopes of the report's a = 1.0 and a = 0.5 tables within their printed rounding. The surface velocities of the 31 tabulated thickness forms at zero lift agree with the v/V columns of Report 824 within 0.011 of the free-stream velocity from 5 % to 95 % of the chord (0.002 root mean square) and within 0.014 at 2.5 %. Ahead of 2.5 % of the chord they do not: the differences reach 0.03 to 0.06 on most sections and 0.11 on the 63-series form of 21 percent. Four stations in the first 1.25 % of the chord and a leading-edge radius do not fix the nose shape closely enough for the velocity there, so the suction peak at an angle of attack, the minimum pressure when it lies at the nose, the critical Mach number derived from it and the ideal angle found from the stagnation point carry that uncertainty for a 6-series section; the lift, the moment and the pressures behind 2.5 % of the chord do not. At the ideal angle of its mean line a cambered 6-series section carries the design load: the velocity difference between the surfaces is within 10 % of the thin-airfoil value from 10 % to 80 % of the chord, and the lift coefficient is 4 % to 17 % above c_li because of the thickness (a real section carries less than the potential-flow value).

Limits

Potential flow has no stall and no drag: the lift keeps rising with angle of attack, and the results stand for attached flow, in practice up to about 8° to 10° on these sections and up to about 10° to 15° of plain-flap deflection. The lift-curve slope is above the measured one, because the boundary layer reduces the effective camber and thickness at the trailing edge. One closed body only: no slots, no multi-element sections, no blunt trailing edge. The compressibility rules are for subcritical flow.