Geometric Parametrization

Warning

The following sections assume you have completed all the sections related to geometry preprocessing and meshing in the airfoil analysis tutorial. Be sure to complete those first before continuing.

Introduction

In order to optimize the shape of a geometry such as an airfoil or a wing, we need some way to translate design variables into actual changes in the shape. We use the Free-form deformation (FFD) technique, popularized by Tom Sederberg in the world of computer-aided graphic design, as our parametrization method. The FFD is a mapping of a region in 2D or 3D that is bounded by a set of B-splines. Every point with the region is mapped to a new location based on the deformation of the bounding B-spline curves. The B-splines themselves are defined by control points, so by adjusting the positions of these control points, we can have a great deal of control over any points embedded in the FFD volume (as shown in the figure below). Since our airfoil coordinates are point-based, we can embed them in the FFD and give the optimizer control over their shape. For more detail on the mathematical principles and implementation, refer to the following article.

The actual implementation of the FFD method is housed in the pygeo repository. The specific class to look for is DVGeometry. Navigate to the DVGeometry repository to explore its capabilities Before diving into the parametrization, however, we need to generate an FFD, which is basically a 3D grid in the PLOT3D format. This file contains the coordinates of the FFD points around the airfoil. These are control points that are fitted to the airfoil using B-splines, which are used to deform the airfoil.

Files

The coordinates for the NACA0012 airfoil that we processed in the last tutorial are in the file n0012_processed.dat. Navigate to the directory airfoilopt/ffd in your tutorial folder. Create the following empty runscript in the current directory.

• run_ffd.py

Import Packages

import numpy as np

Load Airfoil

airfoil = np.loadtxt("../../airfoil/geometry/n0012_processed.dat")
npts = airfoil.shape[0]
offsetTEindex = -8
nmid = (npts + 1) // 2 + offsetTEindex

The following two functions are used to get the upper and lower points of the airfoil. It should be noted that in airfoil geometry processing step in the previous tutorial we created a blunt trailing for the NACA0012. As result, we can’t just simply divide the airfoil’s point set down the middle to get the upper and lower surfaces. Index 0 starts on the upper surface trailing edge however we also want the upper half of the blunt trailing edge to be included in our upper surface. Since our blunt trailing edge is 17 points long we will split at an offset of index -8.

def getupper(xtemp):
    myairfoil = np.ones(npts)
    for i in range(offsetTEindex, nmid):
        myairfoil[i] = abs(airfoil[i, 0] - xtemp)
    myi = np.argmin(myairfoil)
    return airfoil[myi, 1]


def getlower(xtemp):
    myairfoil = np.ones(npts)
    for i in range(nmid, npts + offsetTEindex):
        myairfoil[i] = abs(airfoil[i, 0] - xtemp)
    myi = np.argmin(myairfoil)
    return airfoil[myi, 1]

Clarification on why we split at an offset for the blunt TE airfoil

FFD Box Creation

The FFD box can now be set up.

nffd = 10

FFDbox = np.zeros((nffd, 2, 2, 3))

xslice = np.zeros(nffd)
yupper = np.zeros(nffd)
ylower = np.zeros(nffd)

xmargin = 0.001
ymargin1 = 0.02
ymargin2 = 0.005

for i in range(nffd):
    xtemp = i * 1.0 / (nffd - 1.0)
    xslice[i] = -1.0 * xmargin + (1 + 2.0 * xmargin) * xtemp
    ymargin = ymargin1 + (ymargin2 - ymargin1) * xslice[i]
    yupper[i] = getupper(xslice[i]) + ymargin
    ylower[i] = getlower(xslice[i]) - ymargin

nffd signifies the number of chordwise slices. We pre-allocate an array of generic size (a,b,c,3) to set up an empty FFD box. In this example, a=nffd (number of chordwise sections), b=c=2 (number of spanwise and thickness-wise sections respectively) and the final 3 is “fixed” as we are using 3D coordinates for each point. An empty FFD box is created. xmargin and ymargin specify the closest distance from the airfoil to place the FFD box. xslice, yupper, and ylower store the x- and y- coordinates of the control points for each slice along the chord, taking into account the margins from the airfoil. xslice, yupper, and ylower store the x- and y- coordinates of the control points for each slice along the chord, taking into account the margins from the airfoil.

# X
FFDbox[:, 0, 0, 0] = xslice[:].copy()
FFDbox[:, 1, 0, 0] = xslice[:].copy()
# Y
# lower
FFDbox[:, 0, 0, 1] = ylower[:].copy()
# upper
FFDbox[:, 1, 0, 1] = yupper[:].copy()
# copy
FFDbox[:, :, 1, :] = FFDbox[:, :, 0, :].copy()
# Z
FFDbox[:, :, 0, 2] = 0.0
# Z
FFDbox[:, :, 1, 2] = 1.0

The x- and y- coordinates are transferred to the FFDbox variable. Since the airfoil slices are the same along the z-direction, the x- and y- coordinates are copied over. Since the airfoil slices are the same along the z-direction, the x- and y- coordinates are copied over. The z-coordinates are updated to 0 and 1.

Writing to File

with open("ffd.xyz", "w") as f:
    f.write("1\n")
    f.write(str(nffd) + " 2 2\n")
    for ell in range(3):
        for k in range(2):
            for j in range(2):
                f.write(" ".join(f"{x:.15f}" for x in FFDbox[:, j, k, ell]))
                f.write("\n")

We now have to write our FFD box to PLOT3D format for pygeo. The nested loop above is one typically seen for writing PLOT3D files manually.

Run it yourself!

You can now run the python file with the command:

python run_ffd.py

The above script writes the FFD coordinates to a PLOT3D .xyz file, which will be used for optimization.

3D view of the FFD volume

2D view of the FFD volume, together with the embedded airfoil