Hostname: page-component-74656f9884-8r8nv Total loading time: 0 Render date: 2026-11-08T22:50:45.119Z Has data issue: false hasContentIssue false
  • English
  • Français

Optimal airfoils in the intermediate-Reynolds-number range

Published online by Cambridge University Press:  30 July 2026

Gleb Zhdanko Open the ORCID record for Gleb Zhdanko [Opens in a new window]  and
Dmitry Kolomenskiy Open the ORCID record for Dmitry Kolomenskiy [Opens in a new window]
Gleb Zhdanko
Affiliation:
Center for Materials Technologies, Skolkovo Institute of Science and Technology, Moscow, Russia
Dmitry Kolomenskiy*
Affiliation:
Center for Materials Technologies, Skolkovo Institute of Science and Technology, Moscow, Russia
*
Corresponding author: Dmitry Kolomenskiy, dkolom@gmail.com
Rights & Permissions [Opens in a new window]

Abstract

Content of image described in text.

We revisit a classical airfoil design problem: the search for shapes that maximise aerodynamic performance metrics, targeting the underexplored intermediate Reynolds-number regime between 1 and 3000, relevant to small animals and miniature vehicles. The problem is formally stated as the glide ratio or the endurance factor maximisation for Joukowski airfoil profiles and for more general airfoil shapes with adjustable position of the maximum camber, under steady inflow. It is solved numerically by a hybrid approach combining stochastic search and direct parameter sweep, and using a steady laminar Navier–Stokes solver based on conformal mapping and second-order finite-difference discretisation. Zero-thickness cambered airfoils are found to be globally optimal, within the Joukowski family, across the entire Reynolds-number range considered. The optimal angle of attack decreases monotonically with the Reynolds number ($\textit{Re}$), whereas the optimal camber varies non-monotonically, reaching a pronounced maximum near $\textit{Re} \approx 40{-}50$ before declining at higher $\textit{Re}$. At low Reynolds numbers ($\textit{Re} \lesssim 100$), a broad family of cambered shapes performs within a few per cent of the optimum, indicating weak sensitivity to geometrical parameters. In contrast, for $\textit{Re} \gtrsim 1000$, the performance landscape becomes sharply localised around a single preferred design, for which geometric refinement is critical.

JFM classification

Aerodynamics: Aerodynamics Biological Fluid Dynamics: Swimming/flying

Information

Type
JFM Rapids
Information
Journal of Fluid Mechanics , Volume 1040 , 10 August 2026 , A12
Check for updates
Copyright
© The Author(s), 2026. Published by Cambridge University Press

1. Introduction

Lift and propulsion mechanisms in fluids have been extensively studied across a wide range of the Reynolds numbers ( italic Re $\textit{Re}$ ). At high italic Re $\textit{Re}$ ( greater than 10 Superscript 4 $\gt 10^4$ ), typical of aeronautical applications, the literature is rich with theoretical, computational and experimental insights. Decades of research have established that cambered thin-plate airfoils ensure optimal lift-to-drag ratios for the range of italic Re $\textit{Re}$ up to 10 Superscript 5 $10^5$ . Low- italic Re $\textit{Re}$ regimes ( italic Re much less than 1 $\textit{Re} \ll 1$ ), typical of microbial swimmers, have also received substantial attention. Here, the Stokes flow approximation enables analytical solutions yielding explicit formulas for the fluid-dynamic forces. Intermediate italic Re $\textit{Re}$ values ( 1 less than or equivalent to italic Re less than or equivalent to 10 000 $1 \lesssim \textit{Re} \lesssim 10{\,}000$ ), prevalent in biological locomotion, have historically received less attention. This gap likely stems from analytical intractability (viscous and inertial effects compete without simplifying assumptions) and experimental challenges: neither traditional high- italic Re $\textit{Re}$ wind tunnels nor microfluidic set-ups can reliably cover this range of italic Re $\textit{Re}$ . Recent technological advances in bio-inspired micro air vehicles have spurred practical interest in this flow regime.

Studies on insect flight ( italic Re tilde 10 $\textit{Re} \sim 10$ – 1000 $1000$ ) often employ flat, cambered or even corrugated thin foils as wing models. Yet, systematic efforts to identify optimal airfoil shapes remain scarce, with most work prioritising kinematics over geometry. A notable exception is the shape-optimisation study by Srinath & Mittal (Reference Srinath and Mittal2009), which explored a broad parameter space of airfoil contours at fixed angles of attack left parenthesis alpha equals 4 Superscript ring $(\alpha = 4^\circ$ and 12 Superscript ring Baseline right parenthesis $12^\circ )$ using adjoint optimisation. However, the complex shapes obtained therein reflect the alpha $\alpha$ constraint rather than the necessity of non-zero thickness. Our present study addresses this problem through a low-dimensional parameterisation of airfoil profiles, enabling efficient exploration of the combined effect of main features of the shape (i.e. thickness and camber) and angle of attack across intermediate italic Re $\textit{Re}$ . By coupling a parametric geometrical description with Navier–Stokes simulations, we identify optimal airfoils within the Joukowski family, then more general airfoil shapes with adjustable position of the maximum camber. We relate the observed trends of the optimal parameters as functions of italic Re $\textit{Re}$ to the competing effect of camber-enhanced lift and inertia-driven separation.

2. Mathematical model

2.1. Governing equations

Airfoil optimisation requires flow solutions for many candidate profiles, making repeated mesh generation computationally expensive. We avoid this by using conformal mapping, which simultaneously defines the airfoil, generates an orthogonal grid and yields the metric terms for the transformed Navier–Stokes equations. The map is constructed in two stages, mapping the airfoil exterior to the computational domain via an intermediate annular domain (see figure 1). The starting point is a rectangular grid in left parenthesis r comma s right parenthesis $(r,s)$ coordinates, where r $r$ is the radial coordinate and s $s$ the periodic angular coordinate. The first map, z equals e Superscript a r plus normal i a s Baseline comma $ z = e^{ar + \mathrm{i}as},$ where the parameter a element of double struck upper R $a \in \mathbb{R}$ controls the grid stretching, deforms this rectangular grid into an O-grid in the complex z $z$ -plane: lines of constant r $r$ become concentric circles, and lines of constant s $s$ become radial rays. The second map, phi $\phi$ , from the annulus onto the airfoil exterior, x plus normal i y equals phi left parenthesis z right parenthesis comma $ x + \mathrm{i}y = \phi (z),$ deforms the O-grid into a body-fitted mesh, with the inner boundary conforming to the airfoil surface upper Gamma $\varGamma$ and the outer boundary mapped to the far-field boundary upper Gamma 1 $\varGamma _1$ . We will examine two sets of results using two different maps phi left parenthesis z right parenthesis $\phi (z)$ : Joukowski transform (§ 2.3) and Wegmann transform (§ 2.4). With the composite mapping defined, the incompressible steady Navier–Stokes equations in vorticity–streamfunction form can be rewritten in left parenthesis r comma s right parenthesis $(r,s)$ coordinates. Denoting upper A equals normal upper R normal e left parenthesis phi prime left parenthesis z right parenthesis a z right parenthesis $A = \mathrm{Re}(\phi '(z)az)$ and upper B equals normal upper I normal m left parenthesis phi prime left parenthesis z right parenthesis a z right parenthesis $B = \mathrm{Im}(\phi '(z)az)$ , the Jacobian of the transformation left parenthesis r comma s right parenthesis right arrow left parenthesis x comma y right parenthesis $(r,s)\to (x,y)$ is sans serif bold italic upper J equals Start 2 By 2 Matrix 1st Row 1st Column upper A 2nd Column negative upper B 2nd Row 1st Column upper B 2nd Column upper A EndMatrix comma $ \unicode{x1D645}= \big(\begin{smallmatrix} A & -B \\ B & A \end{smallmatrix}\big),$ and the governing equations are written in computational coordinates

(2.1) StartLayout 1st Row StartFraction 1 Over det sans serif bold italic upper J EndFraction left parenthesis StartFraction partial differential psi Over partial differential s EndFraction StartFraction partial differential omega Over partial differential r EndFraction minus StartFraction partial differential psi Over partial differential r EndFraction StartFraction partial differential omega Over partial differential s EndFraction right parenthesis equals StartFraction 1 Over italic Re EndFraction StartFraction 1 Over det sans serif bold italic upper J EndFraction nabla Subscript r comma s Superscript 2 Baseline omega comma EndLayout \begin{eqnarray} \frac {1}{\det \unicode{x1D645}}\left (\frac {\partial \psi }{\partial s}\frac {\partial \omega }{\partial r} -\frac {\partial \psi }{\partial r}\frac {\partial \omega }{\partial s}\right ) = \frac {1}{\textit{Re}}\frac {1}{\det \unicode{x1D645}}{\nabla} ^2_{r,s}\omega, \end{eqnarray}
(2.2) StartLayout 1st Row omega equals minus StartFraction 1 Over det sans serif bold italic upper J EndFraction nabla Subscript r comma s Superscript 2 Baseline psi comma EndLayout \begin{eqnarray} \omega = -\frac {1}{\det \unicode{x1D645}}{\nabla} ^2_{r,s}\psi, \end{eqnarray}

where omega $\omega$ is the vorticity, psi $\psi$ is the streamfunction, italic Re equals rho upper U Subscript normal infinity Baseline c divided by mu $\textit{Re} = \rho U_\infty c / \mu$ , rho $\rho$ is the fluid density, upper U Subscript normal infinity $U_\infty$ the free-stream velocity, c $c$ the chord length and mu $\mu$ the dynamic viscosity.

Figure 1.

The three-stage conformal mapping: (a) rectangular computational domain in left parenthesis r comma s right parenthesis $(r,s)$ coordinates; (b) intermediate annular domain in the z $z$ -plane; (c) physical domain with the body-fitted grid to the airfoil surface. Inner and outer boundaries are shown in red and blue, respectively.

A diagram illustrating a three-stage conformal mapping process with three panels showing different domains and boundaries.
Figure 1. Long description

Panel A: A rectangular computational domain in coordinates. The domain is divided into a grid with an inner boundary marked in red, an outer boundary marked in blue, and a periodic cut marked in green. Panel B: An intermediate annular domain in the zeta-plane. The domain is circular with a unit circle marked in red and an outer boundary marked in blue. Panel C: A physical domain with a body-fitted grid to the airfoil surface. The domain includes an airfoil marked in red and an outer boundary marked in blue.

2.2. Boundary conditions

On the airfoil surface upper Gamma $\varGamma$ , the no-penetration and no-slip conditions are partial differential psi divided by partial differential bold italic n equals 0 $\partial \psi /\partial \boldsymbol{n} = 0$ and partial differential psi divided by partial differential bold italic tau equals 0 $\partial \psi /\partial \boldsymbol{\tau } = 0$ , where bold italic n $\boldsymbol{n}$ and bold italic tau $\boldsymbol{\tau }$ denote the unit normal and unit tangent vector to upper Gamma $\varGamma$ , respectively. Since conformal maps preserve orthogonality, the conditions retain their form in left parenthesis r comma s right parenthesis $(r,s)$ coordinates. From the boundary condition partial differential psi divided by partial differential s equals 0 $\partial \psi / \partial s = 0$ on upper Gamma $\varGamma$ , it follows that psi $\psi$ is uniform on the airfoil surface. Since psi $\psi$ is defined only up to an additive constant, this value may be chosen as zero, reducing the boundary conditions to psi equals 0 $\psi = 0$ and partial differential psi divided by partial differential r equals 0 $\partial \psi / \partial r = 0$ on upper Gamma $\varGamma$ .

At the outer boundary upper Gamma 1 $\varGamma _1$ , uniform horizontal free-stream conditions are imposed, which in terms of the streamfunction read partial differential psi divided by partial differential y equals 1 comma partial differential psi divided by partial differential x equals 0 $\partial \psi /\partial y = 1,\ \partial \psi /\partial x = 0$ on upper Gamma 1 $\varGamma _1$ . Applying the chain rule for the conformal transformation yields the outer boundary conditions

(2.3a,b) StartLayout 1st Row psi Subscript r Baseline equals upper B comma psi Subscript s Baseline equals upper A on upper Gamma 1 period EndLayout \begin{align} \psi _r = B, \quad \psi _s = A \quad \text{on } \quad \varGamma _1. \end{align}

At the continuous level, both conditions (2.3) are valid and equivalent representations of the imposed far-field velocity. However, in the discrete formulation, prescribing both streamfunction derivatives simultaneously leads to an ill-conditioned system when coupled with the Poisson equation (2.2). Condition (2.3a ) is retained as it directly controls the mass flux through the outer boundary, whereas condition (2.3b ) is redundant given the irrotationality of the undisturbed free stream and therefore dropped in favour of the vorticity condition omega equals 0 $\omega = 0$ . This substitution restores well-posedness of the discrete system without altering the physical model, yielding the system

(2.4) StartLayout 1st Row 1st Column left parenthesis StartFraction partial differential psi Over partial differential s EndFraction StartFraction partial differential omega Over partial differential r EndFraction minus StartFraction partial differential psi Over partial differential r EndFraction StartFraction partial differential omega Over partial differential s EndFraction right parenthesis 2nd Column equals StartFraction 1 Over italic Re EndFraction nabla Subscript r comma s Superscript 2 Baseline omega comma 2nd Row 1st Column omega 2nd Column equals minus StartFraction 1 Over det sans serif bold italic upper J EndFraction nabla Subscript r comma s Superscript 2 Baseline psi comma EndLayout with upper B period upper C period StartLayout 1st Row 1st Column psi 2nd Column equals 0 comma StartFraction partial differential psi Over partial differential r EndFraction equals 0 3rd Column Blank 4th Column on upper Gamma comma 2nd Row 1st Column StartFraction partial differential psi Over partial differential r EndFraction 2nd Column equals upper B comma omega equals 0 3rd Column Blank 4th Column on upper Gamma 1 period EndLayout \begin{equation} \begin{aligned} \left (\frac {\partial \psi }{\partial s}\frac {\partial \omega }{\partial r} - \frac {\partial \psi }{\partial r}\frac {\partial \omega }{\partial s}\right ) & = \frac {1}{\textit{Re}}{\nabla} ^2_{r,s}\omega, \\ \omega & = -\frac {1}{\det \unicode{x1D645}}{\nabla} ^2_{r,s}\psi, \end{aligned} \quad \text{with B.C.}\quad \begin{aligned} \psi & = 0, \quad \frac {\partial \psi }{\partial r} = 0 & & \text{on }\varGamma, \\ \frac {\partial \psi }{\partial r} & = B, \quad \omega = 0 & & \text{on }\varGamma _1. \end{aligned} \end{equation}

2.3. Joukowski conformal map

The Joukowski map is defined on the exterior of the unit disk ( double struck upper C minus upper D 1 $\mathbb{C}\setminus D_1$ ) by analytic formula ModifyingAbove phi left parenthesis z right parenthesis With tilde equals z plus z 0 plus k squared divided by left parenthesis z plus z 0 right parenthesis $\widetilde {\phi (z)} = z + z_0 + k^2/(z + z_0)$ , where three real-valued parameters fully determine the airfoil shape: k greater than or equals 0 $k \geq 0$ controls the thickness, mu Subscript x Baseline equals normal upper R normal e left parenthesis z 0 right parenthesis $\mu _x = \mathrm{Re}(z_0)$ controls the trailing-edge geometry and mu Subscript y Baseline equals normal upper I normal m left parenthesis z 0 right parenthesis $\mu _y = \mathrm{Im}(z_0)$ introduces camber. Despite this low-dimensional parameterisation, Joukowski airfoils represent a wide variety of streamlined shapes. This is an analytical mapping, therefore its numerical evaluation is computationally highly efficient and precise.

Since conformality requires a non-zero derivative throughout the domain, constraints on k $k$ and z 0 $z_0$ must be imposed. For a conformal map from the exterior of the unit disk to the exterior of the airfoil, two conditions must hold: the shifted centre z 0 $z_0$ must lie inside the unit disk, ensuring that the denominator z plus z 0 $z + z_0$ is non-vanishing in double struck upper C minus upper D 1 $\mathbb{C}\setminus D_1$ , and ModifyingAbove phi prime left parenthesis z right parenthesis With tilde not equals 0 $\widetilde {\phi '(z)} \neq 0$ throughout double struck upper C minus upper D 1 $\mathbb{C}\setminus D_1$ .

Preliminary computations over the full parameter space, consistent with Sunada et al. (Reference Sunada, Sakaguchi and Kawachi1997), confirm that an airfoil with a cusped trailing edge ( mu Subscript x Baseline less than or equals 0 $\mu _x \leq 0$ ) and upward camber ( mu Subscript y Baseline greater than or equals 0 $\mu _y \geq 0$ ) achieves higher performance metrics; the parameter space is therefore restricted to

(2.5) mu Subscript x Baseline element of left bracket negative 1 comma 0 right bracket comma mu Subscript y Baseline element of left bracket 0 comma StartRoot 1 minus mu Subscript x Superscript 2 Baseline EndRoot right parenthesis comma k element of left bracket 0 comma StartRoot 1 minus mu Subscript y Superscript 2 Baseline EndRoot minus StartAbsoluteValue mu Subscript x Baseline EndAbsoluteValue right parenthesis period \begin{equation} \mu _x \in [-1,0],\quad \mu _y \in [0,\sqrt {1-\mu _x^2}\ ),\quad k \in [0, \sqrt {1-\mu _y^2} - |\mu _x| \ ). \end{equation}

The airfoil thickness decreases monotonically as k $k$ increases toward its upper bound; in the limit k right arrow StartRoot 1 minus mu Subscript y Superscript 2 Baseline EndRoot minus StartAbsoluteValue mu Subscript x Baseline EndAbsoluteValue $k \to \sqrt {1 - \mu _y^2} - |\mu _x|$ the airfoil degenerates to a flat plate.

Two further operations complete the map definition: normalisation to the unit chord and introduction of the angle of attack alpha $\alpha$ . Defining c equals StartAbsoluteValue z Subscript x Sub Subscript italic max Subscript Baseline minus z Subscript x Sub Subscript italic min Subscript Baseline EndAbsoluteValue $c = |z_{x_{\textit{max}}} - z_{x_{\textit{min}}}|$ as the distance between the leftmost and rightmost points, which coincides with the chord length for the thin airfoils considered here, the normalised map is ModifyingAbove phi left parenthesis z right parenthesis With caret colon equals ModifyingAbove phi left parenthesis z right parenthesis With tilde slash c period $ \widehat {\phi (z)} := \widetilde {\phi (z)} / c.$ The angle of attack alpha $\alpha$ is then introduced by composing ModifyingAbove phi With caret $\widehat {\phi }$ with a rigid rotation: phi left parenthesis z right parenthesis colon equals ModifyingAbove phi left parenthesis z right parenthesis With caret normal e Superscript i alpha Baseline period $ \phi (z) := \widehat {\phi (z)}\mathrm{e}^{{i}\alpha }.$

The composite normalised Joukowski map, with alpha $\alpha$ as an additional parameter, thus takes the form

(2.6) phi left parenthesis z right parenthesis equals StartFraction e Superscript i alpha Baseline Over StartAbsoluteValue z Subscript x Sub Subscript italic max Subscript Baseline minus z Subscript x Sub Subscript italic min Subscript Baseline EndAbsoluteValue EndFraction left parenthesis z plus z 0 plus StartFraction k squared Over z plus z 0 EndFraction right parenthesis comma \begin{equation} \phi (z) = \frac {e^{{i}\alpha }}{|z_{x_{\textit{max}}} - z_{x_{\textit{min}}}|} \left ( z + z_0 + \frac {k^2}{z + z_0} \right )\!, \end{equation}

with the parameter constraints given by (2.5).

2.4. Wegmann transform

Joukowski mapping provides control over the two most important geometric parameters: thickness and camber. However, it does not allow independent variation of the chordwise position of maximum camber p $p$ , which is known to be a third key parameter for thin airfoils. To investigate the effect of p $p$ , we generate the target profile from a four-parameter family and compute the corresponding map numerically following Wegmann (Reference Wegmann1986).

The airfoil is parametrised similarly to NACA profiles. The mean camber line y Subscript c Baseline left parenthesis x right parenthesis $y_c(x)$ , with maximum camber f $f$ at chordwise position p $p$ is defined as

(2.7) y Subscript c Baseline left parenthesis x right parenthesis equals StartLayout Enlarged left brace 1st Row 1st Column StartFraction f Over p squared EndFraction left parenthesis 2 p x minus x squared right parenthesis comma 2nd Column x less than p comma 2nd Row 1st Column StartFraction f Over left parenthesis 1 minus p right parenthesis squared EndFraction left parenthesis 1 minus 2 p plus 2 p x minus x squared right parenthesis comma 2nd Column x greater than or equals p comma EndLayout \begin{equation} y_c(x) = \begin{cases} \dfrac {f}{p^2}\bigl (2px - x^2\bigr ), & x \lt p, \\[9pt] \dfrac {f}{(1-p)^2}\bigl (1 - 2p + 2px - x^2\bigr ), & x \geq p, \end{cases} \end{equation}

and a symmetric thickness distribution y Subscript t Baseline left parenthesis x right parenthesis $y_t(x)$ is added normal to it. The thickness follows the NACA four-digit polynomial, with a factor left parenthesis 1 minus x right parenthesis Superscript q minus 1 $(1-x)^{q-1}$ controlling trailing-edge closure

(2.8) StartLayout 1st Row y Subscript t Baseline left parenthesis x right parenthesis equals StartFraction t Over 0.2 EndFraction left parenthesis 1 minus x right parenthesis Superscript q minus 1 Baseline left parenthesis 0.2969 StartRoot x EndRoot minus 0.1260 x minus 0.3516 x squared plus 0.2843 x cubed minus 0.1015 x Superscript 4 Baseline right parenthesis period EndLayout \begin{align} y_t(x) = \frac {t}{0.2}\,(1-x)^{q-1}\bigl (0.2969\sqrt {x} - 0.1260\,x - 0.3516\,x^2 + 0.2843\,x^3 - 0.1015\,x^4\bigr ). \end{align}

The exponent q $q$ controls the trailing-edge shape: q equals 1.5 $q=1.5$ gives a Joukowski-like cusp, q equals 1 $q=1$ the wedge of the standard NACA series. We fix t equals 0.01 $t = 0.01$ and q equals 1.5 $q = 1.5$ throughout, leaving left parenthesis f comma p comma alpha right parenthesis $(f, p, \alpha )$ as free parameters. The resulting profiles have the leading-edge radius of almost equals 0.011 percent sign ${\approx} 0.011\,\%$ of the chord length and the maximum thickness near the leading edge of almost equals 1.2 percent sign ${\approx} 1.2\,\%$ of the chord length. The map is then computed in two steps. First, the inverse Kármán–Trefftz map is computed

(2.9) StartFraction z minus 1 Over z plus 1 EndFraction equals left parenthesis StartFraction w minus 1 Over w plus 1 EndFraction right parenthesis Superscript 1 divided by n Baseline comma n equals 2 minus StartFraction tau Over pi EndFraction comma \begin{equation} \frac {z-1}{z+1} = \left (\frac {w-1}{w+1}\right )^{1/n}, \quad n = 2 - \frac {\tau }{\pi }, \end{equation}

with tau $\tau$ defining the trailing-edge wedge angle, which the algorithm measures from the surface tangents at the trailing edge of each profile. For the cusped trailing edges considered here ( q equals 1.5 $q = 1.5$ ), tau $\tau$ is small and the exponent is close to 2 $2$ . This map removes the trailing-edge corner and deforms the airfoil upper Gamma $\varGamma$ into a smooth near-circular contour eta left parenthesis s right parenthesis $\eta (s)$ parametrised by arclength s element of left bracket 0 comma 2 pi right parenthesis $s \in [0, 2\pi )$ . Second, Wegmann’s algorithm computes the exterior map g $g$ from the unit circle onto this near circle. The map must project each boundary point normal e Superscript i t $\mathrm{e}^{{i}t}$ to a point on the contour

(2.10) g left parenthesis normal e Superscript i t Baseline right parenthesis equals eta left parenthesis upper S left parenthesis t right parenthesis right parenthesis comma \begin{equation} g(\mathrm{e}^{{i}t}) = \eta (S(t)), \end{equation}

where the boundary correspondence upper S left parenthesis t right parenthesis $S(t)$ is the unknown. Requiring g $g$ to be analytic outside the unit circle turns (2.10) into a nonlinear equation for upper S left parenthesis t right parenthesis $S(t)$ , which Wegmann’s method solves by iteration. The converged map is represented as a Laurent series, g left parenthesis w right parenthesis equals m 0 w plus sigma summation Underscript k equals 0 Overscript normal infinity Endscripts c Subscript k Baseline w Superscript negative k $g(w) = m_0 w + \sum _{k=0}^{\infty } c_k\, w^{-k}$ , of which we retain the first 80 $80$ coefficients. Composing g $g$ with the forward Kármán–Trefftz transformation (2.9) results in the conformal map from the exterior of the unit circle directly onto the exterior of the airfoil, yielding phi left parenthesis z right parenthesis $\phi (z)$ analogous to the Joukowski map of § 2.3.

2.5. Numerical method

The nonlinear system (2.4) is solved by Newton–Raphson iteration with backtracking line search. Because convergence is sensitive to the initial guess, different initialisation procedures are used for low and high Reynolds numbers. For italic Re less than 200 $\textit{Re} \lt 200$ , the potential-flow solution is used as the initial guess. For higher italic Re $\textit{Re}$ , a continuation strategy is employed over an increasing sequence of Reynolds numbers, italic Re Subscript 1 Baseline less than italic Re Subscript 2 Baseline less than midline horizontal ellipsis less than italic Re Subscript italic target $\textit{Re}_1 \lt \textit{Re}_2 \lt {\cdots} \lt \textit{Re}_{\textit{target}}$ , with each converged solution used for the next step. As an example, to reach italic Re equals 1000 $\textit{Re} = 1000$ , the sequence 200 right arrow 300 right arrow 500 right arrow 750 right arrow 1000 $200 \to 300 \to 500 \to 750 \to 1000$ is used.

The linearised system is discretised with second-order finite differences on a uniform upper N Subscript r Baseline times upper N Subscript s Baseline equals 450 times 600 $N_r \times N_s = 450 \times 600$ grid with a periodic boundary condition in the circumferential direction. The outer boundary is placed at the distance of 114 c $114c$ and 78 c $78c$ away from the airfoil for italic Re less than 150 $\textit{Re} \lt 150$ and italic Re greater than or equals 150 $\textit{Re} \geq 150$ , respectively, where the larger domain compensates for the slower decay of velocity perturbations in the low- italic Re $\textit{Re}$ regime (Khalili & Liu Reference Khalili and Liu2017). The conformal stretching concentrates nodes near the wall: the first off-wall layer thickness ranges from 1.3 times 10 Superscript negative 4 Baseline c $1.3 \times 10^{-4}c$ near the trailing-edge cusp to 7.0 times 10 Superscript negative 3 Baseline c $7.0 \times 10^{-3}c$ along the remainder of the chord. Solver accuracy is confirmed in Appendix.

Once the vorticity–streamfunction system is solved, the lift and drag coefficients are recovered from the surface pressure upper P $P$ and wall shear stress tau Subscript w $\tau _w$ . We use the non-dimensionalisation rho Subscript normal infinity Baseline equals 1 $\rho _\infty = 1$ , upper U Subscript normal infinity Baseline equals 1 $U_\infty = 1$ , c equals 1 $c = 1$ , under which the dynamic pressure is q Subscript normal infinity Baseline equals 1 divided by 2 comma $q_\infty = 1/2,$ dynamic viscosity is mu equals 1 slash italic Re $\mu = 1/ \textit{Re}$ and the lift and drag coefficients reduce to upper C Subscript upper L Baseline equals 2 upper F Subscript upper L $C_L = 2F_L$ , upper C Subscript upper D Baseline equals 2 upper F Subscript upper D $C_D = 2F_D$ , where upper F Subscript upper L $F_L$ and upper F Subscript upper D $F_D$ are lift and drag forces.

The surface pressure satisfies normal d upper P divided by normal d s equals minus left parenthesis 1 slash italic Re right parenthesis left parenthesis partial differential omega divided by partial differential r right parenthesis vertical bar Subscript upper Gamma Baseline $\mathrm{d}P/\mathrm{d}s = -(1/ \textit{Re})(\partial \omega /\partial r)|_{\varGamma }$ , which is obtained from the tangential momentum equation at the wall, with gauge condition integral Underscript upper Gamma Endscripts upper P normal d s equals 0 $\int _\varGamma P\,\mathrm{d}s = 0$ . The wall shear stress is tau Subscript w Baseline equals omega StartRoot det bold italic upper J EndRoot slash italic Re $\tau _w = \omega \sqrt {\det \boldsymbol{J}}/ \textit{Re}$ . The lift and drag coefficients follow from integrating the pressure and shear-stress contributions over the airfoil surface, with the Jacobian components upper A $A$ and upper B $B$ entering through the coordinate transformation

(2.11) upper C Subscript upper D Baseline equals minus 2 integral Underscript upper Gamma Endscripts italic PA normal d s minus StartFraction 2 Over italic Re EndFraction integral Underscript upper Gamma Endscripts omega upper B normal d s comma upper C Subscript upper L Baseline equals minus 2 integral Underscript upper Gamma Endscripts italic PB normal d s plus StartFraction 2 Over italic Re EndFraction integral Underscript upper Gamma Endscripts omega upper A normal d s period \begin{equation} C_D = -2\int _\varGamma \textit{PA}\:\mathrm{d}s - \frac {2}{\textit{Re}}\int _\varGamma \omega B \:\mathrm{d}s, \quad C_L = -2\int _\varGamma \textit{PB}\:\mathrm{d}s + \frac {2}{\textit{Re}}\int _\varGamma \omega A \:\mathrm{d}s. \end{equation}

3. Optimisation

We optimise airfoil profiles over italic Re element of left bracket 1 comma 3000 right bracket $\textit{Re} \in [1, 3000]$ in two independent sets of numerical experiments to maximise two different performance metrics: upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ (in the literature termed as lift-to-drag ratio, or glide ratio) and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ (termed as endurance factor, or power factor). These are the two most commonly used objectives when designing airfoils for regimes that typically constitute most of the flight duration, e.g. cruise and loitering flight. Since direct search over the full four-parameter space left parenthesis k comma mu Subscript x Baseline comma mu Subscript y Baseline comma alpha right parenthesis $(k, \mu _x, \mu _y, \alpha )$ , as defined in (2.6), is computationally expensive, we progressively narrow the parameter domain as follows.

We begin by analysing symmetric airfoils: mu Subscript y Baseline equals 0 $\mu _y = 0$ , mu Subscript x Baseline equals 0 $\mu _x = 0$ . This reduction to a two-parameter problem ( k comma alpha $k,\alpha$ ) makes it feasible to construct the objective function over a structured grid and examine its global behaviour. The computed objective surfaces reveal that, within the class of symmetric airfoils, the highest performance is consistently achieved for zero-thickness airfoils, corresponding to a flat plate. To determine whether this tendency persists when camber and trailing-edge geometry are free, we proceed to the full four-parameter space of Joukowski foils, employing three complementary global search strategies in parallel: Sobol sequence sampling, covariance matrix adaptation evolution strategy and Bayesian optimisation. All three strategies consistently identify zero-thickness airfoils with mu Subscript x Baseline almost equals 0 $\mu _x \approx 0$ as optimal across the entire italic Re $\textit{Re}$ range, in agreement with existing results for low-Reynolds-number airfoils (Sunada et al. Reference Sunada, Sakaguchi and Kawachi1997). Accordingly, the parameter mu Subscript x $\mu _x$ is set to 0 $0$ and k $k$ is set to the zero-thickness limit for each mu Subscript y $\mu _y$ , as given by the upper bound in (2.5). This reduction implies k equals StartRoot 1 minus mu Subscript y Superscript 2 Baseline EndRoot $k = \sqrt {1 - \mu _y^2}$ (in practice, k $k$ is set to StartRoot 1 minus mu Subscript y Superscript 2 Baseline EndRoot minus epsilon $\sqrt {1-\mu _y^2} - \varepsilon$ with epsilon equals 10 Superscript negative 7 $\varepsilon = 10^{-7}$ to avoid the map singularity at the boundary of the admissible domain; this has no measurable effect on the aerodynamic coefficients), leaving mu Subscript y $\mu _y$ as the only geometric parameter of the airfoil. For a thin Joukowski profile, mu Subscript y $\mu _y$ is related to the camber f $f$ through the image of the top of the circle ( z equals normal i $z = \mathrm{i}$ )

(3.1) f equals normal upper I normal m left parenthesis normal i plus mu Subscript y Baseline normal i plus left parenthesis 1 minus mu Subscript y Superscript 2 Baseline right parenthesis divided by left parenthesis normal i plus mu Subscript y Baseline normal i right parenthesis right parenthesis divided by c equals mu Subscript y Baseline divided by left parenthesis 2 k right parenthesis comma \begin{equation} f = \mathrm{Im} \big( \mathrm{i} + \mu _y\mathrm{i} + \big(1 - \mu _y^2 \big)/(\mathrm{i} + \mu _y\mathrm{i}) \big )/c = \mu_y/(2k), \end{equation}

which inverts to give mu Subscript y Baseline equals 2 f divided by StartRoot 1 plus 4 f squared EndRoot period $\mu _y = 2f/\sqrt {1 + 4f^2}.$ The optimisation thus reduces to a two-parameter problem in left parenthesis f comma alpha right parenthesis $(f, \alpha )$ space. A structured high-resolution grid search is performed over this plane with steps normal upper Delta alpha equals 0.5 Superscript ring $\Delta \alpha = 0.5^\circ$ and normal upper Delta f equals 0.5 percent sign $\Delta f = 0.5\,\%$ , within windows centred on the regions identified by the global search.

4. Results and discussion

Let us denote as alpha Subscript italic opt comma upper G $\alpha _{\textit{opt},G}$ and f Subscript italic opt comma upper G $f_{\textit{opt},G}$ , respectively, the angle of attack and camber that maximise upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ . Similarly, alpha Subscript italic opt comma upper E $\alpha _{\textit{opt},E}$ and f Subscript italic opt comma upper E $f_{\textit{opt},E}$ will stand for the values that maximise upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ . We will also use shorter notations alpha Subscript italic opt $\alpha _{\textit{opt}}$ and f Subscript italic opt $f_{\textit{opt}}$ in statements that apply to both maximisation objectives. The results of the high-resolution grid search optimisation are presented in figure 2. Solid lines correspond to the Joukowski airfoils. Point markers correspond to the extended optimal shapes obtained with the position of the maximum camber p $p$ as an additional optimisation parameter, by employing a Wegmann transform. For both airfoil families and both objective functions, the results are qualitatively the same and quantitatively close. Let us first discuss the results obtained for Joukowski airfoils. As italic Re $\textit{Re}$ increases from 1 $1$ to 3000 $3000$ , alpha Subscript italic opt $\alpha _{\textit{opt}}$ decreases monotonically from 39 Superscript ring $39^\circ$ to 6 Superscript ring $6^\circ$ , with a maximum difference between alpha Subscript italic opt comma upper E $\alpha _{\textit{opt},E}$ and alpha Subscript italic opt comma upper G $\alpha _{\textit{opt},G}$ of 3.2 Superscript ring $3.2^\circ$ at italic Re equals 50 $\textit{Re} = 50$ . The parameter f Subscript italic opt $f_{\textit{opt}}$ follows a bell-shape curve, with f Subscript italic opt comma upper G $f_{\textit{opt},G}$ and f Subscript italic opt comma upper E $f_{\textit{opt},E}$ behaving nearly identically at the boundaries of the italic Re $\textit{Re}$ range considered: starting and ending at 3.5 percent sign $3.5\,\%$ at italic Re equals 1 $\textit{Re} = 1$ and italic Re equals 3000 $\textit{Re} = 3000$ . Significant differences appear only near the peaks: f Subscript italic opt comma upper G Baseline equals 12.3 percent sign $f_{\textit{opt},G} = 12.3\,\%$ at italic Re equals 40 $\textit{Re} = 40$ and f Subscript italic opt comma upper E Baseline equals 15.2 percent sign $f_{\textit{opt},E} = 15.2\,\%$ at italic Re equals 50 $\textit{Re} = 50$ , and a maximum difference of 2.9 percent sign $2.9\,\%$ at italic Re equals 50 $\textit{Re} = 50$ .

Figure 2.

Results of optimisation of upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ (subscript upper G $G$ ) and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ (subscript upper E $E$ ). (a) Optimal angle of attack alpha Subscript italic opt $\alpha _{\textit{opt}}$ as a function of italic Re $\textit{Re}$ ; (b) optimal camber f Subscript italic opt $f_{\textit{opt}}$ as a function of italic Re $\textit{Re}$ ; (c) optimal chordwise position p Subscript italic opt $p_{\textit{opt}}$ of the maximum camber as a function of italic Re $\textit{Re}$ ; (d) maximum upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ as functions of italic Re $\textit{Re}$ . Solid blue and red lines and the letter J correspond to the results obtained for Joukowski airfoils. Green and yellow triangular markers and the letter W correspond to the results of calculations using the Wegmann method with the chordwise position p $p$ of the maximum camber treated as one of the optimisation parameters. Black dotted lines and letter M represent the asymptotic formula for the maximum upper C Subscript upper L Superscript asterisk Baseline divided by upper C Subscript upper D Superscript asterisk $C_L^*/C_D^*$ value and corresponding alpha Subscript o p t comma asterisk $\alpha _{opt,*}$ for a flat plate. The position of the maximum camber of thin Joukowski foils is constrained at p equals 50 percent sign $p=50\,\%$ and therefore not shown on the plot.

Multiple graphs depict optimization results for aerodynamic parameters.

When the position p $p$ of the maximum camber is treated as an adjustable parameter, for the flow regimes near italic Re equals 10 $\textit{Re} = 10$ , the optimal position of the maximum camber is near 50 % chordwise. However, for larger italic Re $\textit{Re}$ , optimal position of the maximum camber progressively shifts downstream toward the trailing edge. The cambers f $f$ are larger but the angles of attack alpha $\alpha$ are the same as for the Joukowski airfoils. The objective function values are slightly greater than for the Joukowski airfoils. The latter results, however, have been obtained with the fixed maximum thickness near the leading edge of 1.2 % of the chord length as required for the numerical mapping, and this value is considerably larger than the Joukowski airfoil thickness. It can be postulated that trends shown in figure 2 might be generally representative of optimal thin airfoils.

In the low-Reynolds-number regime, the Oseen approximation is valid, and the numerical results show that f Subscript italic opt Baseline right arrow 0 $f_{\textit{opt}} \to 0$ , meaning the optimal shape reduces to a flat plate. It is therefore natural to compare the numerical results against the analytical solution derived by Miyagi (Reference Miyagi1964) for the lift and drag of a flat plate under the Oseen approximation: upper C Subscript upper D Superscript asterisk Baseline equals left parenthesis 16 pi slash italic Re right parenthesis left parenthesis left parenthesis 2 upper S 1 minus 1 minus cosine 2 alpha right parenthesis divided by left parenthesis 4 upper S 1 squared minus 2 left parenthesis 1 plus cosine 2 alpha right parenthesis right parenthesis right parenthesis comma upper C Subscript upper L Superscript asterisk Baseline equals left parenthesis 16 pi slash italic Re right parenthesis left parenthesis left parenthesis sine 2 alpha right parenthesis divided by left parenthesis 4 upper S 1 squared minus 2 left parenthesis 1 plus cosine 2 alpha right parenthesis right parenthesis right parenthesis comma $ C_D^* = ({16\pi }/{\textit{Re}})(({2S_1 - 1 - \cos 2\alpha)}/{(4S_1^2 - 2(1+\cos 2\alpha)})), \quad C_L^* = ({16\pi }/{\textit{Re}}) (({\sin 2\alpha})/{(4S_1^2 - 2(1+\cos 2\alpha )})),$ where upper S 1 equals ln left parenthesis 16 slash italic Re right parenthesis minus gamma $S_1 = \ln (16/ \textit{Re}) - \gamma$ , with gamma $\gamma$ denoting the Euler–Mascheroni constant. Maximising upper C Subscript upper L Superscript asterisk Baseline divided by upper C Subscript upper D Superscript asterisk $C_L^*/C_D^*$ with respect to alpha $\alpha$ yields the optimal angle, alpha Subscript o p t comma asterisk Baseline equals arc tangent StartRoot 1 minus 1 divided by upper S 1 EndRoot $\alpha _{opt,*} = \arctan \sqrt {1 - 1/S_1}$ , which satisfies alpha Subscript o p t comma asterisk Baseline right arrow 45 Superscript ring $\alpha _{opt,*} \to 45^\circ$ as italic Re right arrow 0 $\textit{Re} \to 0$ , explaining the large optimal angles of attack observed at the lower end of the Reynolds-number range. Since upper S 1 $S_1$ is a decreasing function of italic Re $\textit{Re}$ , alpha Subscript o p t comma asterisk $\alpha _{opt,*}$ itself decreases monotonically with italic Re $\textit{Re}$ , consistent with the numerical results.

The non-monotonic behaviour of f Subscript italic opt $f_{\textit{opt}}$ arises from a competition between positive and negative aerodynamic effects. The benefits of non-zero camber are twofold. First, curvature enhances the inertial pressure differential between lower and upper surfaces. Second, these enhanced pressure forces near the leading edge are oriented vertically, contributing to upper C Subscript upper L $C_L$ to a greater extent. On the downside, greater pressure differential on a cambered foil entails larger adverse pressure gradient in the aft portion of it, which leads to flow separation and increased drag. As italic Re $\textit{Re}$ increases from the lower end, increased inertia boosts the curvature-induced excess upper C Subscript upper L $C_L$ . At the same time, viscous diffusion remains strong enough to mitigate flow separation even at large f $f$ and alpha $\alpha$ , permitting highly cambered airfoils to remain near optimal. Together, these effects drive f Subscript italic opt $f_{\textit{opt}}$ upward until italic Re almost equals 40 minus 50 $\textit{Re} \approx 40{-}50$ .

Beyond this peak, the trend reverses: as inertia becomes dominant, viscosity can no longer resist and the flow separates at large f $f$ and alpha $\alpha$ without reattachment. Both f Subscript italic opt $f_{\textit{opt}}$ and alpha Subscript italic opt $\alpha _{\textit{opt}}$ decline toward lower values. Flow detachment simultaneously reduces lift and increases drag, impairs both upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ ; any left parenthesis f comma alpha right parenthesis $(f, \alpha )$ combination at which the bubble opens is suboptimal. The optimal airfoils, therefore, reside in the attached-flow regime, which justifies the use of the steady solver in this study. At italic Re equals 100 $\textit{Re} = 100$ , the optimal cambered airfoil ( f Subscript italic opt comma upper G Baseline equals 11.2 percent sign comma alpha equals 16.6 Superscript ring $f_{\textit{opt},G} = 11.2\,\%, \, \alpha = 16.6^\circ$ ) already carries a trailing-edge recirculation bubble, visible in the streamlines (figure 3 b,c), yet the flow remains steady.  This airfoil achieves 17 percent sign $17\,\%$ higher upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ than the flat plate (figure 3 a,d), confirming that moderate camber remains beneficialeven as a separation bubble begins to form. For the optimal shape found at italic Re equals 3000 $\textit{Re}=3000$ ( f Subscript italic opt comma upper G Baseline equals 3.5 percent sign $f_{\textit{opt},G} = 3.5\,\%$ , alpha equals 6 Superscript ring $\alpha = 6 ^\circ$ ), this was confirmed by supplementary unsteady numerical simulations using OpenFOAM (see https://osf.io/y4pr8), in which the fluctuation amplitude did not exceed 0.2 % of the respective mean values. The Reynolds number based on the perpendicular projection of the chord, italic Re Subscript up tack Baseline equals italic Re sine alpha $\textit{Re}_\perp = \textit{Re} \sin \alpha$ , equals 314 in the optimal case at italic Re equals 3000 $\textit{Re} = 3000$ , and decreases with the decreasing italic Re $\textit{Re}$ , further supporting the steady flow assumption.

Figure 3.

(a) Lift, drag and lift-to-drag ratio curves as functions of alpha $\alpha$ at italic Re equals 100 $\textit{Re} = 100$ for the flat and optimal cambered profiles, with optimal alpha $\alpha$ and values of upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ marked. (b,c) Streamlines for the optimal cambered and flat profiles respectively, shown at levels 0 comma plus or minus 10 Superscript negative 5 Baseline comma plus or minus 5 times 10 Superscript negative 5 Baseline comma plus or minus 10 Superscript negative 4 Baseline comma plus or minus 5 times 10 Superscript negative 4 Baseline comma plus or minus 10 Superscript negative 3 $0,\pm 10^{-5}, \pm 5\times 10^{-5}, \pm 10^{-4}, \pm 5\times 10^{-4}, \pm 10^{-3}$ . (d) Pressure coefficient difference along the chord for both profiles.

Multiple graphs depict aerodynamic coefficients, streamlines, and pressure coefficient differences for flat and cambered profiles.
Figure 3. Long description

Panel A: A line graph shows aerodynamic coefficients as functions of the angle of attack for flat and cambered profiles. The x-axis represents the angle of attack in degrees, and the y-axis represents the aerodynamic coefficients. The graph includes multiple lines representing different coefficients for flat and cambered profiles, with optimal values marked by horizontal and vertical dashed lines. Panel B: A line graph displays streamlines for the optimal cambered profile at specific levels. The x-axis represents the x-coordinate, and the y-axis represents the y-coordinate. Panel C: A line graph shows streamlines for the flat profile at specific levels. The x-axis represents the x-coordinate, and the y-axis represents the y-coordinate. Panel D: A line graph illustrates the pressure coefficient difference along the chord for both flat and cambered profiles. The x-axis represents the x-coordinate, and the y-axis represents the pressure coefficient difference. The graph includes lines for both profiles, with a legend indicating which line corresponds to each profile.

At low italic Re $\textit{Re}$ , strong viscous diffusion spreads the upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ maxima broadly across the left parenthesis f comma alpha right parenthesis $(f, \alpha )$ space; as inertia strengthens with italic Re $\textit{Re}$ , the maxima sharpen, making a geometrically precise optimal airfoil identifiable. This is illustrated in figure 4, which shows contour plots of upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ at italic Re equals 10 comma 100 $\textit{Re} = 10, 100$ and 500 $500$ . All panels span equal ranges delta f times delta alpha equals 7 percent sign times 10 Superscript ring $\delta f \times \delta \alpha = 7\,\% \times 10^\circ$ , centred on left parenthesis f Subscript italic opt Baseline comma alpha Subscript italic opt Baseline right parenthesis $(f_{\textit{opt}}, \alpha _{\textit{opt}})$ at each italic Re $\textit{Re}$ , enabling direct comparison of landscape sharpness. At italic Re equals 10 $\textit{Re} = 10$ , upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ varies by less than 5 percent sign $5\,\%$ across the displayed range, meaning that a wide number of left parenthesis f comma alpha right parenthesis $(f, \alpha )$ combinations yields nearly identical performance. As italic Re $\textit{Re}$ increases, the maximum becomes more distinct: at italic Re equals 100 $\textit{Re} = 100$ the variation reaches 10 percent sign $10\,\%$ , and at italic Re equals 500 $\textit{Re} = 500$ it is 25 percent sign $25\,\%$ . This progressive sharpening has a direct implication: below italic Re almost equals 100 $\textit{Re} \approx 100$ , performance is insensitive to the precise choice of f $f$ and alpha $\alpha$ , whereas above italic Re almost equals 500 $\textit{Re} \approx 500$ , the sharp peak singles out a physically distinguished left parenthesis f comma alpha right parenthesis $(f, \alpha )$ combination, making optimisation both meaningful and effective.

Figure 4.

Counterplots of upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ normalised by its maximum value, in the left parenthesis f comma alpha right parenthesis $(f, \alpha )$ parameter space: (a) for italic Re equals 10 $\textit{Re} = 10$ ; (b) for italic Re equals 100 $\textit{Re} = 100$ ; (c) for italic Re equals 500 $\textit{Re} = 500$ . Stars mark the left parenthesis f Superscript asterisk Baseline comma alpha Superscript asterisk Baseline right parenthesis $(f^*, \alpha ^*)$ in each case. All panels span equal areas of the left parenthesis f comma alpha right parenthesis $(f, \alpha )$ parameter space.

Heat map showing the objective function values across the parameter space.

5. Conclusion

Optimisation of Joukowski airfoils across italic Re element of left bracket 1 comma 3000 right bracket $\textit{Re} \in [1, 3000]$ within the steady, laminar, incompressible Navier–Stokes model has revealed features and trends of optimal airfoils that are common to both the upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ maximisation in the intermediate-Reynolds-number range. Zero-thickness cambered airfoils are optimal for both upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ throughout the entire intermediate- italic Re $\textit{Re}$ range. This result is confirmed by three independent global search strategies. The optimal angle of attack decreases monotonically with italic Re $\textit{Re}$ , while the optimal camber follows a bell-shaped curve. At italic Re less than or equivalent to 10 $\textit{Re} \lesssim 10$ , near-flat plates ( f less than or equivalent to 11 percent sign $f \lesssim 11\,\%$ ) at high angle of attack ( alpha almost equals 35 minus 39 Superscript ring $\alpha \approx 35{-}39^\circ$ ) are optimal. As italic Re $\textit{Re}$ increases, optimal camber attains a peak of 12.3 percent sign $12.3\,\%$ ( upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ ) and 15.2 percent sign $15.2\,\%$ ( upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ ) at italic Re almost equals 40 minus 50 $\textit{Re} \approx 40{-}50$ , marking the transition from viscosity-dominated to inertia-dominated flow. At italic Re equals 3000 $\textit{Re} = 3000$ , both f $f$ and alpha $\alpha$ converge toward the slightly cambered ( tilde 4 percent sign ${\sim }4\,\%$ ) low- alpha $\alpha$ ( tilde 6 Superscript ring ${\sim }6^\circ$ ) geometry, in close agreement with the optimal airfoils identified by Sunada et al. (Reference Sunada, Sakaguchi and Kawachi1997) at italic Re equals 4000 $\textit{Re} = 4000$ and consistent with the slightly cambered geometry characteristic of conventional low-speed airfoils. Another useful property is the sharpness of the upper C Subscript upper L Baseline divided by upper C Subscript upper D $C_L/C_D$ and upper C Subscript upper L Superscript 3 divided by 2 Baseline divided by upper C Subscript upper D $C_L^{3/2}/C_D$ maxima in left parenthesis f comma alpha right parenthesis $(f,\alpha )$ space. Below italic Re almost equals 100 $\textit{Re} \approx 100$ , geometric precision is aerodynamically unimportant: performance varies by less than 10 percent sign $10\,\%$ across a broad neighbourhood of the optimum. Above italic Re almost equals 500 $\textit{Re} \approx 500$ , the objective functions have a distinct optimum, with a 25 percent sign $25\,\%$ change in performance across a neighbourhood of the same size. The optimisation then reveals a distinct optimal shape, clearly superior to neighbouring designs in terms of aerodynamic performance. After extending the parameter space by including the chordwise position p $p$ of the maximum camber as an adjustable parameter, the trends of optimal alpha $\alpha$ and f $f$ have remained the same, although the values of f Subscript italic opt $f_{\textit{opt}}$ have increased by a few per cent. In addition, these optimisation runs have revealed a trend for the optimal chordwise position of the maximum camber to shift, with increasing italic Re $\textit{Re}$ , towards the trailing edge, which has not been reported in earlier studies.

Declaration of interests

The authors report no conflict of interest.

Data availability statement

The solver source code is accessible on GitHub. Datasets and other supplementary information can be downloaded from https://osf.io/y4pr8.

AI statement

AI tools have been used in the writing process.

Appendix. Solver verification

The solver is implemented in Python and is validated against three benchmark cases: flow past a circular cylinder, flow past a thin flat plate and flow past a NACA 0012 airfoil. Figure 5(a) shows drag coefficients for the cylinder at italic Re equals 7 $\textit{Re} = 7$ – 40 $40$ in close agreement with the solutions of Dennis & Chang (Reference Dennis and Chang1970), Nieuwstadt & Keller (Reference Nieuwstadt and Keller1973) and Dennis (Reference Dennis1976). Figure 5(b,e) compares lift and drag coefficients for the flat plate at italic Re equals 5 comma 10 comma 20 $\textit{Re} = 5, 10, 20$ with results of In et al. (Reference In, Choi and Kim1995), showing good agreement across all alpha $\alpha$ . Figure 5(c, f) compares lift and drag coefficients for the NACA 0012 airfoil at italic Re equals 10 comma 100 comma 250 comma 500 $\textit{Re} = 10, 100, 250, 500$ with results of Srinath & Mittal (Reference Srinath and Mittal2009). The coefficients agree closely at italic Re equals 100 comma 250 comma 500 $\textit{Re} = 100, 250, 500$ for all alpha $\alpha$ ; at italic Re equals 10 $\textit{Re} = 10$ the agreement is weaker, which we attribute to the comparatively small computational domain (10 chord lengths) used by Srinath & Mittal (Reference Srinath and Mittal2009), insufficient to fully capture the far-field decay of the flow at this low Reynolds number. Mesh convergence is demonstrated in figure 5(d) for the cylinder at italic Re equals 40 $\textit{Re} = 40$ , confirming second-order spatial accuracy.

Figure 5.

Solver validation. (a) Drag coefficient for flow past a circular cylinder at italic Re equals 7 comma 10 comma 20 comma 40 $\textit{Re} = 7, 10, 20, 40$ , compared with Dennis & Chang (Reference Dennis and Chang1970), Nieuwstadt & Keller (Reference Nieuwstadt and Keller1973) and Dennis (Reference Dennis1976). (b,e) Lift and drag coefficients versus angle of attack alpha $\alpha$ for a thin flat plate at italic Re equals 5 comma 10 comma 20 $\textit{Re} = 5, 10, 20$ , compared with In et al. (Reference In, Choi and Kim1995). (c, f) Lift and drag coefficients versus alpha $\alpha$ for NACA0012 at italic Re equals 10 comma 100 comma 250 comma 500 $\textit{Re} = 10, 100, 250, 500$ , compared with Srinath & Mittal (Reference Srinath and Mittal2009). (d) Mesh convergence at italic Re equals 40 $\textit{Re} = 40$ , showing second-order behaviour.

Multiple graphs depict solver validation for flow past different shapes, comparing drag and lift coefficients across various studies.
Figure 5. Long description

Panel A: A line graph shows the drag coefficient for flow past a circular cylinder at different Reynolds numbers, compared with data from Dennis & Chang, Nieuwstadt & Keller, and Dennis. The x-axis represents the Reynolds number, and the y-axis represents the drag coefficient. Panel B: A line graph displays the lift coefficient versus the angle of attack for a thin flat plate at different Reynolds numbers, compared with data from In et al. The x-axis represents the angle of attack in degrees, and the y-axis represents the lift coefficient. Panel C: A line graph shows the lift coefficient versus the angle of attack for a NACA 0012 airfoil at different Reynolds numbers, compared with data from Srinath & Mittal. The x-axis represents the angle of attack in degrees, and the y-axis represents the lift coefficient. Panel D: A scatter plot with a fitted line shows mesh convergence at a specific Reynolds number, indicating second-order behavior. The x-axis represents the number of mesh points, and the y-axis represents the difference in drag coefficients. Panel E: A line graph displays the drag coefficient versus the angle of attack for a thin flat plate at different Reynolds numbers, compared with data from In et al. The x-axis represents the angle of attack in degrees, and the y-axis represents the drag coefficient. Panel F: A line graph shows the drag coefficient versus the angle of attack for a NACA 0012 airfoil at different Reynolds numbers, compared with data from Srinath & Mittal. The x-axis represents the angle of attack in degrees, and the y-axis represents the drag coefficient.

References

Dennis, S.C.R. 1976 A numerical method for calculating steady flow past a cylinder. In Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics, (ed. A.I. Van De Vooren & P.J. Zandbergen), Lecture Notes in Physics, vol. 59, pp. 165–172. Springer.Google Scholar
Dennis, S.C.R. & Chang, G.Z. 1970 Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J. Fluid Mech. 42, 471–489.10.1017/S0022112070001428CrossRefGoogle Scholar
In, K.M., Choi, D.H. & Kim, M.-U. 1995 Two-dimensional viscous flow past a flat plate. Fluid Dyn. Res. 15 (1), 13–24.10.1016/0169-5983(95)90438-8CrossRefGoogle Scholar
Khalili, A. & Liu, B. 2017 Stokes’ paradox: creeping flow past a two-dimensional cylinder in an infinite domain. J. Fluid Mech. 817, 374–387.10.1017/jfm.2017.127CrossRefGoogle Scholar
Miyagi, T. 1964 Oseen flow past a flat plate inclined to the uniform stream. J. Phys. Soc. Japan 19 (6), 1063–1073.10.1143/JPSJ.19.1063CrossRefGoogle Scholar
Nieuwstadt, F. & Keller, H.B. 1973 Viscous flow past circular cylinders. Comput. Fluids 1 (1), 59–71.10.1016/0045-7930(73)90026-1CrossRefGoogle Scholar
Srinath, D.N. & Mittal, S. 2009 Optimal airfoil shapes for low Reynolds number flows. Intl J. Numer. Meth. Fluids 61 (4), 355–381.10.1002/fld.1960CrossRefGoogle Scholar
Sunada, S., Sakaguchi, A. & Kawachi, K. 1997 Airfoil section characteristics at a low Reynolds number. J. Fluids Engng 119 (1), 129–135.10.1115/1.2819098CrossRefGoogle Scholar
Wegmann, R. 1986 An iterative method for conformal mapping. J. Comput. Appl. Maths 14 (1–2), 7–18.10.1016/0377-0427(86)90128-7CrossRefGoogle Scholar
Figure 0

Figure 1. Figure 1 long description.The three-stage conformal mapping: (a) rectangular computational domain in (r,s)$(r,s)$ coordinates; (b) intermediate annular domain in the z$z$-plane; (c) physical domain with the body-fitted grid to the airfoil surface. Inner and outer boundaries are shown in red and blue, respectively.

Figure 1

Figure 2. Results of optimisation of CL/CD$C_L/C_D$ (subscript G$G$) and CL3/2/CD$C_L^{3/2}/C_D$ (subscript E$E$). (a) Optimal angle of attack αopt$\alpha _{\textit{opt}}$ as a function of Re$\textit{Re}$; (b) optimal camber fopt$f_{\textit{opt}}$ as a function of Re$\textit{Re}$; (c) optimal chordwise position popt$p_{\textit{opt}}$ of the maximum camber as a function of Re$\textit{Re}$; (d) maximum CL/CD$C_L/C_D$ and CL3/2/CD$C_L^{3/2}/C_D$ as functions of Re$\textit{Re}$. Solid blue and red lines and the letter J correspond to the results obtained for Joukowski airfoils. Green and yellow triangular markers and the letter W correspond to the results of calculations using the Wegmann method with the chordwise position p$p$ of the maximum camber treated as one of the optimisation parameters. Black dotted lines and letter M represent the asymptotic formula for the maximum CL∗/CD∗$C_L^*/C_D^*$ value and corresponding αopt,∗$\alpha _{opt,*}$ for a flat plate. The position of the maximum camber of thin Joukowski foils is constrained at p=50%$p=50\,\%$ and therefore not shown on the plot.

Figure 2

Figure 3. Figure 3 long description.(a) Lift, drag and lift-to-drag ratio curves as functions of α$\alpha$ at Re=100$\textit{Re} = 100$ for the flat and optimal cambered profiles, with optimal α$\alpha$ and values of CL/CD$C_L/C_D$ marked. (b,c) Streamlines for the optimal cambered and flat profiles respectively, shown at levels 0,±10−5,±5×10−5,±10−4,±5×10−4,±10−3$0,\pm 10^{-5}, \pm 5\times 10^{-5}, \pm 10^{-4}, \pm 5\times 10^{-4}, \pm 10^{-3}$. (d) Pressure coefficient difference along the chord for both profiles.

Figure 3

Figure 4. Counterplots of CL/CD$C_L/C_D$ normalised by its maximum value, in the (f,α)$(f, \alpha )$ parameter space: (a) for Re=10$\textit{Re} = 10$; (b) for Re=100$\textit{Re} = 100$; (c) for Re=500$\textit{Re} = 500$. Stars mark the (f∗,α∗)$(f^*, \alpha ^*)$ in each case. All panels span equal areas of the (f,α)$(f, \alpha )$ parameter space.

Figure 4

Figure 5. Figure 5 long description.Solver validation. (a) Drag coefficient for flow past a circular cylinder at Re=7,10,20,40$\textit{Re} = 7, 10, 20, 40$, compared with Dennis & Chang (1970), Nieuwstadt & Keller (1973) and Dennis (1976). (b,e) Lift and drag coefficients versus angle of attack α$\alpha$ for a thin flat plate at Re=5,10,20$\textit{Re} = 5, 10, 20$, compared with In et al. (1995). (c, f) Lift and drag coefficients versus α$\alpha$ for NACA0012 at Re=10,100,250,500$\textit{Re} = 10, 100, 250, 500$, compared with Srinath & Mittal (2009). (d) Mesh convergence at Re=40$\textit{Re} = 40$, showing second-order behaviour.