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

回転座標系における着氷現象の数値モデリング
Numerical Modeling of Ice Accretion Phenomena in Rotational Coordinate System○ 松浦大輝, 東理大院, 東京都千代田区九段北 1-14-6, E-mail:matsu.yamalab@gmail.com
鈴木正也, 東理大, 東京都千代田区九段北 1-14-6, E-mail: masaya@rs.kagu.tus.ac.jp
山本誠, 東理大, 東京都千代田区九段北 1-14-6, E-mail:yamamoto@rs.kagu.tus.ac.jp
宍戸進一郎, ＩＨＩ, 東京都西多摩郡瑞穂町殿ヶ谷 229, E-mail: shinichirou_shishido@ihi.co.jp
室岡武, ＩＨＩ, 東京都西多摩郡瑞穂町殿ヶ谷 229, E-mail: takeshi_murooka@ihi.co.jp
宮川浩, ＩＨＩ, 東京都西多摩郡瑞穂町殿ヶ谷 229, E-mail: hiroshi_miyagawa@ihi.co.jp
Taiki MATSUURA, Grad. School of Tokyo Univ. of Science, 1-14-6 Kudankita, Chiyoda-ku, Tokyo
Masaya SUZUKI, Tokyo Univ. of Science, 1-14-6 Kudankita, Chiyoda-ku, Tokyo
Makoto YAMAMOTO, Tokyo Univ. of Science, 1-14-6 Kudankita, Chiyoda-ku, Tokyo
Shishido SHINICHIROU, IHI Corporation, 229, Tonogaya, Mizuho-machi, Nishitama-gun, Tokyo
Takeshi MUROOKA, IHI Corporation, 229, Tonogaya, Mizuho-machi, Nishitama-gun, TokyoHiroshi MIYAGAWA, IHI Corporation, 229, Tonogaya, Mizuho-machi, Nishitama-gun, TokyoIce accretion is the phenomenon that super-cooled water droplets impinge and accrete on a body. It is well known thatice accretion on blades and airfoils leads to performance degradation and severe accidents. For this reason,
experimental investigations have been carried out using flight tests or icing tunnels. However, it is too expensive,dangerous, and difficult to set actual icing conditions. Hence, computational fluid dynamics is useful to predict iceaccretion. A rotor blade is one of jet engine components where ice accretes and the centrifugal force sheds the ice.
Therefore, the purpose of this study is developing the model of the ice shedding in rotational coordinate system forturbo-machine, and clarifying the icing phenomena on a rotor blade. The authors carried out icing calculation on therotor blade of a commercial axial blower, and ice accretion on the rotor blade is numerically investigated.

１． 緒言
着氷とは，大気中の過冷却水滴が物体に衝突することによって，固体表面上に氷層を形成する現象である．
航空機における着氷は，離着陸時や雲を通過する際に生じ，航空機の性能を低下させる原因となる．たとえばジェットエンジンインテークにおける着氷の発生は，エンジンに流入する空気の流れを乱すほか，剥離した氷片が飛散した際にそれがエンジンに吸い込まれ，エンジン内部に機械的損傷を与える恐れがある．このため，エンジンの型式承認時に着氷に対する安全性を確認するためのテストが課されており，この際に大規模な着氷が発生すると，改良設計への対応が非常に困難である．それゆえ，着氷形状とそれによる流れの変化を，設計の初期段階で予測することは，設計上重要であると考えられる．
しかし，着氷は流体力学と熱力学の側面を併せ持つ複雑な現象であり，また実験によって着氷現象を繰り返し再現する事は容易ではなく，地上での模擬試験や実機を用いた飛行試験には多額の費用を要し，危険を伴うことすらある．このため，数値計算による着氷シミュレーションへの関心が高まっている．
エンジンファンなどの動翼は着氷が生じやすい要素の一つである．しかしながら主翼や尾翼に比べ，研究例は少ない．著者らは過去の研究(1)でエンジンファンを想定した市販の軸流送風機の動翼に対する着氷現象を予測し，着氷形状に関する知見を得た．
動翼に特有の現象として，氷の飛散・脱落（シェディング）の危険性があることは知られているが，そのメカニズムには未解明の点が多い．簑谷ら(2)は，軸流送風機の動翼に対する着氷現象の実験的研究を行い，ファンブレードのスタガ角による着氷形状とシェディングへの影響を明らかにした．また，Papadakis ら(3)は，氷片を模擬した長方形平板にかかる力とモーメント係数を風洞設備で測定し，得られた空力データを用い，主翼から飛散する氷片の追跡計算を行った．Baruzzi ら(4)は氷片を模擬した立方体に対する流体構造連成問題を解き，飛散する氷片の運動を解くことに成功した．彼らの研究では飛散した氷の追跡は行われているが，現状では成長した氷片のシェディングを模擬できる数値モデルは存在せず，着氷とシェディングの双方を考慮した数値計算は行われていない．
したがって，本研究では，氷成長に伴うシェディングを模擬できる新たなモデルの開発を行い，軸流送風機の動翼に対するシェディングを伴う着氷現象を予測する．
２． 数値解法
２.１． 概要
本研究では着氷現象を以下のように弱練成解析する．
(1) 流れ場計算を行い，収束解を得る．
(2) 水滴を投入し，軌道を追跡する．このとき，壁面との衝突判定を行い，水滴衝突特性を得る．
(3) 水滴の収集効率を求め，熱力学計算を行う．これより，動翼表面上における氷層の成長率を算出し，着氷形状を決定する．
(4) 着氷形状に対し氷のシェディング判定を行い，条件を満たす場合，氷形状の変更を行う．
２.２． 流れ場計算
本研究では，流れ場を圧縮性乱流とし，支配方程式としてファーブル平均された連続の式，ナビエ・ストークス方程式，エネルギ式を用いる．乱流モデルとして Launder-Spalding (1974)の高レイノルズ数型 k-εモデル(5)を適用する．非粘性項は Yee-Harten (1987)の二次精度風上型 TVD スキーム(6)，粘性項は二次精度中心差分でそれぞれ離散化し，時間進行には四段階ルンゲクッタ法を用いる．
２.３． 水滴軌道計算
水滴軌道計算は，水滴の衝突質量分布（収集効率）を得るため
1

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C10-3
北海道工業大学にて行われた氷付着力試験により，物体と氷の間の付着力が得られており，これを氷のシェディング判定に用いた．氷がシェディングする条件は以下の式で表される．に行われる．計算はラグランジュ的手法により行われ，水滴にかかる力は抗力・遠心力・コリオリ力とする．また，水滴は完全な球形であり，衝突・変形・合体・分裂は無視し，流れ場に影響を及ぼさない(one-way coupling)と仮定する．

∫ ρ rω dV ≥ τ2

ice

(1)

Eim + Ein + E ac + q f = Eva + Eou + qc

(2)

× S all(3)

ここでρice は翼の各表面セルにおける氷の密度， rω2はそのセルに働く遠心力，dV はそのセル上の氷の体積，τlim は氷付着力試験により求められた付着力，Sall は，着氷している部分における翼の面積を表している．つまり，全セルの氷にかかる遠心力の合計が氷の付着力を超えたときにシェディングが生じると判定する．

２.４． 熱力学・氷堆積計算
本計算では，1 セルを検査体積とし，そこでの質量バランス及びエネルギバランスを計算することにより，検査体積内に堆積する着氷質量を決定する．図 1 に概念図を示す．質量バランスとエネルギバランスは，以下の式で表される(7)．

mim + min = mva + mou + maclim

３． 計算条件
計算対象は市販の軸流送風機(昭和電気製，快流シリーズ
A2D6H-411)であり，モータ出力は11kW，定格回転数は2920 rpm，ケーシング内径は 630 mm，設計流量は 8.70 kg/s である．計算条
件の詳細を表 1 に示す．計算格子は図 2 に示すように，重合格子法に基づき，流路に沿ったメイングリッド(図中青)と翼周りのサブグリッド(図中赤)で構成される．総格子点数は，約 150 万点であり，このうちメイングリッドが約 90 万点，サブグリッドが約60 万点である．本研究では，スイープによる着氷・シェディングへの影響を明らかにするために，図 3(b)に示す通常翼と図 3(c)に示すスイープ翼の 2 ケースについて計算を行った．計算は回転を考慮した相対座標系で行われ，境界条件は次のように設定する．
流入境界では全温・全圧・流入角度・乱れの諸量を固定し，流出境界では静圧を固定する．壁面は壁関数を用いた滑りなし・断熱条件を課す．周方向には周期条件を適用する．なお，今回は簡単のために，チップクリアランスを無視しており，ケーシング壁面のみ滑りあり・断熱条件を課している．

ここで，m は質量流量，E はエネルギを表している．添え字は imが対象表面に衝突する水滴，in が周囲のセルから検査体積への流入（ランバック・イン），ac が検査体積内に堆積する氷の潜熱，va が蒸発・昇華，ou が検査体積からの流出（ランバック・アウト）である．また qf は摩擦によって空気から検査体積に与えられる熱エネルギ，そして qc は空気の対流により損失する熱エネルギである．ここで本研究では回転体を取り扱うため min および mou の輸送方向は，流体力および遠心力の合力により判定される．
２.５． シェディング判定
本研究では氷のシェディングを模擬するため以下の仮定をおこなった．
(1) シェディングは遠心力のみに依存して起こり，流体力等の他の外力は無視する．
(2) 翼面上のすべての氷は同時にシェディングする．

zxy(a) Bird viewFig. 2 Computational gridFig. 1 Schematic energy balance of ice accretion(a) Overview(b) Blade surface(b) Normal blade(c) Swept blade

Fig. 3 Computational target
2Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C10-3(a) Rotor passage080(b) Suction Side near hub[m/s]Fig. 4 Streamline around rotor blade(a) Top view(b) Front view
065

[m/s]

Fig. 5 Droplet trajectory
Table 1 Calculation Condition
Inlet total pressure[kPa]
Inlet total temperature[K]
Rotational speed[rpm]
Mass flow rate[kg/s]
Liquid water content[g/m3]
Mean volume diameter of droplet[mm]
Exposure time of droplet
[s]

101.3
268.15
1800
5.21
1.54
30300

４． 結果と考察
４.１． 軸流送風機の流れ場と水滴軌道
本節では，軸流送風機内部の流れとそれに基づく水滴の軌道の様子を示す．通常翼とスイープ翼に関して同様の傾向が見られるため，通常翼に焦点をあて説明する．
図 4 は着氷前の動翼まわりの流線である．色は流速を意味し，相対座標系で表示されている．このため，ハブ側では低速に，チップに向かうにつれて高速になるといった傾向が見て取れる．この送風機は亜音速翼列であり，計算条件の範囲では，流入マッハ数で 0.1 以下，局所的なマッハ数の最大値で 0.3 程度に収まってい

T.E.

P.S.
0

L.E.

S.S.

T.E.
1.0

Fig. 6 Collection efficiency on normal blade wall
3

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C10-3
いる．すなわち，値が大きくなるにつれて，水滴の衝突数が多いことになる．図の左から正圧面，前縁，そして負圧面を表している．これより前縁で水滴が顕著に衝突している様子が伺える．しかし，ハブ近傍では水滴の衝突がないこともわかる．これは前述の通り，流れの剥離と遠心力が軌道に影響を及ぼしたためであると考えられる．また，正圧面にも水滴の衝突が見られ，多くの水滴が正圧面に衝突している様子が見てとれる．正圧面の一部は水滴の衝突が少ない場所が見受けられるが，これはねじれ等の翼形状の影響により，水滴が翼面に衝突せず，流れに追従するためだと考えられる．

る．図 4(b)に示す通り，負圧面のハブ近傍で流れの剥離が見られる．これは半径方向の圧力差や遠心力等の外力により，流れが引き寄せられるため起こると考えられる．
図 5 に代表的な水滴の軌道のパターンを示す．なお，流れ場同様，相対座標系で示されている．図 5(a)は上面から，図 5(b)は前面から見た軌道を示す．
流入した水滴は，流れの影響を強く受け，ハブ付近では低速であり，チップに近づくにつれて速度を増していることがわかる．ハブ近傍では前述の剥離と遠心力の影響により半径方向に水滴が向かっている様子が伺える．またそれ以外の場所では，水滴は直進する傾向にある．これは流路方向の速度が上昇しているため，遠心力よりも慣性力の方が勝るからであると考えられる．
図 6 に壁面における水分の収集効率を示す．これはその翼面位置における単位面積あたりの水分量を最大値で正規化して表して４.２． 着氷形状
本節では，通常翼とスイープ翼の着氷履歴を示し，着氷とシェディングの傾向の違いについて議論をおこなっていく．

1s

11s

21s

31s

41s

51s

61s

71s

81s

91s

101s

111s

121s

131s

141s

(a) Normal blade
1s

11s

21s

31s

41s

51s

61s

71s

81s

91s

101s

111s

121s

131s

141s

(b) Swept blade
0.0

9.5[mm]Fig. 7 History of ice geometry and thickness
4Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C10-3
図 9 に式(3)中の左辺と右辺の時間変化を示す．ここで実線は左辺，点線は右辺を示し，青が通常翼，赤がスイープ翼である．式の類推から左辺は着氷体積の成長傾向を表し，右辺は着氷面積の成長傾向を表すと言える．
従って，左辺は氷の高さに依存するが，右辺は氷の高さには依存しない．図 9 より，右辺（着氷面積）の増加量において，通常翼よりもスイープ翼のほうが早い段階で頭打ちになっているのが見てとれる．図 7 の着氷形状履歴からもこの傾向は見られる．これよりスイープ翼が早期にシェディングすると考えられる．
スイープ翼の方が着氷面積が増大しない理由を考える．図 10，11 に通常翼とスイープ翼の正圧面における収集効率と温度分布を示す．
図 11 のコンター上限は水の凍結点(273.15K)としてある．
図 10 より，水滴は正圧面にまんべんなく衝突しており，両者に決定的な差はない．このため，水滴軌道が原因とは考えにくい．一方，図 11 を見てみると，通常翼よりもスイープ翼の方が凍結限界（スパン方向での凍結点を超える位置）が低いことがわかる．これにより，スイープ翼の方が着氷できない領域が広く，着氷面積が増大しづらい結果となったと考えられる．以上より，翼面温度が原因であることが明らかとなった．

図 7(a)に通常翼の着氷の様子を示し，図 7(b)にスイープ翼の着氷の様子を示す．履歴は 1s から 141s まで 10s おきにまとめてあり，色は氷の厚みを示す．それぞれに共通の傾向を述べると，着氷は前縁のハブから少し離れた部分に多く発生することがわかる．
これは前述の通り，前縁に水滴が多く衝突しているからである．
また，チップに近づくにつれて，着氷量は少なくなり，ケーシン
グ近傍では全く着氷が発生していないことも見てとれる．これはチップ付近では静温が水の凍結点(273.15K)よりも高くなるためである．
履歴に着目すると，通常翼とスイープ翼の双方である時間まで着氷量が増大し，その後着氷が減少し，また氷が成長するという傾向が繰り返されるのが見てとれる．これは，本研究で導入したシェディングモデルにより，氷のシェディングが起こったためである．
図 8 に時間に対する着氷体積の変化を示す．これからも同様の傾向が見てとれる．ここで，通常翼が初めてシェディングする時間は 79s であり，スイープ翼では 73s であった．つまり，通常翼よりもスイープ翼のほうがわずかに早くシェディングすることがわかる．
また，着氷の絶対量についてもスイープ翼の方が少ない．
この原因に関して以下に考察する．

4500

1.6

4000

1.4

Value of LHS&RHS

Amount of accreted ice [10-5m3]

5000

1.8

1.2
10.8
0.6
0.4

3500
3000
2500
2000

0.2

1000

0500

0100

200

Exposure time[s]0Normal RHS
Swept LHS
Swept RHS

0300
Normal
Swept

020

4060

80

Exposure time [s]Fig. 8 History of amount of accreted ice(a) Normal

Normal LHS

1500

Fig. 9 Balance of LHS and RHS in Eq. (3)(a) Normal(b) Swept
1.0(b) Swept
273.15 [K]268.15

Fig. 10 Collection efficiency on blade wall of P.S.

Fig. 11 Temperature on blade wall of P.S.
5Copyright © 2011 by JSFM

５． 結言回転座標系において働く遠心力によるシェディング現象をモデル化し，軸流送風機の動翼における着氷・シェディングを再現できた．また，本研究により以下の知見を得た．
(1) 水滴はハブ近傍では半径方向に移動する．これにより，ハブ近傍では着氷が発生しない．
(2) 着氷はハブから少し離れた，前縁に多く発生し，チップに向かうにつれて着氷量は減少する．
(3) 通常翼よりもスイープ翼の方が翼面温度のため，わずかに早くシェディングする．また，着氷量の絶対値も小さくなる．

(2)

(3)

(4)

謝辞
本研究を進めるにあたりまして，氷付着力計測試験の実験デー
タをご提供いただいた北海道工業大学 平元先生に感謝いたします．

(5)

参考文献
(1) Matsuura, T., Suzuki, M., Yamamoto, M., Shishido, S., Murooka, T.
and Miyagawa, H., “Numerical Simulation of Ice Accretion

(7)

(6)

6

第 25 回数値流体力学シンポジウム
C10-3
Phenomena on Rotor Blade of Axial Blower,” Proceedings of the10th International Symposium on Experimental and Computational
Aerothermodynamics of Internal Flows, (2011), ISAIF10-103.
簑谷，平元，室岡，宍戸, “軸流ファンブレードへの着氷におけるスタガ角の影響,”日本機械学会 2011 年次大会講演論
文集, (2011), G050054.
Papadakis, M., Yeong, H., Suares, I. G. and Jacob, J., “Experimentaland Computational Investigation of Ice Shedding from Aircraft
Surfaces,” (2006), AIAA 2006-1010.
Baruzzi, G. S., Lagace, P., Aube, M.S. and Habashi, W. G.,“Development of a Shed-Ice Trajectory Simulation inFENSAP-ICE,” (2007), SAE technical paper series 2007-01-3360.
Launder, B. E., Spalding, D. B., “The Numerical Computation ofTurbulent Flows,” Computer Methods in Applied Mechanics andEngineering, Vol.3 No.2, (1974), pp. 269-289.
Yee, H. C., “Upwind and Symmetric Shock-capturing Schemes,”(1987), NASA-TM-89464.
Messinger, B. L. “Equilibrium temperature of an unheated icing
surface as a function of air speed,” Journal of Aeronautical Sciences,January, (1953), pp. 29-42

Copyright © 2011 by JSFM