・軸対称ラバルノズルの設計方法

・まえがき

　以前，二次元ラバルノズルの計算でFoelschの方法によるノズル設計法を紹介した．NASAから提供されている技術資料(NACA Technical note 1651)は以下．

http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19930082268_1993082268.pdf

しかし，こいつを読みこなして実際にノズルを設計するのは容易なことではない(と思う)．そこで，Technical note 1651で言及されているFoelschの方法によるノズル設計法について詳しく述べたい．また，前回はスリット状ノズルだったが，今回は軸対称ノズルの設計を行いたい．違いは断面積がrに比例するか，r^2に比例するかで，その思想は同じである．しかし結果として得られる形状は思ったほど直感的でない(スリットノズル形状の平方根を取るだけではだめ)．この辺の理論的背景について述べた資料(AE-TM-17-77)を以下に示す．前回はMathematicaで作ったスクリプトを説明抜きで公開したが，今回はその設計手法を詳細に解説したいと思う．ただし，理論的背景は抜きで，計算アルゴリズムのみの解説にとどめたい．

http://en.scientificcommons.org/15385460

・設計思想

図1はTechnical note 1651の図12を書き直したものである．設計の目的は，ノズルの「始まり部」で生じた膨張波を「終わり部」の曲線で相殺し，結果的にx軸に平行な超音速流を得ることである．「始まり部」においてダクトの断面積変化率(∂^2A/∂x^2)は正であり，「終わり部」では負となる．つまり，「始まり部」と「終わり部」を分けるのはノズル形状関数の変曲点である．変曲点の手前，「始まり部」で発生した膨張波が中心軸を超えて反対側のノズル壁面にぶつかるとき，ノズル壁面で発生するの圧縮波と相殺，消滅するように「終わり部」の曲線を設計するのが「ノズル設計」というわけである．


図1: 基本的設計思想の説明

・設計手順

●設計に当たってまず行うことは，ノズルの設計マッハ数を決定することである．これはアプリケーションからの要求で決まることで，設計上の任意パラメータではない．

●続いて，変曲点におけるノズルの傾きθ1を決定する．これはある程度の自由度のあるパラメータで，θ1を大きくするとノズル全長が短くなるが，ある制限が課せられる．

スリットのズルの場合，ノズル出口のプラントル・マイヤー関数νt(下式に設計マッハ数を代入)の値

　(1)

に対してνt > 2θ1でなくてはいけない[1]．聞きかじりの知識では，M<3位の超音速ノズルにおいてはθ1～νt/4くらいが良いらしい．Technical note 1651によると

　(2)

なる経験式があるらしい．

一方，これが軸対称の場合どうなるかというと，文献[1]のプラントル・マイヤー関数と流線の角度に関する保存則が

　　(3)

になるため，計算するとνt > 4θ1を得る．AE-TM-17-77によると，「θ1は15度を超えてはならず，通常は12度が用いられる」とある．結局，経験則ということか．

●θ1が決まったら，ノズルの「始まり部」の形状を決定する．Foelschの設計では，流線は変曲点より上流では放射状と近似するため，上流の形状に厳密な設計手法はない．x=0におけるノズル半径(y)と変化率(dy/dx)，x=x1におけるy，dy/dx，d^2y/dx^2を与えられるとそれを満足する三次曲線が一意に決定される．Technical note 1651でも三次曲線による設計方法が紹介されているのでこれに従おう．

表1: 始まり部の三次関数

x y dy/dx d^2y/dx2
0 y0 0 0
x1 y1 tanθ1 0

ここで問題は，x1とy1をどのように決定するか，である．普通，Foelschの設計方法ではまず終わり部を計算し，θ=θ1における(x，y)からx1とy1を求める，というのが伝統的な方法である．しかし，これは直感的に分かりにくいし，我々はMathematicaという強力な武器を持っている．ここはまず始まり部を決定してしまおう．前提として以下の仮定を置く．


図2: y1を決定する計算

すると図2の計算によりy1が定まる．ただしここで，A1/A*を知る必要があるが，これは以下の方法で求める．

  1. 軸対称ノズルの場合，θ1とνtは以下の関係で結びつけられる．ここでν1は(x1,y1)点におけるプラントル・マイヤー関数である．
    　　(5)
  2. ここから未知数ν1を解析的に求めることができ，式(1)を陰的に解けば(x1,y1)点のマッハ数M1が求まる．
  3. M1が知れれば，式(4)に代入することでA1/A*を求められる．
  4. x1が未知のまま，表1の条件を満たす三次関数を書き下すと以下のようになる．
    (6)
  5. y=y1を代入し，x1について解けば
    　(7)
    によりx1が定まる．

以上の手続きにより始まり部の三次関数が決定された．

●つづいて，終わり部の設計に移る．Foelschの方法は，終わり部において以下の仮定を置き，プラントル・マイヤー関数と連続の式からノズル形状を決定する．


図3: 終わり部の計算

  1. 変曲点を含む半径r1の円弧において流線は円弧に垂直，かつ一様．この点におけるマッハ数をM1とする．
  2. 変曲点を出発する右向きマッハ線(1)上で，流線はr1と同じ原点を発する円弧状．ただしマッハ数は場所によって異なる．
  3. 上記マッハ線(1)上の各点Pにおけるマッハ数は，ベクトルrの長さだけで決まる．

マッハ線(1)上の点P(x2,y2)から発する左向きマッハ線(2)を考える．このマッハ線がノズル壁と交差する点をQ(x,y)とすると，いま興味はQ(x,y)をθまたはP点のマッハ数Mの関数として求めることになる．ヒントとなる材料は以下．

計算手順は以下のようになる．予め，マッハ数Mにおけるダクト断面積Aを計算するため，放射流領域のM=1におけるベクトル半径r0を求めておく．

　(9)

  1. マッハ数Mを決定．これはM=M1から設計マッハ数Mtまで微小ステップでくり返す独立変数．
  2. マッハ数から半径rを以下の関係から決定．
    　(10)
  3. 式(8)にマッハ数を代入，P点の角度θを求める．
  4. 連続の式から以下の関係が導かれる．
    　(11)
    ちなみに，左辺は弧P'Pが張る，球面を切り取った面積でこいつは簡単．右辺はQQ'が張る傘の面積だが，円錐台の側面積の公式[2]を使う．各部寸法は図4を参照．
  5. 式(11)は独立変数がlのみなので，lについて解くことができる．
    　(12)
  6. lが決定されれば(x,y)は以下のように求められる．
    　(13)
    　(14)
  7. Mを増加させ，手順1-6をくり返す．


図4: QQ'によって張られる円錐台の側面積を求める．

[1]松尾一泰「圧縮性流体力学」，(理工学社, 1994)，p. 257
[2]<美しい数学の話>　第27話「回転体の表面積　http://www2.ocn.ne.jp/~mizuryu/toukou/toukou27.html

・自動設計スクリプト

　以上のように，軸対称ノズルの場合は計算は多少面倒になる．しかし手順さえきちんと定義してしまえばこれを自動化するのは易しい．Mathematicaによるノズル自動設計スクリプトをここに示す．

・テスト計算

　設計マッハ数2のノズルを上記のスクリプトで設計し，Phoenicsで計算させてみた．計算は軸対称2次元，壁面はスリップとする．その他の条件は以下の通り．

表2: テスト計算の条件

設計マッハ数 2.0
モデル流体 分子量28，理想気体
比熱比=1.4
粘性係数=0
メッシュ分割 45×150
スロート半径 1.0mm
質量流量 5.76e-3kg/s(206mmol/s)
外部圧力 100000Pa

メッシュ分割図
ノズル出口のマッハ数分布数ノズル出口の圧力分布
NEP面の圧力，マッハ数分布

綺麗な一様流を得ることができた．計算結果のマッハ数は設計マッハ数に微妙に足りないが，こいつはFoelschの設計手法の近似誤差か離散化誤差のどちらかだろう．

・まとめ