Airfoil optimization for Vertical Axis Hydrokinetic Turbine using genetic algorithm(2)---Modified PARSEC Parametrization

1 Modified PARSEC Parametrization

As mentioned before, the PARSEC parametrization popularized by Sobieczky actually defines a set of shape functions, and the upper and lower curve of airfoil can be determined by a combination of these functions.

$$
Z_{k}=\sum_{n=1}^{6} a_{n,k}X_{k}^{n-\frac{1}{2}}
\tag{1-1}$$

The original PARSEC parametrization has these 11 parameters, each affecting a physical characteristic of airfoil. However, there is a large range of airfoil profiles that set the trailing edge thickness to zero. Then the physical meaning of trailing edge direction becomes abstract. To avoid this circumstances, several adaptations have been proposed to the original PARSEC, which directly use the trailing edge angle of upper and lower curve( $\theta_{TE,up}$ , $\theta_{TE,low}$) instead of $\alpha_{TE}$ and $\beta_{TE}$. Besides, the x coordinate of trailing edge is set to be 1 for convenience, and the difference of leading edge radius between upper and lower curve is also considered.
Then the 11 parameters become:

• Leading edge radius $R_{LE,UP} , R_{LE,LOW}$
• upper crest location $X_{UP},Z_{UP}$
• lower crest location $X_{LO},Z_{LO}$
• upper and lower curvature $Z_{XXUP},Z_{XXLO}$
• trailing edge coordinate $Z_{TE}$
• upper trailing edge angle $\theta_{TE,up}$
• lower trailing edge angle $\theta_{TE,low}$

The coefficients of shape functions can be found by imposing some constraints, take the upper curve as an example:

$$
\begin{aligned}
Z_{up}=\sum_{n=1}^{6} a_{n}x_{k}^{n-\frac{1}{2}}=a_{1} x^{0.5}+a_{2} x^{1.5}+a_{3} x^{2.5}+a_{4} x^{3.5}+a_{5} x^{4.5}+a_{6} x^{5.5}
\end{aligned}
\tag{1-2}$$

The first and second order derivative can be expressed as:

$$
\begin{aligned}
\frac{dZ_{u p}}{dx}=0.5 a_{1} x^{-0.5}+1.5 a_{2} x^{0.5}+2.5 a_{3} x^{1.5}+3.5 a_{4} x^{2.5}+4.5 a_{5} x^{3.5}+5.5 a_{6} x^{4.5}
\end{aligned}
\tag{1-3}$$

$$
\begin{aligned}
\frac{d^{2}Z_{u p}}{dx^{2}}=&-0.25 a_{1} x^{-1.5}+0.75 a_{2} x^{-0.5}+3.75 a_{3} x^{0.5}+8.75 a_{4} x^{1.5}\&+15.75 a_{5} x^{2.5}+24.75 a_{6} x^{3.5} \end{aligned}
\tag{1-4}$$

1.1 Leading edge radius can be calculated by the following equation:

$$
\begin{aligned}
R_{LE}&=\left| \frac{\left( 1+\left( \frac{dZ}{dx}\right)^{2} \right) ^{\frac{3}{2}}}{\frac{d^{2}Z}{dx^{2}}}\right| {x=x{te}} =\left| \frac{\left[1+\left(\sum_{n=1}^{6}\left(n-\frac{1}{2}\right) a_{n} \cdot x^{n-\frac{3}{2}}\right)^{2}\right]^{\frac{3}{2}}}{\sum_{n=1}^{6}\left(n-\frac{1}{2}\right)\left(n-\frac{3}{2}\right) a_{n} \cdot x^{n-\frac{5}{2}}} \right|\ &=\left|\frac{\left[1+\left(\frac{a_{1}}{2 x^{1 / 2}}+\sum_{n=2}^{6}\left(n-\frac{1}{2}\right) a_{n} \cdot x^{n-\frac{3}{2}}\right)^{2}\right]^{3 / 2} }{-\frac{a}{4 x^{3 / 2}}+\sum_{n=2}^{6}\left(n-\frac{1}{2}\right)\left(n-\frac{3}{2}\right) a_{n} \cdot x^{n-\frac{5}{2}}}\right|\&=\left|\frac{\left(\frac{a_{1}}{2 x^{1 / 2}}\right)^{3}\left[\left(\frac{2 x^{1 / 2}}{a_{1}}\right)^{2}+\left(1+\frac{2 x^{1 / 2}}{a_{1}} \sum_{n=2}^{6}\left(n-\frac{1}{2}\right) a_{n} \cdot x^{n-\frac{3}{2}}\right)^{2}\right]^{3 / 2}}{\left(-\frac{a_{1}}{4 x^{3 / 2}}\right)\left(1+\left(-\frac{4 x^{3 / 2}}{a_{1}}\right) \sum_{n=2}^{6}\left(n-\frac{1}{2}\right)\left(n-\frac{3}{2}\right) a_{n} \cdot x^{n-\frac{5}{2}}\right)}\right|
\end{aligned}
\tag{1-5}$$
When x tends to be zero, then:

$$
R_{LE} \stackrel{x \rightarrow 0}{\longrightarrow}\left|\frac{a_{1}^{2}}{2}\right|
\tag{1-6}$$
Thus the value $a_1$ for the upper curve $a_{1,up}=\sqrt{2R_{LE,up}}$ , and for the lower curve $a_{1,low}=-\sqrt{2R_{LE,low}}$ .

1.2 At the upper crest position:

$$
z\left( x_{up}\right)=z_{up}=\sum_{n=1}^{6} a_{n}x_{up}^{n-\frac{1}{2}}
\tag{1-7}$$

1.3 The first order derivative at the upper crest position is supposed to be zero according to PARSEC parametrization:

$$
\left[\frac{dZ_{u p}}{dx}\right]{x=x{up}}=\sum_{n=1}^{6}\left(n-\frac{1}{2}\right) a_{n} \cdot x_{up}^{n-\frac{3}{2}}=0
\tag{1-8}$$

1.4 The second order derivative at upper crest position is defined as a control parameter $Z_{XXUP}$ :

$$
\left[\frac{d^{2}Z_{u p}}{dx^{2}}\right]{x=x{up}}=Z_{XXUP}=\sum_{n=2}^{6}\left(n-\frac{1}{2}\right)\left(n-\frac{3}{2}\right) a_{n} \cdot x_{up}^{n-\frac{5}{2}}
\tag{1-9}$$

1.5 The trailing edge angle can be expressed as:

$$
tan\left(\theta_{te,up}\right)=\left[\frac{dy}{dx}\right]{x=x{te}}=\sum_{n=1}^{6}\left(n-\frac{1}{2}\right) a_{n} \cdot x_{te}^{n-\frac{3}{2}}
\tag{1-10}$$

1.6 The trailing edge coordinate $\left(X_{te}, Z_{te}\right)$ should also satisfy the function:

$$
Z\left(X_{te}\right)=Z_{te}=\sum_{n=1}^{6} a_{n,k}X_{te}^{n-\frac{1}{2}}
\tag{1-11}$$

The coefficients $a_{n}$ can therefore be found by solving (1-6)—(1-11). The matrix form of these equations:

$$
\begin{aligned}
\begin{bmatrix}
1& 0 &0 &0 &0 &0 \
x_{te}^{\frac{1}{2}} & x_{te}^{\frac{3}{2}}& x_{te}^{\frac{5}{2}}& x_{te}^{\frac{7}{2}}& x_{te}^{\frac{9}{2}}& x_{te}^{\frac{11}{2}}\
x_{up}^{\frac{1}{2}} & x_{up}^{\frac{3}{2}}& x_{up}^{\frac{5}{2}}& x_{up}^{\frac{7}{2}}& x_{up}^{\frac{9}{2}}& x_{up}^{\frac{11}{2}}\
\frac{1}{2}x_{te}^{-\frac{1}{2}}& \frac{3}{2}x_{te}^{\frac{1}{2}} & \frac{5}{2}x_{te}^{\frac{3}{2}} & \frac{7}{2}x_{te}^{\frac{5}{2}} & \frac{9}{2}x_{te}^{\frac{7}{2}} & \frac{11}{2}x_{te}^{\frac{9}{2}}\
\frac{1}{2}x_{up}^{-\frac{1}{2}}& \frac{3}{2}x_{up}^{\frac{1}{2}} & \frac{5}{2}x_{up}^{\frac{3}{2}} & \frac{7}{2}x_{up}^{\frac{5}{2}} & \frac{9}{2}x_{up}^{\frac{7}{2}} & \frac{11}{2}x_{up}^{\frac{9}{2}}\
-\frac{1}{4}x_{up}^{-\frac{3}{2}} & \frac{3}{4}x_{up}^{-\frac{1}{2}} & \frac{15}{4}x_{up}^{\frac{1}{2}} & \frac{35}{4}x_{up}^{\frac{3}{2}} & \frac{63}{4}x_{up}^{\frac{5}{2}} & \frac{99}{4}x_{up}^{\frac{7}{2}}
\end{bmatrix}\begin{bmatrix}
a_1\
a_2\
a_3\
a_4\
a_5\
a_6\
\end{bmatrix}=\begin{bmatrix}
\sqrt{2R_{LE,up}} \
Z_{TE}\
Z_{up}\
tan\left(\theta_{te,up}\right)\
0\
Z_{XXUP}\
\end{bmatrix}
\end{aligned}
\tag{1-12}$$
Similar equations apply for the lower curve.

2 Python Code

It’s not difficult to find some scripts about PARSEC parametrization, and actually This person‘s code helped me a lot(I just need to make some small changes).

from __future__ import division
from math import *
import numpy as np
import matplotlib.pyplot as plt


# ========================================================================================
def pcoef(
        xte, yte, rle,
        x_cre, y_cre, d2ydx2_cre, th_cre,
        surface):
    """evaluate the PARSEC coefficients"""
    # Initialize coefficients
    coef = np.zeros(6)

    # 1st coefficient depends on surface (pressure or suction)
    if surface.startswith('p'):
        coef[0] = -sqrt(2 * rle)
    else:
        coef[0] = sqrt(2 * rle)

    # Form system of equations
    A = np.array([
        [xte ** 1.5, xte ** 2.5, xte ** 3.5, xte ** 4.5, xte ** 5.5],
        [x_cre ** 1.5, x_cre ** 2.5, x_cre ** 3.5, x_cre ** 4.5,
         x_cre ** 5.5],
        [1.5 * sqrt(xte), 2.5 * xte ** 1.5, 3.5 * xte ** 2.5,
         4.5 * xte ** 3.5, 5.5 * xte ** 4.5],
        [1.5 * sqrt(x_cre), 2.5 * x_cre ** 1.5, 3.5 * x_cre ** 2.5,
         4.5 * x_cre ** 3.5, 5.5 * x_cre ** 4.5],
        [0.75 * (1 / sqrt(x_cre)), 3.75 * sqrt(x_cre), 8.75 * x_cre ** 1.5,
         15.75 * x_cre ** 2.5, 24.75 * x_cre ** 3.5]
    ])

    B = np.array([
        [yte - coef[0] * sqrt(xte)],
        [y_cre - coef[0] * sqrt(x_cre)],
        [tan(th_cre * pi / 180) - 0.5 * coef[0] * (1 / sqrt(xte))],
        [-0.5 * coef[0] * (1 / sqrt(x_cre))],
        [d2ydx2_cre + 0.25 * coef[0] * x_cre ** (-1.5)]
    ])

    # Solve system of linear equations
    X = np.linalg.solve(A, B)

    # Gather all coefficients
    coef[1:6] = X[0:5, 0]

    # Return coefficients
    return coef


# ========================================================================================
# Sample coefficients  NACA0012
'''
leading edge radius (rle_suc rle_pre)
pressure and suction surface crest locations (x_pre, y_pre, x_suc, y_suc)
curvatures at the pressure and suction surface crest locations (d2y/dx2_pre, d2y/dx2_suc)
trailing edge coordinates (x_TE, y_TE)
trailing edge angles between the pressure and suction surface and the horizontal axis (th_pre, th_suc)
'''
rle_suc= .014927
rle_pre= .014181
x_suc = .29866
y_suc = .059404
x_pre = .29962
y_pre = -0.059632
d2ydx2_suc = -.42399
d2ydx2_pre = 0.445281
yte = 0
th_suc = -7.672047
th_pre = 7.59506

xte = 1.0 # Trailing edge x position
# Evaluate pressure (lower) surface coefficients
cf_pre = pcoef(xte, yte, rle_pre,
               x_pre, y_pre, d2ydx2_pre, th_pre,
               'pre')
# Evaluate suction (upper) surface coefficients
cf_suc = pcoef(xte, yte, rle_suc,
               x_suc, y_suc, d2ydx2_suc, th_suc,
               'suc')
x = (1 - np.cos(np.linspace(0, 1, np.ceil(1e3)) * np.pi)) / 2
uppery = np.array([0] * len(x))
lowery = np.array([0] * len(x))
for i in range(1, 7):
    uppery = uppery + cf_suc[i - 1] * x ** (i - 1 / 2)
    lowery = lowery + cf_pre[i - 1] * x ** (i - 1 / 2)
# ========================================================================================
# Plot this airfoil
def Naca0012(x):
    y=0.594689181*(0.298222773*sqrt(x) - 0.127125232*x - 0.357907906*x**2 + 0.291984971*x**3 - 0.105174606*x**4)
    return y
f2 = np.vectorize(Naca0012)
Nacay=f2(x)
Error_up=Nacay-uppery
Error_lo=-Nacay-lowery
RMS_up=0
RMS_lo=0
for item in range(1,len(x)):
    RMS_up=RMS_up+(uppery[item-1]-Nacay[item-1])**2
    RMS_lo=RMS_lo+(lowery[item-1]+Nacay[item-1])**2
RMS_up=sqrt(RMS_up/len(x))
RMS_lo=sqrt(RMS_lo/len(x))
RMS=sqrt(0.5*RMS_up**2+0.5*RMS_lo**2)
#plt.plot(x, uppery, x, lowery)
#plt.gca().axis('equal')
lines=plt.plot(x, uppery,'o--',x,lowery, 'o--')
plt.setp(lines, color='b', linewidth=2.0)
lines2=plt.plot(x,Nacay,x,-Nacay)
plt.setp(lines2, color='r', linewidth=2.0)
plt.show()
with open('airfoil_data.txt', 'w') as f:
    for i in range(1,len(x)):
        f.write("%s\t\t%s\t\t%s\n" % (x[i-1], Error_lo[i-1],Error_up[i-1]))