模型変形計測データを反映した CFD 表面格子修正法の開発

保江 かな子

^*1

　口石 茂

^*1

　橋本 敦

^*2

　村上 桂一

^*2

加藤 裕之

^*3

　中北 和之

^*3

　渡辺 重哉

^*1

　菱田 学

^*4

Modification of CFD Surface Mesh Based on Model Deformation Measurement

^*

Kanako YASUE

^*1

, Shigeru KUCHI-ISHI

^*1

, Atsushi HASHIMOTO

^*2

, Keiichi MURAKAMI

^*2

Hiroyuki KATO

^*3

, Kazuyuki NAKAKITA

^*3

, Shigeya WATANABE

^*1

and Manabu HISHIDA

^*4

Abstruct

　

Modification approach of a computational fluid dynamics

（

CFD

）

surface mesh are examined using measurement data of wind tunnel model deformation. Model deformation measurements (MDM) were successfully made at Japan aerospace exploration agency (JAXA) using stereo photogrammetry with markers. The model’s deformation was then approximated using a deformation law which is identified by a few parameters with markers attached on the wind tunnel model. In this paper, the approach is used to rearrange the CFD surface mesh of a CFD model to conform its configura- tion with that of the wind tunnel test model. We duplicate the deformed configuration of the wind tunnel model using three deformation laws and examine the resulting CFD models. A high-fidelity RANS simulation of the DLR-F6 FX2B wind tunnel model is then performed considering the static deformation of the model and the effct of model deformation on the aerodynamic characteristics are also explored by performing the CFD simulation for the deformed configuration as well as the original (non-deformed) configuration.

Keywords: static aeroelasticity, wind-tunnel model deformation, transonic flow simulation.

概　要

　風洞試験で計測した模型変形データを用いて，

Computational Fluid Dynamics

（

CFD

）表面格子を修正する 方法を検討する．

JAXA

ではマーカーを使ったステレオ写真法により風試模型の模型変形量を計測しており，

計測したマーカー座標値を使うことで主翼の変形則を同定し，変形後の形状を定義することができる．本報 告では，計測したマーカー座標値を使って変形後の形状を同定し，

CFD

の表面格子が通風時の風試模型形状 と一致するように修正する方法を検討する．ここでは三種類の変形手法を検討する．そして，実際に模型の 変形量を計測した

DLR-F6FX2B

モデルに対して本手法を適用し，変形モデルの検証をおこなう．また，変 形前後の形状に対して

Reynolds-averaged Navier-Stokes

（

RANS

）解析を実施することで，変形を考慮してい ない場合と考慮した場合とで空力特性にどの程度影響を及ぼすかを検討する．

*

平成

24

年

11

月

26

日受付 （

Received 26 November, 2012

）

*1

研究開発本部　流体グループ

（

Fluid Dynamics Group, Aerospace Research and Development and Directorate

）

*2

研究開発本部　数値解析グループ

（

Numerical Analysis Group, Aerospace Research and Development and Directorate

）

*3

研究開発本部　風洞技術開発センター

（

Wind Tunnel Technnology Centor, Aerospace Research and Development and Directorate

）

*4

（株）菱友システムズ

（

Ryoyu Systems Co. Ltd.

）

1　緒言

　実機と風洞試験との間には，レイノルズ数効果や壁

/

支持装置の影響，空力弾性による模型変形効果など，

様々な違いが存在する

¹⁾

．これらの違いの中でも模型変形効果は実機空力特性を正確に予測するための重要 な課題の一つであり，レイノルズ数効果や壁・支持干渉などと同様，風洞試験で得られたデータに対して何 らかの補正が必要である．近年では，

National Transonic Facility

（

NTF

）

²⁾

や

European Transonic Wind-tunnel

（

ETW

）

³⁾

のような極低温加圧型風洞の開発により， 気流温度を変化させることで

q/E

一定のまま（

q:

動圧，

E:

ヤング率）

レイノルズ数を変化させることで模型変形の影響を考慮することがが可能となったが，実機と同程度の高い レイノルズ数を実現するためにはやはり気流圧力を上げる必要があるため，このような最新の風洞であって も模型変形の影響は少なからず補正する必要がある．

　通常，風試模型は剛体であり変形しないと仮定している．しかし，模型変形が生じることで，通風時の風 試模型形状が設計した

1G

形状と異なってしまうため，得られる空力特性に違いが生じてしまう．また，模 型が変形することで，

PSP

（

Pressure Sensitive Paint

）計測の精度が低下してしまうといった問題もある

⁴⁾

．風 洞試験で表面圧力分布を取得する

PSP

計測を実施する際，感圧塗料を塗布した模型表面に光を照射し，通風 前と通風後に撮影した画像の発光強度の比から表面圧力分布を算出する．その際，

2

枚の画像の模型位置を 合わせるために，マーカーを使った位置合わせを実施するが，模型変形が生じると模型の位置が動いてしま い正確に表面圧力分布の値を算出できない部分が生じ，

PSP

の誤差要因の一つとなってしまう．その他にも，

実機の設計開発段階で広く使われるようになっている

CFD

（

Computational Fluid Dynamics

）では模型を剛体 と仮定しているため，

CFD

解析結果と風試結果との比較において模型変形の影響により両者の形状に差異が 生じてしまい，結果比較や

CFD

コードの妥当性検証が厳密ではなくなるといった問題もある．

　これまでに，模型変形効果を把握することを目的として，変形量の計測や予測に関する多くの技術が開発 されてきた．風洞試験においては，風試模型の模型変形量を知ることを目的として，ステレオ写真法

⁵⁾

やモ アレ干渉法

⁶⁾

など，様々な模型変形計測法が開発されてきた．また，

PSP

計測画像の模型変形効果を補正す るためにもこれらの技術が適応されている

⁴⁾

．一方，模型変形量を予測するために，

FEM

と

CFD

を用いた 風試模型の静的空力弾性解析が行われており，模型の変形量予測だけでなく模型変形が空力特性に及ぼす影 響も検証されている

^7,8,9)

．しかし，風試模型が平衡形状に達するまでに何度か連成解析を実施する必要があ るため，複雑形状の静的空力弾性解析を実施するには高い計算負荷がかかる．さらに，風試模型には，圧力 配管用の溝などが存在することにより，内部構造が複雑となり，それを模擬するためにはさらなる計算負荷 の増大に繋がる．また，簡単のために内部構造を無視した中実形状での流体構造連成解析が行われているが，

その妥当性を検証するには，風試で取得した模型変形変位計測データや空力特性と十分に比較検討する必要 があると考えられる．そこで我々は，模型変形計測データを使って

CFD

表面格子あるいは

CAD

データを 変形形状に合わせて修正できれば，流体構造連成解析をすることなく短時間で変形が空力特性へ及ぼす影響 を検証できると考えた．また，これにより風試模型を忠実に再現した形状での

CFD

解析が可能となるため，

風試データと

CFD

解析結果との比較においてもより詳細な検証が可能となる．

　そこで本報告では風試模型の変形計測結果を

CFD

表面格子へ反映する方法について検討する．本研究で は，

PSP

画像の模型変形効果を補正するために

Le Sant

により開発された多項式近似による簡易モデル

⁴⁾

（以 下，簡易多項式モデルと呼ぶ）を用いて，変形計測データを元にした

CFD

表面格子の変形を試みる．この 方法では主翼の垂直方向（

z

方向）変位のみを考慮しているため，本報告では前後方向（

x

方向）および横方 向（

y

方向）変位にも対応できるよう，簡易多項式モデルの改良を行い，

x

，

y

方向の変位を考慮する必要が あるかどうかの検証も併せて行う．最後に，本手法で修正した変形後の

CFD

表面格子形状に対して

Reynold

Averaged Navier-Stokes (RANS)

解析を実施し，初期形状の

RANS

解析結果と比較することで，変形が空力特

性へ及ぼす影響を検証する．

2　表面格子修正法

　風洞試験では，通風時に模型だけでなく支持装置のたわみも生じる．そのため，本手法では，始めに胴体

に取り付けられたマーカーの変位量計測結果から支持装置のたわみによる模型位置の移動・回転を算出し，

その後，翼上面に取り付けられたマーカーの変位量計測データを使って空力荷重による翼の変形を算出す

る．

Fig. 1

および

Fig.2

にそれぞれ模型移動補正および模型変形補正の流れを示す．

　本報告では翼の変形を求める際に，

Le Sant

によって開発された方法（簡易多項式モデル）

⁴⁾

を適用する．

この方法は，通風前後のマーカーの位置を計測することでマーカーの変位量を算出し，

PSP

画像の模型変形 効果の補正を行うために開発されたものであり，本研究ではその方法を

CFD

表面格子の変形に適用する．

ここでは主翼の変形をはりのたわみおよびねじりと仮定し，多項式で近似する．そして計測データを使って 多項式の係数を同定することで，主翼の変形形状を定義し，

CFD

表面格子を修正する．しかし，簡易多項 式モデルは

z

方向変位のみを考慮している手法であるため，このモデルを改良して

x

，

y

，

z

全方向の変位を 考慮した手法の開発も行う．

2.1　座標系の定義

　移動・回転パラメータおよび変形パラメータを算出するにあたり，複数の座標系を導入するため，まず始 めに座標系の定義を示す．迎角

0

度の模型初期取り付け位置での模型先端を原点とする座標系を模型座標系 と呼ぶこととし，

x = (x, y, z)

で表わす．また，通風時の支持装置のたわみによる模型の移動・回転を考慮する 必要があるため，

Fig. 3

のように通風時に支持装置がたわみ移動・回転した後の座標系を移動後模型座標系と 呼び

x’= (x’ , y’ , z’ )

で表わす．このとき移動模型座標系の原点は移動後の模型先端とし，機軸方向を

x’

と している．また，主翼の変形を考慮するために，翼根付近を原点とした主軸方向の座標系を ξ

= (

ξ

,

η

,

ζ

)

と定義し，これを変形主軸座標と呼ぶ．なお，変形主軸座標は設定に任意性があるため，変形主軸の設定も 検討が必要である．検討した結果については後述する．

2.2　模型の移動・回転パラメータ算出方法

　支持装置のたわみによる模型の移動および回転を算出する方法を以下に示す．本手法では，胴体を剛体と 仮定し，胴体マーカーを用いて模型の移動・回転パラメータを算出する．初期胴体マーカーを移動・回転さ せ た時の位置と通風時に変形計測した胴体マーカー位置とを比較することで，通風時のたわみによる移動・

回転量を算出することが出来る．ここでは，模型座標系上でのマーカー位置の移動・回転を考える．無風時

Start Correction of motion

Calculate Motion parameters

Modify CFD mesh/STL data

End Re-input

Initial value Input

Body marker data

Input CFD mesh/STL data

Output Motion parameters Output

Error

YES NO

Output CFD mesh/STL data Converge?

Start Correction of

deformation

Calculate Deformation parameters

Modify CFD mesh/STL data

End Input

Order of accuracy

Input Wing marker data

Input CFD mesh/STL data

Output CFD mesh/STL data Output

Error Input Motion parameters

YES NO Converge?

Fig. 1: Flowchart for correction of wind

tunnel model motion. Fig. 2: Flowchart for correction of wind tunnel model deformation.

の胴体マーカー初期位置

x0 = (x0 , y0 , z0 )

を原点周りに

ZYX

系

Euler

角

(

φ

,

θ

,

ψ

)

で回転移動し，

∆ x = (∆ x, ∆y,

∆z)

平行移動した時の模型座標系における移動後の胴体マーカ－位置は次式で表される．

ここで回転行列

T^m

は，

xrot

は回転中心である．ここでは座標原点を回転中心とした．また，

(

φ

,

θ

,

ψ

)

はそれぞれロール角， ピッチ 角，

ヨー角に対応している．

　本手法では最小二乗法により平行移動量

∆ x

および回転角

(

φ

,

θ

,

ψ

)

を算出することを考え，次式で表され る残差二乗和

S ( p)

を最小にするパラメータ

p = ( ∆ x, ∆y, ∆z,

φ

,

θ

,

ψ

)

を

Newton

反復法により求める．

n

はマーカー点数，

ωi

は重みを表す．測定値の誤差の分散から重み

ωi

を算出することで，測定誤差を考慮 することが出来る．フィッティング関数

Fi

は次式で定義する．

ここで

x^ei = (x^ei , y^ei , z^ei )

は通風時に計測されたマーカー座標値である．式

(3)

を最小にするパラメータ

p

を求め るためには，次式を解けばよい．

Newton

反復法では，

を

∆ p

について解くことにより

,

次の近似値

Fig. 3: Coordinates of model for calculating motion parameters and deformation parameters.

䃖 䃝

䃝

䃖 䃕 䃕

X’ Y’

Z’

が得られる．ここで，

である．なお，

(x, y, z)

の一次導関数，二次導関数については，付録

A

に示した．

2.3　模型の変形パラメータ算出方法

　次に，空力荷重による主翼の変形形状を算出する方法を以下に示す．まず模型座標系

x = (x, y, z)

から変形 主軸座標系 ξ

= (

ξ

,

η

,

ζ

)

への座標変換を考える．このときの変換式は以下のように表される．

xw = (xw , yw , zw )

は変形則適用開始位置である．また回転行列

T^t

は

であり，

(

Γ

,

Θ

,

Λ

)

はそれぞれ翼の上半角，取付角，後退角である．

　また，通風時に計測されたマーカ座標値には支持装置のたわみによる変位が含まれているため，前節で求 めた移動・回転パラメータを用いて支持装置のたわみによる変位分を考慮して変形主軸座標系へ変換する必 要がある．そこでまず模型座標系

x = (x, y, z)

から移動後模型座標

x’ = (x’ , y’ , z’ )

へ座標変換し，さらに移動 後模型座標系から変形主軸座標系 ξ

= (

ξ

,

η

,

ζ

)

へ座標変換する．模型座標系から移動後模型座標系への変 換式は，

で表される．ここで

T^mv

は変換行列

T^m

の逆行列（転置行列）である．

2.3.1　z 方向の変位のみ考慮する簡易多項式モデル

　簡易多項式モデルでは，まず始めに模型に適した変形モデルを考えなければならない．ここでは片持ち

はりのたわみとねじりを変形モデルとして適用し，変形主軸座標 ξ

= (

ξ

,

η

,

ζ

)

における ζ 方向変位を多

項式 で近似する．本報告では，曲げを二次，捻りを一次を基本とした

Z21

モデルおよび，曲げを

4

次，捻り を

4

次で近似した

Z44

モデルの

2

種類を検討するが，ここでは曲げを

4

次，捻りを

4

次で近似した

Z44

モデ ルの場合を例に説明する．

　変形モデルを適用した際の変形主軸座標系における変形後の模型の位置

(

ξ

_d ,

η

_d ,

ζ

_d )

は以下のように表され る．

ここで

o2

−

o4, t1

−

t4

は変形パラメータであり，次式で表される残差二乗和が最小となるパラメータ

p = (t1, o2, t2, o3, t3, o4, t4)

を求めることで，主翼の変形を再現する．なお，

Z21

モデルの場合には，

p = (t1, o2)

となる．

n

はマーカー点数，ω

_i

は重みを表す．またフィッティング関数

Fi

は次式で定義する．

ここで ζ

^e_a,i

は変形主軸座標系に変換された通風時のマーカー位置である．移動・回転パラメータの算出方法 と同様に最小二乗法を適用して変形パラメータを同定する．

Newton

法を適用する際必要となる，

Fi

の一次導 関数および二次導関数を以下に示す．

2.3.2 x，y，z 方向の全変位を考慮する改良モデル

　ここでは簡易多項式モデルを改良し，

x, y, z

方向全ての変位を考慮することを考える．簡易多項式モデルを

PSP

画像の模型変形補正に適用した際，結果に大きく影響を及ぼすパラメータは曲げ二次のパラメータ

o2

， 捻り一次のパラメータ

t1

および変形則開始位置

yw

と報告されている

⁴⁾

．そこで，ここでは曲げを二次， 捻りを 一次と仮定する．

z

方向変位は上述した簡易多項式モデルを適用する．一方

x

，

y

方向変位に関しては，

x

方 向あるいは

y

方向の変位を考慮していない場合の近似式と

x

あるいは

y

方向の変位を考慮した場合の近似式と の差 δ

x

，δ

y

をとることで以下のような変形モデルを適用する．

ここでは

p = (t1, o2, dt1, do2)

が変形パラメータであり，残差二乗和を最小にするような各パラメータを _求め

ることで，変形形状を近似する．

また，フィッティング関数

Fi

は次式で定義する．

なお，

Newton

反復法を適用する際に必要となる

(

ξ

_d ,

η

_d ,

ζ

_d )

の一次導関数，二次導関数については，付録

B

に示した．

2.4　CFD 表面格子修正方法

　模型の変形パラメータが求まれば，変形主軸座標系における変形後の位置は式

(17)

あるいは式

(21)

で定義 できるため，

CFD

表面格子についても同様に各格子点にこれらの式を当てはめることで表面格子を変形させ ることが出来る．

CFD

の表面格子修正の際にも，まず

CFD

格子を模型座標系から変形主軸座標系へ変換し，

変形主軸座標系において式

(17)

あるいは式

(21)

を適用し，変形後の格子点座標を求める．その後，変形主軸 座標系から模型座標系へ逆変換し，変形を反映した

CFD

格子を作成する．なお，ここでは主翼のみの変形を 考慮しているため，主翼の表面格子点のみ変形させている．

3　結果および考察 3.1　模型変形計測を反映した表面格子修正

　本報告では，標準模型である

DLR-F6

模型

¹⁰⁾

を用いて本手法を検証する．マーカーを用いたステレオ写真 法により計測された模型の変形量

¹¹⁾

を元に，上述した手法にて表面格子修正を行う．胴体および翼上面に設 置されたマーカー位置を

Fig. 4

に示す．ステレオ写真法は，

2

台の高解像度

CCD

カメラを用いて模型に貼り 付けられたマーカーの三次元位置情報を計測することで，模型の変形量を算出する方法である．解析対象と した主流条件を

Table 1

に示す．前述した通り，結果に大きく影響を及ぼすパラメータは曲げ二次のパラメー タ

o2

，捻り一次のパラメータ

t1

および変形則開始位置

yw

と報告されている．そこで，ここでは曲げを二次，

捻りを一次を基本とし，簡易多項式モデルを用いた

Z21

モデルおよび簡易多項式モデルを改良して

x, y, z

方 向全変位を考慮した

XYZ

モデルを使って検証する．また，確認のために簡易多項式モデルを用いて曲げを

4

次，捻りを

4

次で近似したもの（

Z44

モデル）についても併せて検討する．また，本計算では，式

(3), (18)

お よび

(22)

の重みは ω

_i = 1

としている．

Fig. 4: Schematic illustration of marker positions.

Wing markers Body markers

3.1.1　変形主軸設定方法の検討

　それぞれの変形モデルの比較を行う前に，変形主軸を設定する際に考慮すべきパラメータである変形モ デル適用開始位置

(xw , yw , zw )

が結果におよぼす影響を

Z21

モデルで検証する．また，検証には

Case 5

の 主流条件を用いた．ここでは，翼根前縁位置

(Position-1)

，最も翼根側に近いマーカーの前縁マーカー位置

(Position-2)

および後縁キンク位置前縁

(Position-3)

の

3

つの位置で検討した．また，η 軸は翼前縁の後退角と

平行になるよう設定した．

　それぞれの変形モデル適用開始位置で本手法を適用した際のマーカーのたわみ量およびねじり量を

Fig. 5

および

6

に示す．たわみ量を見ると，後縁キンク位置前縁を変形主軸原点とした場合には，計測データから 大きく離れてしまう．また，ねじり量を見ると，翼根前縁位置に主軸原点を設定した場合には過小評価，後 縁キンク位置前縁に設定した場合には過大評価している．次に，得られた結果の誤差を

Table 2

に示す．こ こでは，次式で表されるように，本手法により推算された変形後のマーカー位置

xm

と計測により取得したマー カー位置

x^expm

との距離の

L2

ノルムを誤差と定義した．

Table 2

より，最翼根側にある前縁マーカーを変形主軸原点とした場合が，最も誤差が少ないことがわかる．

これらの結果に基づいて，以降の計算については変形モデル開始位置を最も翼根にある前縁側のマーカー 位置と設定する．今後， さまざまな模型形状に適用することを考えた場合， 変形モデル開始位置も変形パラメー タとして，他のパラメータと同様に推算する必要があると考えられる．

3.1.2　変形モデルの影響

　次に，各変形モデルの妥当性検証を行う．

Z21

モデル，

XYZ

モデル，

Z44

モデルの

3

つのモデルを適用し て得られた表面格子を初期形状の表面格子と共に

Fig. 7

に示す．ここでは見やすさのために変位量を

10

倍に して表示している．どのモデルにおいても定性的に妥当な変形形状を示している．変形形状を定量的に評価 するために，表面格子修正法により得られたたわみ量およびねじり量のスパン方向分布を，マーカーの計測

Table 1: Freestream conditions.

Case Mach Number Chord Reynolds number Angle of attack

[Deg] Total pressure

[kPa] Total temperature [K]

1 0.75 1.5 × 10⁶ -3 100 322

2 0.75 1.5 × 10⁶ -2 100 322

3 0.75 1.5 × 10⁶ -1 100 322

4 0.75 1.5 × 10⁶ 0 100 322

5 0.75 1.5 × 10⁶ 1 100 322

Table 2: Error between marker position estimated by deformation models and that measured in wind tunnel test.

変形モデル適用開始位置 誤差 ε

(

無次元距離

)

Position-1

翼根前縁位置

7.245 × 10⁻⁴

Position-2

最翼根側の前縁マーカー位置

6.615 × 10⁻⁴

Position-3

キンク前縁位置

1.130 × 10⁻³

結果と共に

Fig. 8

および

Fig.9

に示す．たわみ量に関しては，どのモデルも計測データを良く再現している．

一方，ねじり量を見ると，

Z21

モデルと

XYZ

モデルの間に翼端で最大で

0.07deg

の差が生じている．また，

Z44

モデルを見ると，全てのケースに関して，計測データを最も良く再現しているが，

Case-1

や

Case-3

の翼 端部分を見ると，計測データを忠実に再現しようとするためにオーバーシュートしてしまっている部分があ る．

Z44

モデルは高次にねじり量の分布を再現できるが，計測誤差などが入った場合にも忠実に再現してし まうため，

Z44

モデルを適用する際には，計測誤差を考慮する必要があると考えられる．

　

Table 3

に本手法により得られたマーカー変位量と計測により得られたマーカ変位量との誤差を示す．ここ

では，本手法で得られたマーカーの位置（ねじり量）と計測したマーカー位置（ねじり量）との差の

L2

ノ ルムを誤差とした．

Z21

モデルと

XYZ

モデルを比較すると，たわみ量の誤差，ねじり量の誤差共に

XYZ

モ デルの方が小さく，

x

，

y

，

z

全方向の変位を考慮することで，より計測値に近づき，変形形状をよく再現で きることがわかる．

Z44

モデルを見ると，たわみ量に関しては

Z21

モデルと同程度だが，ねじり量の誤差は 最も小さくなっている．しかし，高次の多項式近似であるが故に計測点位置での誤差が小さくなるように近 似されていると考えられるため，前述したとおり，

Z44

モデルを適用する際には誤差を考慮した手法を適用 する必要がある．

Fig. 5: Obtained wing deflection in deformation

coordinates. Fig. 6: Obtained wing twist change in deformation

coordinates.

Table 3: Error between marker position estimated by deformation models and that measured in wind tunnel test.

Case

たわみ量の誤差 ε

×10⁻⁴ (

無次元距離

)

ねじり量の誤差

×10⁻² [deg]

Z21-model XYZ-model Z44-model Z21-model XYZ-model Z44-model

1 5.00 3.64 4.89 7.70 4.50 3.20

2 4.72 3.93 4.59 9.80 8.00 3.30

3 4.50 3.90 4.28 12.3 1.15 3.00

4 4.54 3.94 4.41 7.70 6.20 3.00

5 6.62 4.84 6.47 9.20 6.70 4.10

3.2　模型変形計測を反映した CFD 解析

　ここでは，前節で算出した変形後の表面格子を用いて，

DLR-F6

模型の

RANS

解析を実施する．初期形 状および

3

種類の変形モデルを適用した形状それぞれに対して，

RANS

解析を行い，実験結果と比較する ことで，変形による流れ場への影響を検証する．計算には，

JAXA

で開発された高速

CFD

ソルバーである

FaSTAR¹²⁾

を用いた．離散化には，セル中心有限体積法を用い，対流流束の計算には

HLLEW

，乱流モデルは

Fig. 7: Obtained wing twist change in deformation coordinates.

Gray : Original configuration Blue : Deformed configuration (Z21-model) Red : Deformed confiuration (XYZ-model) Green : Deformed configuration (Z44-model)

Zoom up view of wing tip

(a) Case-1

(d) Case-4 (e) Case-5

(b) Case-2 (c) Case-3

Fig. 8: Obtained wing bending in deformation coordinates.

Spalart-Allmaras

一方程式モデル，時間発展には

LU-SGS

陰解法を適用した．

　ここで，変形計測を実施した試験ではラフネスなしの自然遷移で計測されたものであるが，本計算では全 面乱流を仮定しているため，

CFD

結果との比較検証のために，別の試験で取得した同条件でラフネスありの 場合の力および圧力計測結果と比較することとした．ラフネスがある場合とない場合とで，揚力係数の値が 変わり，模型の変形量も異なると考えられるが，本報告で用いた迎角

1

度以下の条件下では，ラフネスの有 無に関わらず空力係数結果が概ね一致していたため，ラフネスの有無による変形量の差はないものとした．

　計算格子は初期形状に対しては

JAXA

で開発された自動格子生成ソルバー

HexaGrid¹³⁾

を用いて作成した．

また，変形後の形状の計算格子を初期形状の格子のトポロジーと一致させるために，

Surface Infruence

法を 適 用して初期形状の空間格子点を表面格子の変位量を基に修正することで，変形後の形状に対する空間格 子を

作成した

14, 15, 16)

．

　

Fig. 10

に初期形状の表面格子および模型近傍の対称面の格子を示す．六面体セル，ピラミッドセル，プリ

ズムセル，四面体セルの

4

種類のセル形状から成り，セル数は約

850

万，接点数は

732

万点である．また模 型表面のセル数は約

22

万である．

　まず始めに

Table 1

の

Case 5

の条件で得られた初期形状に対する模型表面圧力係数分布を

Fig. 11

に，また 各スパン断面における初期形状および変形形状の圧力係数分布を実験結果と共に

Fig.12

に示す．内翼側 η

=

0.150

および η

= 0.331

に関しては，変形が小さいため初期形状に対する分布と変形を考慮した場合の 分布

とで大きな差はなく，実験値とも概ね良い一致を示している．一方，η

= 0.514

に関しては，計算の方が実 験と比べて衝撃波位置が多少後退しているが，変形を考慮することによって衝撃波位置は前方に移動し，実 験結果に近づいていることがわかる．また外翼側の η

= 0.638

および η

= 0.847

でも，計算で得られた衝撃 波位置が実験よりも後退しているが，変形を考慮することで，大幅に衝撃波位置が前進している．さらに変 形を考慮していない初期形状での計算結果では，サクションピークが実験より高くなっているが，変形を考 慮することで実験値に近づいていることがわかる．翼端に近づくにつれ衝撃波位置の前方への移動量が大き くなっているが，これは

Fig. 9

にあるように，模型変形の影響により翼端に近づくにつれてねじり量が増加

(a) Case-1

(d) Case-4 (e) Case-5

(b) Case-2 (c) Case-3

Fig. 9: Obtained wing twist change in deformation coordinates.

し翼の局所迎角が減少するためと考えられる．一方，η

= 0.847

断面では，変形モデルの違いにより各モデ ルで衝撃波位置が多少異なっており

Z44

モデルが最も実験に近づいているが，ねじり量を過大評価したこと による影響である可能性も考えられる．

　次に計算により得られた空力係数を

Table 4

および

Fig. 13

に

, Drag polar

カーブを

Fig.14

に示す．実験値と 計算値を比較すると，

CL

および

CD

は実験値よりも大きく，

CM

に関しては実験値よりも小さく見積もられ ている．変形を考慮した場合でも，計算結果と実験値との間に差が残るが，これは壁や支持装置の影響，遷移 位置の影響によるものと考えられる．変形を考慮していない場合に比べて変形を考慮した場合には，

CL

は

Z21

モデルで

0.0135

，

XYZ

モデルで

0.0173

，

Z44

モデルで

0.0163

減少した．これは，

Fig. 9

に示した よう に変形によりねじり下げが生じ，局所的な有効迎角が減少したためと考えられる．この

CL

の変化は

CL

− α カーブの傾きから算出すると迎角

0.1

度程度の減少に相当する．

CD

は

Z21

モデルで

8cnts

，

XYZ

モデルで

11cnts

，

Z44

モデルで

10cnts

減少した（

1cnt=0.0001

） ．これら

CD

の減少は，

Table 4

に示すように， 圧力抵抗の 減少によるものである．一方，

CM

は

Z21

モデルで

0.0034

，

XYZ

モデルで

0.0044

，

Z44

モデルで

0.0057

増加 した．また，

Fig.14

を見ると，変形を考慮しても，

polar

カーブの形状は変化せず，カーブ上 に沿って変化し ていることがわかる．

　次に各スパン断面の局所揚力係数と局所コード長の積で表されるスパン方向の局所揚力係数分布を

Fig.15

に示す．ここでも変形によるねじり下げの影響で実験値に近づいていることが分かる．また変形の大きい翼 端だけでなく，変形量の小さな翼根付近でも揚力の減少が見られるが，これは，翼根付近でも衝撃波位置が 多少前方に移動しているためと考えられる．衝撃波位置が前方に移動した理由としては，翼根付近でもねじ り下げが生じていることと，外翼側に引きずられて翼根付近でも移動したことが考えられるが，現状ではそ の切り分けは困難である．

4　結言

　模型変形による空力特性への影響を簡単に検証するために，風試模型の変形計測データを反映した

CFD

表 面格子修正法の検討を行った．通風時の風試模型と同形状の

CFD

表面格子を作成するために，

Z21

モデ ル，

XYZ

モデル，

Z44

モデルの

3

種類の多項式モデルを用いたところ，たわみ量に関しては，低次多項式近似で あっても

x

，

y

，

z

全方向の変位を考慮した

XYZ

モデルを適用することで，誤差を大幅に低減することができ，

より正確に形状を再現できた．一方，ねじり量においては，高次多項式近似を用いた

Z44

モデ ルが実験デー

Fig. 10: Computational mesh for RANS simulation

over DLR-F6 FX2B wind tunnel model. Fig. 11: Obtained pressure coefficient distributions of Case 5 for original (non-deformed) configuration.

タを最も良く再現できており，このことからねじりについては

4

次まで考慮する必要があると考えられる．

しかし高次多項式を使って分布を近似しているため計測誤差がある場合でも忠実に再現している可能性があ り，

Z44

モデルの有用性を示すためにはさらなる検証が必要である．

　次に，修正した変形形状に対して

RANS

解析を実施し，変形が空力特性へ及ぼす影響を調べた．変形を考 慮することにより，衝撃波位置や空力特性に関して変形を考慮しない

CFD

解析結果よりも風試結果に近づく 結果が得られた．このことから，本手法を適用することで，風洞試験により忠実な

CFD

解析を実施すること が可能となった．一方で，

3

種類の変形モデルの影響を検討したところ，横方向の変位を考慮することによる 空力特性での大きな改善は見られなかった．また高次多項式近似モデルを用いた場合には他のモデルよりも 実験に近づく結果が得られたことからも，ねじり量は高次多項式で近似する必要があると考えられる．しか し高次多項式を用いることでねじり量を過大評価している部分があり，それにより空力特性も過大に変化し た可能性も考えられるため，今後は，計測誤差を考慮した検討を進める．

謝辞

　簡易多項式モデルを適用するにあたり，

JAXA

研究開発本部 風洞技術開発センター栗田充氏に有益な情報 を頂きました．ここに感謝の意を表します．

参考文献

1) Bushnell D. M., “Scaling: Wind tunnel to flight,” Annu. Rev. Fluid Mech., 38(2006), 111-128.

2) “National Transonic Facility,” http://www.aeronautics.nasa.gov/atp/facilities/ntf/index.html.

3) “European Transonic Windtunnel,” http://www.etw.de/index.html.

4) Le Sant Y., “A Model Deformation Method Applied to PSP Measurement,” Proceedings of 20th International Congress on Instrumentation in Aerospace Simulation Facilities, (2003).

5) Burner A. W., Goad W. K., Massey E. A., Goad L. R., Goodliff S. L., and Bissett O. W., “Wing deformation measurements of the DLR-F6 transport configuration in National Transonic Facility,” AIAA paper 2008-6921, (2008).

Table 4: Obtained aerodynamic coefficient of Case 5.

CL CD CD, pres CD, f ric CM

実験値

0.5545 0.0320 — — -0.1111

初期形状

0.5902 0.0356 0.0217 0.0139 -0.1261

Z21 0.5767 0.0348 0.0208 0.0140 -0.1227

XYZ 0.5726 0.0346 0.0206 0.0140 -0.1204

Z44 0.5739 0.0347 0.0207 0.0140 -0.1218

6) Pallek D., Bütefisch K. A., Quest J., and Strudthoff W., “Model deformation measurement in ETW using the Moiré technique,” Proceedings of 20th International Congress on Instrumentation in Aerospace Simulation Facilities, (2003).

7) Yasue K. and Sawada K., “Static Aeroelasticity Analysis of Wind Tunnel Model Using Discontinuous Galerkin CFD Solver,” AIAA Paper 2009-604, (2009).

8) Yasue K. and Sawada K., “Effect of Model Deformation on Aerodynamic Coefficients for AGARD-B Wind Tunnel Model,” Transactions of the Japan Society for Aeronautical and Space Sciences, (2011).

9) Keye S. and Rudnik R., “Aero-elastic Simulation of DLR’s F6 Transport Aircraft Configuration and Comparison to Experimental Data,” AIAA Paper 2009-580, (2009).

10) “3rd AIAA Drag Presiction Workshop,” http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-dpw/Workshop3/

11)

加藤，中北，栗田，中島，山谷

, “風洞試験におけるマーカーを用いた写真測量法による模型変形量計測,”

第

48

回飛行機シンポジウム

, (2009).

12) Hashimoto A., Murakami K. and Aoyama T., “JAXA Digital/Analog Hybrid Wind Tunnel: Decelopment of Digital Wind Tunnel,” Proceedings of the 2nd Workshop on Integration of EFD and CFD, JAXA-SP-09-003 (2003).

13) Hashimoto A., Murakami K., Aoyama T. and Lahur P., “Lift and Drag Prediction Using Automatic Hexahedra Grid Generation Method,” AIAA Paper 2009-1365, (2009).

14) Morton S. A., Melville R. B. and Visbal M. R., “Accuracy and Coupling Issues of Aeroelastic Navier-Stokes Solutions on Deforming Meshes,” Journal of Aircraft, 35(1998), 798-805.

15) C. B. Allen, “Parallel Universal Approach to Mesh Motion and Application to Rotors in Forward Flight,” Int. J.

Numer. Meth. Engng., 69(2007), 2126-2149.

16)

菱田， 橋本， 保江， 村上

, “高速 CFD

空間格子変形法の検討

,”

第

43

回流体力学講演会

/

航空宇宙数値シミュ

レーション技術シンポジウム

2011., (2011).

(a)

η

= 0.150

(a) CL

(d)

η

= 0.638 (e)

η

= 0.847

(b)

η

= 0.331

(b) CD

(c)

η

= 0.514

(c) CM

Fig. 12: Obtained pressure coefficient (Cp) profiles of Case 5.

Fig. 14: Obtained drag polar curve. Fig. 15: Spanwise loading profiles of Case 5.

Fig. 13: Obtained aerodynamic coefficients.

付録 A.

　x

の一次導関数および二次導関数を以下に示す．

変換行列

T^m

の一次導関数および二次導関数を以下に示す．

付録 B.

　ξ

_d

の一次導関数および二次導関数は以下のようになる．