第 25 回数値流体力学シンポジウム
C01-2

圧縮性 Building-Cube Method を用いた翼周りの非定常非粘性流れ解析
Unsteady and inviscid flow analysis around the wing
using the compressible Building-Cube Method○ 西村康孝,東北大院, 宮城県仙台市青葉区荒巻字青葉 6-6-01, yasutaka@ad.mech.tohoku.ac.jp
佐々木大輔, 東北大工, 宮城県仙台市青葉区荒巻字青葉 6-6-01, sasaki@ad.mech.tohoku.ac.jp
中橋和博,
東北大工, 宮城県仙台市青葉区荒巻字青葉 6-6-01, naka@ad.mech.tohoku.ac.jp
Yasutaka NISHIMURA, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Sendai, 980-8759, Japan
Daisuke SASAKI,
Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Sendai, 980-8759, Japan
Kazuhiro NAKAHASHI, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Sendai, 980-8759, JapanThe goal of this study is to establish a practical and efficient flutter simulation tool based on Immersed Boundary
Method-based Cartesian mesh solver , which achieves low computational costs and high accuracy to simulate the wing
flutter phenomenon. The unsteady Building Cube Method(BCM) is firstly validated by the pitching NACA0012 airfoil.
The tool is then extended for flutter simulation by integrating structural analysis. AGARD 445.6 weakened model isused for the validation.
1. 緒言
近年の航空機設計においては, 数値流体力学(CFD)技術の進歩により, 巡航時の性能推算のみならず, フラッタのような動的不安定現象解析も行われ, 遷音速域の衝撃波の影響による遷音速デ
ィップのような線形解析等では推算が難しかった現象も捉える事が出来るようになってきている(1). しかしながら, 従来の構造格
子法や非構造格子法でのフラッタ解析では物体の移動に伴う以下の処理が多大な計算負荷となっており, 実用性を損なっている.

・複雑形状に対しても自動で高速に格子生成可能
・Cube 細分化が容易なので流れの局所特性に格子を適合可能
・全 Cube は同数の Cell を含むので並列計算に容易に対応可能
・空間高次精度を容易に確保可能
・大規模データ後処理が容易
・アルゴリズム, 計算コード構造が簡略化可能等の長所を持つ.

フラッタ解析における処理①表面格子変形②表面格子変形速度に基づく壁面境界条件の更新③空間格子変形④空間格子の面積, 体積の再計算⑤空間格子の変形速度に伴う付加流束項計算
フラッタ解析への直交格子法の適用例は構造格子法,非構造格子法と比較すると尐ないが, この手法は物体壁面の移動を, その移動速度を加味した境界条件の付加だけで扱えるため, 上記の③～⑤のステップを省略でき, 大幅な解析時間の短縮が期待できる.
以上より, 中橋により提案されたブロック構造を用いた直交格子法に基づく Building-Cube Method(BCM) (2)のフラッタ解析への適用を検討する.
直交格子法で移動物体を取り扱う際には物体近傍での流体 /壁フラグの効率的な設定及び壁面境界条件の取り扱いが重要となる. 本研究ではBCMについてまず, 壁面のフラグ切替手法の構築,
Immersed Boundary Method の導入, 物体移動速度を加味した壁面
境界条件の付加を施した非粘性非定常流解析ツールを開発した.
続いて, ピッチング運動を行うNACA0012 翼型を用いて移動物体計算検証を行い, その後 3 次元に拡張し, モーダル解析法による構造計算ルーチンを組み込み, フラッタ計算の代表的な検証例である AGARD weakened 445.6 model に対してフラッタ検証計算を行った.

Fig.1 Cubes and cells around the circle (2D)
本研究で開発したソルバーの支配方程式には以下に示す, 圧縮
性 Euler 方程式を用いた. ここで Q は保存量ベクトル, F は非粘性
流束ベクトルである.

Q F j0
t x j(2.1)

この式をセル中心有限体積法で離散化すると以下の式(2.2)が得ら
れる.


Q i
1 
   S ij f (Q ij , Q ij , n ij ) (2.2)
tVi  j (i )

Vi は cell の体積, S ij は cell 境界の面積, f は非粘性流束ベクトル

2. 流体計算手法
2-1. Building-Cube Method (BCM)
BCM は Fig.1 のように, 流れ場を”Cube”と呼ばれる正方形(2 次元)または立方体(3 次元)領域に分割し, 各 Cube 内に”Cell”と呼ばれる等間隔直交格子を生成して流体計算を行う. この手法は,でQij とQij は cell 境界を挟んだ両側の物理量を表す. 各段階における計算法には Table1 の手法を用いる.
1

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C01-2
Fig.2 に本研究で使用する BCM の計算結果と, 非構造格子法(TAS-code)による非粘性計算結果(機体表面 Cp 分布)を示す. 両者の比較よりBCMとTASは非常によく一致していることが分かる.

Table1 BCM compressible Euler solver
支配方程式
圧縮性 Euler 方程式
離散化手法
cell 中心有限体積法
HLLEW (3)非粘性流束評価
高次精度化手法
3 次精度 MUSCL
Differentiable limiter(4)
制限関数
時間積分法
LU-SGS 陰解法(5)2-2. 壁面近傍流体 / 壁フラグ切替手法
直交格子による移動物体解析においては物体移動に伴い物体近傍における流体 / 壁フラグ切替, 及び n-1 時刻においては壁であったがn時刻では流体となるようなcell (Fresh Cell, FC)の取り扱いが重要となる.
Fig.3 のように破線部で表される n-1 時刻における物体表面が移動して, n 時刻には実線部に移動した場合を考える. まず, n-1 時刻
における物体と n 時刻における物体表面を用い bounding box を作成する. 続いて bounding box に含まれるすべての cell に対して各軸方向に Tomas Möller の交差判定法(6)を応用して物体を検索する.
そして, 最終的に物体の法線情報を用いて cell の内外判定を行い,
Fig.4 のように流体 / 壁フラグ切替が終了する. なお図中の Ghost
cell は壁面境界条件構築に必要な cell で, 各軸方向において流体
cell と接している物体 cell と定義される.

LUSGS陰解法は式(2.3)のように右辺項に3点後退差分と内部反復を用いて時間 2 次精度を維持する. 本研究では内部反復数を 5回とした.

RHS j ,k ,l 

3Q mj ,k,1l  4Q nj ,k ,l  Q nj ,k1,l2t  S ij f (Q ij , Q ij , n ij ) (2.3)
j (i )

ここで m は内部反復数を表す. m-1 時刻において流体, n 時刻や
n-1 時刻に壁のような場合には局所的に時間 1 次精度とすることで対処した.
本計算法で DLR-F6 の定常計算を実施した. 以下の Table2 に計算条件を示す.
Table2 calculation condition
Mach number
0.7
Angle of attack [deg]0.5
Number of cube
4,662
cells in a cube
16x16x16
Total number of cells
19,095,552
Minimum cell size
0.00092 / length of fuselage
dt
1.0

Fig.3 fluid / wall flag change (1)(a)BCM

Fig.4 fluid / wall flag change (2)続いて, FC の取り扱いを説明する. FC へは Image Point (IP)を用いる Immersed Boundary Method(7)を応用して適切な値を内挿する.
注目する FC から 2 次元では 1.5dx, 3 次元では 1.75dx (dx は最小格
子幅, 1.5 や 1.75 は cell 対角長以上を意図する.)の物体法線方向距
離に IP を作成する. 続いてこの IP の周囲 9 点(2 次元の場合, 3 次
元では 27 点)に着目し, 流体 cell のみから距離で重み付け(8)を行って物理量を内挿する. 最後に IP の値を用いて物体壁面で slip 条件を満たすように FC の値を定める. 速度は式(2.4)のように定める.
圧力は Neumann 条件(式(2.5)), 密度は式(2.6)のように一様流総エ
ンタルピより求める.

(b)TAS
Fig.2 Surface Cp distribution of the DLR-F62

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C01-2

Fig.5 Interpolation to the fresh cell




 
VFC  VIP  (dip /( dip  db))(VIP・n)n (2.4)

Fig.7 Schematic of interpolation to GC from IP

pFC  pIP (2.5)
2FC  pFC /(  )( H  0.5VFC) (2.6)2-3. 壁面境界条件
2-2節のFCへの物理量の内挿で用いたように壁面境界条件を構築する手法は Ghost Cell(GC)と Image Point(IP)を使用する手法をベースとする Immersed Boundary Method を適用する. まず, Fig. 6 に示すように流体 cell と接しているcell を Ghost Cell と名付け, 物体
から法線方向に最小cell サイズの1.5 倍の長さにある地点にImage
Point と呼ばれる点を定義する. このときに, GC を含むIP の周囲9点(3 次元では 27 点)の cell に注目する.この cell のうち流体の cellのみから距離に基づいて重み付けを行い, IP へ周囲の流体 cell より値を内挿する.

Fig.8 Interpolation to the boundary velocity




 
VGC  VIP  (1  (dIP / dGC ))((VIP  Vboundary・) n)n (2.7)

P    IP (aboundary n )
n
(2.8)
T
 0
n
(2.9)

3. 移動物体検証計算 (NACA0012 pitching airfoil)
3-1. 計算対象及び計算条件
2-2 節, 2-3 節で述べた壁面における流体 / 壁フラグ切替手法を検証すべく, ピッチング運動をする NACA0012 airfoil を取り上げる. この流れ場に対する1982年に行われたLandonらによる実験(9)及び, 非構造格子法(TAS-code)による結果(10)と比較することで,本研究で開発した BCM 移動物体計算コードの精度を検証する.
計算条件は Landon の実験に従った. 以下の Table 3 に計算条件
及び解析格子についてまとめた. また Fig.9 に使用した格子を示す. 解析に使用した計算機は本研究室保有のワークステーションである(CPU : Xeon 2.93[GHz] x 8, Memory : 4[GB] x 16).

Fig.6 Schematic of GC and IP
続いて IP から GC への内挿であるが Fig.7 に示すように速度に関しては物体の移動速度を加味して slip 条件を課した(式(2.7)). このときの物体の移動速度は Fig.8 のように各 node の速度から内挿して求める. また, 式(2.8)のように壁面圧力勾配には物体の加速度を加味し, 断熱壁条件も式(2.9)のように課す.
翼の後縁等の薄い部分で, wall cell が 1 層のみとなるような場合でも計算できるように, 多価の GC を定義する. 例えば翼型後縁近くの点では GC 上に上面と下面それぞれの値を持たせる. そして, 最終的な壁面の境界条件はその値をもとに壁面第1層目のcell境界で空間 1 次精度の流束として与える.

Table3 calculation condition
Mach number
0.755
Angle of attack [deg]0.0Number of cube
226
cells in a cube
32x32
Total number of cells
231,424
Minimum cell size
0.0018/chorddt
0.01
非定常計算は定常計算の結果を初期値として用いる. 非定常計
算において NACA0012 airfoil が行うピッチング運動は以下の式で表される. ピッチング運動の中心は 1/4 コード長である.

 (t )     am sin(t ) (3.1)
3Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C01-2(2) 迎え角-2.0deg]で上向きにピッチング運動

ここで 0 は平均迎え角,  m は振幅,  は角振動数である. それ
ぞれの値は 0  0.016[deg] ,  m  2.51[deg] であり, 無次元振動数 k  0.0814 とした. 無次元振動数は角振動数 , 一様流マッ
ハ数 M  , コード長 c を用い, k  c / M  と定義される.

Fig.12 Cp contour (angle of attack=-2.0[deg])

Fig.9 Cubes and cells around the NACA0012 airfoil
3-2. 計算結果及び考察
非定常性が顕著になる 3 周期目における, (1)迎え角 0.52 度で下向きにピッチング運動している場合, (2) 迎え角-2.0 度で上向きにピッチング運動している場合の 2 つについて, Cp コンターを示し(Fig.10, Fig.12), Cp分布を非構造格子(TAS-code)を用いたEuler計算
結果及び実験値と比較する.
(1) 迎え角 0.52[deg]で下向きにピッチング運動

Fig.13 Cp distribution (angle of attack = -2.0[deg])
ピッチング運動によりにより見かけの迎え角が変化するので流れ場の様子は非定常性が強くなる. Fig.11 及び Fig.13 の Cp 分布を見ると, まず BCM と TAS を比較すると 2 者はよく一致している.
しかしながら BCM は計算値が振動している. これは壁面の流体/ 壁フラグ切替とそれに伴う FC への物理量内挿により不安定になっているものと考えられる.
また BCM と実験値を比較するとよく一致しているが, BCM の方がより強く衝撃波を捉えていることが分かる. これは BCM が非粘性計算で行われているために, 衝撃波と境界層の干渉を考慮に入れていないためであると考えられる. 実際の流れでは翼表面に境界層が発達し, 衝撃波と干渉を起こし, 急激な圧力勾配が緩和される.
続いて迎え角と揚力係数cl の履歴線図をBCM とTAS それぞれ
Fig.14 と Fig.15 に示す.ピッチング運動を開始してから 3 周期目までを示した. 対称翼型なので迎え角 0.0[deg], 揚力係数 0.0 を始点にピッチング運動に伴うヒステリシスループを描く. BCM と TASを比較すると両者はよく一致している. しかしながら, いずれの計算も実験値よりも大きなヒステリシスループを描いており, 非定常性をより強く捉え, 全体的に揚力係数を小さく評価している.
この原因は非粘性計算にあると考えられ, 物体の滑り速度の許容と境界層と衝撃波の干渉を考慮していないことによる, より強い
衝撃波の生成が原因であると考えられる.

Fig.10 Cp contour (angle of attack=0.52[deg])

Fig.11 Cp distribution (angle of attack = 0.52[deg])
4

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C01-2

i2 は i 番目振動モードの固有角振動数を表す.
最終的には式(4.3)は 3 点後退差分を用いて式(4.4)のように時間
2 次精度で離散化される.
1

m 1
n
n 1

2

 2   3S i  4 S i  S i 
m 1
m 1
S im   I  tAi t    Ai S i  f i 32t
 3  



(4.4)
5. フラッタ検証計算 (AGARD weakened 445.6 model)
5-1. 計算対象
AGARD weakened 445.6 model(11) (Fig.16)を対象とするフラッタ
解析を実施した. この翼はアスペクト比が 1.6525, テーパー比が0.6576, 1/4 コード長における後退角が 45.0[deg], 翼断面に使用している翼型は NACA65A004 である. 解析で考慮したモードは計6モードであり, その固有振動数を Table4 に示す. また, それぞれの振動モード形を Fig.17 に示す. この 6 つのモード形についてス
プライン関数を用いて表面格子にマッピングした.

Fig.14 Cl hysteresis loop (BCM)

Fig.15 Cl hysteresis loop (TAS)
4. 構造計算手法
前節において, 実験値, 非構造格子法との比較によって直交格
子法に基づく移動問題解析の妥当性を確認した. そこで, 解析ツールを 3 次元に拡張して構造との連成解析機能を導入し, フラッ
タ解析を実施する. 構造の弾性振動方程式は次式で表される.

Fig.16 AGARD weakened 445.6 model
Table4 List of vibration modes for AGARD weakened 445.6 model
Mode No.
Mode
Frequency [Hz]1
Wing 1st Bending
9.6
2
Wing 1st Torsion
38.1
3
Wing 2nd Bending
50.7
4
Wing 2nd Torsion
98.5
5
Wing 3rd Bending
118.1
6
Wing In-Plane Bending
140.2[ M ]{d}  [ K ]{d}  {F } (4.1)
ここで [M ] は質量マトリクス, [K ] は剛性マトリクス, {d }は変位, {F } は空気力を表す. 式(4.1)に Rayleigh-Ritz 法を適用して運動方程式を構築し, 変位{d } に対して複数の固有振動モードの重ね合わせで振動を表現するモーダル解析法を用いた. モーダ
ル解析法により, 変位は式(4.2)のように表される.

{d}  []{q} (4.2)[] は固有振動モードを表す固有ベクトル(モーダルベクトル),{q} は一般化変位である. 次に, 固有ベクトルの直交性を利用して式(4.1)を変形すると式(4.3)が求まる.

S  Ai S i  f i (4.3)
ここで S i  q i

T0q i  Ai   2
 i

 1
0 

0


fi   Ti F M i 

i2  []Ti [ K ][]i M i  []Ti [ M ][]i
5Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C01-2
MP で並列化した.
Fig.19~Fig.21 に等価対気速度(Vf[m/s]EAS)で整理した, 安定, 不安定, 調和振動それぞれのモード変位時歴(mode1~mode4)を示す.
(1)M=0.678(a) mode 1(b) mode 2(c) mode 3(d) mode 4(e) mode 5(f) mode 6(a)Stable (Vf=104.9 [m/s EAS])

Fig.17 AGARD weakened 445.6 model with mapped mode shapes
5-2. 計算条件
定常計算の結果を初期値としてフラッタ計算を行った. 使用した格子を Fig.18 に示し, 計算条件等を Table5 に示す.

(b)Low. dump (Vf=105.7 [m/s EAS])

Fig.18 cubes and cells around the AGARD weakened 445.6 model
Table5 calculation condition
Mach number
0.678 / 0.96 / 1.141
Angle of attack [deg]0.0Number of cube
3,028
cells in a cube
8x8x8
Total number of cells
1,550,336
Minimum cell size
0.0073/MAC
dt0.06(c)Unstable (Vf=106.3 [m/s EAS])
5-3. 計算結果及び考察
フラッタ解析では, 亜音速域(M=0.678), 遷音速域(M=0.96), 超音速域(M=1.141)それぞれにおいて動圧を変化させた解析を実施,
各モードの一般化変位を用いて系の安定不安定を判断する. 最終
的に安定限界であるフラッタバウンダリを同定する. 解析には本研究室保有のワークステーション(CPU : Xeon 2.93[GHz] x 8,
Memory : 4[GB] x 16)を使用し, 流体解析部, 構造解析部を Open

Fig.19 modal displacement (M=0.678)

6

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C01-2(2)M=0.96(3)M=1.141(a)Stable (Vf=63.8 [m/s EAS])(a)Stable (Vf= 81.1 [m/s EAS])(b)Low. dump (Vf=65.7 [m/s EAS])(b)Low. dump (Vf=93.0 [m/s EAS])(c)Unstable (Vf=69.1 [m/s EAS])(c)Unstable (Vf=95.2 [m/s EAS])

Fig.20 modal displacement (M=0.96)

Fig.21 modal displacement (M=1.141)

7

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C01-2
結果と比較し妥当な結果を得た.
AGARD weakened 445.6 model のフラッタ解析では遷音速ディップを適切に捉え，フラッタバウンダリの推定として概ね妥当な結果を得たが, 特に亜音速域について精度改善の余地がある事が明らかとなった．一方，計算コストの面では従来手法である非構造格子法に対する優位性が示された．更なる大規模並列計算に容易に対応できる BCM の従来手法に対する優位性，現状でも大型計算機を用いず実用的な効率である点は設計者にとって魅力であり, 高速な格子生成と組み合わせての大幅なスループットの短縮が可能である. 今後は, 解析時間の更なる効率化を図りつつも,フラッタバウンダリ予測精度の向上を達成し, 実用的なフラッタ設計解析ツールを構築したい.

また Fig.22 に本計算より求められたフラッタバウンダリを, 風洞試験結果(12), NASTRAN による結果, 非構造格子法 Euler コード
である TAS-code による結果(1)と比較して示す.

謝辞
本研究は科研費（21226018）の助成を受けたものである. ここに深く感謝の意を表す.
本研究を行う上で, 三菱重工業株式会社 航空宇宙事業本部
研究部 機体システム課 竹中啓三博士, 空力統合課 米本聖氏には多大なご協力を頂いた. ここに深く感謝の意を表す.
Fig.22 flutter boundary
参考文献(1) 熊野孝保, 森野裕行, 大林茂, “航空機実用設計に向けた非構
造 CFD 空弾解析ツールの開発”, 第 21 回数値流体力学シンポジウム, 東京, 2007 年 12 月.
(2) K. Nakahashi, and L. S. Kim, “Building-Cube Method forLarge-Scale, High Resolution Flow Computations,” AIAA Paper
2004-0423, 2004.
(3) S. Obayashi and Y. Wada, "Practical Formulation of a Positively
Conservative Scheme, "AIAA Journal, Vol. 32, pp. 1093-1095,
1994(4) V. Venkatakrishnan, "Preconditioned Conjugate Gradient Methods
for the Compressible Navier - stokes Equations", AIAA journal, Vol.
29, No.7, pp 1092-1099, 1991
(5) S. Yoon, and A. Jameson, “Lower-Upper Symmetric-Gauss Seidel
Method for the Euler and Navier-Stokes Equations”, AIAA Journal,
Vol. 26 No. 9, pp. 1025-1026.
(6) T. Möller, B. Trumbore, “Fast, Minimum Storage Ray-Triangle
Intersection,” J. Graphics Tools 2(1), 1997.
(7) R. Mittal et al.,”A versatile sharp interface immersed boundary
method for incompressible flows with complex boundaries”,
Journal of Computational. Physics. 227(2008)4827-4825(8) R. Franke, Scattered data interpolation: tests of some methods, Math.
Compt. 38(1082) 181-200(9) R. H. Landon, "Compendium of Unsteady Aerodynamic
Measurements", AGARD Report No. 702, Aug. 1982.
(10) 高橋俊, 文殊川一恵, 中橋和博, "重合非構造格子法を用いた
非定常流れの数値計算", 第 19 回数値流体力学シンポジウム,
東京, 2005 年 12 月

Fig.22 より遷音速域(M=0.96)においては, 遷音速ディップを適切に捉えられていることが分かる. しかしながら, 亜音速域
(M=0.678)については TAS-code や風洞試験結果よりも危険側へフラッタバウンダリを予測している. また, 超音速域(M=1.141)においては風洞試験とよく一致しているものの, TAS-code とは大きく結果が異なるが, 本翼は超音速域において衝撃波剥離が顕著である事が風洞試験において報告されており, 本領域で非粘性計算である TAS と BCM の両者の優劣を議論することは適切でないと考える. 概して, 現状のBCMによるフラッタ解析では特に亜音速領域での精度改善の必要があり, 今後の課題である.
計算効率については, 流体解析における直交格子法の効率化も寄与してはいるが, 冒頭で述べた移動変形問題に対する効率的な対処によるところが大きい. 計算時間は普及型のPCにて5000stepの非定常計算が6hであった. 従来のTAS-codeでは大型計算機による大規模並列計算でも同等の解析時間を要する．同一計算機で比較すれば10倍のオーダーでの飛躍的な計算効率の向上と言える．
近未来のトレンドである，更なる大規模並列計算においてもBCMの原理であるCube間のロードバランス均等である点から容易に対応可能であるが，大型計算機を用いずとも普及型のPCで十分設計解析に耐えうる効率である点は設計者にとって大いに魅力であると考える．
BCMではいかなる複雑形状に対しても，CFDや格子生成を専門としないユーザー，設計者がGUIを利用して解析格子を数十秒という効率で自動生成できる. 格子生成の高速化, 流体計算の高速化といったBCMの長所を生かした, 前処理から後処理までのスループットの大幅な短縮化によって，本ツールが現状のフラッタ設計解析ツールに対する一つのブレークスルーとなる事が期待できる.

(11) E. Yates, Carson, Jr., ”AGARD Standard Aeroelastic
Configurations for Dynamic Response. Candidate
Configuration I. - Wing 445.6,” NASA TM-100492, 1987.
(12) E. C. Yates, jr, "AGARD standard aeroelastic configurations
for dynamic response", NASA technical memorandum
100492, August, 1987.

6. 結言
本研究では, 圧縮性Euler Building-Cube Method を移動物体問題
へ拡張し, 翼型周りの非定常計算を行った. また, 構造方程式との連成解析により, 三次元翼のフラッタ解析を行った.
ピッチング運動をするNACA0012 airfoilの解析では,直交格子法に固有である壁面の取り扱いに起因する空気力の振動は見られたものの断面 Cp 分布，揚力係数のヒステリシスループともに,
非構造格子法である TAS-cod e による非粘性解析結果や風洞試験
8Copyright © 2011 by JSFM