JAXA Repository AIREX: ソニックブーム推算に関する委託研究: 平成2123年(20092011年)度研究成果報告書

目

次

第

1

部

総括

... 1

第1章 はじめに ...2

第2章 本委託研究の内容 ...5

1) ソニックブーム推算のための高精度近傍場圧力波形推算技術に関する研究...5

2) 低ソニックブーム設計技術の近傍場実証に関する実験的研究 ...6

3) ソニックブーム長距離非線形伝播CFD解析ツールの開発 ...6

4) ソニックブーム長距離非線形伝播音響解析ツール開発 およびソニックブーム波形に対する大気乱流効果に関する研究 ...7

5) ソニックブーム波形に対する大気乱流効果モデル構築のための実験的研究...9

第3章 総括と今後の展望 ...9

参考文献 ... 11

第

2

部

各研究成果報告（詳細）

... 13

ソニックブーム推算のための高精度近傍場圧力波形推算技術に関する研究

... 14

第1章 研究概要と目的 ...14

第2章 解適合格子細分化手法 ...15

2.1 誤差評価法 ...15

2.2 Hessian行列を用いた座標変換 ...17

2.3 Shock Function ...20

2.4 超音速域でのShock Functionの取扱い ...21

2.5 Edge-based Error Estimationによるフィルタ指標 ...23

2.6 細分化手法 ...23

2.6.1 Bisection Algorithm ...23

2.6.2 非等方性の考慮 ...25

2.7 疎化手法 ...29

2.7.1 細分化レベルの設定 ...29

2.7.2 分割情報の保持 ...30

2.8 細分化領域 ...34

2.9 解適合細分化手順 ...35

第3章 解適合格子細分化結果（静粛超音速研究機モデル） ...37

3.1 解析条件 ...37

3.1.1 解析対象および解析条件 ...37

3.1.2 解適合格子法設定 ...39

3.2 解析結果 ...40

3.2.1 モデル直下方向(Case1) ...41

3.2.2 周方向の解析(Case2) ...42

参考文献 ...47

低ソニックブーム設計技術の近傍場実証に関する実験的研究

... 48

第1章 研究目的及び内容 ...48

第2章 実験装置 ...48

2.1 バリスティックレンジ ...48

2.2 計測装置 ...48

第3章 近傍圧力場取得システムに対する評価実験 ...49

3.1 試験模型 ...49

3.2 実験結果 ...49

3.3 本章のまとめ ...51

第4章 新サボ分離方法の開発 ...52

4.1 従来方式(空気力学的分離) ...52

4.2 分割サボ方式 ...52

4.2.1 設計概要 ...52

4.2.2 性能評価試験 ...53

4.3 空力分離／固体衝突併用方式(ハイブリッドサボ分離法) ...54

4.3.1 ハイブリット分離法概要 ...54

4.3.2 空気力学的分離評価実験 ...55

4.4 本章のまとめ ...60

第5章 D–SEND#1の近傍圧力場計測 ...60

5.1 試験模型とサボ ...60

5.2 近傍圧力場計測 ...61

5.3 近傍場圧力同時4点計測 ...62

5.4 迎角に対する考察 ...64

5.5 試験模型フィンオリエンテーションに対する考察 ...65

5.6 本章のまとめ ...66

第6章 まとめ ...67

参考文献 ...67

ソニックブーム長距離非線形伝播

CFD

解析ツールの開発

... 68

記号表 ...68

第1章 はじめに ...68

第2章 基礎方程式および数値計算法 ...69

2.1 振動エネルギーの保存式を解く場合の基礎方程式(以下，2009モデル) ...69

2.2 振動励起第1準位の化学種保存式を解く場合の基礎方程式(以下，2010モデル) ...70

2.3 数値計算法 ...71

第3章 計算結果 ...71

3.1 平面衝撃波伝播問題 ...71

3.3 2010モデルの平面衝撃波伝播問題での計算結果 ...75

3.3.1 分散波形成の概要 ...75

3.3.2 テスト問題による2010モデルの検証 ...76

3.3.3 H2O-O2 振動-振動エネルギー交換反応の分散波厚さへの影響 ...77

3.3.4 N2-O2 振動-振動エネルギー交換反応の分散波厚さへの影響 ...78

3.4 2009モデルと2010モデル間の分散波厚さの比較 ...79

3.5 有限立ち上がり時間予測精度の検証 ...80

3.6 平面衝撃波が長距離伝播して収束した際の立ち上がり時間データ ...83

第4章 本解析法の応用 ...85

第5章 まとめ ...86

謝辞 ...86

参考文献 ...86

ソニックブーム長距離非線形伝播音響解析ツール開発

およびソニックブーム波形に対する大気乱流効果に関する研究

... 88

第1章 はじめに ...88

第2章 大気減衰効果を考慮した長距離非線形伝播音響解析ツール開発 ...89

2.1 背景 ...89

2.2 推算コードの開発指針 ...89

2.3 計算手法 ...90

2.4 結果と考察...92

第3章 大気乱流モデリングの高精度化に関する検討 ...94

3.1 背景 ...94

3.2 数値計算手法 ...95

3.2.1 大気擾乱モデル ...95

3.2.2 ソニックブーム伝播解析 ...96

3.3 結果と考察...97

3.3.1 大気境界層を模擬した擾乱場におけるブーム強度のばらつき ...97

3.3.2 近傍場圧力値の抽出位置によるブーム強度のばらつき ...98

第4章 低ソニックブーム設計技術実証試験へ向けた軸対称供試体のブーム解析 ... 100

4.1 背景 ... 100

4.2 計算手法と計算条件 ... 100

4.3 結果と考察... 103

4.3.1 近傍場圧力波形 ... 103

4.3.1.1 近傍場圧力係数分布の比較 ... 103

4.3.1.2 飛行マッハ数変化による影響の評価 ... 104

4.3.2 遠方場圧力波形 ... 105

4.3.2.1 飛行マッハ数変化による影響の評価 ... 105

4.3.2.2 計測位置変化による影響の評価 ... 105

謝辞 ... 108

参考文献 ... 108

ソニックブーム波形に対する大気乱流効果モデル構築のための実験的研究

... 111

第1章 背景 ... 111

第2章 実験装置と実験条件 ... 111

第3章 実験結果 ... 113

3.1 一点計測 ... 113

3.2 多点計測 ... 115

ソニックブーム推算に関する委託研究

平成

21～23

年（2009～2011

年）度

研究成果報告書

第

1

部

総括

第１章 はじめに

宇宙航空研究開発機構 航空プログラムグループ 超音速機チームでは，環境適合性に優れた超音速旅客機の設計技 術を確立することを目的とした静粛超音速機技術の研究開発プログラム[1] を推進している．この研究開発においては 第一の目標として，超音速機特有の騒音である超音速飛行時のソニックブーム（図１）を半減させる技術の獲得を目 指しており，ソニックブーム低減化のための機体形状設計技術の研究開発を行っている．言うまでも無く低ソニック ブーム設計の質を左右する最大の要因はソニックブームの推算精度であるため，設計技術と並行して高精度ソニック ブーム推算技術の研究開発も進めている．また，ソニックブーム低減化の目標とすべきソニックブーム基準は未だ決 められておらず（現在はコンコルド時代のままの基準でありソニックブーム強度にかかわらず陸地上空での超音速飛 行は禁止されている），近年国際民間航空機関（ICAO）において進められている将来的なソニックブーム基準策定プ ロセスにおいては，ソニックブームが地上の人や建物に与える影響評価が重要とされているが，その評価においても 正確なソニックブーム推算技術が求められている．これら「低ソニックブーム機体設計」「ソニックブーム推算」「ソ ニックブーム評価」の 3 者で，図２に示す様なフィードバックサイクルを形成することにより，全体としてソニック ブーム低減技術が向上するものと考えられる

図１．ソニックブーム 図２．ソニックブーム低減技術サイクル

ソニックブーム推算においては超音速機から発生する衝撃波が地上に届くまでを推算するわけであるが，機体から の距離Hと機体長Lの比に応じて下記の3つの領域に分けて扱われることが多い．

（１） 近傍場（H/L<10程度） ： 衝撃波の発生（圧縮性空気力学）

（２） 中間場（H/L=10~数100） ： 機体近傍から地上付近までの音の伝播（長距離非線形音響伝播）

（３） 遠方場（大気境界層内） ： 地上付近や建築物周りの伝播（大気乱流，地面効果，音響／構造連成）

これらを図に表すと図３の様になり，まず機体近傍場において機体から発生する衝撃波が形成する機体近傍場圧力波 形を推算し，これを入力として中間場において長距離非線形音響伝播解析を行い，地上におけるソニックブーム圧力 波形を推算する．ただし，地表から1~2kmの高度までは大気境界層と呼ばれる大気乱れが強い領域が存在し，ソニッ クブーム圧力波形にも大きく影響することが知られているため（図４），大気境界層内の伝播における大気乱流の効果 は別途扱う必要がある．また，ソニックブーム評価において重要な要因となる室内におけるソニックブームを推算す るには，地上屋外におけるソニックブーム波形を入力として音響／構造連成解析を行い，室内への音響透過及び室内 反響，あるいは構造振動による 2 次的な効果などを推算する必要もあるため，建築物との干渉も別に取り扱われるこ とが多い．

以下に各領域における推算手法の現状と課題について紹介する．

低ソニックブーム 機体設計

ソニックブーム推算

ソニックブーム 評価 低減目標・

図３．ソニックブーム推算

図４．大気乱流によるソニックブーム音圧波形の変形[2]

＜機体近傍場圧力波形推算＞

ソニックブーム推算における機体近傍場圧力波形推算には，Euler方程式やNS方程式など衝撃波の非線形性を正し く評価できる流体解析を用いて機体近傍場の流れ場を解く必要があるが，機体から発生する衝撃波の伝播を正確に捕 らえるには計算格子を密にすることによって数値粘性による減衰を抑えなくてはならず，計算コストが増加する傾向 にある．一方で最適設計等を用いた低ソニックブーム設計のためには，実機の複雑形状に対応可能でかつ計算コスト を抑えた実用的な解析ツールが求められる．これまで JAXAでは機体近傍場圧力波形推算法として，線形パネル法の 推算近傍場波形に衝撃波に起因する非線形波形変形効果としてAging補正を加えた手法[3]や，3次元機体形状をマッ ハ面で切断した断面を機軸方向に投影した面積と，機軸方向の揚力分布と等価な影響を有する断面積との和として定 義される等価断面積分布を有する軸対称物体（等価軸対称物体）に置き換えた上で軸対称CFD解析を行う手法，ある いはシングルブロック構造格子を用いて円錐状の格子を作り衝撃波をシャープにとらえる手法や，そこに重合格子法 を適用して複雑形状に対応させる手法[4]等を用いて低ソニックブーム最適設計を実施してきており，更に複雑形状に

大気境界層 （ 地表から上空1～2km）

（ 大気乱流あり）

①近傍場

機 体 形状

近 傍 場圧力波形

②中間場

近 傍 場圧力波形

大 気 境 界層上方圧力波形

③遠方場

大 気 境 界層上空圧力波形

地 上 圧力波形

大 気 乱流による変形 長 距 離非線形伝播

高 精 度CFD

④室内ブーム

地 上 屋外圧力波形

室 内 圧力波形

対してはマルチブロック構造格子による UPACS 解析で対応してきたが，近傍場波形推算精度の高さと，複雑形状へ

の対応，そして計算コストの全ての要求を満たすことは困難であった．一方米国においては，NASAが2008年にNASA

内部でソニックブーム推算ワークショップを開催し，非構造格子法や直交格子法にAdjoint感度解析を利用した解適合

格子法を適用した推算法の開発を進めており，我が国においても同様な取組みが求められる． ＜長距離非線形音響伝播解析＞

ソニックブーム伝播解析技術では，従来よりWhithamのF関数法[5]を数値解析的に解く手法である波形パラメー

タ法（Thomas 法）[6]が用いられてきた．Thomas 法は幾何音響理論をベースとしており，近傍場圧力波形を入力と

して 4 本の波線で構成される波管に沿って計算を進めるため，大気の層状効果を取り入れることが可能であるが，衝

撃波の非線形性による波形変形にWhitham のAging補正を適用しているため，衝撃波は常に不連続な圧力上昇とし

て扱われる．一方，地上におけるソニックブームの影響評価においては，圧力上昇の立ち上がり時間が重要となるた め，音圧変化量と立ち上がり時間の双方を含めたソニックブームの圧力時刻歴波形の正確な予測が求められる．そこ で，音響学的なアプローチによるソニックブーム伝播解析として，大気の非均一性による回折効果を考慮した伝播経

路を音線法を用いて求め，その伝達経路に沿ってKZK方程式（あるいは回折が生じない平面波に対する特別なケース

としてBurgers方程式）[7]を用いて大気の散逸効果も考慮した有限振幅音波の1次元非線形伝播解析が欧米で行われ

ている．特に欧州では音線法を高度化し，ブームカーペットの分布やセカンダリブームに関する研究が行われている．

米国においてはソニックブームの 2 次元的な伝播解析への取り組みも行われ始めているが（2008 年段階），計算負荷

が大きいために解析可能な伝播距離には大きな制限がある． ＜近傍場／中間場マッチング＞

上述のソニックブーム推算法においては，通常，近傍場で得られた圧力波形を入力波形として，直接長距離非線形 音響伝播解析を実施するが，より正確に近傍場と中間場のマッチングを行う手法として，米国において多重極分布法

（Multipole法）[8]が提案されており，CFD解析結果から得られる近傍場圧力波形分布を多重極音源近似することに

よって，より正確に長距離非線形伝播解析に受け渡す手法が用いられている． ＜大気乱流効果推算＞

大気乱流がソニックブーム波形に与える影響は，先に見た通り波形がスパイク状になって圧力値が高くなったり，

逆に鈍って低くなったりと，気象条件等によりその結果は大きく変化する．従って決定論的に推算すると言うよりは，

統計的手法により変動の影響を評価する傾向にある．Crow[9]は1次散乱方程式から，地上付近の慣性散乱と熱散乱が

ソニックブーム波形に与える影響について検討し，波形の変動範囲を示した．また近年，米国では数多くの実フライ ト試験において大気境界層より上空の空中で計測されたブーム波形と，地上で観測されたブーム波形のデータを収集

し，大気乱流フィルタを構成する試み[10]が行われている．また欧州では，実験室において大気乱流を模擬した媒質中

を伝播する衝撃波を観察する実験的取り組み[11]も行われている．

＜音響／構造連成解析＞

ソニックブームの構造物への影響評価に関しては，FAAとNASAが2006年に実施したPARTNERプログラムにお

けるLow-Boom/No-Boom飛行試験（F-18）の結果[12]から，室内において観測されるソニックブームと建屋の振動が

重要であることが指摘されており，NASA では計測データを元にソニックブームが引き起す建屋の振動現象の解明を

進めており，NASA 資金による音響／構造振動連成の数値解析も近年行われてきている．音響と構造振動の数値解析

は，低周波帯域に対しては有限要素法などが実用段階にあるが，ソニックブームの急激な圧力上昇による影響を考慮 するにはより高い周波数帯域に対する解析が必要であり，計算効率を向上させる定式化や連成方法が検討されている．

以上の様な状況の中で，JAXAは平成21年度から3年に渡って国内の大学に計5件の委託研究を依頼し，ソニック

第２章 本委託研究の内容

本委託研究においては，JAXA 内で実施が困難，あるいは大学において長期的に実施すべき項目として，実験ベー

スの研究2件と数値解析手法開発3件の計5件の研究を実施した．以下にそれぞれの委託研究の概要と主要な成果に

ついてまとめる．各々の詳細な報告については第2部に報告書を載せる．

1) ソニックブーム推算のための高精度近傍場圧力波形推算技術に関する研究

研究分担者：

東北大学・教授・中橋和博（代表者） 東北大学・助教・佐々木大輔

本研究では，数値流体力学(CFD：Computational Fluid Dynamics )を用いた機体近傍場における圧力波形推算精度

の向上を目的とし，非構造格子CFD（TASコード）に対して解適合格子細分化法の開発を行う．解適合格子細分化を

用いることで，適切な指標に基づき物理量変化の大きい領域の細分化を行い，局所的に解の精度を高めることが可能 である．衝撃波波面の同定に用いられる衝撃波関数を指標として用いることにより，衝撃波を精度良くとらえること が可能であるが，適切なフィルタを用いることによる数値誤差の影響低減も必要である．また，格子細分化を繰り返 すことによる格子数の大幅な増大は解析速度を遅くしてしまうため格子数の不必要な増大は避けるべきである．

そこで本研究では，解析コストの増加を抑えた精度の改善のために，新たな細分化指標の提案および導入，等方的 な分割法である二分割アルゴリズムの非等方細分化への拡張，疎化による不要な格子点の除去を取り入れた．

細分化指標は空間格子の格子解像度不足による，圧力変動の減衰の影響を受けにくいものとして，衝撃波関数を細

分化指標として適用した．超音速域での数値誤差の影響を取り除くフィルタ指標にを圧力勾配絶対値から，圧力の 2

階微分値と格子の辺の長さから算出されるEdge-based Error Estimation(Edge Error)を取り入れた．本研究で取り入

れた細分化指標により，格子解像度が不足する格子において細分化判定が可能な領域の改善を実現した．また，フィ

ルタ指標の特性により，圧力変動に応じた格子サイズ(最小格子幅)が与えられ過度の細分化を抑制するため大幅な格子

点数の低減効果を得た．

細分化手法では，特定方向への圧力変動が大きい(非等方性が強い)衝撃波や膨張波に対して効率的に細分化を行うた

めに，等方的な分割法である二分割アルゴリズムに対して流れ場の圧力変動に応じた座標変換を用いて非等方細分化 へ拡張した．非等方細分化の導入により非等方性が強い領域の最小格子幅が局所的に小さく設定され精度の改善効果

が得られた．また，格子点数に対しては3割強の低減効果を得た．

疎化手法では，細分化を進める過程で生じた必要以上に追加された格子点を取り除くDerefinementを導入し，格子

点数が2割弱低減する効果を得た．その一方で，疎化を用いることで波形のピーク値が減衰する結果となった．

本手法での三次元形状に対する検証は，翼付き単純形状であるDeltaWing形状と実機複雑形状である静粛超音速研

究機モデル(S3TD)を用いて行った．DeltaWing形状での検証では，既存の近傍場解析手法であるCart3Dとの結果と

よく一致しており，三次元性に対して良好な結果を得た．その一方で，最遠方の圧力波形取得位置において格子点数

増加の抑制のために行ったフィルタ閾値引き上げの影響により波形が減衰する結果となった．S3TD形状での検証では，

HexaGridにより生成された計算格子との比較を行った．波形形状は良く一致するもののフィルタ閾値の影響と考えら

れる減衰が生じる結果となった．複雑形状での検証結果よりフィルタ閾値の影響は圧力変動が小さい領域で強く現れ る傾向が得られた．本解析では，格子点数の制御にフィルタ閾値の引き上げで対応したが，幅広い強さの圧力変動を 生じる実機形状で適切に解析を行うにはフィルタ閾値の過度の引き上げは波形推算精度への影響を及ぼすといえる． 今後，推算精度への影響が小さい格子点数の制御方法の検討が必要であるといえる．

とで，境界層や後流のように異方性の強い現象に対応できると考えられる．

2) 低ソニックブーム設計技術の近傍場実証に関する実験的研究

研究分担者：

名古屋大学・教授・佐宗章弘（代表者）

本研究では，低ソニックブーム設計技術を機体近傍場において実証するための実験システムの獲得，及び機体近傍 場圧力波形推算技術検証のための実験データ取得を目的とした，バリスティックレンジを用いた自由飛行体近傍場の 圧力波形計測技術の研究を行い，以下の成果が得られた．

バリスティックレンジを用いたモデル射出試験において，自由飛行する飛行体に対し，近傍場の圧力波形（実験室 定点で測定した圧力履歴）計測を可能とする手法を開発した．飛行体のロール運動を抑えるため，矩形断面バリステ ィックレンジによる三次元物体射出のためのサボを設計した．低ブーム軸対称物体の飛行体モデルを製作し，自由飛 行状態における近傍場圧力波形計測及び，高速度カメラを用いたフレーミングシュリーレン可視化計測を行った．

D-SEND#1の射出実験をめざし，新しいサボ分離方式を検討した．具体的には，①これまでの管内空力的サボ分離

法では，サボ分離部でのモデル姿勢が問題となっていたのに対して，サボを駆動ガスシール部とモデル支持部に分割 し，サボ分離部内でもモデル姿勢が変化しない方式，および②サボを分割せず部分的な分離管内空力分離とマズルで

の固体衝突を併用した方式(ハイブリットサボ分離)を検討した．①については，マズルでの衝突時にモデル自体の姿勢

が変わってしまうことがわかった．②については， ハイブリッドサボ分離を行う上で，どのような条件で模型及び

サボを設計し試験を行えばよいかの指針を得るために，サボ分離距離評価実験を行った．その結果，まず試験模型の 飛行マッハ数は駆動ガス圧力及び試験模型とサボの質量の和でから求めることができることが分かった．求まった飛 行マッハ数と駆動ガス圧力からサボ前圧力，飛行マッハ数と試験部圧力からサボ背面圧力が求まり，その差圧から分 離距離が予測できることが確認できた．この経験則に基づきサボ及び試験模型を設計し試験条件を決定する事で，サ

ボ分離部下流端で試験模型はサボから分離し，ハイブリッドサボ分離が有効に行える．これを用いて迎角1deg.以内で

試験模型を射出でき，また高い確率で試験模型を姿勢よく射出できることが可能になった．

D-SEND#1の試験模型を射出し，その近傍圧力場を取得することに成功した．そして以下の知見が得られた．

・近傍圧力場について，迎角が付くことにより後端の圧力波形のプロファイルが大きく変化する． ・近傍圧力場について，オリエンテーションにより概ね一致はするものの多少の相違が発生する．

3) ソニックブーム長距離非線形伝播CFD解析ツールの開発

研究分担者：

名古屋大学・准教授・酒井武治（代表者）

ソニックブームが大気中を伝播する際には，大気による様々な効果を受けて波形が変形されるため，ソニックブー ム推算技術においてはその効果を統合的にモデル化できることが望ましい．このうち，いわゆる分子振動緩和による 波の変形や減衰は，大気の主要構成分子である酸素と窒素の振動エネルギーモードが励起されることでおこる．励起 過程に重要な役割を果たすのは相対的には微量な水や二酸化炭素である．特に，分子振動緩和は，地上付近の大気乱 流層に入る前までのブーム波面形成において重要な役割を荷なう物理メカニズムの一つであり，分子振動緩和による

立ち上がり時間増加に強く影響する．従って，そのメカニズムの理解に加え，大気乱流に起因するものと切り分けて，

な状況にある．

このような背景を踏まえ本研究では， CFD法をベースとした弱い衝撃波背後での分子振動緩和過程に関する非平衡

解析手法の開発を行う．先述したように，大気の音の吸収に関する遷移交換反応速度係数が比較的よく整備されてい

る．加えて，極超音速流れの高温反応効果のCFD解析でこれまで開発されてきた手法を適用することで，メカニズム

に関する詳細な情報が得られるのではないかと期待できる．これまでにも，流体力学方程式に Landau-Teller 型の分

子振動エネルギー緩和方程式を加えて解析する手法の開発が行われ，弱い衝撃波の分散過程について解析が行われて いる．本研究では，このような解析手法の開発に加え，詳細な振動励起反応モデルを組みこんだ非平衡計算コードの 開発を行い，湿度に関する影響をより詳細に理解することをターゲットとした．そして，弱い衝撃波が長距離伝播す る際の変形過程を主にシミュレーションし，波の空間分布変化の特性および定常状態における分散波の厚さや立ち上 がり時間を定量的に評価する．

開発された CFD コードを用いて，静止一様大気中を長距離伝播する弱い衝撃波の波動-実在気体干渉効果を予測し

た．本解析手法では，伝播する波の座標系で基礎方程式を解くことで，固定座標系での解法に比べて，衝撃波の分散

過程が 10～20km 程度で定常状態に収束するような条件についても，大幅に計算時間を低減しかつ現実的な時間内で

流れ場のシミュレーションが可能である．また，本解析により得られた，振動緩和によって変形した分散波の構造は，

厳密解と非常によく一致した．さらに，本解析法を球状衝撃波の伝播問題に適用し，有限立ち上がり時間を程よく予 測できることがわかった．

振動エネルギー保存式にLandau-Teller型の振動エネルギーソース項を含めた形で熱的非平衡を解く解析法（2009

モデル）に加え，本研究では，音波の強さが分子振動緩和により弱められる際の詳細な分子間の振動励起反応機構を

組み込んだモデル（2010モデル）の開発も行った．同様な平面波の長距離伝播の数値シミュレーションを行うことで

わかったモデル間の違いは，非常に湿度の低い空気中を伝播するケースについてのみ現れることである．特に本解析

条件においては，湿度が0%に近づくにつれて，定常状態の分散された波厚さは最大で10倍程度異なる：2009モデル

で計算した分散された波厚さの方が薄い．これは，湿度が非常に低い場合に限り，Bassらのモデルで評価した実効緩

和時間が短くなることを示唆する．

大気状態で観測される典型的な相対湿度に対しては，2010 モデルによる分散波厚さの予測値は，2009 モデルより

も約 30%高い．この傾向は，ソニックブームが完全に分散された波の構造を持つと仮定した場合，ソニックブームの

立ち上がり時間が増大することを意味する．従って， Bass らの実効緩和時間モデルをベースとしたブーム伝播予測

手法との比較を想定すると，2010年モデルの結果は，立ち上がり時間の上限値を与える．

4) ソニックブーム長距離非線形伝播音響解析ツール開発

およびソニックブーム波形に対する大気乱流効果に関する研究

研究分担者：

東北大学 流体科学研究所・教授・大林茂（代表者）

超音速旅客機の飛行高度は，英仏が開発したコンコルドの飛行高度を参考にすれば約18 [km] (≈ 60,000 [ft]) であり，

超音速機から発生した衝撃波はこの厚みの大気層を伝播することになるため，ソニックブーム推算における長距離非 線形音響伝播解析では，実在大気伝播中に大気の状態がソニックブーム波形に及ぼす影響を考慮することが重要であ

る．また，大気の状態によってソニックブーム強度ΔP (過剰圧) や波形の立ち上がり時間 (Rise time) が変動するこ

本研究ではソニックブーム高精度推算に関する委託研究に関して次に挙げる 3 項目の研究を実施し，それぞれ以下 に示す成果が得られた．

①大気減衰効果を考慮した長距離非線形伝播音響解析ツール開発

波形パラメータ法において，周波数領域で大気減衰効果を考慮した推算コードの開発を試みた．数値実験の結果か ら，立ち上がり時間ゼロであったソニックブーム波形に対して大気の減衰係数を利用することで，有限な立ち上がり 時間をもつブーム波形が得られた．ブーム波形の形は全体が丸みを帯びた波形となり，ソニックブーム強度も大きく 低減した．結果から，低湿度領域をソニックブームが伝播することで，騒音値は減少すると予想される．今後の課題

としては，他に開発されている実在大気効果を考慮可能なソニックブーム伝播解析手法 (SHOCKN，ZEPHYRUS，

THOR など) や文献と推算結果の詳細な検証が必要である．また，本手法を元にブーム強度や立ち上がり時間の変化

が，ソニックブームの感覚量評価へどのように影響を及ぼすか，新たな知見の発見につながることが期待できる．ま た，大気減衰項の計算周期についてパラメトリックスタディを行った結果から，本モデルには大気減衰項の計算周期

に対する依存性があると考えられる．この依存性と適切な計算周期の設定についても引き続き検討する必要がある．

②大気乱流モデリングの高精度化に関する検討

開発中の幾何音響理論をベースとした大気擾乱を考慮するソニックブーム伝播解析コードに関して数値実験を実施 した．本モデルは大気擾乱によるブーム強度のばらつきを定性的に予測可能であるといえる．大気乱流モデリングの 高精度化については二つの方法が考えられる．一つはソニックブームの波形変形を考慮した高精度推算手法を構築す ることである．これには任意の計測点に影響を及ぼす複数の波形を，少なくとも大気境界層中は伝播解析を行う必要 があり，計算コストが問題点としてあげられる．二つめは，乱流と衝撃波の干渉に関する実験結果から一般的な渦ス ケール，乱流速度変動，そして伝搬距離とソニックブーム強度や立ち上がり時間との関係を明らかにし，モデルを構 築することである．問題点としては，風洞実験で発生させる乱流スケールと実際の大気境界層に存在する渦スケール との差である．この取り扱いについては今後議論を深める必要がある．両手法とも実験結果を統計的に評価し，大気 乱流効果として検証が可能となるよう十分な試験回数を確保する必要がある．とくにモデルの検証が可能な十分な量 の屋外実験でのデータ獲得が望ましい．

また本研究では，一様大気擾乱場を乱流の統計理論に基づき，Von KarmanとPaoが示したエネルギースペクトル

から求めている．空間に一様に不規則な擾乱速度を作り出すためにランダムフーリエモードを用い，波数ベクトルや 位相などに擬似乱流を用いている．このランダムフーリエモードによって作り出された擾乱速度は，完全に不均一な

ものではなく，局所的に偏りが生じやすい．本研究では一様大気擾乱場を，ランダムフーリエモードを用いて生成後，

そのままの状態でブーム推算に使用した．従って，擾乱場には局所的に擾乱速度が大きい箇所や小さい箇所が存在す る．このためブームが伝播中に擾乱の影響を受ける量が，ブーム伝播距離に比例するとは必ずしも保証できず，今回 の結果が得られた一因と考えられる．擾乱速度の偏りを解消し，空間に一様に擾乱速度をもつ擾乱場を生成するため には，ランダムフーリエモードで得た擾乱場に対し，一度流体計算を行うことが必要となる．しかし推算に用いる擾 乱場空間は非常に大きいため，その計算コストが課題となる．

③低ソニックブーム設計技術実証試験へ向けた軸対称供試体のブーム解析

JAXA D-SEND#1で使用される軸対称供試体についてブーム解析を行った．まず，近傍場圧力係数分布の比較を行

ったところ，LBMはNWMに比べ，機体胴体側面部での圧力上昇が小さく，ソニックブームの低減効果が見てとれた．

また，各供試体の近傍場圧力波形を比較した結果，NWMは典型的なN型波形を発生したのに対して，LBMは低ブー

ム波形に特徴的な，いわゆるFlat-top型の波形を発生した．さらに，1.4 ≦ M

∞ ≦ 1.7の範囲において，各供試体の

飛行マッハ数の変化がそれぞれの近傍場圧力波形に与える影響について評価を行ったところ，供試体の飛行マッハ数 の変化が近傍場圧力波形に与える影響は小さいことがわかった．

体は定常落下するものと仮定した．まず，飛行マッハ数がM

∞ = 1.4およびM∞ = 1.7の2通りで変化した場合に遠方

場圧力波形に与える影響を評価したところ，M

∞ = 1.7の場合はM∞ = 1.4の場合に比べて，わずかに圧力ピーク値が

大きくなるものの，飛行マッハ数の変化が遠方場圧力波形に与える影響は小さいことがわかった．しかし，低ブーム

波形に特徴的なFlat-top部の持続時間は，M

∞ = 1.4の場合の方がわずかに長かった．このことからD-SEND#1試験

で重要となる低ブーム波形の評価を行う際には，この低ブーム波形に特徴的なFlat-top部の持続時間の長いM

∞ = 1.4

の場合の方が，試験条件として適していることがいえる．さらに，計測地点の供試体落下地点からの距離を変化が遠

方場圧力波形に与える影響を評価したところ，N型波形と低ブーム波形を識別するためには，落下地点からの距離R =

10 [km] 以内の範囲で計測を行うことが望ましいことがわかった．

なお遠方場圧力波形を推算する際に，供試体は定常落下するものと仮定した．しかしながら，実際の試験時には供 試体の落下は加減速運動である．今後は，この供試体の加減速運動が遠方場圧力波形に与える影響について評価を行 うことで，より実際の試験条件を反映した解析を行うことが課題である．さらに，バリスティックレンジ試験等で得 られる実験結果と比較し，本研究で得られた解析結果の妥当性の検討を行う必要がある．

5) ソニックブーム波形に対する大気乱流効果モデル構築のための実験的研究

研究分担者：

名古屋大学・教授・佐宗章弘（代表者）

現在ソニックブームの強度を評価する際には，擾乱の無い標準大気を用いて計算されることが多いが，正しくソニ ックブームの強度を評価するには大気擾乱の影響を考慮する必要がある．何故なら機体近傍から生じた圧力波が，大 気中を伝播する際に大気乱流や温度変化の影響を受け，圧力波形が変形することが実機を用いた試験より報告されて いるからである．

そこで本研究では，等方性格子乱流と球状衝撃波を干渉させ，干渉した衝撃波を一点計測することで衝撃波背後の 過剰圧ピーク値の変調について統計的実験データを得た．等方性格子乱流と球状衝撃波を干渉させ，干渉した衝撃波 を一点計測することで衝撃波背後の過剰圧ピーク値の変調について統計的実験データを得た．この実験により，圧力 波が大気中を伝播することで圧力波形に変化が生じ，過剰圧のピーク値が大きく変化することが確認された．これは 地上で観測されるソニックブームの強度が擾乱の無い条件で予測した値と異なることを示唆しており，低ブーム機体 の設計において圧力波が大気中を伝播する際の大気の影響も考慮する必要があると言える．そこでソニックブーム強 度をより正しく予測するには，大気の影響により過剰圧がどれほど変化するかを知っておく必要があることから，本 実験では多点計測を行い，場所による過剰圧ピーク値の違いや到達時間の違いを調査した．

本研究で実施した等方性格子乱流と球状衝撃波の干渉実験により，以下の知見が得られた．

・等方性格子乱流と干渉した衝撃波の過剰圧ピーク値の標準偏差と速度変動r.m.s値には比例の関係がある．

・等方性乱流と干渉した衝撃波の過剰圧ピーク値と到達時間の間に弱い正の相関があり，格子乱流中を伝播した 衝撃波の到達が遅れるほど，過剰圧のピーク値が増加する傾向にあることが確認できた．

第３章 総括と今後の展望

本委託研究により，ソニックブーム推算技術向上に寄与する機体近傍場圧力波形推算ツール，長距離非線形音響伝 播解析ツール，及び大気乱流がソニックブーム波形に及ぼす影響を評価する手法が開発・高精度化されるとともに， それらを検証し得る機体近傍場圧力波形計測データ，空気分子の振動励起反応モデルを組みこんだ非平衡計算コード， 及び実験室レベルでの衝撃波／乱流干渉データが得られた．

た非構造格子法の衝撃波捕獲性能が向上されたことで，実用的な推算ツールが整備できたものと考える．今後は，JAXA

で開発を進めているハイブリッド格子CFDツールHexagrid/FaSTARや構造（UPACS）/非構造（TAS）重合格子法

CFD解析等，他の推算手法や，風洞試験データ及びバリスティックレンジデータ等の実験データとの比較を通して，

その推算精度について検証を進め，実用性を高めるとともにTASの通常オプションとして実装されてゆくことが期待

される．

バリスティックレンジによる機体近傍場波形取得技術においては，空気力学的分離と個体衝突分離を組み合わせた ハイブリッドサボ分離手法の開発により，軸対称物体に関しては高い確率で姿勢よく射出できることが可能となった

ことで，検証データ取得手法としての価値が高まったと言える．現在（2012年），本研究で得られたD-SEND#1供試

体形状（NWM, LBM）の近傍場圧力波形を用いて，ワークショップ形式で近傍場圧力波形推算手法を検証する活動が

JAXA を中心に進められており，本研究で獲得した試験技術がソニックブーム推算技術の高精度化にますます活用さ

れてゆくものと思われる．今後は，本手法を 3 次元機体形状模型に適用した上で，軸対称模型同様の姿勢角制御精度

を達成するための技術開発や，M=1.2~1.6 程度の低マッハ数において先行衝撃波との干渉を避ける技術，あるいは圧

力センサーの物理的サイズと模型サイズの制約の中で計測圧力波形の解像度を向上する技術などが求められる．

CFDベースのソニックブーム長距離非線形伝播解析ツール開発では，極超音速域での高温反応効果の CFD 解析手

法を，弱い衝撃波背後での分子振動緩和過程に適用し，詳細な振動励起反応モデルを組み込んだ非平衡解析を行うこ

とにより，弱い衝撃波の立ち上がり時間を定量的に評価し，特に湿度に関する影響を調べた結果，Bassらの提唱する

半経験的な湿度効果モデルでは，低湿度条件において誤差を大きく含む可能性を示した．通常の長距離非線形音響伝 播解析に用いられる音響学的手法の結果を，より詳細なモデルを適用した流体力学的手法により検証したことは非常 に意義深いと言える．今後とも引き続き，他の音響学的伝播解析ツールの検証コードとして活用されることが期待さ れる．

Thomas法に周波数領域での大気減衰効果を加えた長距離非線形音響伝播解析ツールの開発では，現在最も普及して

おり非常に実用的な Thomas 法に，その欠点である立ち上がり時間推算機能を付加するツールを開発したことで，米

国の ZEPHYRUSに相当する時間領域／周波数領域組み合わせ解析手法を獲得できたものと考える．今後は，本ツー

ルの実用性（計算負荷の小ささ）についても念頭に入れつつ，時間領域解析と周波数領域解析の頻度や，波形の位相

に関する大気減衰の影響等について検討を進め高精度化を図るとともに，D-SEND#1等の実験データに対する適用を

通して推算精度検証を行い， JAXAが開発しているBurgers方程式ベースの時間領域非線形音響伝播解析ツールとの

使い分けが可能な実用ツールとして位置づけられることが期待される．

Thomas法に大気擾乱効果のモデルを組み入れた伝播解析ツールの開発では，国内唯一と言ってよいソニックブーム

に及ぼす大気乱流効果推算ツールを開発して数値実験を行い，圧力ピーク値の変動を統計的に調べた結果フライト試 験での地上測定結果と定性的に良い一致を示したことは重要な意味を持つと思われる．今後は，更に現実に近い大気 擾乱場モデルを構築するとともに，上述の周波数領域での大気減衰効果推算と組み合わせることにより，圧力ピーク 値だけでなく立ち上がり時間も含めた詳細な波形変形効果についても調べてゆくことが期待される．

球状衝撃波と格子乱流の干渉に関する実験的研究においては，クイックピストン式の無隔膜衝撃波管による衝撃波 発生システムを構築し，統計データ取得に十分な実験回数を現実的な時間で得られる様になったことは，実験データ からより明確に干渉効果を抽出する上で非常に重要な成果であった．そして実験結果から，過剰圧ピーク値の標準偏

差が速度変動r.m.s値に比例するという結果が得られたことは，本実験手法により大気乱流とソニックブームの干渉効

参考文献

[1] 村上哲，「静粛超音速研究機技術の研究開発について」第 40 回流体力学講演会/航空宇宙数値シュミレーション技

術シンポジウム論文集，pp.79-82 , 2008

[2] Plotkin, K. J. and Maglieri, D. , “ Sonic Boom Research: History and Future.” AIAA 2003-3575, 2003.

[3] Y.Makino and I.Kroo, “Robust objective functions for sonic-boom minimization”, Journal of Aircraft, Vol.43,

No.5, pp.1301-1306, 2005.

[4] Y.Makino, et al., “Nonaxisymmetrical fuselage shape modification for drag reduction of low-sonic-boom airplane”, AIAA Journal, Vol.41, No.8, pp.1413-1420, 2003.

[5] Whitham, G.B., “The Flow Pattern of a Supersonic Projectile,” Communications in Pure and Applied

Mathematics, Vol.5, 1952, pp.301-348.

[6] Thomas, C.L., “Extrapolation of Sonic Boom Pressure Signatures by the Waveform Parameter Method,” NASA TN D-6832, June 1972.

[7] 鎌倉友男, 「超音波における非線形現象」数理解析研究所講究録, 993巻, pp.1-10, 1997.

[8] Page, J. A. and Plotkin, K. J., “An efficient method for incorporating computational fluid dynamics into sonic boom prediction,” Proceedings of the AIAA 9th Applied Aerodynamics Conference, AIAA Paper No. 91-3275, 1991. [9] S. C. Crow, ‘‘Distortion of sonic bangs by atmospheric turbulence,’’ J.Fluid Mech. 37, 529–563, 1969.

[10] Locey, L. L. and Sparrow, V. W., “Modeling atmospheric turbulence as a filter for sonic boom propagation,” Noise Control Engineering Journal, Vol. 55, No. 6, 2007, pp. 495–503.

[11] P.V.Yuldashev, et al., “Nonlinear propagation of spark-generated N-waves in atmosphere: theoretical and experimental assessment of the shock front structure”, Acoustics’08 Paris, pp.1813-1818, 2008.

ソニックブーム推算に関する委託研究

平成

21～23

年（2009～2011

年）度

研究成果報告書

高精度近傍場圧力波形推算技術に関する研究

中橋和博

，佐々木大輔

，大木裕介

第１章 研究概要と目的

将来の超音速旅客機開発には様々な技術障壁が存在しており，その一つとして地上でのソニックブームによる騒音 問題が挙げられる．ソニックブームは超音速で飛行する機体周りの衝撃波が地上へ伝播する過程で統合され，地上到

達時に急激な圧力変化を伴う波形になる現象である（図1）．ソニックブームにより地上の人や物がさまざまな被害を

受けるため，機体の低ブーム設計が非常に重要となる．超音速機設計の評価の指標として地上におけるソニックブー

ムの圧力波形（以下，地上ブーム波形）を評価する必要があり，地上ブーム波形推算精度の向上が必要不可欠である．

数値計算による地上ブーム波形推算手法の一つとして，機体の近傍場での圧力波形を元に地上波形の推算を行う方法 がある．この手法では，機体近傍の圧力波形を元に伝播解析を行い地上での圧力波形を求めるため，近傍場の圧力波 形を精度よく解析する必要がある．

本受託研究では，数値流体力学(CFD：Computational Fluid Dynamics )を用いた近傍場での圧力波形推算精度の向

上を目的とし，解適合格子細分化法の開発を行う．解適合格子細分化を用いることで，適切な指標に基づき物理量変

化の大きい領域の細分化を行い，局所的に解の精度を高めることが可能である．衝撃波波面の同定に用いられるShock

Functionを指標として用いることにより，衝撃波を精度良くとらえることが可能であるが，適切なフィルタを用いる

ことによる数値誤差の影響低減も必要である．また，格子細分化を繰り返すことによる格子数の大幅な増大は解析速 度を遅くしてしまうため格子数の不必要な増大は避けるべきである．

本研究では，非構造格子CFD（TASコード，[1]）を用いて精度よく解適合格子細分化を行うために，以下の点につ

いての研究を行う．その際，計算時間の短縮を目指した格子疎化法の開発及び衝撃波の特性に沿って細分化を行う非 等方細分化手法の開発を行う．

・ 近傍場圧力波形を精度良く捕らえるための解適合格子細分化手法の開発

・ 衝撃波に適した細分化指標及びフィルタの検討

・ 格子疎化法の開発

・ 非等方細分化法の開発

・ 軸対称物体及び翼付き単純形状であるDeltaWing形状による妥当性の検証

・ 静粛超音速研究機モデル(S3TD)形状による解析

*1 東北大学 大学院工学研究科 航空宇宙工学専攻 教授

*2 東北大学 大学院工学研究科 航空宇宙工学専攻 助教

ソニックブーム推算のための

高精度近傍場圧力波形推算技術に関する研究

中橋和博

，佐々木大輔

，大木裕介

第１章 研究概要と目的

将来の超音速旅客機開発には様々な技術障壁が存在しており，その一つとして地上でのソニックブームによる騒音 問題が挙げられる．ソニックブームは超音速で飛行する機体周りの衝撃波が地上へ伝播する過程で統合され，地上到

達時に急激な圧力変化を伴う波形になる現象である（図 ）．ソニックブームにより地上の人や物がさまざまな被害を

受けるため，機体の低ブーム設計が非常に重要となる．超音速機設計の評価の指標として地上におけるソニックブー ムの圧力波形（以下，地上ブーム波形）を評価する必要があり，地上ブーム波形推算精度の向上が必要不可欠である． 数値計算による地上ブーム波形推算手法の一つとして，機体の近傍場での圧力波形を元に地上波形の推算を行う方法 がある．この手法では，機体近傍の圧力波形を元に伝播解析を行い地上での圧力波形を求めるため，近傍場の圧力波 形を精度よく解析する必要がある．

本受託研究では，数値流体力学 ： を用いた近傍場での圧力波形推算精度の向

上を目的とし，解適合格子細分化法の開発を行う．解適合格子細分化を用いることで，適切な指標に基づき物理量変 化の大きい領域の細分化を行い，局所的に解の精度を高めることが可能である．衝撃波波面の同定に用いられる

を指標として用いることにより，衝撃波を精度良くとらえることが可能であるが，適切なフィルタを用いる ことによる数値誤差の影響低減も必要である．また，格子細分化を繰り返すことによる格子数の大幅な増大は解析速 度を遅くしてしまうため格子数の不必要な増大は避けるべきである．

本研究では，非構造格子 （ コード， ）を用いて精度よく解適合格子細分化を行うために，以下の点につ

いての研究を行う．その際，計算時間の短縮を目指した格子疎化法の開発及び衝撃波の特性に沿って細分化を行う非 等方細分化手法の開発を行う．

・ 近傍場圧力波形を精度良く捕らえるための解適合格子細分化手法の開発 ・ 衝撃波に適した細分化指標及びフィルタの検討

・ 格子疎化法の開発 ・ 非等方細分化法の開発

・ 軸対称物体及び翼付き単純形状である 形状による妥当性の検証

・ 静粛超音速研究機モデル 形状による解析

東北大学 大学院工学研究科 航空宇宙工学専攻 教授

東北大学 大学院工学研究科 航空宇宙工学専攻 助教

東北大学 大学院工学研究科 航空宇宙工学専攻 修士 年

図1 超音速飛行で生じるソニックブーム

第２章 解適合格子細分化手法

数値解析において，解析コストと精度の両立を図るには流れ場の変化に応じた格子点の配置が求められる．解適合 格子法では，格子及び流れ場から細分化(疎化)が必要な領域を抽出し，格子操作を行う．特に，近傍場においては圧力 波形取得位置まで伝播空間における格子の改善が必要なため，モデル長に対し数倍離れた位置までの細分化しなけれ ばならない．本手法では，格子の改善と，格子操作後の流れ場の更新の反復処理を行う．解析コストを低減するには 効率良く格子の改善を行い，反復数を減らす工夫が必要である．格子解像度が不足する格子では，物理量の変動の減 衰により細分化判定が困難となるため，1回の細分化処理で改善できる領域は制限される．解適合格子法を効率良く行 うには，細分化指標が減衰の影響を受けにくいものを選択し，判定可能領域をできる限り確保することが重要である． 本研究では，2種の減衰の影響を受けにくい指標を組合せる事により細分化判定可能領域の改善を行う．

また，超音速域での圧力変動は衝撃波や膨張波が支配的で，これらの現象は特定方向への圧力勾配が大きく(異方性 が強く)，方向により必要とする格子解像度が異なる．変動方向を考慮した格子点の追加(非等方細分化)は，考慮しな い場合(等方細分化)に比べ少ない格子点の追加で効果的に精度の改善を行うことが期待できるため，非等方細分化を取 り入れる．

2.1 誤差評価法

評価指標にはAit-Ali-YahiaとHabashiによって提案されたEdge-based Error Estimation [2] (以下，Edge Error) を用いる．ある連続した関数を1次元上で離散化するとき格子点間の値は補間によって表され，理論値gと補間値gh の間に生じる誤差は

( )

x

g

( )

x

g

x

e

(

)

=

−

として表せる(図2.1)．ここで�̅は離散化後の要素内における局所座標系で[0,hE]の値を取る．ここで2次の補間値と，

1次の補間値の間に生じる差を考えると，

2 2 2

2

2

)

(

dx

g

d

h

x

x

x

e

E









₋

=

近傍場

中間場

と表せ，最大値は E h E E h E

dx

g

d

h

dx

g

d

h

e

8

∝

=

と表せる．離散化した格子の評価は2点間の距離と2階微分の値から式(2.4)のように行える．なお，本研究では圧力

波形の推算精度の改善を目的とするため，圧力の2階微分値を用いる．

E h E E

dx

g

d

h

e

=

この評価方法で3次元格子の各辺で算出する場合には，ベクトルを用いた拡張を行う．単位ベクトルV方向における

ghの2階微分値は

HV

V

V

g

=

∂

∂

として与えられる．ここでHはghのHessian行列であり以下のように表記される．





































∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

=

z

g

y

z

g

x

z

g

z

y

g

y

g

x

y

g

z

x

g

y

x

g

x

g

H

このとき辺における最大誤差は辺の単位接線ベクトルτと辺の長さlを用いて

E h E h

d

g

d

l

d

g

d

l

e

8

τ

∝

τ

=

と表せ，指標値は

E h edge

d

g

d

l

e

τ

=

と与えられる．各辺の値はHessian行列H(x)を用いて座標変換後の辺の長さと考えられる．Hを変換行列として扱う

ために，対称行列であるHを対角化し，その固有値の絶対値を得ることで変換行列Mを得る．

T

R

R

H

=

Λ

T

R

R

M

=

Λ





















=

Λ





















=

0

0

0

0

0

0

,

λ

λ

λ

e

e

e

R







ここでΛは行列Hの固有値λ1,λ2,λ3からなる対角行列で，Rは�⃗1,�⃗2,�⃗3(行ベクトル)からなる直交行列である．このと

と表せ，最大値は

∝

=

と表せる．離散化した格子の評価は 点間の距離と 階微分の値から式 のように行える．なお，本研究では圧力

波形の推算精度の改善を目的とするため，圧力の 階微分値を用いる．

=

この評価方法で 次元格子の各辺で算出する場合には，ベクトルを用いた拡張を行う．単位ベクトル 方向における

の 階微分値は

=

∂

∂

として与えられる．ここで は の 行列であり以下のように表記される．





































∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

=

このとき辺における最大誤差は辺の単位接線ベクトルτと辺の長さ を用いて

τ

τ

∝

=

と表せ，指標値は

τ

=

と与えられる．各辺の値は 行列 を用いて座標変換後の辺の長さと考えられる． を変換行列として扱う

ために，対称行列である を対角化し，その固有値の絶対値を得ることで変換行列 を得る．

Λ

=

Λ

=





















=

Λ





















=

λ

λ

λ







ここでΛは行列 の固有値λ1 λ2 λ3からなる対角行列で， は�⃗1 �⃗2 �⃗3 行ベクトルからなる直交行列である．このと

き変換前の座標系における ベクトルの座標変換後における長さは

MV

V

V

=

として表せる．行列Mは空間座標の関数であるから2点間における指標値は

∫

−

−

=

0

(

x

x

)

M

(

t

)(

x

x

)

dt

E

(2.13)

となる．

図2.1 離散化による近似解(線形補間)

2.2 Hessian行列を用いた座標変換

行列Mを変換行列とするとき，行列Mにより表現される二次曲面�= [�1,�2,�3]はベクトル座標系において

1

)

(

P

X

=

XM

として表現でき，図(a)のような楕円体を描く．この曲面は点Pを起点とする線分の最大誤差が1となる位置の等値面

を示している．ここで点Pから曲面までの距離をℎ(�,�⃗)として単位ベクトル�⃗を用いて表すと

1

)

(

P

X

=

XM

1

)

,

(

)

(

)

,

(

P

e

e

M

P

h

P

e

e

=

h









であるから

T

e

P

M

e

e

P

h



_

_

)

(

1

)

,

(

=

と変換行列Mと方向ベクトル�⃗によって表せる．M(P)は実対称行列であるとき以下のように分解できる．

(distortion) (rotation)

ここで，�(�)は点Pの座標における関数で









































=

0

0

0

0

0

0

)

(

λ

λ

λ

e

e

e

P

Q







で表される．回転成分と変形成分から構成され，この行列による変換後の座標ベクトル�′は

)

(

P

XQ

X

′

=

1

)

(

)

(

)

(

))

(

)(

(

=

=

=

=

′

′

X

P

XT

X

P

Q

P

XQ

P

XQ

P

XQ

X

X

と与えられる．楕円体表面はQ(P)により半径が1の球表面へと変換される(図2.2 (b))．これを格子が配置している座 標系に適用すると，流れ場の変動に応じて座標系が変換される．四面体に適用すると，T(P)によって定まる最適な辺 の長さを持つ四面体は辺の長さ1の正四面体に変換される図2.2 (c))．変換後の辺の長さは式(2.13)により算出される 値に等しく

∫

∫

=

−

−

=

=

′

dx

e

x

M

e

dt

x

x

t

M

x

x

M

AB

L





)

(

)

)(

(

)

(

)

(

0 (2.22)

のように，変換行列Mと方向ベクトル�⃗で表すことができる．T(P)は空間座標の関数であるため，辺の2点間で異な る計量場を持つ(図2.3)．T(P)が2点間の間で線型に変化すると仮定したとき式(2.22)は













−

−

=

′

otherwise

)

,

(

)

,

(

ln

)

,

(

)

,

(

)

,

(

)

,

(

if

)

,

(

)

(

e

B

h

e

A

h

e

A

h

e

A

h

L

e

A

h

e

A

h

e

A

h

L

AB

L















ここで，� � は点 の座標における関数で









































=

λ

λ

λ







で表される．回転成分と変形成分から構成され，この行列による変換後の座標ベクトル�′は

=

′

=

=

=

=

′

′

と与えられる．楕円体表面は により半径が の球表面へと変換される図 ．これを格子が配置している座

標系に適用すると，流れ場の変動に応じて座標系が変換される．四面体に適用すると， によって定まる最適な辺

の長さを持つ四面体は辺の長さ の正四面体に変換される図 ．変換後の辺の長さは式 により算出される

値に等しく

∫

∫

=

−

−

=

=

′





のように，変換行列 と方向ベクトル�⃗で表すことができる． は空間座標の関数であるため，辺の 点間で異な

る計量場を持つ図 ． が 点間の間で線型に変化すると仮定したとき式 は













−

−

=

′















と変換し算出することができる ．

(a)M(P)により表現される楕円体 (b)楕円体に対する座標変換

(c)四面体への適用

図2.2 変換行列による座標変換

図2.3 辺AB上に展開される計量場

P

x

x

x

)

,

(

P

e

h





1 ) (_{P X}T =

XM

P

x

x

x

1

)

,

(

=

′

P

e

h



)

(

P

Q

e



1

=

′

′

X

X

)

(

P

XQ

X

′

=