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

DES による軸流圧縮機動翼列の失速発現メカニズムの解析
Detached Eddy Simulation for Onset Mechanisms of Rotating Stallin an Axial Flow Compressor Rotor○ 喜久田 啓明, 九大院, 〒819-0395 福岡市西区元岡 744, E-mail: kikuta@haira.mech.kyushu-u.ac.jp
山田 和豊, 九大, 〒819-0395 福岡市西区元岡 744
郡司嶋 智, 九大院, 〒819-0395 福岡市西区元岡 744
原 靖典, 九大院, 〒819-0395 福岡市西区元岡 744
石飛 悠希, 九大, 〒819-0395 福岡市西区元岡 744
古川 雅人, 九大, 〒819-0395 福岡市西区元岡 744
Hiroaki KIKUTA, Kyushu Univ. Graduate School, 744 Motooka, Nishi-ku, Fukuoka, 819-0395
Kazutoyo YAMADA, Kyushu Univ., 744 Motooka, Nishi-ku, Fukuoka, 819-0395
Satoshi GUNJISHIMA, Kyushu Univ. Graduate School, 744 Motooka, Nishi-ku, Fukuoka, 819-0395
Yasunori HARA, Kyushu Univ. Graduate School, 744 Motooka, Nishi-ku, Fukuoka, 819-0395
Yuki ISHITOBI, Kyushu Univ., 744 Motooka, Nishi-ku, Fukuoka, 819-0395
Masato FURUKAWA, Kyushu Univ., 744 Motooka, Nishi-ku, Fukuoka, 819-0395
DES(Detached Eddy Simulation), a hybrid simulation of RANS(Reynolds Averaged Navier-Stokes equations) andLES(Large Eddy Simulation), becomes increasingly important in CFD analyses for unsteady flow fields inturbomachineries, as a method to simulate high Reynolds number flow with feasible computational cost. In this paper,from a standpoint of concrete and practical application, DES, based on k-ω two equation turbulence model, has been
carried out for flow fields at stall onset in a low-speed axial flow compressor rotor with variable tip clearance. By
comparing the results to unsteady data acquired in measurements, availability of DES is shown for unsteady flow fields
and also the details of the flow fields are shown for each tip clearance.

１．緒言
軸流圧縮機において発生する旋回失速は，性能の急激な低下や構造体の疲労破壊等を引き起こす異常流動現象である．この失速
発現メカニズムの解明は，失速検知や失速マージンの拡大等の問題に関して有用な知見を与える．失速発現時の流れ場は過渡的な非定常内部流れであり，計測のみによって現象を十分に把握するのは不可能なため，CFD による解析が不可欠である．これまでに筆者らは計算コストを抑えつつも精度よくターボ機械内の非定常流れ場を予測する CFD の手法として DES（Detached Eddy
Simulation）が有効であることを示してきた(1)．本報では，DES の
具体的および応用的な利用の立場から，可変の翼端隙間を有する
軸流圧縮機における失速発現時の流れ場に対し非定常計測およびDES を行い，両者の結果を比較して DES の非定常流れの解析における有効性を確認するとともに，失速発現時の流れ場の詳細を示す．

ーボ機械において問題となる非定常流動現象では，はく離渦流れ
構造が空間的に大きく変動することから， LES領域とRANS領域の決定にはこの渦流れ構造の非定常挙動を考慮する必要がある．
この観点から，ターボ機械に対するDES解析では，k-ω二方程式乱流モデルに基づいたStreletsら(3)によるDES手法の適用が適切である．本モデルでは，長さスケールとしてk-ω乱流モデルから算出される渦の乱れ長さスケールを用いており，流れ場の時間変化に応じて各計算セル毎にRANS計算とLES計算の選択が逐次なされる結果，大規模なはく離を伴う翼の失速現象を高精度に予測可能であり，DES解析は航空宇宙分野を中心にして様々な流れ場に適用されている．以下にk-ω二方程式乱流モデルに基づくDES解析手法(3)の概要について述べる．
k   二方程式モデルにおいて，乱流運動エネルギー k の散逸項である Dk   *k は渦の乱れ長さスケール lk  を用いて，次のように表わされる．
Dk   *k 

２．解析手法
２．１ DES 解析手法
DES用乱流モデルでは，既存のRANS用乱流モデルを基礎方程式として渦粘性が評価されるとともに，LES解析が可能な領域においては，RANS用乱流モデル中の長さスケールを格子幅スケー
ルに置き換えることで評価するSGS（Sub Grid Scale）粘性が用いられる．その際，一般にRANS用乱流モデルには一方程式形のSpalart-Allmarasモデルが用いられることが多い．DES 用のSpalart-Allmarasモデル(2)では，渦粘性の生成と崩壊には壁面からの距離が重要な役割を果たすとの観点から，長さスケールとして壁面からの距離が用いられる．このため，同モデルではLES領域とRANS領域は与えられた格子点分布のみによって一意に決まり，
流れ場の時間的変化に対しては無関係となる．しかしながら，タk 3/ 2lk ,

lk  k 1/ 2(2)

 *

この k   乱流モデルをベースにして DES 計算を行う際には，上式中の乱れ長さスケール lk  を，次式から算出される lDES に置き換える．
lDES  min  lk  , lSGS (3)

上式中の lSGS は格子幅スケールであり，次式で与えられる．

lSGS  CDES i , j , k , i , j , k  max  i ,  j , k 

1

(4)

ここに， i ,  j , k は各計算格子点における３方向の格子幅ある．
また CDES は DES でのモデル定数であり，流れ場が定常な一様等方性乱流である場合，その値は 0.78 がとられる．この DES モデ
Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A01-2
数法を用いていない．なお，DESによる非定常解析の時間刻み幅
は，１ブレードパッシング当たりτ=1.0%翼端コードで212.5ステップ，τ=3.0%翼端コードで425.0ステップに設定した．

ルでは，各計算セルでの乱れ長さスケール lk  がその点における格子幅スケール lDES で捉えられているかが式(2)によって評価されることにより，LES 領域と RANS 領域の決定が計算結果に応じてなされる．すなわち，格子幅スケールが乱れ長さスケールより大きな計算セル（ lk   lDES ）では，LES 解析ができないものとして lDES には乱れ長さスケール lk  が設定されて，RANS（ 乱流モデル）に相当する渦粘性が算出される．一方，格子幅スケール
が乱れ長さスケールより小さな計算セル（ lk   lDES ）では，LES解析が可能であるとして lDES には格子幅スケール lSGS が設定されて，LES の Smagorinsky モデル(4)に相当する SGS 粘性が算出される．
DESモデルとしては，RANSモデルに基づいたDESモデルのほかに，RANS領域とLES領域でそれぞれ異なる乱流モデル（Baldwin-LomaxモデルとSmagorinskyモデルの組合せ等）を用いる手法(5)(6)(7)も数多く提案されている．しかしながら，これらの手法ではLES領域とRANS領域間で使用している乱流モデルが異なるために，両領域の接続境界において平均速度勾配が過大となる
結果，速度分布に不連続が生じてしまう問題がある．このため，
この接続境界においてLESとRANSの解をいかに滑らかに接続するかということに様々な試行錯誤がなされている．一方，上述のk   二方程式乱流モデルに基づくDESでは，LES領域とRANS
領域において乱流運動エネルギー k および比散逸率  に関する
同一の輸送方程式が解かれるために，両領域の接続境界で渦粘性は滑らかに分布し，特別な処理を行うことなく，接続境界における不連続の発生を防ぐことができる．これは，lDESがlk - ωからlSGSに切り替わった直後では渦粘性が通常のLESよりも大きく評価され，一種のVLES（Veri Large Eddy Simulation）となることによる．
２．２ 実験解析手法
ケーシング壁面上に高応答圧力センサ（Kulite XCS-062）を図 2（図中の黒丸がセンサ位置）のように配置して同時面計測(8)（Synchronous Field Measurement Technique）を行った．
３．解析対象および計算格子
本解析対象である供試低速軸流圧縮機の概略図を図1に示す．
ケーシング内径を450mmで，試験部は入口案内羽根，静翼1，動翼，
静翼2，出口案内羽根で構成され，設計回転数は1800rpmである．
動静翼を組み合わせた段落は，設計流量係数φ（主流の軸方向速
度を動翼先端周速で無次元化）0.5，設計圧力上昇係数ψ（全圧上昇を動翼先端周速に相当する動圧で無次元化）0.43，ハブ比0.7の軸流圧縮機中間段を対象としたものである．動翼入口における設計渦形式は平均半径で反動度50％の一定旋回角形式である．動翼
枚数24枚，静翼枚数22枚，翼弦長50.1mmであり，翼端隙間は，動翼部のケーシングを交換することにより可変であり，τ=1%翼端コード（0.5mm）とτ=3%翼端コード（1.5mm）である．以下，それぞれの翼端隙間を有する場合の対象をτ=1%翼端コードおよびτ=3%翼端コードと呼称する．
DES解析においては，図3に示すとおり，動翼列の上流と下流の静翼列（前置静翼および後置静翼）まで含めた1.5段を計算対象とし，τ=1%翼端コードでは1/3周（8ピッチ），τ=3%翼端コードでは全周(24ピッチ)を計算領域とした．本解析では静翼および動翼
枚数ともに同数の24枚とした．計算格子はH型の構造格子を用い，
τ=1%翼端コードでは計算領域全体（1/3周）で3,086万セル（動翼：208万セル/ピッチ，前置静翼：84万セル/ピッチ，後置静翼：94万
セル/ピッチ），τ=3%翼端コードでは，計算領域全体（24ピッチ）で約12,109万セル（動翼：273万セル/ピッチ，前静翼：109万セル/
ピッチ，後置静翼：122万セル/ピッチ）を設定した．壁面に対する最小格子幅は3.0×10-5 [m] (y+<1)とし，壁面境界条件として壁関
2

４．解析結果４．１ 性能曲線
図 5(a)および(b)に，それぞれτ=1%翼端コードとτ=3%翼端コードの性能曲線を示す．どちらの場合も，実験と DES の結果はよく一致している．特に，τ=1%翼端コードの実験の失速初生点は
φ=0.338 であり，DES の初生点と一致する．ここで，φ=0.337 における DES の結果は，実験結果から大きく外れている．これは，実験においては失速に至ると自然に流量が減少するが，本 DES においては失速の初生を予測するという観点から流量を固定して計算を行っており，これらの条件の違いからくるものである．
４．２ τ=1%翼端コードとτ=3%翼端コードの擾乱の比較
図 4 に実験により得られたτ=1%翼端コードとτ=3%翼端コードの動翼 4.1mm 上流におけるケーシング壁面の静圧変動波形（青）とそのローパスフィルタ波形（赤：翼通過周波数の 50%以上をカット）を示す．τ=1%翼端コードでは，図 4(a)に示す通り，
283 回転目から発生した擾乱が10～20回転のうちに発達して旋回失速に至り，その途中で擾乱が消失することはない．一方，図 4(b)に示す通り，τ=3%翼端コードでは 216 回転目から擾乱が発生す
るが，261 回転目には擾乱が消失する．このようにτ=3%翼端コ
ードでは，失速点近傍で擾乱が発生と消失を繰り返す．
本供試軸流圧縮機では，上記のように，翼端隙間の大きさによって失速に関わる擾乱の発現の仕方が異なる．同時面計測および
DES の結果から，それぞれの擾乱発生時の流れ場を詳しく解析することで，失速発現メカニズムを明らかにする．
４．３ τ=1%翼端コードの場合
同時面計測結果より，失速初生時の擾乱において段階的な発達が見られた．擾乱発生直後には弱い低圧域が翼間を移動し，ある程度擾乱が発達した後には強い低圧域が翼間を移動する様子が見られた．
同時面計測結果である図 6 の T=0.32~1.28（T は無次元時間で1.00 が翼通過周期に対応する．
）
における黒点線で囲まれた領域は翼間を移動する弱い低圧域であり，それに対応する DES 結果（相対系）が図 7 である．図 7 の上段は動翼列における無次元ヘリシティ（Hn）で色づけされた渦コアを示している．下段は翼端付近の渦コアおよびケーシング壁面静圧を示している．翼端側で前縁剥離が発生する（T=0.12，上段）と，翼端もれ渦が小規模に崩壊
し，翼端もれ渦の一部が分離する．その分離した渦は隣接翼に向かって移動する（T=0.33，上段）．この渦はケーシング壁面と干渉し，干渉した領域において弱い低圧域を形成する．そのため，ケーシング壁面上に翼間を移動する弱い低圧域が現れる（T=0.33，
下段）．
同時面計測結果である図 8 の T=0.00~1.28 における黒点線で囲まれた領域は翼間を移動する強い低圧域であり，それに対応する
DES 結果（相対系）が図 9（図 7 と同様に表示）である．Blade3の負圧面から竜巻状のはく離渦がケーシング壁面に向かって立ち上がり（T=1.33，上段），はく離渦がケーシング壁面と干渉している領域において強い低圧域を形成する．この低圧域は，隣接翼に向かって移動する（T=1.33，下段）．隣接翼に近づくと，隣接翼の翼端側で前縁はく離が発生し，前縁はく離渦が翼端もれ渦と連結する．連結した渦は翼面から離れ（T=1.67，上段），やがて竜巻状のはく離渦に発達する（T=1.76，上段）．このようにして次々に竜巻状のはく離渦が発生し翼間を移動することで，ケーシング壁面
上に翼間を移動する強い低圧域が現れる（T=1.33 および 1.76，下段）．
Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A01-2
４．４ τ=3%翼端コードの場合
図 10 に失速点近傍における擾乱発生時のケーシング壁面静圧とその変動値の分布（同時面計測結果）を示す．ここで，変動値とは，注目する時刻における壁面静圧の瞬間値と，その時刻と同じ動翼回転位相時における擾乱発生が見られない状態での壁面静圧のアンサンブル平均値との差として定義している．図中の黒点
線 z1～z4 は変動値が特徴的な領域を示している．z1 は負の変動
値が大きな領域であり，Blade1 から Blade2 へと移動し，隣接翼（Blade2）に干渉する．また，Blade1 の動翼前縁付近の領域 z2 は
z1 が Blade1 から離れるに従って，徐々に変動値が高くなる（T=0.00～0.48）．隣接翼 Blade2 の動翼前縁付近の領域 z3 は z1 がBlade2 に近づくにつれて，変動値が高い値から 0 付近まで変化している（T=0.48～0.96）．更に，z4 のような低圧のスポット（変動
値が負の小規模な領域）が現れる．
図 11(a)～(d)は DES 結果における相対系上の瞬時の流れ場を表示しており，上段がケーシング壁面静圧分布，中段がその変動値
分布（定義は前述のとおり），下段が無次元ヘリシティで色づけされた渦コアの構造である．図中の黒点線 z1*～z4*は変動値が特徴的な領域を示しており，図 10 の同時面計測結果における領域 z1～z4 に対応させてある．大きな負の変動値を示す領域 z1*はBlade1*を出発して動翼の回転と逆方向に移動し，Blade2*に干渉する(T=0.00～4.12)．このとき，翼端もれ渦はスパイラル型の強い崩壊を起こし，崩壊した渦の一部が隣接翼（Blade2*）前縁側へ移動する．z1*の移動は崩壊した渦の移動に対応している．z2*は z1*が Blade1*から離れるに従って変動値が高くなる（T=0.00～2.94）．
z3*は z1*が Blade2*に近づくにつれて，徐々に変動値が高くなり（T=1.65～2.94），z1*が Blade2*に到達すると変動値が急激に 0 付
近まで変化する（T=2.94～4.12）．z2*や z3*の変動値の変化は，翼端もれ渦が崩壊して流路を部分的に塞いだり，崩壊した渦が隣接翼に干渉したりすることで起こる翼負荷の変化に伴うものである．
また，崩壊した渦の一部が，ケーシング壁面と干渉することで，z4*のように低圧のスポット（変動値が負の小規模な領域）が現れる（T=4.12）．以上より，DES の結果（図 11）における領域 z1*～z4*のケーシング壁面上の変動値分布の非定常変化は，同時面計測結果（図 10）における領域 z1～z4 の非定常変化とよく一致している．DES の結果における渦流れ場を考慮することで，τ=3%
翼端コードの失速点近傍において発生する擾乱は，翼端もれ渦の崩壊により引き起こされていることがわかった．

５．結言
DES の具体的および応用的利用の立場から，軸流圧縮機の失速
発現メカニズムの解析に対して DES を適用した．
翼端隙間がτ=1%翼端コードおよびτ=3%翼端コードの場合にそれぞれ現れる 2 種類の失速発現メカニズムの解析を行い，以下の知見を得た．
（１） 同時面計測と DES の非定常流動解析結果は，互いによく一致しており，ターボ機械における DES の有効性を非定
常の流れ場において確認することができた．
（２） 翼端隙間がτ=1.0%翼端コードの場合，失速発現プロセスにおいて，失速に関わる擾乱の段階的な発達が見られた．
擾乱発生直後はケーシング壁面上に翼間を移動する弱い
低圧域が見られ，
ある程度擾乱が発達すると翼間を移動す
る強い低圧域が見られた．前者の低圧域は，翼端側の前縁剥離が起こるとともに翼端もれ渦が小規模に崩壊して一部が分離し，
分離した渦がケーシング壁面と干渉しながら翼間を移動することで現れる．後者の低圧域は，動翼負圧
面から立ち上がった竜巻状のはく離渦がケーシング壁面
上に干渉し，
その干渉領域が翼間を移動することで現れる．
（３） 翼端隙間がτ=3%翼端コードの場合，失速点近傍で発生と消失を繰り返す擾乱が現れる．
この擾乱発生時には動翼間を移動する低圧域や，低圧のスポットが現れる．これらの低圧域は，翼端もれ渦の崩壊（スパイラル型）によって引き起こされている．
参考文献(1) 古川, 岩切, 喜久田, "Detached Eddy Simulation による翼列の非定常流動解析," ターボ機械協会第 63 回総会講演会, 5(2010)(2) Spalart, P. R., Jou, W. H., Sterlets, M., and Allmaras, S. R.,
Proceedings of 1st AFOSR International Conference on DNS/LES,(1997), pp.137-147.
(3) Strelets, M., AIAA Paper, No.2001-0879 (2001).
(4) Smagorinsky, J., Monthly Weather Review, Vol.91 (1963),
pp.99-164.
(5) Camelli, F. E., and Lohner, R., AIAA Paper, No.2002-0426 (2002).
(6) Hamba, F., International Journal of Theoretical and Computational
Fluid Dynamics, No.16 (2003), pp.387-403.
(7) Kawai, S., and Fujii, K., AIAA Paper, No.2004-0068 (2004).
(8) Iwakiri, K., Furukawa, M., Tomita, I., Kameda, T., Kuroumaru, M.,
and Inoue, M., Proceedings of the International Gas Turbine
Congress 2007, Paper No. TS-046 (2007).
41.0

39.6

33.6

35.6

41.0

[mm]

CasingFlowτ =1% chord
3% chord
Stator-1

Rotor

Stator-2

67.5

Hub

(a) A schematic of test compressor(b) Detail of test section
Fig.1 Test compressor
3

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A01-2(a) Perspective view of compressor stage

Fig.2 Pressure holes on casing wall(b) Close-up view of computational grid system
Fig. 3 Computational grid system(a)τ=1% chord (stall inception)(b)τ=3% chord (near-stall)Fig.4 Pressure coefficient trace on casing wall (blue: raw waveform, red: LPF waveform, experiment)

Experiment
φ=0.338(stall inception)
DES

0.55

0.50
0.45
φ=0.337
0.40
0.35

0.60

Total Pressure Rise Coefficient ψ

Total Pressure Rise Coefficient ψ

0.60

Experiment
DES

0.55

0.50
0.45
0.40
0.35

0.30

0.30
0.20

0.30
0.40
Flow Rate φ(a) τ=1% chord

0.50: Design Point

0.20

Fig.5 Characteristics of compressor stage
4

0.30
0.40
Flow Rate φ(b) τ=3% chord

0.50: Design Point

Copyright © 2011 by JSFM

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

Blade1

rotation

0.24
Ps
-0.72

LE

Blade2TE(a)T=0.00(b)T=0.32(c)T=0.64(d)T=0.96(e)T=1.28(f)T=1.60
Fig.6 Time valiation of casing wall pressure distribution, where weak moving low pressure region (LPR) appears (τ=1% chord, experiment)TE

tip leakage vortex

LE

Blade1

Blade2LE separation
at the tip

1.0

fragmented vortex
Hn

LE separation vortex
rotationrotationTE
Blade1
Blade2-1.0

fragmented vortex

0.24
Ps

-0.72

LE

weak moving LPR

tip leakage vortex(a)T=0.00(b)T=0.12(c)T=0.33
Fig.7 Time valiation of vortex cores colored with normalized helicity (upper) and casing wall pressure distribution (lower)with vortex cores near blade tips at stall inception, where weak moving low pressure region (LPR) appears on casing wall (τ=1% chord, DES)

0.24
Ps
-0.72

LE
Blade1
rotation

Blade2TE(a)T=0.00(b)T=0.32(c)T=0.64(d)T=0.96(e)T=1.28(f)T=1.60
Fig.8 Time valiation of casing wall pressure distribution, where strong moving low pressure region (LPR) appears (τ=1% chord, experiment)TELE
tornado-type
separation vortex
Blade3rotation

Blade4

1.0

LE separation
tornado-type
separation vortex

rotationTE
Blade3
Blade4

Hn

-1.0
0.24
Ps

-0.72

LE
strong moving LPR

strong moving LPR(a)T=1.33(b)T=1.67(c)T=1.76
Fig.9 Time valiation of vortex cores colored with normalized helicity (upper) and casing wall pressure distribution (lower)with vortex cores near blade tips at stall inception, where strong moving low pressure region (LPR) appears on casing wall (τ=1% chord, DES)
5
Copyright © 2011 by JSFM

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

Blade1

Blade2

Rotation

0.24
Ps
-0.72

z1

0.12
ΔPs
-0.12

z2
z1

z2 z3

z3
z1

z1

z4

z3
z1

(a)T=0.00(b)T=0.24(c)T=0.48(d)T=0.72(e)T=0.96(f)T=1.20
Fig.10 Time variation of casing wall pressure distribution (upper) and pressure deviation field (lower) on casing wall (τ=3% chord, experiment)

Blade2* LE Blade1*
TE Rotation
z3*
z1*

z2*

z1*

z3*

0.24
Ps
-0.72

z2*

0.12
ΔPs
-0.12
1.0
Hn
-1.0

Spiral type vortex breakdown(a) T= 0.00(b) T=1.65

0.24
Ps
-0.72

z3*

z1*

z2*

z3*

z2*
z1*

0.12
ΔPs
-0.12

z4*
1.0
Hn
-1.0

(c) T= 2.94
(d) T=4.12
Fig.11 Time variation of casing wall pressure distribution (upper), pressure deviation field (middle) and tip leakage vortex cores (lower)
at near-stall, where breakdown occurs on tip leakage vortex (τ=3% chord, DES)

6

Copyright © 2011 by JSFM