先駆的科学計算に関するフォーラム2013(Forum on Advanced Scientific Computing 2013)〜GPGPU & 先端的計算科学研究プロジェクト成果報告〜

ターボ機械における不安定流動現象の発現メカニズムの解明
古川雅⼈（九州⼤学）
橋本⻯⼀（九州⼤学）
⻄澤 崇（九州⼤学）
⽯⾶悠希（九州⼤学）

発表内容
研究背景・⽬的
供試圧縮機
実験⽅法

1.
2.
3.

◦

瞬時ケーシング圧⼒場計測◦

時空間補間法

数値計算⽅法
結果

4.
5.

6.
7.

◦

圧縮機性能特性◦

スパイク失速初⽣における⾮定常ケーシング圧⼒挙動◦

スパイク初⽣の流れ構造およびメカニズム◦

クリアランス流れの“spillage”

まとめ
今後の予定

研究背景
 旋回失速やサージといった⾮定常流動現象は騒⾳発⽣や翼振動励起の原因◦ 翼の疲労破壊，圧縮機の故障現象の理解 失速マージン予測の改善

 圧縮機の⾼信頼性化，失速制御技術の開発
 スパイク形失速とは…

◦ ⼩スケール擾乱（数ピッチ程度）◦ 擾乱は数回転のうちに発達し失速◦ 翼端側での局所的な境界層剥離◦ 擾乱成⻑後，⻯巻状渦へと発達して失速状態に⾄る（Inoueら, 2001）

Flow model of stall cell

研究⽬的
 本研究の⽬的◦ 低速軸流圧縮機におけるスパイク型失速の初⽣メカニズムを明らかにすること

 “瞬時ケーシング圧⼒場計測”の実施◦ ⾼応答圧⼒センサをケーシング壁⾯の動翼１流路の領域に30個埋設◦ すべての圧⼒センサを同時に動作させることで瞬時ケーシング圧⼒場を計測◦ 時空間補間の実装
 Detached eddy simulation (DES)◦ スパイク失速初⽣の⾮定常流れ場の詳細を解析

供試圧縮機 低速軸流圧縮機 （九州⼤学）
IGV

Front
StatorRear
Rotor Stator

OGV

Rotor Blade Number

24

Design Condition

Stator Blade Number

22

Rotational Speed: 1,800 rpm

Casing Inner Diameter

450 mmHub/Tip Ratio

0.7

Flow Rate Coefficient φ: 0.500

Rotor Chord Length50.1 mm

Total Pressure‐rise Coefficient ψ : 0.43

Tip Clearance(Chord/Tip Ratio)

0.5 mm(1%)( total pressure‐rise divided by the equivalent
dynamic pressure for rotor tip speed )( mean velocity divided by rotor tip speed )実験⽅法
 瞬時ケーシング圧⼒場計測

Rotor LE

Rotor TE

Pressure coefficient :at wallin inlet chamber

dynamic pressure based on tip speed●Instantaneous Casing Pressure Field Measurement●
・30 high-response pressure sensors（Kulite XCS-062）
・All the sensors operate concurrently
・Spatial-temporal interpolation

瞬時ケーシング圧⼒場計測
 時空間補間法
Rotatione.g.) The sensor colored in red can obtain 25 data for one rotor pitch,while the rotor passes through the measurement window
Rotation

The blue circle point is measured three times, passing throughthe measurement window-> Equal-time data can be obtained for 25 points by temporallyinterpolating the data from the three different sensors

瞬時ケーシング圧⼒場計測
 時空間補間法
Point N

Temporal interpolation curve for point N

空間補間のみtime

時空間補間数値計算⽅法
 離散化 : 有限体積法

 時間積分 : Gauss‐Seidel法（完全陰解法）
 数値流束 : Roe法 + MUSCL (3rd‐order)
 乱流モデル : LRN k‐omega based DES

 計算格⼦ : 1.2億セル（全周），3,100万セル（8流路）◦ 動翼： 200万セル/流路，静翼： 100万セル/流路数値計算⼿法
 k‐omega based DES ‐ Strelets（2001） k  t ui*( k ) (  u j k )  (    t )   ijt
x j
x j 
x j 
x j (  ) (  u j ) (   t )t
x j
x j 
x jDiffusion

Convection  *  k･･･ k‐equation

 t ui  2    ij･･･ omega‐equation
k x jProduction Dissipation

k‐⽅程式の散逸項 Dk の修正
Dk    k *

k3 2lDES; wherelDES  min(lk ω, lSGS )Separation Region
Wake Region
Leakage Flow

PSTip Leakage
VortexSSPS

SS
Rotation

Rotation
LES REGIONRANS REGION

数値計算⽅法
 並列化：MPI並列（領域分割）
 並列化率：α = 0.999707◦ 実測マシン：PRIMERGY RX200 S6α◦ 問題サイズ：4,522,800セル◦ 2点での実⾏時間
 Tm = 1.88146 sec/step （ m = 24 ）
 Tn = 0.486288 sec/step （ n = 96 ）

 飽和台数：約34万 ベンチマーク結果
並列数

計算時間

28750 sec/step

3,384

4 sec/step

並列数

計算時間

3,384

4.98 sec/step

33,8400.32 sec/stepTn  Tmm 1n 1Tn Tm
mn100α
飽和台数 1-α200
175
150
125
100
75
50
25
0
0

10,00020,00030,00040,000圧縮機性能特性
 全圧上昇特性
Stall point: =0.338

mild stallcondition

 流量係数の時間履歴IIIケーシング圧⼒波形
 流量係数の時間履歴I ケーシング圧⼒波形

II瞬時ケーシング圧⼒場
 Comparison between EXP and DES
ISurrounding point A (spike emergence point)

EXPDES
Surrounding point B (flow-rate drop point)

EXPDES

II

スパイク失速初⽣スパイク失速初⽣スパイク失速初⽣
 スパイク失速初⽣における⻯巻状剥離渦クリアランス流れの”Spillage”

まとめ
1.

瞬時ケーシング圧⼒場計測◦ スパイク失速初⽣において，低圧領域が動翼前縁から離れ，周⽅向に伝播するのが確認される
2.

Detached Eddy Simulation◦ スパイク失速初⽣において，動翼先端で前縁剥離が⽣じ，それが⻯巻状剥離渦となる◦ 動翼負圧⾯に発⽣した⻯巻状剥離渦は，ケーシング壁⾯に⾜を持つ結果，ケーシング壁⾯圧⼒に低圧領域が発⽣する◦ ⻯巻状剥離渦の成⻑とともに，低圧領域は動翼前縁付近を周⽅向に移動◦ ⻯巻状剥離渦は⼤きなブロッケージ効果を有するため，隣接翼の前縁剥離を誘起する．したがって，前縁剥離に伴う⻯巻状剥離渦がスパイク型擾乱となって動翼列を伝播する

今後の予定 実機モデルにおける失速初⽣現象の解明◦ ⾼速圧縮機（翼端周速が遷⾳速，衝撃波の発⽣）◦ 多段圧縮機（14段，IGVおよびOGVを含めると30翼列）◦ 実機試験結果との⽐較
14‐stage Axial Compressor

今後の予定 計算格⼦（IGV〜7段動翼：20億セル）
IGV

1stRotor
1stStatorS
2.6
0.0