“With four parameters I can fit an elephant, and with five I can make him wiggle his trunk.” — John von Neumann.

Turbomachinery design is, ultimately, a choice of geometry. Although we reason in terms of fluid dynamics, it is the shape of the turbine blades which dictates the flow field.

Sensible aerofoils have smooth surfaces, to avoid unnecessary accelerations and an associated increase in boundary-layer loss. In an automated, computational fluid dynamics based design process, the optimiser should have sufficient control to select the global optimum geometry, but use as few independent variables as possible. These goals are not straightforward to balance.

Fortunately, geometric parametrisation of aerofoil sections is a solved problem, thanks to the work of Brenda Kulfan ¹ at Boeing. In this post, I will give an outline of the mathematics of the method, and also walk through a Python implementation. All code is available on sourcehut ².

Prerequisite perpendicular thickness

The key to an efficient parametrisation is to reduce curvature of the function to be fitted. For example, linear interpolation requires many samples around a curved leading edge to accurately capture the geometry. Defining a curving camber line between the upper and lower surfaces of an aerofoil and fitting a perpendicular thickness with respect to this camber line is a first step to removing curvature from the problem. This is particularly effective with highly cambered turbine blades.

We begin with two sets of Cartesian coordinates \(\newcommand{\L}{\mathrm{L}}\newcommand{\U}{\mathrm{U}}x_\U,y_\U\) and \(x_\L,y_\L\) which describe the upper and lower surfaces of a datum aerofoil, which we want to run through an optimiser. The first step is to choose a camber line, and calculate the resulting perpendicular thickness.

A quadratic curve is the simplest that allows direct specification of the inlet and exit angles, \(\chi\). Define the camber line, and slope, $$ \newcommand{\C}{\mathrm{C}} \begin{align} y_\C(x) &= \tfrac{1}{2}x^2 \big(\tan \chi_2 - \tan \chi_1 \big)+ x \tan \chi_1 , \\ y_\C’(x) &= x \big(\tan \chi_2 - \tan \chi_1 \big)+ \tan \chi_1 \end{align} $$ or as code,

import numpy as np

def evaluate_camber(x, chi):
    """Camber line as a function of x, given inlet and exit angles."""
    tanchi = np.tan(np.radians(chi))
    return tanchi[0] * x + 0.5 * (tanchi[1] - tanchi[0]) * x ** 2.0


def evaluate_camber_slope(x, chi):
    """Camber line slope as a function of x, given inlet and exit angles."""
    tanchi = np.tan(np.radians(chi))
    return tanchi[0] + (tanchi[1] - tanchi[0]) * x

In the general case, the intersection of a line perpendicular to the camber line with a point on the upper surface is the solution \(x_\C\) of, $$ y_\U - y_\C(x_\C) = -\frac{ x_\U - x_\C }{y_\C’(x_\C)} \tag{I} $$ The thickness \(t_\U(x)\) is then the distance between \(x_\U,y_\U\) and \(x_\C,y_\C(x_\C)\). An identical relation applies for the lower surface, subscript \(\L\). We can use the SciPy newton solver on Eqn. (I),

from scipy.optimize import newton

def coord_to_thickness(xy, chi):
    """Perpendicular thickness distribution given camber line angles."""
    xu, yu = xy
    tanchi = np.tan(np.radians(chi))
    # Find intersections of xu, yu with camber line perpendicular
    def iterate(x):
        return yu - evaluate_camber(x, chi) + (xu - x) / evaluate_camber_slope(x, chi)
    xc = newton(iterate, xu)
    yc = evaluate_camber(xc, chi)
    # Now thickness
    t = np.sqrt(np.sum(np.stack((xu - xc, yu - yc)) ** 2.0, 0))
    return np.stack((xc, yc, t))

Let’s plot what we have so far.

Shape space transformation

The clever part of the shape space method is a transformation between real coordinates \(t(x)\) and a “shape space” \(s(x)\) where the aerofoil surface becomes a benign function, $$ t(x) = \sqrt{x}(1-x) s(x) + x \Delta , \quad \Rightarrow \quad s(x) = \frac{t(x) - x \Delta}{\sqrt{x}(1-x)} $$ where \(\Delta\) is the (blunt) trailing edge thickness. The \(\sqrt{x}\) multiplier brings the nose of the aerofoil to a rounded leading edge, reducing the curvature that must be captured in shape space; the multiplier \((1-x)\) together with the \(x \Delta\) term prescribes the trailing edge thickness. The value of \(s(0)\) sets the leading edge radius.

The Python code for these transformations is,

def from_shape_space(x, s, zte):
    """Transform shape space to real coordinates."""
    return np.sqrt(x) * (1.0 - x) * s + x * zte


def to_shape_space(x, z, zte):
    """Transform real thickness to shape space."""
    # Avoid singularities at leading and trailing edges
    eps = 1e-6
    ii = np.abs(x - 0.5) < (0.5 - eps)
    s = np.ones(x.shape) * np.nan
    s[ii] = (z[ii] - x[ii] * zte) / (np.sqrt(x[ii]) * (1.0 - x[ii]))
    return s

The reverse transformation, to shape space, has numerical singularities at \(x=0\) and \(x=1\). In practice, an optimiser will work in shape space and only use the forward transformation to produce geometry for analysis.

Fitting the elephant

Now we have wrestled the aerofoil surface into a smooth function, amenable to fitting, in the shape space. Here we use a sum of \(n+1\) Bernstein polynomials of order \(n\) to approximate \(s(x)\) (although there are many other possible ways). Let, $$ s(x) \approx \sum_{i=0}^n A_i \begin{pmatrix}n\\i\end{pmatrix} x^i(1-x)^{n-i} $$ where \(A_i\) are unknown coefficients to be determined. Bernstein polynomials have a number of attractive properties for this purpose.

1. Bernstein polynomials are infinitely differentiable, so their sum is always a smooth function with continuous curvature;
2. The maxima of each successive polynomial are equispaced along the blade chord, so each coefficient is a measure of local thickness;
3. The value of \(s(0)\) depends only on \(A_0\), so allows the leading edge radius to be directly controlled.

Functions to evaluate \(s(x)\) for known \(A_i\) go like,

from scipy.special import binom

def bernstein(x, n, i):
    """Evaluate ith Bernstein polynomial of order n at some x-coordinates."""
    return binom(n, i) * x ** i * (1.0 - x) ** (n - i)


def evaluate_coefficients(x, A):
    """Evaluate a set of Bernstein polynomial coefficients at some x-coords."""
    n = len(A) - 1
    return np.sum(
        np.stack([A[i] * bernstein(x, n, i) for i in range(0, n + 1)]), axis=0
    )

To determine the coefficients that best approximate a given \(s(x)\) we can use a least-squares fit \(B\mathbf{a} = \mathbf{s}\), where the \(n+1\) columns of \(B\) are evaluations of the Bernstein polynomials, \(\mathbf{a}\) is a column vector of \(n+1\) unknown coefficients, and \(\mathbf{s}\) is a column vector of the desired shape space distribution to approximate. To work around any singularities at the leading and trailing edges, the values of \(s(x<\delta)\) and \(s(x>1-\delta)\) for some displacement \(\delta\) may be safely ignored in the fitting process. The implementation looks like,

from scipy.linalg import lstsq

def fit_coefficients(x, s, order, dx=0.0):
    """Fit shape-space distribution with Bernstein polynomial coefficients."""
    n = order - 1
    # When converting from real coordinates to shape space, we end up with
    # singularities and numerical instability at leading and trailing edges.
    # So in these cases, ignore within dx at LE and TE
    itrim = np.abs(x - 0.5) < (0.5 - dx)
    xtrim = x[itrim]
    strim = s[itrim]
    # Evaluate all polynomials
    X = np.stack([bernstein(xtrim, n, i) for i in range(0, n + 1)]).T
    # Fit, return coefficients and residual
    return lstsq(X, strim, rcond=None)[:2]

Putting this all together, the figure shows our upper surface thickness, \(t(x)\), the transformed shape space \(s(x)\), and a Bernstein polynomial fit with \(n=3\). The transformation removes the leading edge radius and creates a singularity at \(s(1)\) which is ignored during fitting.

Wiggling his trunk

It might seem like we are home, but there is one more complication to ensure that our surface remains curvature continuous everywhere. We have treated the upper and lower halves of the aerofoil separately, so when they come together at the leading edge we must intervene to match radii of curvature. This is equivalent to matching \(s(0)\) and hence \(A_{0\U} = A_{0\L}\).

To match \(A_0\) requires fitting the upper and lower surfaces together in a bigger least-squares problem, $$ \begin{bmatrix} \uparrow & \uparrow & & \uparrow & \uparrow \\ B_{0\U} & B_{1\U} & \dots & 0 & 0 \\ \downarrow & \downarrow & &\downarrow & \downarrow \\ \uparrow & \uparrow & \uparrow & \uparrow & \\ B_{0\L} & 0 & 0& B_{1\L} & \dots \\ \downarrow & \downarrow &\downarrow &\downarrow & \\ \end{bmatrix} \begin{bmatrix} A_0 \\ A_{1\U} \\ \vdots \\ A_{1\L} \\ \vdots \end{bmatrix} = \begin{bmatrix}\uparrow \\ \mathbf{s}_\U \\ \downarrow \\ \uparrow \\ \mathbf{s}_\L \\ \downarrow \end{bmatrix} $$ Here, \(B_{i\U}\) and \(B_{i\L}\) are the \(i\)th Bernstein polynomials evaluated on the camber line \(x\)-coordinates for the upper and lower surfaces. The code to set up this least-squares problem is a bit lengthy, but we can still use scipy.linalg.lstsq as before.

Fitted in this way, and reversing the process described above to transform our approximated \(s(x)\) back into real coordinates, the aerofoil is shown in the Figure. Not bad for ten parameters (seven Bernstein coefficients plus trailing edge thickness and two camber line angles.)

Conclusion

I gave an outline of the shape space method of aerofoil geometry parametrisation. The key is to define a transform that reduces curvature of the function to be fitted. The method captures a section geometry with ten parameters, few enough to be input to an optimiser. This is a compression ratio of twelve compared to the 62 pairs of Cartesian coordinates provided as input.

With some tidying up, the Python implementation will form part of a larger project I am working on. The full code, including the script to generate the Figures, is on my sourcehut ².

I think von Neumann’s point was that number of parameters is not a useful metric, because you can “cheat” using prior knowledge of what you are trying to fit. This is exactly what we want to exploit here: because we know what aerofoils should look like, we can reduce the number of variables an optimiser needs to explore the feasible design space. Of course, we would need a different transformation if we wanted to fit elephants.

References