翼型点列データを出力するpythonスクリプト

NACA 4 digit

※ スクリプト中で4桁を指定しています。

# Import libraries
from __future__ import division
import os
import sys
from sys import exit
from math import *

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

import time

# Command line args.
argvs = sys.argv
# Number of args.
argc = len(argvs)

print argc, argvs

## Setup constants (Note that chord length is normalized and scaled at final stage)
# For example, the NACA 4412 airfoils (our tip airfoil) has
#   a maximum camber of 4% located 40% (0.4 chords)
#   from the leading edge with a maximum thickness of 12 % of the chord.
#   4-digit series airfoils by default have maximum thickness at 30% of the chord
#   (0.3 chords) from the leading edge.

# NACA4420
m = 0.04
p = 0.40
t = 0.20

N = 50
x   = [0]*(N+1)
xu = [0]*(N+1)
xl = [0]*(N+1)
yu = [0]*(N+1)
yl = [0]*(N+1)
ycc = [0]*(N+1)
dycdx = [0]*(N+1)
theta = [0]*(N+1)
yt    = [0]*(N+1)

# make X coordinate
dx = float(1)/N
for i in range(N+1):
    # X
    x[i] = 0.5*(1-cos(dx * i * pi))

    # ycc and thete (thickness to be applied perpendicular to camber line)
    if(x[i] <= p):
        ycc[i] = m/p**2 * (2*p*x[i]-x[i]**2 )
        dycdx[i] = 2*m/p**2*(p-x[i])
    else:
        ycc[i] = m/(1-p)**2 * ( (1.0-2*p)+2*p*x[i]-x[i]**2 )
        dycdx[i] = 2*m/(1-p)**2*(p-x[i])
    theta[i] = atan(dycdx[i])

    # yt ---- set last digit 0.10150 to 0.10360 if closed T.E.
    yt[i] = t/0.2*(0.29690*sqrt(x[i])-0.12600*x[i]-0.35160*x[i]**2+0.28430*x[i]**3-0.10150*x[i]**4 )

    # xu, xl, yu, yl
    xu[i] = x[i] - yt[i] * sin(theta[i])
    xl[i] = x[i] + yt[i] * sin(theta[i])
    yu[i] = ycc[i] + yt[i] * cos(theta[i])
    yl[i] = ycc[i] - yt[i] * cos(theta[i])

# Output
print i,x,ycc,yu,yl

plt.figure()
plt.plot(x,ycc,"o-",label="Camber")
plt.plot(xu,yu,"o-",label="Upper")
plt.plot(xl,yl,"o-",label="Lower")
plt.grid(True)
plt.xlim( (0, 1) )
plt.legend()
plt.xticks(np.arange(0,1.1,0.1))
plt.gca().axis('equal')
plt.title("NACA-4 airfoil upper/lower/camber profiles")

plt.show()

修正PARSEC翼型（modPARSEC.py, parsecCoeff.py）

※ データファイルから読み込む形です。

# Import libraries
from __future__ import division
import os

from sys import exit
from math import *

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

import parsecCoeff as pc

import time

### Extermal (12) profile design parameters (to be optimized)
### Note that profiles are normalizes by chord length C

# Read parsec parameters (user input) & assign to array
# I/O files path
path = './'
pp = 'parsec_parameters.dat'
pparray = np.genfromtxt(os.path.join(path,pp), delimiter=",", dtype = float, skiprows=0)

# TE & LE of airfoil (normalized, chord = 1)
xle = 0.0
zle = 0.0
xte = 1.0 # zte is given as parameter

# Number of points of profile
N = 101

# thickness profile 6 parameters (symm. to camber)
rle      = pparray[0]          # L.E.thickness carvature radius (positive)
xt       = pparray[1]          # location of thickness crest(x)
zt       = pparray[2]          # location of thickness crest(y)
zxxt     = pparray[3]          # curvature of thickness crest
betate   = pparray[4]          # T.E. Wedge angle(positive in deg.)
dzte     = pparray[5]          # thickness at T.E.

# camber profile 6 parameters
rc       = pparray[6]          # L.E. camber line carvature radius (positive)
xc       = pparray[7]          # location of camber crest(x)
zc       = pparray[8]          # location of camber crest(y)
zxxc     = pparray[9]          # curvature of camber crest
alphate = pparray[10]         # T.E. camber line angle(in deg.)
zte      = pparray[11]         # T.E. camber vertical displacement

th_coef, cm_coef = pc.pcoef(xte,
                            rle,xt,zt,zxxt,betate,dzte,
                            rc ,xc,zc,zxxc,alphate,zte)

# Evaluate suction (upper) surface points
xx_cm = np.zeros(N,dtype='float')
for i in range(N):
x = i * (xte-xle)/N
xx_cm[i] = 0.5*(1.0-cos(x*pi))

zz_cm = (
          cm_coef[0]*xx_cm**0.5 +
          cm_coef[1]*xx_cm      +
          cm_coef[2]*xx_cm**2.0 +
          cm_coef[3]*xx_cm**3.0 +
          cm_coef[4]*xx_cm**4.0 +
          cm_coef[5]*xx_cm**5.0
         )
dzzdx_cm = (
             0.5*cm_coef[0]*xx_cm**(-0.5) +
                 cm_coef[1]               +
             2.0*cm_coef[2]*xx_cm         +
             3.0*cm_coef[3]*xx_cm**2.0    +
             4.0*cm_coef[4]*xx_cm**3.0    +
             5.0*cm_coef[5]*xx_cm**4.0
           )

theta = np.zeros(N,dtype="float")
for i in range(N):
theta[i] = atan(dzzdx_cm[i])

print theta

# Evaluate thickness profile points (note that thickness profile is symmetrical)

zz_th = (
         th_coef[0]*xx_cm**0.5 +
         th_coef[1]*xx_cm**1.5 +
         th_coef[2]*xx_cm**2.5 +
         th_coef[3]*xx_cm**3.5 +
         th_coef[4]*xx_cm**4.5 +
         th_coef[5]*xx_cm**5.5
        )

# Evaluate upper/lower points
xx_u = np.zeros(N,dtype="float")
zz_u = np.zeros(N,dtype="float")

xx_l = np.zeros(N,dtype="float")
zz_l = np.zeros(N,dtype="float")

for i in range(N):
   xx_u[i] = xx_cm[i] - zz_th[i] * sin(theta[i])
   zz_u[i] = zz_cm[i] + zz_th[i] * cos(theta[i])
   xx_l[i] = xx_cm[i] + zz_th[i] * sin(theta[i])
   zz_l[i] = zz_cm[i] - zz_th[i] * cos(theta[i])

##### make Plot ###################################################

plt.figure()

plt.subplot(211)
plt.plot(xx_u,zz_u,'b',label='upper(suction)')
plt.plot(xx_l,zz_l,'r',label='lower(pressure)')
plt.grid(True)
plt.xlim([0,1])
plt.xticks(np.arange(0,1.1,0.1))
plt.legend()
plt.gca().axis('equal')
plt.title("modified PARSEC airfoil upper/lower profiles")

plt.subplot(212)
plt.plot(xx_cm,zz_cm,'y',label='camber')
plt.plot(xx_cm,zz_th,'g',label='thickness')
plt.grid(True)
plt.xlim([0,1])
plt.xticks(np.arange(0,1.1,0.1))
plt.legend()
plt.gca().axis('equal')
plt.title("modified PARSEC airfoil camber/thickness profiles")

plt.savefig(os.path.join('modified_parsec_airfoil.pdf'))
plt.savefig(os.path.join('modified_parsec_airfoil.png'))
plt.show()

インポートしているライブラリ（？）parsecCoeff.py

# Import libraries
from math import sqrt, tan, pi
import numpy as np

# User function pcoef
def pcoef(xte,
          rle,xt,zt,zxxt,betate,dzte,
          rc ,xc,zc,zxxc,alphate,zte):

    # Docstrings
    """evaluate the modefied PARSEC thickness/camber coefficients"""

    # Initialize coefficients
    acoef = np.zeros(6)
    bcoef = np.zeros(6)

    ### Thickness profile (as basic PARSEC symmetrical aifoil) ###

    # Form system of equations
    acoef[0]   = sqrt(2*rle)
    A = np.array([
                 [xte**1.5, xte**2.5, xte**3.5, xte**4.5, xte**5.5],
                 [xt **1.5, xt **2.5, xt **3.5, xt **4.5, xt **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(xt ), 2.5*xt **1.5, 3.5*xt **2.5, 4.5*xt **3.5, 5.5*xt **4.5],
                 [0.75*(1/sqrt(xt)), 3.75*sqrt(xt), 8.75*xt**1.5, 15.75*xt**2.5, 24.75*xt**3.5]
                ])

    B = np.array([
                 [0.5*dzte - acoef[0]*sqrt(xte)],
                 [zt       - acoef[0]*sqrt(xt )],
                 [-tan(0.5*betate*pi/180) - 0.5*acoef[0]*(1/sqrt(xte))],
                 [                        - 0.5*acoef[0]*(1/sqrt(xt ))],
                 [zxxt + 0.25*acoef[0]*xt**(-1.5)]
                ])

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

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

    ### Camber profile ###

    # Form system of equations
    bcoef[0]   = sqrt(2*rc)

    A = np.array([
                 [xte, xte**2, xte**3, xte**4, xte**5],
                 [xc , xc **2, xc **3, xc **4, xc **5],
                 [1.0, 2*xte, 3*xte**2, 4*xte**3, 5*xte**4],
                 [1.0, 2*xc , 3*xc **2, 4*xc **3, 5*xc **4],
                 [0.0, 2.0 , 6*xc    ,12*xc **2,20*xc **3]
                ])

    B = np.array([
                 [zte - bcoef[0]*sqrt(xte)],
                 [zc - bcoef[0]*sqrt(xc )],
                 [-tan(alphate*pi/180) - 0.5*bcoef[0]*(1/sqrt(xte))],
                 [                     - 0.5*bcoef[0]*(1/sqrt(xc ))],
                 [zxxc + 0.25*bcoef[0]*xc**(-1.5)]
                ])

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

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

    # Return coefficients
    return acoef, bcoef