PythonでNACA4桁翼の断面形状を作図する方法

NACA4桁翼とは、航空機に使われている代表的な翼型の事です。空力性能が良い翼の設計は航空機業界だけの仕事ではなく、翼の断面形状を自由に作図できるようになる事で、色々な解析が捗ります。この記事はPythonでNACA4桁翼を作図する方法を紹介します。

こんにちは。wat(@watlablog)です。ここではNACA4桁翼の断面形状をPythonで作図してみます！

目次（項目クリックでジャンプできます）

NACA4桁翼とは？

航空力学分野で代表的な翼型

NACA4桁翼型(NACA 4-digit series)とは、前縁(Leading Edge)が丸く、後縁(Trailing Edge)が先細っている一般的な航空機用の翼を、4桁の数字で表したものです。他にも5桁型、6桁型があり、4桁型は1930年代に開発された古い翼型です。

NACA4桁翼型のパラメータ

下図がNACA4桁翼型のパラメータです。

翼は先頭が前縁、最後端を後縁と呼びますが、前縁から後縁にまっすぐ引いた線をコード線(Chord line)と呼びます。

コード線を上方向に反らせた線を反り曲線、またはキャンバー線(Camber line)と呼び、このように翼を反らせる事で空力性能を向上させます。

翼は図中\(y\)軸方向に厚みを持っており、横軸\(x\)で一様ではなく分布を持ちます。この分布の事を翼厚分布と呼びます。

作図のためには最大厚さ比\(t\)、反りの最大位置\(p\)、反りの最大値\(m\)をパラメータにとります。

\(t\)は翼のコード線の長さを1とした時の厚み比で、%表示の2桁の数値で表します（20[%]等)。\(p\)もコード線の長さを1とした時の反り位置で指定しますが、4と指定したら0.4を意味します。\(m\)も同様に比ですが、4を指定すると0.04を意味します。

これだけ知っていれば4桁の数字を聞いただけで翼の形状がわかりますね。

NACA4桁翼型の式\(^{[1]}\)

翼厚分布

翼厚分布\(y_{t}\)は式(1)で定義します。

\[ \pm y_{c}=\frac{t}{0.20}(0.29690\sqrt{x}-0.12600x-0.35160x^{2}+0.28430x^{3}-0.10150x^{4}) \tag{1} \]

キャンバー線

キャンバー線は式(2)で定義します。

\[ \left\{\begin{eqnarray} y_{c}&=&\frac{m}{p^{2}}(2px-x^{2})&[0\leq x \leq p]\\ y_{c}&=&\frac{m}{(1-p)^{2}}\left\{(1-2p)+2px-x^{2}\right \}&[p＜x\leq 1] \end{eqnarray}\right. \tag{2} \]

NACA4桁翼型の例

式(1)で定義した翼厚分布の上面と下面に対し、上方向に式(2)のキャンバー線を加算するとNACA翼が完成します。

以下にNACA4桁翼型のいくつかのパターンを示します。

翼型をプログラムで作図する理由

今回はPythonで翼型を作図するプログラムを作る事を目標としています。ここでは翼型をプログラムで作図する理由をお気持ち表明しておきます。

作図の自動化

プログラムで作図する最も大きなメリットは作図の自動化ができるようになるという事です。

作図自動化ができると、CADソフトを使っていちいち形状を手で描かなくても良くなります。

作図自動化はCADソフト上でマクロを組んでしまう方法が一般的だと思いますが、CADソフトのためにマクロ記法を覚えたりする習得コストもそれなりにあると思います。

また、工学系のソフトは意外とPythonコードを使えるものもあります（某電磁界解析ソフトや某最適化ソフト等…多数）。Pythonで色々できるようになっておくとこういうソフトとも連携ができるため便利です。

シミュレーションの自動化

先ほどの続きになりますが、作図の自動化ができると様々なソフトと繋げる事ができます。

そして、翼型は主に流体シミュレーションで様々なパラメータによる揚力や抗力、流れの剥離、再付着位置、騒音…等の多数の検討に使われると思います。

Pythonでパラメータを変更→作図→シミュレーション→評価…と自動で回せると人の工数は激減するはずです。

PythonでNACA4桁翼の断面を作図するコード

動作環境

本プログラムは以下の環境で動作を確認しています。

Windows	OS	Windows10 64bit
	CPU	2.4[GHz]
	メモリ	4[GB]

Mac	OS	macOS Catalina 10.15.7
	CPU	1.4[GHz]
	メモリ	8[GB]

Python	Python 3.9.6
PyCharm (IDE)	PyCharm CE 2020.1
Numpy	1.21.1
matplotlib	3.4.3

全コード（コピペ動作可能）

今回はそんなに行数がありません。以下の全コード内に示したコメントで説明は十分かと思いますが、仮にNACA3415を作図する場合は、
「NACA4(m=3, p=4, t=15, x_resolution=300)」
と対応色文字部分を変更して実行します。

x_resolutionを増やすと断面が滑らかになります。おそらく3Dプリンタで印刷する時にも滑らかさが変わると思います。

「# 後でサーフェスを閉じるために末尾を0にする」や「# z座標を設定」の部分は、次回記事で3Dモデルとする時に使います。

import numpy as np
from matplotlib import pyplot as plt

# 最大反り量m, 最大反り位置p, 最大厚み比tからNACA4桁翼形状を作成する関数
def NACA4(m, p, t, x_resolution):
    # 無次元横軸(0<=x<=1)を作成
    x = np.linspace(0, 1, x_resolution)

    # 公称パラメータを実際の桁に合わせる
    t = t / 100
    m = m / 100
    p = p / 10
    print('t=', t, 'm=', m, 'p=', p)

    # 翼厚分布
    yt = (t / 0.20) * \
         (0.29690 * np.sqrt(x)
          - 0.12600 * x
          - 0.35160 * (x ** 2)
          + 0.28430 * (x ** 3)
          - 0.10150 * (x ** 4))

    # キャンバー線(反り)
    yc_front = (m / (p ** 2)) * ((2 * p * x) - (x ** 2))
    yc_back = (m / (1 - p) ** 2) * ((1 - 2 * p) + (2 * p * x) - (x ** 2))

    # 翼厚分布とキャンバー線から翼形状を作成
    yc = []
    for i in range(len(x)):
        if x[i] < p:
            yc.append(yc_front[i])
        else:
            yc.append(yc_back[i])

    # 翼厚分布と反りを形状に反映
    y_top = yt + yc
    y_bottom = - yt + yc

    # 後でサーフェスを閉じるために末尾を0にする
    y_top[-1] = 0
    y_bottom[-1] = 0

    # Z座標を設定
    z = np.zeros_like(x)

    return x, y_top, y_bottom, yc, z

# NACA4桁翼の断面3D点列座標を計算する関数を実行
x, y_top, y_bottom, yc, z = NACA4(m=3, p=4, t=15, x_resolution=300)


# ここからグラフ描画-------------------------------------
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。
ax1.set_xlabel('x [-]')
ax1.set_ylabel('y [-]')

# スケールの設定をする。
ax1.set_xlim(0, 1)
ax1.set_ylim(-1, 1)

# データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。
ax1.plot(x, y_top, label='Top', lw=1, color='green')
ax1.plot(x, y_bottom, label='Bottom', lw=1, color='blue')
ax1.plot(x, yc, label='Camber', lw=1, color='red')
ax1.grid()
ax1.legend()

# レイアウト設定
fig.tight_layout()

# グラフを表示する。
plt.show()
plt.close()
# ---------------------------------------------------

import numpy as np

from matplotlib import pyplot as plt

# 最大反り量m, 最大反り位置p, 最大厚み比tからNACA4桁翼形状を作成する関数

def NACA4(m, p, t, x_resolution):

# 無次元横軸(0<=x<=1)を作成

x = np.linspace(0, 1, x_resolution)

# 公称パラメータを実際の桁に合わせる

t = t / 100

m = m / 100

p = p / 10

print('t=', t, 'm=', m, 'p=', p)

# 翼厚分布

yt = (t / 0.20) * \

(0.29690 * np.sqrt(x)

- 0.12600 * x

- 0.35160 * (x ** 2)

+ 0.28430 * (x ** 3)

- 0.10150 * (x ** 4))

# キャンバー線(反り)

yc_front = (m / (p ** 2)) * ((2 * p * x) - (x ** 2))

yc_back = (m / (1 - p) ** 2) * ((1 - 2 * p) + (2 * p * x) - (x ** 2))

# 翼厚分布とキャンバー線から翼形状を作成

yc = []

for i in range(len(x)):

if x[i] < p:

yc.append(yc_front[i])

else:

yc.append(yc_back[i])

# 翼厚分布と反りを形状に反映

y_top = yt + yc

y_bottom = - yt + yc

# 後でサーフェスを閉じるために末尾を0にする

y_top[-1] = 0

y_bottom[-1] = 0

# Z座標を設定

z = np.zeros_like(x)

return x, y_top, y_bottom, yc, z

# NACA4桁翼の断面3D点列座標を計算する関数を実行

x, y_top, y_bottom, yc, z = NACA4(m=3, p=4, t=15, x_resolution=300)

# ここからグラフ描画-------------------------------------

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(111)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel('x [-]')

ax1.set_ylabel('y [-]')

# スケールの設定をする。

ax1.set_xlim(0, 1)

ax1.set_ylim(-1, 1)

# データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。

ax1.plot(x, y_top, label='Top', lw=1, color='green')

ax1.plot(x, y_bottom, label='Bottom', lw=1, color='blue')

ax1.plot(x, yc, label='Camber', lw=1, color='red')

ax1.grid()

ax1.legend()

# レイアウト設定

fig.tight_layout()

# グラフを表示する。

plt.show()

plt.close()

# ---------------------------------------------------

実行結果

下図が実行結果で得られるmatplotlibのプロット結果です。緑線が上面、青線が下面、赤線がキャンバー線です。

まとめ

本記事ではNACA4桁翼型の概要と作図に必要なパラメータ、翼厚分布とキャンバー線の式を紹介しました。

最後はPythonによる自動作図関数を紹介し、実際にmatplotlibによるプロットをしてみました。

翼の自動作図ができると航空機業界以外にも風等の流体を扱う人の工数を少し削減できるかも知れません。

次回はこちらの翼をメッシュデータにしていきたいと思います（記事の構成の都合上分けました）。

参考文献

NACA翼については以下の参考文献に詳しく記載があります。本書は航空力学に興味のある方であれば楽しく読めると思うので是非書店で見かけたら手に取ってみて下さい。

航空機設計法-軽飛行機から超音速旅客機の概念設計まで

航空機設計とは，まったくゼロの状態から作業を開始して，いままでにない新しい機能や性能を有する航空機の形を決定していく作業である。本書では，航空機開発のプロセスの一部である航空機概念設計の流れが把握できるよう解説した。

Amazon

楽天

[1] 李家賢一, 航空機設計法-軽飛行機から超音速旅客機の概念設計まで-, コロナ社, pp.109-110, (2017)

今回は作図会！みんなでNACA翼作図しましょう！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

後輩ちゃんより:

2022年1月27日 11:16 PM

ここに書くことではないかもしれませんが…
翼厚の定義は，翼弦に対して垂直にとるべきでしょうか，それともキャンバーラインに対して垂直にとるべきでしょうか．

ps.
y軸のスケールがx軸の2倍になっているので,プロットされる翼型がゆがんだ見えてしまっているようです．

返信
1. wat より:
  
  2022年1月29日 6:53 AM
  
  ご訪問ありがとうございます。
  翼厚の定義ですが、
  本記事で紹介している参考書「航空機の設計法」には特に反り曲線に対する厚みの設定に対する言及はありませんでした。
  仮に厚みを反り曲線に垂直にとる場合は各xで接線や法線方向を求めて幾何計算をしなければならなくなり、複雑となると思います。
  航空機の専門家ではないので、参考書やその他ネットの情報で複雑な式になっていないから、単純にy方向の加算で考えれば良いかな…という感覚です。
  
  スケールのご指摘ありがとうございました（考えていませんでした！）。
  翼形状は見た目も重要と思いますので、記事修正したいと思います。
  
  返信
  1. 後輩ちゃんより:
    
    2022年2月18日 4:22 AM
    
    ご返信ありがとうございます。
    確かに参考書等を見ても翼弦に垂直に取られているものも多いですね
    
    航空機の力学に関する記事は日本語で探しても数が多くないので、参考にさせていただきました！
    ぜひ他の記事も拝見しようと思います。
    
    ご対応下さりありがとうございました。
    
    返信
    1. wat より:
      
      2022年2月18日 5:18 AM
      
      こちらも記事改善のきっかけになりました。
      ありがとうございます。
      他の記事も是非ご覧ください！
      
      返信