• 検索結果がありません。

2.1 基礎方程式

基礎方程式として,保存形の2次元圧縮性Navier-Stokes 方程式を用いる.これを基準となる長さL*,音速c0*,密度ı0*, 温度T0*,粘性係数Ò0* を用いて式(2.1)のように無次元化を行う.

(2.1)

ここで,u, vx, y方向の速度,tは時間,pは圧力,Tは温度であり,上付きの添え字*は有次元量を表す.このような無 次元化を施した基礎方程式は式(2.2)のようになる.

(2.2)

Reは基準音速によるReynolds 数であり式(2.3)のようにかける.

(2.3)

また,Îは比熱比で1.4,PrはPrandtl数で定数として扱う.すると,単位体積あたりの全エネルギーeは式(2.4),

圧力pは式(2.5)のように表される.

(2.4)

(2.5)

本研究では,対流項をChakravarthy-Osher の3次精度TVD スキーム,粘性項を2次精度の中心差分で評価した.また,時 間積分は陽解法である2次精度のRunge-Kutta 法を用いた.

2.2 検証計算

作成したプログラムを検証するために,Fig. 2.1のように,一様な超音速流れと平行に置かれた平板に発達する層流境界 層を計算し,理論解と比較した.ただし,理論解はIllingworth-Stewartson 変換による自己相似解を用いた.

2.2.1 計算条件

検証計算を行う際には,標準大気(海面上)を貯気槽状態として,これをMach 数M=1.8に等エントロピー的に加速し た一様な流れの中に,平板を平行に置くものとする.具体的な計算条件は,以下の通りとする.

>一様流条件

−Mach 数M0=1.8

−Prandtl数Pr=1.0

−圧力p0*=1.763×104[Pa]

−密度ı0*=0.3514[kg/m3

−温度T0*=174.8[K]

>無次元化

−状態量:一様流の値を基準として無次元化

−代表長さ:L*=7.230×105[m]

−Reynolds 数:Re=620.2

>計算領域

−格子形状:Fig. 2.2を参照

−大きさ150.5×5.0(無次元)

−格子数185×160

>境界条件

−上流:主流条件

−下流・上方:零次外挿

−下方(平板上):滑りなし,断熱条件

−下方(平板なし):滑り条件

>初期条件

−全領域を一様流状態とする

プラントル数を1.0としたのは理論解との比較のためであり,粘性係数の評価には式(2.6)を用いた.

(2.6)

2.2.2 計算結果

あるxにおける,y方向の座標を式(2.7)を用いて新たな座標≈に変換する.そして,得られた結果のx=150.5におけ る≈方向の境界層の速度,温度,密度分布をプロットしたものが,Fig. 2.3である.

(2.7)

Fig. 2.3には,理論解も同時に示している.図より,作成したプログラムが概ね理論と一致していることが分かる.次に,

理論解との相対誤差を境界層内すべてにおいて調べた結果をFig. 2.4に示す.縦軸は各xにおける最大誤差(%)を対数表 Fig. 2.1 計算領域

においてはı, u, T ともに最大誤差が1%を下回っていることが確認できる.平板の先端付近の大きな誤差は衝撃波による 影響と,境界層内に入っている格子点数が少ないことによると考えられる.また,x=14付近に生じている振動は,Fig.

2.5に示すように,平板先端から生じた衝撃波が計算領域上方で反射し,再び平板に到達しているため生じた.計算領域上 方に課した零次外挿の条件は,計算領域上端の状態が,そのままy方向に続いているというものであり,斜め衝撃波の波 面の傾きまで考慮されていない.このために,計算領域の上端で衝撃波の反射のような現象が生じた.よって,斜めに衝 撃波が入射する境界に零次外挿の条件を課す場合は,その反射が計算の対象に悪影響を及ぼさないよう考慮しなくてはな らない.最後に,計算領域の後端で特にuに関して誤差が大きくなっている原因も,零次外挿の影響であると考えられる.

計算領域後端では大部分が超音速領域になっているので,この境界条件は特に問題はない.しかし,境界層内では亜音速 領域が存在し,零次外挿の条件が問題となる.というのも,このような領域では情報が上流側に伝えられなくてはならな いが,この境界条件のために,情報が一方的に下流側に伝えられているからである.そのため,計算格子を作成する際に は,計算対象と計算領域の後端には空間を大き目にとる配慮が必要となる.

Fig. 2.2 格子形状(平板)

Fig. 2.3: 自己相似解との比較(x= 150.5)

2.3 非粘性圧縮性流れの理論

非粘性圧縮性流れは第2.1節に示した,基礎方程式(2.2)の右辺を0とすることで,一般的に表現することが出来きる.

この式をEuler方程式と呼ぶこともある.ここでは凹角または凸角を流れる場合の理論を示しておく.

まず,超音速流がFig. 2.6(a)のような凹角に沿って曲がる時には斜め衝撃波が発生する.この時流れの偏向角をÙ,衝 撃波角をÇ,衝撃波の上流側の状態量には下付添え字1,下流側は2を用いると,衝撃波前後の関係式は式(2.8)〜式

(2.10)のように表すことができる.よって,上流側の値,M1, p1およびÙが既知である場合には,式(2.8)より,Çを求 め,さらに,式(2.9)からp2,式(2.10)からM2と順に求めることができる.さらに,ここで,断熱流れを仮定すると,

全温が保存されるので,衝撃波前後の温度が決定され,状態方程式を用いることで,密度も決定することができる.

(2.8)

(2.9)

(2.10)

次に,超音速流がFig. 2.6(b)のような凸角に沿って曲がる場合であるが,これは凹角の時とは異なり,等エントロピー 的に膨張し,Plandtl -Meyer 膨張扇を形成する.この流れにおける偏向角ÙとMach 数Mの関係は,Plandtl -Meyer 関数式

(2.11)を導入することで求めることができる.それは流れが等エントロピー的に曲げられる場合には,式(2.3)が成立し,

Fig. 2.4 平板境界層内の最大誤差

Fig. 2.5 平板先端付近の圧力分布

るからである.このようにして,偏向後のMach 数M2が求まれば,等エントロピー変化を考慮して,その他の状態量の導 出が可能となる.

(2.11)

(2.12)

関連したドキュメント