02 Sep

DHMTU Aerofoil Definition

Without doubt, the nation that has progressed ground-effect-vehicles the furthest is Russia. More specifically, this was the USSR from the 1960s through to the late 1980s with their Ekranoplan aircraft. These GEVs (Ground Effect Vehicles) were designed by the Central Hydrofoil Design Bureau, at the heart of the aerodynamic design was the Department of Hydromechanics of the Marine Technical University (DHMTU) in Saint Petersburg.

Undoubtedly the most famous Ekranoplan, the KM, or in the west, the Caspian Sea Monster.

The DHMTU established an aerofoil series whose geometries are defined analytically, which to my knowledge is the only series of parametrised aerofoils designed specifically for ground effect flight. The definition varies quite significantly from the NACA 4-digit aerofoils, already in GEVfoil. The aerofoil has a partly or fully flat base with an s-shaped mean camber line.

DHMTU Notation

The DHMTU aerofoils notation has eight digits, given as follows:

DHMTU Y1-X1-Y2-X2-Y3-X3up-RLE

Y1maximum ordinate of the upper surface (%c)
X1abscissa of Y1 (%c)
Y2ordinate of the start of the flat section (%c, positive value represents negative y)
X2abscissa of Y2 (%c)
Y3ordinate of the end of the flat section (%c, positive value represents negative y)
X3absissa of Y3 (%c)
δupslope parameter of the upper trailing edge
RLEleading edge radius parameter

Some DHMTU sections from the literature include DHMTU 12-35-3-10-2-80-12-2 [1] and DHMTU 10-40-2-10-2-60-21-5 [2]. These names don’t exactly register in memory as readily as the likes of the NACA 0012! The JavaFoil manual has a very nice figure of the DHMTU section:

DHMTU aerofoil section, illustration from the manual of JavaFoil. Source: JavaFoil [3]

The above information is readily available in some western papers such as Moore’s [1], theses such as a Smuts [2] and manuals such as JavaFoil by Hepperle [3]. However, no further information is given. Moore, from 2002, references a South Korean website, which is now defunct. Hepperle references a Korean document from 1996; ‘Development of S-shaped Section (DHMTU Family)‘ by Chun Ho-Hwan, Chang Chong-Hee of Pusan University. No further information is given about this document, for example whether it’s a conference paper, a journal paper or perhaps a thesis, and despite my best efforts I was initially unable to find this document online.

Lacking further information, the notation is only sufficient to define these basic points and shapes. Clearly more information is required than the notation alone.

With thanks to my friends in the Wing In Ground Effect group on groups.io I have managed to obtain the ‘Development of S-shaped Section’ paper in Korean [4] from a member’s personal archive (thank you, Marc!). I cannot recommend this group enough for those with an interest in ground effect craft, it is the successor to the long standing group of the same name on Yahoo Groups, whose service is moribund.

I have absolutely no command of the Korean language, however the equations and diagrams are mutually intelligible, and the online translators are adequate enough. I have done my best to capture this information below. At a future date I will explore options for more formally capturing my findings following some further analysis, in the hope that this fascinating aerofoil family can be recorded for prosperity in western literature.

The following is my best attempt (along with online translators) to understand this document.

Defining the Upper-Surface

The upper-surface is split into a fore and aft section, the fore starting at the leading edge and terminating at the maximum y ordinate, Y1. The aft section somewhat intuitively starts at Y1 and ends at the trailing edge.

Upper-surface – Fore Section

\[ y_{upper} = a_0 \sqrt{x} + a_1 x + a_2 x^2 +a_3 x^3 \]

subject to \[0 < x < x_1\]

And secondly the upper-surface aft of the maximum ordinate:

\[ y_{upper} = d_0 + d_1(1-x) + d_2(1-x)^2 + d_3(1-x)^3 \]

subject to \[ x_1 < x < 1\]

The parameters a0 to a4 are obtained via a fairly convoluted process.

Leading edge radius \[r_t = a_0/2\]

Leading edge parameter \[ r_{nose} = r / t+u^2\]

Radius at the point X1 is \[R = \frac{1}{(2d_2+6d_3(1-x_1)}\]

The paper then works through some derivation – and it is unclear throughout the paper whether this is the interpretation of the author, or whether this is the author translating from an original and unreferenced Russian document. Four conditions are given as a set of linear equations:

Condition 1 (for a0):

\[t_u = a_0 \sqrt{x_t} + a_1 x_t+a_2 x_t^2+a_3 x_t^3\]
\[t_u – a_0 \sqrt{x_t} = a_1 x_t+a_2 x_t^2+a_3 x_t^3\]

Condition 2 (for a1):

\[y’_{upper}(x=x_t) = \frac{a_0}{2\sqrt{x_t}} + a_1 + 2 a_2 x_t + 3a_3 x_t^3= 0\] \[-\frac{a_0}{2\sqrt{x_t}} = a_1 + 2a_2x_t+3a_3x_t^3\]

Condition 3(for a2):

\[r=kt_u^2\] \[k t_u^2 = a_0^2 / 2\] \[\therefore a_0 = \sqrt{2k}t_u\]

Condition 4 (for a3):

\[y_{upper} = \frac{1}{R} \] \[2a_2 + 6a_3x_t = 2d_2 + 6d_3(1-x_t) + a_0/4x_t^{\frac{3}{2}}\]

Upper-surface – Aft Section

For the aerofoil aft of the maximum ordinate, the four parameters d0 to d3 are obtained via similarly convoluted process.

Condition 1:

\[t_u = d_0 + d_1(1-x_t) + d_2(1-x_t)^2 + d_3(1-x_t)^3\] \[t_u = -\delta_b(1-x_t) + d_2(1-x_t)^2 + d_3(1-x_t)^3\]

Condition 2:

\[y_{up} ^ {\prime} = – d_1 + 2d_2(1-x_t) – 3d_3(1-x_t)^2 = \delta_b – 2d_2(1-x_t) -3d_3(1-x_t)^2\] \[2d_2(1-x_t)=\delta_b -3d_3(1-x_t)^2\] \[d_2 = \frac{\delta_b – 3 d_3( 1-x_t)^2}{2(1 – x_t)}\] \[t_u=-\delta_b(1-x_t)+d_2(1-x_t)^2+d_3(1-x_t)^3\] \[=-\delta_b(1-x_t)+\frac{\delta_b-3d_3(1-x_t)^2}{2(1-x_t}(1-x_t)^2+d_3(1-x_t)^3\] \[=-\delta_b(1-x_t)+\frac{1}{2}\delta_b(1-x_t)-\frac{3}{2}d_3(1-x_t)^3+d_3(1-x_t)^3\] \[=\frac{1}{2}\delta_b(1-x_t)-\frac{1}{2}\delta_3(1-x_t)^3\] \[\therefore \frac{1}{2}d_3(1-x_t)^3=–\frac{1}{2}\delta_b(1-x_t)-t_u\] \[d_3 = \frac{-\delta_b(1-x_t)-t_u}{(1-x_t)^2}=\frac{-\delta_b-\frac{t_u}{(1-x_t)}}{(1-x_t)^2}\] \[d_2 = \frac{t_u+\delta_b(1-x_t)-d_3(1-11x_t)^3}{(1-x_t)^2}\]

Condition 3:

Mercifully, the third condition simple:

\[ d_0 = 0 \]

Condition 4:

\[y_{upper}\prime=-d_1-2d_2(1-x)-3d_3(1-x)^2\] \[\delta_b = -d_1\hspace{0.5cm}at\hspace{0.5cm} x = 1\] \[\therefore d_1 = – \delta_b \]

Defining the lower-surface

The lower surface is split into three sections, the curved fore section, similar to the upper-section, a flat mid section and a curved aft section.

Lower-surface – Fore Section

As with the upper-surface:

\[ y_{lower} = b_0 \sqrt{x} + b_1 x + b_2 x^2 +b_3 x^3 \]

subject to: \[0 < x < x_1\]

The four parameters a0 to a4 are obtained via a fairly convoluted process.

tlower is the ordinate of the lower surface at point x1.

y” = 0 at point x1.

Leading edge radius:

\[r_t=\frac{b_0^2}{2}\]

At the point x1:

\[\frac{dy}{dx}fore =\frac{dy}{dx}mid\]

That is to say, there is continuity at the interface of the sections.

Condition 1

\[b_1x_1+b_2x_1^2+b_3x_1^3=t_{lower}-b_0\sqrt{x_1}\]

Condition 2

\[2b_2+6b_3x_1=\frac{b_0}{4x_1^{\frac{3}{2}}}\]

Condition 3

\[k=\frac{r}{t_u^2}\] \[r=t_u^2k\] \[b_0\sqrt{2kt_u}=a_0\]

Condition 4

\[b_1+2b_2x_1+3b_3x_1^2=c_1-b_0 / 2\sqrt{x_1}\]

Lower Surface – Mid-Section

The flat middle section is somewhat simply defined, between x1<x<x2:

\[y_{lower} = -c_0 – c_1 (x-x_1)\] \[c_0 =t_{lower}\] \[c_1 = \frac{t_{lower2}-t_{lower}}{x_2-x_1} \]

Lower Surface – Aft-Section

Finally, the aft-section is defined in a manner not entirely dissimilar to the fore, between x2 < x < (c=1):

\[y_{lower} = -(e_0+e_1(1-x_t)+e_2(1-x_t)^2+e_3(1-x_t)^3)\]

tlower2 is the ordinate of the lower surface at point x2.

y” (x2) = 0 at point x1.

At x=1 y(1)=0 – that is to say the trailing edge is (1, 0), i.e. an abscissa of 1 with an ordinate of zero. Noting that we defining this aerofoil with a normalised chord at unity, as is customary with aerofoils.

At the point x2:

\[\frac{dy}{dx}mid =\frac{dy}{dx}aft\]

Once again I believe the paper is just telling us there is continuity and no gradient change at the point where the mid and aft sections meet.

The pattern continues once one with linear expressions and coefficient definition:

Condition 1

\[e_0+e_1(1-x_2)+e_2(1-x_2)^2+e_3(1-x_2)^3 = t_{lower2}\]

Condition 2

\[2e_2+6e_3(1-x_2)=0\]

Condition 3

\[e_0 = d_0 = 0\]

Condition 4

\[e_1-2e_2(1-x_2)-3e_3(1-x_2)^2 = – c_1\]

Thoughts and Conclusions

The mathematics behind the DHMTU sections are certainly more complex than other analytically defined geometries such as the well known and used NACA 4-series and 6-series.

It cannot be stressed strongly enough – I have made my best attempt to understand a paper in a foreign language (Korean), I have no understanding as to the validity of said reference. Consequently, I cannot have much confidence in these equations as they are presented.

I do have access to some images of DHMTU sections. My next step will be to implement this series into GEVfoil and attempt to reproduce the shapes, to increase confidence.


[1] Moore, N., Wilson, P.A. and Peters, A.J. (2002) An investigation into wing in ground effect aerofoil geometry. In RTO-MP-095. NATO RTO.
[2] Smuts, E. (2009) A Computational Study of a Lifting Wing in Close Proximity to a Moving Ground Plane. Masters Dissertation.
[3] Hepperle, M (2017). JavaFoil Theory Document. User Manual.
[4] Ho-Hwan, C., Chong-Hee C. (1996). Development of S-shaped Section (DHMTU Family). Unknown publication type.

06 Mar

Introducing GEVfoil

I have decided to rename my little aerofoil project. I found no fewer than four other projects using the previous name, which I’ll avoid stating to reduce confusing search engines any further.

I’ve settled on GEVfoil (Ground Effect Vehicle Aerofoil) which may seem a little obvious, but clarity should at least be in my favour. The objectives remain the same;

  • To create a set of tools to aid in the analysis of aerofoils in ground effect.
  • To primarily focus on 2D effects, but not exclude 3D effects where appropriate and relevant.
  • To utilise open-source tools as much as possible, for example links into the OpenFOAM CFD (Computational Fluid Dynamics) project.
  • To maintain a blog of the development, and open-source the whole project when it has content enough to be worthy.

I have renamed previous blog posts with this in mind:

18 Feb

Ramer–Douglas–Peucker algorithm

Generally speaking, for low-speed flow around an aerofoil, the largest gradients in the flow-properties coincide with the greatest curvature of the aerofoil, and thus are the region that requires the greatest concentration of cells (or nodes) in the computational domain.

Many of the simplified treatments for spacing points along an aerofoil use a logarithmic spacing biased towards the leading edge. This of course can create issues if spaced in the chord-wise direction (x/c). The surface-normal at the leading edge of a symmetric aerofoil is perpendicular to the stream-wise direction, and thus can become under-resolved, and even faceted, with a simple linear or logarithmic chord-wise spacing.

One of the more significant issues with linear or logarithmic spacing is that it simply does not discriminate on the basis of curvature, thus the growth factor becomes a qualitative parameter controlled by the aerodynamicist. Any other curvature, for example the upper-surface curvature caused by the s-shaped chord of a DHTMU aerofoil, is neglected.

S-shaped chord of the DHMTU aerofoils. Source: JavaFoil Manual

Some computational treatments, though able to nicely capture curvature, are a little too complex for this simplified tool intended for 2D aerofoil analyses. I believe that the Ramer-Douglas-Peucker (RDP) algorithm is a simple yet powerful tool ideally suited to be applied to spacing points/nodes along a computational grid for for 2D aerofoil analyses.

RDP is a a decimation (down-sampling) algorithm, which although cannot increase resolution (create points) in areas of high-curvature, it can remove remove points in areas where curvature is low. This approach is quite simple and elegant. The algorithm works by creating a line between the start and end points, then recursively divides this line whilst each time finding the furthest point, marking these points to be kept if they are greater than a user-specified distance (epsilon).

Ramer-Douglas-Peucker Animation. Source: Mysid

As this requires a set of points, rather than a spline or other analytical curve, there are upsides and downsides. In the case of of analytic aerofoils, such as the DHTMU, or the NACA 4 that I originally scripted into GEVgoil, these need to be discretised before being down-sampled. However, for aerofoils that are already being supplied as a set of points, for example using the CSV method that is now also in GEVfoil, RDP can be applied directly.

To simplify matters further, RDP is available as a Python module via PIP: https://pypi.org/project/rdp/ Using this excellent ready-made library, I have constructed perhaps the simplest Python class, implementing this RDP module into GEVfoil:

from rdp import rdp
import numpy as np

class RDP(object):
    """A class to decimate aerofoil points"""
    def __init__(self, x, y, e=0.0001):
        xy = np.vstack((x, y)).T
        self.reducedXY = rdp(xy, epsilon=e)
        self.reducedX = self.reducedXY[:,0]
        self.reducedY = self.reducedXY[:,1]

And of course, the proof is in the pudding. Some examples of the RDP module in action in GEVfoil with varying epsilon:

Original NACA 2412 as generated by GEVfoil with 800 points each on the upper and lower surfaces
Decimated with RDP using an epsilon of 1e-05
Decimated with RDP using an epsilon of 1e-04
Decimated with RDP using an epsilon of 1e-03
21 Oct

GEVfoil – NACA 4 Series Part 2 and CSV Reading

Previously I had made a simple function that will generate a number of points/coordinates for any specified NACA 4 series aerofoil. The limitation with this approach was that it did a fairly poor job at capturing the leading-edge. The leading edge is of critical importance to the aerodynamic performance and characteristics of the aerofoil. After some consideration, I have decided not to modify this technique. I am designing the GEVfoil tools to be as modular, generic and extendable as possible. This isn’t to say I am neglecting the importance of the leading edge, or that the method is flawed.

Initially I had considered returning the aerofoil as a mathematical function, rather than a set of points. This concept may well be revisited in the future. However, most of the other input aerofoils will arrive as a set of coordinates of varying fidelity, from sources such as a the UIUC Airfoil Coordinates Database or from Comma Separated Variable (CSV) files. Therefore, the classes which will receive the aerofoil coordinate objects will deal with interpolating and appropriately converting coordinates into functions anyway. Consequently, as long as the linear-spacing is arbitrarily high enough a simple piece-wise cubic interpolation will suffice.

Note that the CSV reader itself is an abstraction of the AerofoilReader class, which I plan to later extend to read the DAT file from UIUC,

To that end, the CSV reader currently looks as follows:

import os

class AerofoilReader(object):
    """ An abstract class to be inherited by classes for specific
        aerofoil file types """
    def __init__(self, filename):
        super(AerofoilReader, self).__init__()
        self.coord = []
        self.title = []  # Column titles
        self.label = []  # Aerofoil Label

        foil = open(filename, 'r')
        self.lines = foil.readlines()


class CSVReader(AerofoilReader):
    """ Read in CSV. Assumes a header is present."""
    def __init__(self, filepath, xcol, ycol, zcol=False):
        super(CSVReader, self).__init__(filepath)

        self.label = os.path.basename(filepath).split('.')[0]

        inpcols = 3        # Two or three columns
        if not zcol:
            inpcols = 2

        for i, line in enumerate(self.lines):
            parts = line.rstrip().replace('"', '').split(',')
            additional = []

            if i > 0:
                # Not the header
                if len(parts) > inpcols:
                    # More than x, y &amp; z - store other columns
                    for j in range(len(parts)):
                        if j != xcol and j != ycol and j != zcol:
                            additional.append(parts[j])

                xyz = [parts[xcol], parts[ycol], parts[zcol]]
                self.coord.append(xyz + additional)
            else:
                # It's the header
                self.title = parts
13 Feb

GEVfoil – NACA 4 Series

The intention behind GEVfoil is to establish an open-source set of Python classes to handle aerofoils (American English: airfoils). The intention is that these can be used for anything from generating simple point clouds, analytical representations all the way to STL files for CFD meshing and OpenFOAM BlockMesh output.

The first step is starting with something basic, in this case a simple class to generate NACA 4 series aerofoils:

import numpy as np
import matplotlib.pyplot as plt


class NACA4(object):
    """A class for generating  NACA 4 series aerofoils
        TO-DO - linspace not good at LE capture """
    def __init__(self, foil, points=200, chord=1):
        super(NACA4, self).__init__()
        m = float(foil[0])/100      # max camber
        p = float(foil[1])/10       # chordwise position of max camber
        t = float(foil[2:])/100     # thickness
        x = np.linspace(0, chord, points)
        self.x = x
        yt = self.thickness(x, t)
        yc = []
        for coord in x:
            if m:
                yc.append(self.camber(m, p, coord))
            else:
                # No camber
                yc.append(0)
        y1 = yc + yt
        y2 = yc - yt

        self.y = y1

        plt.plot(self.x, yc)
        plt.plot(self.x, y1, '.')
        plt.plot(self.x, y2, '.')
        plt.axis('equal')
        plt.show()

    def thickness(self, x, t):
        # Y thickness at given x point
        return t / 0.2 * \
            (0.2969*np.sqrt(x)-0.126*x-0.3516*x**2+0.2843*x**3-0.1015*x**4)

    def camber(self, m, p, x):
        # Return the camber of the aerofoil
        if x <= p:
            return m/p**2 * (2*p*x-x**2)
        return m/(1-p)**2*((1-2*p)+2*p*x-x**2)

if __name__ == "__main__":
    test = NACA4('2412')

For now, Matplotlib has been built into the class for visual reference. For example, 2412 generates:

And 0010 creates:

At this stage there is a clear issue around the leading edge, the function was interrogated with a linearly spaced x/chord-wise interval, this must be improved for the future to cluster points around the greatest curvature.