■技術論文■「京」を用いたプロペラ周りの流場解析
Estimation of Flow Field around Ship Propeller by using "K computer"
技術研究所 流体研究グループ

藤澤

概竹春

FUJISAWA Takeharu

要

スーパーコンピュータ「京」を用いてプロペラ単独試験状態の大規模 LES 計算を実施し、流場計測や従来の計算手法では解像できなかった、詳細な流場情報を得た。

Summary
Large scale LES calculation of the propeller open test state was carried out using "K computer", and detailed
flow information which could not be resolved by Experimental measurement or RANS simulation was obtained.
が小さいなどレイノルズ数が低い場合には、乱れの無1.

緒言

い規則正しい流れである層流と呼ばれる状態になるが、船舶の性能開発においては、設計した船舶の性能を流速が早くなるなどしてレイノルズ数が高くなると、建造開始前に正確に推定することが重要となる。当社

乱流と呼ばれる、流れの中で非常に微細な渦が複雑にを含め造船各社では、縮尺模型を用いた実験により性絡まり合い不規則に変動する状態になる。

能推定を行っている。この模型実験の結果から実際の造船分野においては、模型実験での流れは層流と乱船舶の性能を知るためには、尺度影響と呼ばれるサイ流が併存する遷移状態であり、実際の船舶周りにおいズの違いによる性能の違いを考慮した外挿が必要となては流場全体が乱流となっている。そこで、模型実験

るが、尺度影響を理論的に解明することは難しい。そおよび実船における性能を精度よく推定するためにはのため、実用的には蓄積した海上公試結果と模型実験乱流を推定することが重要となる。

結果の相関などの経験的な知見が用いられており、近CFD 手法の内、DNS（Direct Numerical Simulation ）年活発に開発が行われている省エネ技術を新しく適用

は Navier - Stokes 方程式をそのまま、直接解く方法で、した場合に模型実験結果から実船性能を正確に推定で原理的に正確な推定が可能であるが、乱流中に存在すきない可能性がある。

る微細な渦を渦の中心部まで含めて計算しなければい

他方、CFD（Computational Fluid Dynamics、数値流けないため、非常に大きなメモリ・CPU 性能が必要体力学）に基づく数値計算手法は、コンピュータを用となり、現時点では「京」のようなスーパーコンピュいて流体運動の基礎方程式である Navier - Stokes 方程ータでも実施困難である。近年の急速な計算機性能の式を質量保存則である連続の式と連立させて解くこと発達を考慮しても計算可能となるのは数十年以上先ので、流れの動きや流れの中の物体に働く力を推定することと思われる。そのため Navier - Stokes 方程式にモものである。現在、数値計算による推定結果を参考すデル化と呼ばれる操作を施して計算コストを小さくすることで模型実験の実施回数を削減させることが可能る工夫がされている。
Navier-Stokes 方程式のモデル化は一般に乱流モデルとなっているが、模型実験を代替するまでには至っていない。と呼ばれさまざまなものが提案されているが、いずれ一般に物体周りの流れは、流速が遅い或いはサイズの手法でも精度を上げるためには実験的に計測された-1-JMU テクニカル・レビュー No.5

（2019 年 1 月）

流場との比較により推定と実際の差を把握することが重要となる。現在当社では RANS (Reynolds Averaged
Navier - Stokes equation) と呼ばれる、乱流を時間的に平均化してモデル化する乱流モデルを用いており、新規に開発された船舶の推定を行う場合には、開発された船舶に対する数値計算に加えて、過去に実験を実施したことのある類似船の水槽試験と数値計算の比較を参考にすることで、模型実験状態での性能を推定し、そこから実際の船舶の性能を推定している。

図 1 模型プロペラ

将来的に実船の性能を直接推定することを目標とし

Fig.1 Model Propeller.

た場合、乱流モデルを用いた推定手法では実際の船舶
表 1 プロペラの要目の性能を計測した結果との比較検証が必要となる。しTable.1 Particulars of Propeller
SymbolValue
Diameter
D0.240 m
Number of Blades
Z4
Pitch Ratio
H/D0.700(Const.)
Expanded Area Ratio
EAR0.500
Boss Ratio
0.170
Angle of rake
9.000 deg
Chord Length at 0.7R
L0.7R
0.0669 m
Section Type
MAU

かし、現状では乱流モデルの検証に十分な詳細で高精度な実船周囲流場のデータを得ることは難しい。そのため、乱流モデルに依らない数値計算手法が必要とされるが、先に述べたように DNS では計算コストが実用的ではないため何らかの簡略化が必要となる。
LES ( Large Eddy Simulation ) は、流場を空間的にモデル化する方法で、一定のサイズ以下の流れのみをモデル化し、それより大きい流れはモデル化せずに直接ラ単独試験で標準的な条件であるレイノルズ数解く方法である。LES では計算格子と呼ばれる空間の分割を細かくしていくことで、モデル化の影響が小(=⁄ ) = 1.1 × 10 とした。通常本プロペラでプロペラ単独試験を行う場合 前進係数さくなっていき、推定結果が DNS での推定結果に漸⁄1）近する 。DNS と同等の結果を得るためには RANS(=)が 0.1 から 0.8 程度の範囲で実施しているが、による計算よりも数千倍の大規模な計算が必要となる大規模計算では J = 0.5 の状態を計算した。ここで Nが、DNS で必要とされる計算規模よりは小さく、現はプロペラ回転数、ν は動粘性係数、VA はプロペラ前時点でもスーパーコンピュータなら実施可能である。

進速度である。

これらの状況を踏まえ、実船尺度における乱流現象

また、本プロペラに対し 3D レーザ計測により形状を考慮した数値計算実用化の第一段階として、本研究を忠実に再現した CAD を作成し、計算に使用した。

では、船舶の性能推定に必要な模型実験の内で比較的2.2 解析方法

単純な条件であるプロペラ単独試験の状態を対象とし本研究では大規模 LES 計算を実施するとともに、て、大規模な LES 計算を実施し、その計算結果の検RANS による計算も比較のために実施した。

証を行った。

非圧縮流れの基礎方程式は次の Navier - Stokes 方程本研究の大規模 LES 計算は HPCI（革新的ハイパ式と連続の式である。

フォーマンス・コンピューティング・インフラ）産業+利用課題の採択を経て、スーパーコンピュータ「京」2)を利用して実施した（課題番号：hp150039 ） 。

=0

HPCI 産業利用については http://www.hpci-office.jp/materials/sangyo_jp_pp を参照願う。

2.

=−1+(1)(2)

ここで u は速度ベクトル、p は圧力、t は時刻、ρ は密度、ν は動粘性係数である。

解析条件および解析方法

LES に対する基礎方程式は式(1),(2)に空間フィルタ操作を施すことで以下のように得られる。

2.1 解析条件

〈 〉本研究で対象とした模型プロペラの形状と要目をFig.1、Table.1 に示す。本プロペラは典型的な MAU+〈 〉翼型のプロペラであり、計算で実施する条件もプロペ-2-〈 〉〈 〉

=0=−

1 〈 〉

+

JMU テクニカル・レビュー No.5

〈 〉

−(3)

(4)

（2019 年 1 月）

ここで、〈 〉は空間平均、 は Sub Grid Scale (SGS)応力を表す。一方、RANS の基礎方程式は式(1) , (2)のu , p を時間平均成分と変動成分に分けた上でのアンサンブル平均をとることで以下のように得られる。
+=0=−

̅+− ′ ′

(5)

(6)

ここで、 ̅ は時間平均、′ は平均からの変動、 ′ ′はレイノルズ応力である。

以下に本研究で実施した計算方法について示す。
2.2.1 LES

大規模 LES 計算にはオープンソースソフトウェア(a) for LES

FFR（FrontFlow/Red3) ）を「京」用にチューニングしたものを用いた。FFR は LES による大規模流場解析を目的として開発されたソルバである。FFR を用いた計算においては、SGS 応力項のモデル化には
Dynamic Smagorinsky モデル 5) を採用するともに、収束性を向上させて計算負荷を低減するため疑似圧縮条件 4) を導入した。対流項の離散化は 2 次中心差分を基本として、数値振動を抑制するために 1 次風上差分を 2 %の割合で混合した。時間発展の解法として(b) for RANS
図 2 計算領域
Fig.2 Calculation Regions.

は Euler 陰解法を、行列解法としては ICCG 法を、それぞれ採用している。なお、基礎方程式は静止系で与 え 、 回 転 は ALE 法 (Arbitrary LagrangianEulerian Method) 6) を用いて格子を回転させている。となる長さを作動条件と Schöenherr の摩擦抵抗式か2.2.2 RANS

ら推定すると約 0.0025 mm となる。

RANS 計算には汎用流体解析ソフトウェアである

RANS 計算格子は当社で標準的に実施している計算

ANSYS FLUENT 17.0 を用いた。離散化は有限体積法と同程度で計算格子は 100 万セル程度である。大規模に基づいており、レイノルズ応力項のモデル化には低LES 計算では計算実施可能な範囲で出来る限りプロレイノルズ数型 SST k-ω モデルを用いた。計算は定常ペラ翼面の解像度を細かくするように格子を作成して計算とし、空間離散化には QUICK スキーム 、圧力おり、計算格子は 8 億セル程度である。
格子作成には RANS 用には GAMBIT 2.4.6 及び解法には SIMPLE スキームを採用した。
基礎方程式は回転系で与え、格子は固定して回転のTGrid 17.0 を組み合わせて用いて作成した。 LES 用効果を外力として与えた。また、回転方向に周期境界には Pointwise V17 を用いて、まず格子解像度を条件を用いて１翼分の領域のみを対象とした。

Table2 の 2 倍とした約 1 億セルの格子を作成し、この2.2.3 計算格子

表 2 格子解像度

計算領域の概要を Fig.2 に示す。計算領域はプロペ

Table.2 Grid Resolutions.

ラの回転軸を中心とする円筒形であり、プロペラ直径

Location

D を基準としてプロペラ前方を 4D、後方を 10D 直径を 12D とした。Fig.2(a)のように LES ではプロペラ全Outer Boundary [mm]Propeller Edge [mm]23.5000.029

96.0000.480

体を計算しているが、Fig.2(b)の RANS 計算において

Propeller Blade [mm]0.200( abt. 80 y +)1.800示す。有次元値とともに格子解像度の評価基準としてPropeller Boss [mm]First Prism
Layer Height [mm]

0.725
0.005
( abt. 2 y +)

3.600
0.002
( abt. 1 y +)重要とされる無次元長さ y+で表した値も示す。1.0y+

Total Cells

8.07×1080.01×108は周期境界を用いた１翼分の領域のみを計算対象としている。Table.2 に各位置での代表的な格子解像度を-3-for LES

JMU テクニカル・レビュー No.5for RANS

（2019 年 1 月）

Back

1 億格子から個々の格子を各方向に 2 等分する処理

Face（Refine 処理）をすることで、8 億格子を作成した。

3.

計算結果と実験結果の比較

プロペラ単独試験は当社 技術研究所の船型試験水槽(寸法: 長さ 240m、幅 18m、深さ 8m)で行った。
条件は 2 節に示したとおりである。Fig.3 に試験結果および計算結果を示す。RANS、LES とも推力係数KT(=Thrust/ρN2D4 )、トルク係数 KQ(=Torque/ρN2D5 )を試験値より低く推定している。KT は LES の方が試(a) ReD=1.1×106

験値に近いが、KQ は RANS の方が試験値に近くなっており、どちらの推定精度が良いともいえない。大規模計算を 1 条件しか実施できていないため、推定精度についてはこれ以上の議論は差し控える。
限界流線の観察実験は当社 技術研究所のキャビテーション水槽（計測胴寸法: 長さ 2,600mm、幅600mm、高さ 600mm）で行った。油膜を塗布した状態で表面に流線が形成されるまで作動させてから、取出して写真を撮影した。キャビテーション水槽での(b) ReD=0.6×106
図 4 実験による限界流線
Fig. 4 Experimental Restricted Flow Patterns

計測では、ReD を変えた場合についても計測した。
以下の実験結果および CFD 結果可視化ではプロペラ表面に直径の 20 (CFD のみ), 30, 40, 50, 60, 70, 80,

していないが、Face 面側では翼先端側で差が出てい90, 95 %の円を表示する。

ることがわかる。模型プロペラ表面の流れには層流域と乱流域が有り 7) 、Fig.4 (a) の Face 面側の翼先端側3.1 限界流線の比較
Fig.4 に油膜法で可視化した模型プロペラ翼面の限に見られるような流線が回転方向を向く部分では、流界流線を、また Fig.5 に計算によって得られた摩擦応れが乱流状態となっているとされる。また、レイノル力ベクトルの包絡線として可視化した限界流線を、そズ数の増大により乱流領域は翼根側に広がるとされ 7)、

れぞれ示す。LES 計算の結果は時間平均している。

Fig.4 (a) (b)の Face 面ではその傾向が見られる。 限界

Fig.4 (b) に ReD が低い場合も併せて示している。

流線の比較から RANS 計算では同じレイノルズ数条

実験による限界流線は、前縁側と後端側で回転方向件の実験の流場よりも乱流となりやすく、LES 計算から少し翼先端側に向いた流線となり、その間で大きでは逆に乱流になりにくい傾向がわかる。大規模

く先端方向に曲がっている。この傾向は LES 計算結

LES 計算では、Table.3 に示したように今回用いた格果でも同様だが、RANS 計算では翼根側を除いて流線子の翼面での解像度は 8 億格子でも y+ = 80 程度であが回転方向を向いている。また、Fig.4 (a) と(b) を比り、壁面乱流の LES 計算で DNS 相当の推定を行うた

較すると、Back 面側はレイノルズ数に対しほぼ変化めに必要とされる 10 y+ の格子解像度と比べて粗いため、渦構造が実際より太くなっていると考えられる。
また、渦構造はレイノルズ数が高い程細くなることから、今回の LES 計算結果(Fig.5(a)) が実際より低いレイノルズ数の実験結果(Fig.4(b)) に近くなったものと推定される。

3.2 渦構造の比較
Fig.6 に速度勾配テンソルの第二不変量 Q の等値面により翼面近傍での渦構造を可視化した結果を示す。
図 3 試験結果および計算結果(ReD=1.1×106,J=0.5)

Fig.6(a)は瞬間流場での可視化であり、Fig.6(b)は LES

Fig.3 Measured and Computed Characteristics

計算での時間平均した流場、および Fig6(c)は RANS-4-JMU テクニカル・レビュー No.5

（2019 年 1 月）

Back

計算による時間平均流場での可視化である。Fig.6(a)Faceには部分的に拡大した結果も併せて示している。
Fig.6(a)の LES 瞬間流場を見ると翼面上に多数の細長い渦構造が確認できるが、Fig.6(b)の平均流場ではこれらの微細な渦は現れておらず、渦が時間的に変動しており平均的には存在しないことがわかる。瞬間流場を見ると翼の前縁で変動の少ない層状の渦が形成され、後縁側へ進むに従って、渦層から流線に垂直な横渦が流れ方向に並ぶ Tollmiean - Schlichiting 波が形成される。Back 面側では翼中央付近で Face 面付近では(a) LES (800Mcells) Result

前縁に近い位置で渦の向きが変わり、後端側で流線に平行な縦渦に変化している。縦渦に変化した後はBack 面側では管状の渦構造が維持されず複雑化しているが、Face 面側では後端まで縦渦が維持され比較的単純な渦構造となっている。
これらの渦構造の変遷は平板境界層の発達過程 8)と概ね同様ながら、平板境界層で横渦が縦渦に遷移する過程において見られるヘアピン渦が存在しない。
Fig.6 (a)に示した A-A’の断面におけるプロペラ相対(b) RANS Result
図 5 計算による限界流線( ReD=1.1×106 )速度から求めたヘリシティ分布を Fig.7 に示す。ヘリシティは Back 面側では負値、Face 面側では正値にな

Fig. 5 Restricted Flow Patterns of Averaged Flowっており渦の回転方向がそれぞれの面で同方向に揃っ3.3 乱流エネルギーと表面摩擦の比較ていることがわかる。ヘアピン渦を介して縦渦が生成される場合、回転方向の異なる渦が作られるので、縦Fig.8 にプロペラ近傍での単位体積当たりの乱流運渦の回転方向が揃っていることはプロペラ翼面上での動エネルギーk の等値面を示す。乱流運動エネルギー

境界層の発達が通常の平板境界層の発達と異なることの可視化には、RANS ではモデル変数として計算されを示している。片方の回転方向の渦だけが生成されるる値を用い、LES ではプロペラが 5 回転する時間の過程には回転によるコリオリ力が影響していると考え流場から算出した値を用いた。乱流運動エネルギーがられ、プロペラにおける乱流遷移を検討する上で重要増加する位置を示すため、いくつかの等値面を半透明であると思われる。にして重ねて示している。
流入場は一様流であるためプロペラ翼前縁において乱流エネルギーは 0 であり、翼面の途中で急速に増大し翼面を通過後は減衰する。Fig8(a)で乱流エネルギーの増大する領域は、Fig.6(a)で渦構造が複雑になっている領域とも対応している。RANS 計算では LES とA
A’(a)LES(instantaneous flow)AA’(b)LES(averaged flow)(c)RANS

図 6 第二不変量 Q による渦の可視化
Fig.6 Iso-surfaces of Q/N2 (=8000)-5-JMU テクニカル・レビュー No.5

（2019 年 1 月）

BackFace

Back Side
AA’Face Side(a) LES Result

図 7 図 6(a) A-A’断面におけるヘリシティ分布
Fig.7 Distributions of Helicity on A-A’ Slice on Fig.6(a).
比較して、翼先端に近いほど前縁に近い位置で乱流遷移しており、限界流線と同様に乱流に遷移しやすい傾向がわかる。

4.

結 言

本研究では、プロペラ周りの流れについてスーパーコンピュータ「京」を用いた大規模 LES 計算を実施(b) RANS Result
図 8 乱流エネルギの可視化し、得られた流場の解析から以下の知見を得た。
・大規模 LES では従来の実験的な計測や RANS

Fig. 8 Iso-surfaces of k / ρN2D2 =0.01, 0.02, 0.04, 0.06, 0.08.

による CFD では解像できていなかったプロペラ翼面上の微細な渦構造を解像することが可能であ学会論文集 Vol.27 ,pp1-8,2018

り、プロペラ周りの流場を概ね再現できている。

3) 畝村 毅, 張 会来: 次世代流体解析ソフトウェア

・翼面上の渦構造の変化と限界流線が対応しており、FrontFlow/Red の開発, 生産研究 56 巻, pp40-43, 2004.

横渦が縦渦に変化する付近で限界流線は大きく半4) M.Germano, U.Piomelli, P.Moin, W.H.Cabot: A dynamic

径方向を向く。

subgrid-scale eddy viscosity model, Physics of Fluids A3-7,
pp.1760-1765, 1991.

・翼面上での境界層の発達過程では横渦から縦渦に変化する際にヘアピン渦を介しておらず、平板上

5) A.J.Chorin: A Numerical Method for Solvingでの境界層の発達過程と異なっている。

Incompressible Viscous Flow Problems, J.Comp.Physics
Vol. 2, pp12-26, 1967.

LES 計算の限界流線がより低いレイノルズ数条件の実験と一致する傾向を持つことの原因として、解像

6) R.K.Chan: A Generalized Arbitrary Lagrangian-Eulerian

度が十分でないことが考えられるため、今後さらに格Method for Incompressible Flows with Sharp Interfaces,
J.Comp.Physics Vol.17, pp311-331, 1975.

子解像度を高めた DNS 相当の精度を持つ計算を実施し検討する予定である。このような検討を通じて

7) 武井 幸男, 角川 明, 右近 良孝: 舶用プロペラ翼

LES 大規模計算が乱流遷移を含む流場を正確に捉え面上の流れの可視化, 流れの可視化 Vol.7 ,No.26, 1987,

られることを明らかにできれば、将来的には実船尺度pp.141-144において船体・付加物を含めた船舶の性能を正確に推8) 小橋 安次郎: 境界層遷移の機構と乱流の構造, 水定することが可能になると期待される。

工学に関する夏期研修会講義集 Vol.17, pp.A.7.1-A.7.16,
1981.

参考文献
1) 沢田 龍作: 数値流体力学の工業利用における標準解の整備に関わる研究, HPCI システム利用研究課題
利用報告書, hp140117, 2015
2) 藤澤 竹春, 坪倉 誠, 田中 寿夫: 大規模 LES 計算
藤澤 竹春

による舶用プロペラ周りの流場推定,日本船舶海洋工-6-JMU テクニカル・レビュー No.5

（2019 年 1 月）