格子ボルツマン法を用いたファンの空力騒音の直接解析
山田和豊（九大），岡 信仁（九大院），岡 明宏（九大院），古川雅人（九大）
緒

1.

q言

空力騒音の問題で流れが低マッハ数である場合は多く，空力騒音解析手法には音響アナロジーを用いた分離解法がよく用いられる．しかし，コンパクトの仮定を用いる分離解法では，乱流音のような高周波数の空力音に関して，予測精度が十分でない．これについては流体音の直接計算が望ましいが，高レイノルズ数かつ低マッハ数流れでは，計算規模および計算精度の面で実現が難しい．

x2x2

x1x02q − 1 Wall(b) q ≥ 1/2

Fig. 1 Illustration of interpolated Bounce-Back scheme

2.2

マルチスケールモデル

上述のとおり，直交等間隔格子を使用するため，物体近傍での計算精度低下および計算格子数増加の問題が生じる．そこで，マルチスケールモデル (6) を導入し，解像度の異なる格子を組み合わせて計算を行うことで，物体近傍の格子解像度を局所的に高めた．
．粗い格子と細かい格子との境界では，以下の式により分布関数を変換後，受け渡しを行う．
′′

fi = F̃ieq + (F̃i − F̃ieq )/Ω′′
Fi = fieq + (fi − fieq )ΩΩ=nτc − 1
τf − 1(2)(3)(4)

ここで，大文字 F は粗い格子ブロック，小文字 f は細か′い格子ブロックでの分布関数を表し，上付き添え字 は衝突後の分布関数であることを表す．また，τc および τf
それぞれ粗い格子および細かい格子における緩和係数を表す．n はブロック間での格子幅の比 ∆xc /∆xf を表し，
本研究では n = 2 とした．なお，解像度の違いから対応する格子点が存在しない場合は，周辺の格子点の分布関数を用いて補間する．

LBM では，本来，粒子速度モデルにあわせて等間隔規則格子を用いる必要がある．そのため，境界適合格子を使用する有限差分格子ボルツマン法（FDLBM）や有限体積格子ボルツマン法（FVLBM）が提案されているが，これらの手法では LBM のメリットである単純なアルゴリズム構造が失われる．そこで本研究では，従来どおり直交等間隔格子を使用し，物体境界は Immersed boundary法により取り扱った．また，格子点数の増加を抑えるために，マルチスケールモデルを使用した．

2.3

境界条件

物体表面の境界条件として，Interpolated BounceBack(7) を適用した．Interpolated Bounce-Back は比較的単純なスキームで任意物体形状を通常の Bounce-Backスキームよりも高精度に表現できる．図 1 に壁面境界条件の説明を示す．x0 は物体表面に最も近い物体内部の格子点であり，x1 は流体内の物体表面に最近傍の格子点，x2 はその二番目の格子点である．q は壁面と点 x1 との距離を表す．次のステップで壁面の方向から x1 点に到達する粒子は，壁面での跳ね返りを考慮すると現時刻では赤い丸印の位置にある．その位置の分布関数を補間し，その分布関数を境界条件として x0 に与える．したがって，x0 には次式で計算される分布関数が与えられる．

格子ボルツマン法

本計算では，次式の格子 BGK モデルを用いた．

fi (x + ci ∆t, t + ∆t)= fi (x, t) −x0(a) q < 1/2

数値計算手法2.1

x12q 1 − 2q Wall

格子ボルツマン法 (LBM: Lattice Boltzmann Method)
は，Navier-Stokes 方程式を離散化して解く従来の計算手法とは異なる新しい流れの計算手法である (1) ．LBM は，計算アルゴリズムが単純なことから，従来の計算手法に比べて計算負荷が小さく高速に計算が行えるともに，並列計算に適している．また極めて低いマッハ数の流れを高精度に計算できるといわれる．実際，Navier-Stokes 解
析では空間 6 次精度・時間 4 次精度のスキームが必要な流体音の直接計算について，LBM では空間 3 次精度・時間 2 次精度でも音圧に関して同程度の結果が得られている (2) ．したがって，LBM は低マッハ数流れの空力騒音の予測に適した計算手法であるといえ，近年空力騒音解析の有力な手法として応用されつつある (3)(4) ．そこで本
研究では，まず空力音の検証として，単独翼から発生する乱流音について LBM の予測精度を調査した．次に，半開放型プロペラファンまわりの流れについて，LBM による大規模直接計算を実施し，実験結果との比較を行った．

2.

q]
1[fi (x, t) − fieq (x, t) (1)τ

ここで，f は仮想粒子の分布関数，c は粒子速度，τ は単一時間緩和係数，f eq は局所平衡分布関数を表す．また，添え字 i は離散化された粒子速度方向に対応した整数である．本研究では，粒子速度モデルに非熱流体モデルの3 次元 15 速度（D3Q15）モデルを採用した．非熱流体モ
デルは，熱流体モデルに比べ粒子数が大幅に少なく，使用メモリや演算量が削減できるため大規模計算に適する．
また，低マッハ数の場合には，本モデルであっても十分な精度で空力音を予測できることが報告されている (5) ．

12
12q − 11fi (x0 , t) =fi (x1 , t) +fi (x1 , t), q ≥2q2q
2fi (x0 , t) = 2qfi (x1 , t) + (1 − 2q)fi (x2 , t), q <ここで，i は i と反対向きの速度を持つ粒子を指す．

1(5)

′
urms
0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.6

0.8

1.0

1.2

0.3
EXP (Miyazawa)
0.2

LBM

0.1

yy/Cyxx0.0

Measurementline
line
Measurement-0.1

Fig. 2 Computational mesh for NACA0012 (block linesare only shown)

-0.2
-0.2

3.

単独翼

0.00.20.4

umean

Fig. 3 Time-averaged velocity and RMS velocity distributions at the downstream of airfoil

空力音の検証として，NACA0012 翼から発生する広帯域空力騒音について解析を実施した．流れ場および空力音の評価においては，宮澤らの実験結果 (10) と比較した．
実験が行われた条件は，レイノルズ数 Re = 2.0 × 105 ，迎え角 9.5◦ ，主流流速 20m/s，コード長 C = 0.15m である．この条件では，流れは翼負圧面側前縁付近で一端はく離し，すぐに乱流へと遷移して再付着する．計算では，領域全体で温度が 20◦ C で一様であると仮定し，音速一定として主流マッハ数 0.058 で計算を行った．
計算領域は，流れ方向，垂直方向，スパン方向の長さがそれぞれ 24C ，24C ，0.125C の直方体領域とした．図2 に翼周りの計算格子を示す．本研究では，計算領域の分割に Building-Cube Method（BCM）(8) を用いて，計算格子を作成した．同図にはその計算ブロックが示されている．各ブロックには，17 × 17 × 17 の格子点を設定した．翼表面近傍の最小格子幅は壁座標で 2.3 である．主流部の最大格子幅は 0.0078C であり，これは 15kHz の音波についてその波長を 20 格子点で解像することができる．総格子点数は約 12 億点である．
スパン方向の境界条件には周期条件を用いた．計算領域の境界は翼から十分に離れていると仮定し，遠方境界条件として一様な密度および速度を設定してそれらのマクロな物理量から計算される局所分布関数を固定して与えた．このとき，計算領域外縁での音波の反射を防ぐため，計算領域の外側には Damping region を設け，音波を減衰させた．Damping region では，式 (1) に変えて次式を用いた (9) ．

Observation point

6.7C

3.5 ×10−4ps′

0.5 ρ∞U ∞ 2−3.5 ×10−4Fig. 4 Instantaneous pressure ﬂuctuation ﬁeld

速度をみると，実験結果には主流中でもわずかながら速度変動がみられ，主流乱れの存在を確認できる．LBM では，主流乱れは導入していないため，主流の RMS 速度はゼロとなっている．後流中について実験と比較すると，速度変動は過大評価されている．翼上面側の境界層は，後縁ではかなり厚くなっており，そこでは境界層外縁付近の格子解像度が十分でないことが考えられる．
図 4 に等スパン断面における瞬時の圧力変動分布を示す．前縁および後縁付近から不規則な圧力波が伝播している様子がみられる．前縁付近から発生している圧力波は，前縁近傍で生じた剥離の非定常的な変動による．前縁剥離したせん断層は周期的に巻き上がり三次元的に崩壊することで再付着するが，その過程で微細な乱流渦を生成するとともに前縁近傍の再付着点付近に圧力変動を発生させている．また後縁では，翼上面の発達した乱流境界層内の渦が干渉して騒音を引き起こしている．
図 5 は，翼後縁から流れ方向に対して垂直上方に 6.7Cだけ離れた位置（図 4 の観測点）における空力音のスペクトルに関して，実験結果 (11) との比較を示している．同図において，赤色の線はスパン長に関する音圧レベル補1
fi (t+∆t, x+ci ∆t) = fi (t, x)− [fi (t, x) − fieq (t, x)]
τ( )2[ eq]xd−αfi (x, t) − fitarget (x, t) (6)Ld
ここで，f target は分布関数の目標値であり，主流の密度および速度から計算される平衡分布関数を与えた．Ld はDamping region 入口の位置から計算領域外縁までの距離であり，xd は Damping region 入口からの距離である．
係数 α は 0.05 とした．本計算では，計算領域中心から直径 16C の円の外側を Damping region とした．
図 3 は，翼下流（図 2 の赤破線）における時間平均速度分布および RMS 速度分布の比較を示している．時間平均速度分布をみると，LBM 結果は後流の幅および後流の深さともに実験結果と非常によく一致している．RMS

2

SpectralSpectral Density
Density [dB/Hz][dB/Hz]50

Inlet
ChamberExp.(Kato)(Kato)
EXP

40

LBM
LBM

Inflow

LBM(correction)(correction)
LBM

30

OutletChamber
OutflowRotor

7.8 Dt20104 Dt

012 Dt(a) Overview (Cubes in the whole domain)-10100100010000

f [Hz]

Fig. 5 Frequency spectra of far-ﬁeld soundShroudRotor

ShroudSSTEPSSS

Tip

LEFlowRotation(b) Close-up view around fan
Fig. 7 Computational mesh for propeller fan
ここで，Q は体積流量，Utip は動翼先端の周速度，∆Psは羽根車前後の静圧差である．本研究では，設計流量点の流れ場を解析した．動翼先端半径および動翼先端周速に基づくレイノルズ数は 2.6 × 105 であり，動翼先端周速のマッハ数は 0.05 である．
計算格子を図 7 に示す．図 7(a) は計算領域全体のブロック分布を示している．実験と同様に羽根車前後に入口チャンバーおよび出口チャンバーを設置し，入口境界および出口境界を同図に示すとおり設定した．同図において，赤色の領域は羽根車とともに回転する相対座標系で計算し，青色の領域は静止座標系で計算を行った．なお，本計算では両計算領域の周方向相対位置の時間変化は考慮せず，フローズンローターとした．図 7(b) には羽根車まわりの計算格子を示している．ただし，各方向 7 点おきに格子を間引いて表示している．羽根車およびシュラウドの壁面近傍に細かい格子が配置されていることがわかる．翼先端半径で無次元化した翼面近傍の最小格子
幅は 9.8 × 10−4 (∆xmin = 12) である．ブロックの総数は46,580 個であり，各ブロックは 333 の格子点を有している．総格子点数は約 16.7 億点である．
図 8 に，ファンの内部流れ場について，時間平均および瞬時の渦流れ構造を示す．同図 (a) に示した時間平均流れ場では，翼端渦をはっきりと確認でき，ファン下流まで発達している様子がわかる．また，ハブ側後縁付近にはコーナー剥離がみられる．一方，瞬時の流れ場では，同図 (b) に示されるように，翼端渦は非定常的に変動しているため羽根車下流では明確ではない．ハブ面の上流側角部で剥離した流れが，翼間内の二次流れによって運ばれ，ファン動翼負圧面上を駆け上がっている．同図 (c)に翼先端前縁部の拡大図を示している．同図に示される
ように，翼先端部では前縁剥離が発生しており，これにFig. 6 Half-ducted propeller fan

正前の計算結果であり，緑色の線は補正後の結果である．
補正前の計算結果は，実験結果よりも全体的大きく評価されている．補正後の結果は，700Hz 以下の低周波数域
では，実験結果に比べて過大評価されているが，700Hz から 2kHz の高周波数の領域については，LBM 結果は実験結果と良い一致が得られており，空力騒音を定量的に予測できていることがわかる．本計算のスパン長は 0.125Cと実験のスパン長 3.3C に比べて短いため，スパン方向の大きな流れ構造を再現できていないことが，低周波数域の差異の原因としてが考えられる．

4.

プロペラファン

解析対象は図 6 に示すエアコン室外機用の半開放形プロペラファンである．翼枚数は 3 枚であり，動翼先端の半径 rtip は 235mm, 回転数は 683rpm である．ハブ比 ν は0.3404, 動翼先端とシュラウド間の隙間は 7.76mm（翼先端半径の 3.3%）である．翼先端のコード長は 353mm である．シュラウドの軸方向長さは動翼先端の軸方向コード長の 30%程度であり，シュラウド後端部が動翼後縁の軸方向位置と一致するように取り付けてある．したがって，動翼先端の前方側 70%は上流に開放されている．シュラウドの前縁部はベルマウス形状となっている．設計点の流量係数 ϕ および静圧上昇係数 ψ はそれぞれ 0.291，0.130 である．なお，ϕ および ψ は次式で定義される．

ϕ=Q,π(1 − ν 2 )rtip 2 Utip

ψ=∆Ps
ρUtip 2 /2(7)3Tip Vortex
TIPTip Vortex

Rotation

Corner
Separation
VortexLE Separation
VortexLE1.0Hn−1.0(a) time-average(b) instant(c) instant (Close-up)

Fig. 8 Vortical ﬂow structures inside fan
0.58590.36970.3732cz / U tip

cr / U tipcθ / U tip−0.2156−0.4274

0.02490.58590.36970.3732cz / U tip

cr / U tipcθ / U tip−0.2156−0.4274

0.0249(a) Experiment(b) LBM
Fig. 9 Comparisons of ﬂow ﬁeld at fan exit

伴う渦構造を確認できる．この前縁剥離は小規模ですぐに再付着している．翼端渦のすぐ下では多数の縦渦構造が非定常的に発生し移流する様子が確認される．
図 9 に羽根車出口における流れ場の比較を示す．図中の白丸は翼端渦の位置を示している．翼端渦の周辺には，その旋回流れに伴い，渦中心を対称とした速度場が形成される．LBM 結果をみると，翼端渦の位置およびその周辺の速度場は実験とよく一致していることがわかる．また，破線で囲んだ領域はコーナー剥離を表しているが，このコーナー剥離についても LBM はよく再現していることがわかる．しかし，旋回速度分布において，LBM 結果には翼端渦と翼圧力面との間に低旋回速度域がみられるのに対し，実験結果にはそのような分布は確認されない．
この点を除けば，LBM 結果は実験結果とよく一致しているといえる．
図 10 は羽根車出口における乱れ度分布の比較を示している．同図に示される乱れ度は羽根車の先端周速に基づいて定義されている．LBM 結果は実験結果と定性的に一致しているが，定量的には乱れ度が過大に評価されている．本計算では乱流モデルを用いていないが，翼端渦やコーナー剥離渦等の周辺では壁面に比べて格子解像度がかなり低く，格子解像度が不足したことがこの要因として考えられる．上述したように，LBM 結果では図 9(b)に示す旋回速度分布において，翼端渦と隣接翼圧力面との間に旋回速度の低い領域がみられるが，この低旋回速度域と同じ位置に高い乱れ度の領域を確認できる．この高乱れ度域は，LBM 結果だけでなく実験結果にも確認される．この周辺では実際に流れ場に非定常性が誘起されており，この非定常流動現象を LBM は再現している．
最後に，図 11 に翼負圧面上の圧力変動分布の比較を示す．同図 (a) の実験結果をみると，最も顕著な圧力変動は翼先端の前縁付近にみられる．これは，同図 (b) の LBM
解析結果においても同様の箇所に確認される．この周辺は翼端渦の巻き上がり位置と対応するが，翼端渦の非定常性に起因したものではない．この箇所は図 8(c) に示した前縁剥離に対応しており，この圧力変動は前縁剥離に伴う流れの非定常挙動に起因することがわかる．翼端側では全般にコード方向に圧力変動が高くなっており，前縁剥離によって放出された渦が下流に拡散しながら移流するためである．またハブ側には，コーナー剥離に起因

4

Tip vortex

Wake0.083

単独翼の計算では，後流速度分布に関する流れ場の比較において，LBM 結果は実験結果と良い一致を示した．
また，遠方場の空力音スペクトルにおいても，700Hz から2kHz の高周波数帯域で実験結果と良い一致が得られた．
プロペラファンの計算では，時間平均流れ場の比較において，計算結果は実験結果とよく一致した．流れの非定常性に関して，乱れ度分布および壁面圧力変動分布を比較した．計算結果は，定量的に過大評価したものの，定性的には実験結果とよく一致した．また，前縁剥離に起因した大きな圧力変動が生じていることが明らかとなった．

T .L.

SS

0.005

PSCornerseparation(a) Experiment

参考文献(1) Chen, S. and Doolen, G. D., “Lattice Boltzmann
Method for Fluid Flows”, Annu Rev Fluid Mech(1998), pp. 329–364

0.083

Tip vortex

T .L.

SS

0.005

PS

(2) Tsutahara, M., et al., “Direct Simulation of Acoustic Waves by the Finite Volume Lattice Boltzmann
Method”, Transactions of the Japan Society of Mechanical Engineering, Series B, Vol. 72, No. 724,(2006), pp. 2999–3006(b) LBM(3) Goﬀ, V. L., et al., “Flow-Induced Noise Predictionsof an Automotive Alternator Using a Lattice Boltzmann Method”, Proc. ASME Turbo Expo 2014,
GT2014-26386 (2014).

Fig. 10 Turbulence intensity distributions at rotor exit

Tip

0.08LE(4) Zhu, T., et al., “Experimental and Numerical Investigation of Tip Clearance Noise of an Axial Fan
Using a Lattice Boltzmann Method”, Proc. the21th International Congress on Sound and Vibration, (2014).

TEPf
0.02

(5) 田村, 蔦原, “格子ボルツマン非熱流体モデルの空力音シミュレーションへの適用性について”, 第 21 回
数値流体力学シンポジウム, A4-2 (2007).

Hub(a) Experiment(6) Filippova, O., et al., “Multiscale Lattice Boltzmann
Schemes with Turbulence Modeling”, Journal ofComputational Physics, Vol. 170, (2001), pp. 812–829.

Tip

0.08

LE(7) Bouzidi, M., et al., “Momentum transfer of a
Boltzmann-lattice ﬂuid with boundaries,” Phys.
Fluids, 13 (2001), pp.3452-3459.

TEPf(8) 石川, 他 2 名, “大規模並列計算に向けた Building
Cube Method の拡張に関する研究,” 第 24 回数値流体力学シンポジウム, E4-2.

Hub

0.00

(b) LBM(9) Kam, E. W. S., et al., “Lattice Boltzmann Method
Simulation of Aeroacoustics and Nonreﬂecting
Boundary Condition,” AIAA Journal, Vol.45, No.7(2007), pp.1703-1712.

Fig. 11 Pressure ﬂuctuation distributions on blade suction surface
した圧力変動が，領域の大きさに差異はあるものの，実験結果および LBM 結果ともにみられる．LBM 結果では，翼端の後縁部においても高い圧力変動がみられる．これはシュラウド内部で発生する翼端漏れ渦によるものと考えられるが，更なる調査が必要である．

(10) 宮澤, 他 3 名, “二次元翼周り非定常流の空力音響解析（第 1 報，翼周りはく離遷移流れの LES 解析と精度検証）,” 日本機械学会論文集 B 編, Vol. 72, No.
721 (2006), pp. 2140-2147.

結 言
本研究では，格子ボルツマン法を単独翼および半開放型プロペラファンまわりの流れ解析に適用した．10 億点超の計算格子を用いて大規模数値解析を実施し，流れ場および空力音に関して，実験結果と比較検証を行った．

5.

(11) Kato, C., et al., “Numerical prediction of soundgenerated from ﬂows with a low Mach number,”Computer and Fluids, 36 (2007), pp. 53-68.

5