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

熱流体格子ボルツマン法を用いた乱流騒音の数値解析
山田和豊，草野和也，木戸康太，古川雅人
九州大学大学院 工学研究院 機械工学部門

先駆的科学計算に関するフォーラム2014

目次
1.
2.
3.
4.
5.

研究背景・目的
格子ボルツマン法（LBM）
非熱流体モデルと熱流体モデルの比較
単独翼から発生する乱流音
まとめ1

先駆的科学計算に関するフォーラム2014

研究背景

2

エアコン室外機用プロペラファン

電子機器冷却用小型軸流ファン

Jang et al. [2001]

高効率化＋空力騒音の低減化の要求•低マッハ数流れ（ Ma ≤ 0.2 ）•小型（Re = 10⁵ 程度）•翼周りの複雑な流れ（はく離・遷移・大規模な渦構造）
空力性能・空力騒音を支配

自動車HVACシステム用多翼ファン

先駆的科学計算に関するフォーラム2014

研究背景

3

低速ファンの空力騒音の低減化•低速ファンの翼周りの複雑なはく離渦流れ場の数値解析•低速ファンから発生する乱流音の高精度な予測
高周波数音について音響解析（反射，散乱，放射，透過）
Reflecting, Scattering

Transition

Laminarseparation

Turbulent boundary layer

Eddies先駆的科学計算に関するフォーラム2014

研究の背景および目的

4

低マッハ数流れの空力音の数値解析手法•非圧縮性流れ解析＋Lighthillの音響アナロジー（Kato, et al., 2007）•非圧縮性解法の拡張（稲垣他，2000）•前処理付圧縮性解法＋全速度スキームSLAU（嶋，2008）•格子ボルツマン法（Tsutahara, 2012）
低速ファンの乱流音の予測・低減化の課題
乱流のLES/DNS等の大規模計算が必要
低マッハ数流れの圧縮性解析の困難さ•計算規模が膨大•計算精度・収束性の悪化
研究の目的
低速ファンから発生する乱流音の高精度予測の実現を目指して，低マッハ数の流れ場と音響場の直接計算に対する格子ボルツマン法の有効性を検証する

先駆的科学計算に関するフォーラム2014

格子ボルツマン法（LBM）
乱流解析・空力音解析・・・高次精度の数値解法が必要
計算機コア数の増大 CFDにおける並列化技術の重要性
高精度化と並列化の両立の課題
格子ボルツマン法
長所
アルゴリズムが単純
計算が高速並列化に優れている

短所
直交等間隔格子を用いる制約
高Re数流れの数値不安定性

有限差分格子ボルツマン法（Tsutahara, 2012）
空力音の直接解析（高精度に計算）
LBM本来の単純なアルゴリズム構造が失われる

5

先駆的科学計算に関するフォーラム2014

格子ボルツマン法（LBM）

6

格子ボルツマン法（Lattice Boltzmann Method）

Start
Initial condition•流体粒子の速度分布関数の時間発展を計算•計算アルゴリズムが単純高速な計算が可能，並列化に優れている

Collision
Boundary condition
Streaming

基礎方程式：格子BGK方程式
f i (t   t , x  c i t )  f i (t , x ) マクロな変数1

a[ f i (t , x)  f i eq (t , x)]1   f i ,  u   ci f i , p   c 23ii

D2Q993fi (t , x)   wi 1  3ci  u  (ci  u) 2  u  u  wi 22

eq

D3Q15
8

75109

粒子速度モデル
D2Q9, D3Q15
局所平衡分布関数

Continue ?
NoEndYes

23

14z

9 / 4 (i  0)
1/ 9 (i  1 ~ 4)
1/ 36 (i  5 ~ 8)y12x 1311614

先駆的科学計算に関するフォーラム2014

格子ボルツマン法（LBM）

7

非熱流体モデル•質量保存則，運動量保存則を満足する•エネルギー保存則を考慮しない•計算負荷が小さい
f i (t   t , x  c i t )  f i (t , x ) 1a

   f i ,  u   ci f ii

[ f i (t , x)  f i eq (t , x)]i

1p   c2
3

熱流体モデル•質量保存則，運動量保存則，エネルギー保存則を満足する•計算負荷が大きい
f i (t   t , x  c i  t )  f i ( t , x ) 1

a[ f i (t , x)  fi eq (t , x)]

1 1
g i (t  t , x  ci t )  g i (t , x)  [ g i (t , x)  g i (t , x)]     Z i [ f i (t , x)  f i eq (t , x)]bb a    f i ,  u   ci f i ,  E   g i
Z i  ci  u  u 2 / 21ip   RTiieq

先駆的科学計算に関するフォーラム2014

格子ボルツマン法（LBM）
マルチスケールモデル•局所的に格子解像度の変更を可能にするモデル•分布関数のスケール変換則  fˆ  f
fˆi f  fi eq,c  fˆi c  fi eq,cfˆi c  fi eq,ffieq, fi c 1
x c
n f, n f 2
 1x8

先駆的科学計算に関するフォーラム2014

格子ボルツマン法（LBM）9Building-Cube Method (Nakahashi, 2004)•計算領域をCubeに分割•Cube内は等方等間隔直交格子
Cubeの生成手順
1.
2.
3.
4.
5.

初期設定（物体形状，計算領域，最小Cube幅，Cube内の格子点数）
Cubeの分割
スムージング
各格子点の内外判定（物体表面の検出）
Cubeのグループ化

物体表面の検出法
基準線

基準線に垂直な面

100×100×100
任意の格子点

物体断面形状

基準線上の点
基準点

500×500×500

先駆的科学計算に関するフォーラム2014

格子ボルツマン法（LBM）

10

壁面境界条件
Interpolated Bounce-Backcix2

Re  100 U c  0.1

1.5

cix1xS2qi x 2qi Wallqi
ci ci(x  2q)
 2q
f(t,x)fi (t, x2 ),  q  x / 2 x i 1x
fi (t, xs )   1 fi (t, x1 )  (2q x) fi (t, x1 ),  q  x / 2
2q
 2q

球まわりの流れqix2x1xS

2qi x Wall

球表面の圧力係数の分布

Multi-scale LBM

LBM(bounce back)
LBM(interpolated bouce back)

1

Kim
0.5

Cp0

N-S(Johnson)‐0.5‐10306090

 [deg]120150180

先駆的科学計算に関するフォーラム2014

非熱流体モデルおよび熱流体モデルの比較

11Spectral Density [dB/Hz]単独翼から発生する離散周波数音 音圧スペクトル(Arcondoulis, et. al, 2009）
計算条件
1422 [Hz]Re  1.0  105 , Ma  0.2,   0[deg]計算領域：40C×40C
音波の吸収領域：半径20Cの円の外側
1009080706050403020100100100010000

f [Hz]

Paterson’s formula0.011U 1.5f
C

Cube内の格子数

Cube数 最小格子幅

17×17199000.00012C

総格子点数580万先駆的科学計算に関するフォーラム2014

非熱流体モデルおよび熱流体モデルの比較

音圧スペクトル(Arcondoulis, et. al, 2009）
Spectral Density [dB/Hz]

マッハ数によるトーン周波数の変化
LBM (Isothermal)
LBM (Thermal)
EXP (Arcondoulis)
1012

Paterson's formulaSt

100901422 [Hz]80706050403020100100100010000

f [Hz]
10.04

Ma0.40Paterson’s formula

両結果ともトーン周波数が予測できている

0.011U 1.5f
C先駆的科学計算に関するフォーラム2014

非熱流体モデルおよび熱流体モデルの比較

13

音圧レベルの指向性分布（後縁から垂直方向に10Cの位置）
90

180

9080 [dB]706050403020270

90

110[dB]100
90800

Isothermal
Thermal70

180060

Isothermal
Thermal270

熱流体モデルおよび非熱流体モデルの結果に有意な差は認められない

先駆的科学計算に関するフォーラム2014

単独翼から発生する乱流音

14
音圧スペクトル(Kato, et. al, 2007）
Spectral Density [dB/Hz]

計算条件
NACA0012翼型
Re  2.0  105 , Ma  0.058,   9.5[deg]50403020100

乱流遷移・乱流境界層のDNS
広帯域スペクトルにおいて高周波数域(数kHz以上）まで予測‐101001000f [Hz]10000

スパン方向長さ： 0.125CDamping region

16C

24C

Cube内の格子数

Cube数

最小格子幅

総格子点数

17×17×1724万y   2.312億

先駆的科学計算に関するフォーラム2014

単独翼から発生する乱流音

15

翼周りの渦度分布

三次元渦構造（Q値の等値面）
翼前縁近傍における時間平均の流線および流れ方向速度分布

1.0Hn1.0先駆的科学計算に関するフォーラム2014

単独翼から発生する乱流音

16

翼面上の静圧係数の分布
時間平均

分散値

‐4.0

0.5

‐3.5

EXP (Miyazawa)‐3.0LBM

EXP (Miyazawa)
0.4

LBM

‐2.5
‐2.0C pC p ‐1.5
‐1.0

0.3
0.2

‐0.5
0.0

0.1

0.5
1.0

0.00.0

0.1

0.2

0.3

0.4

0.5

x/C

0.6

0.7

0.8

0.91.00.0

0.1

0.2

0.3

0.4

0.5

x/C

0.6

0.7

0.8

0.9

1.0

先駆的科学計算に関するフォーラム2014

単独翼から発生する乱流音
後流の速度分布0.00.10.2

 /U
urms

0.3

0.4

0.5

17

0.6

0.7

1.0

1.2

0.3

EXP (Miyazawa)
0.2

LBM

0.1

y/C

S.S.
RMS

Mean

0.0

P.S.
‐0.1
‐0.2
‐0.2

0.0

0.2

0.4

0.6

umean / U 0.8

先駆的科学計算に関するフォーラム2014

単独翼から発生する乱流音

18

翼面上の圧力変動のパワースペクトルx / C  0.90

1.0E‐01

1.0E‐01

1.0E‐02

1.0E‐02

1.0E‐03

1.0E‐03

1.0E‐04

1.0E‐04PSD of Cp

PSD of Cp

x / C  0.075

1.0E‐05
1.0E‐06

1.0E‐06

EXP (Miyazawa)

1.0E‐08

LBM (x/c=0.075)

1.0E‐08

1.0E‐09

LBM (x/c=0.100)

1.0E‐09

1.0E‐07

1.0E‐10110x / C  0.25

100

St

1.0E‐02

EXP (Miyazawa)

Slope ‐5/3

LBM

1.0E‐03
1.0E‐04
1.0E‐05
1.0E‐06
1.0E‐07
1.0E‐08
1.0E‐09
1.0E‐10110100St110100St

1.0E‐01

PSD of Cp

LBM

1.0E‐05

1.0E‐07

1.0E‐10

EXP (Miyazawa)

Slope ‐5/3

翼面上の圧力変動分布先駆的科学計算に関するフォーラム2014

単独翼から発生する乱流音

19

圧力変動量分布音圧スペクトル（観測点：後縁から垂直方向に6.7Cの位置）
Spectral Density [dB/Hz]

50

EXP (Kato)
40

LBM
LBM (correction)

30
20
100Airfoil‐10100

1000

10000

f [Hz]

700Hzより高周波数側は定量的に良い一致低周波数域が若干大きく評価

3.5 104p3.5 104

先駆的科学計算に関するフォーラム2014

単独翼から発生する乱流音

20

スパン長さに関する音圧レベルの補正•周期境界を使用・・・スパン長無限大•スパン長無限大→解析領域のスパン長（0.125C）に補正•球面波のエネルギーは音源からの距離の2乗に逆比例
微小長さ x から観測点に放射される
音響エネルギー I  I xr2
無限長から観測点に放射される
音響エネルギーLII   lim dxI L R2  x2RL

YObservation point

Point sound source row R(Uniform acoustic energy I )LoLsr  R2  x2xxL

解析領域のスパン長 Lc から観測点に放射される音響エネルギー
LII c  c2
R補正式

I 
LSPLc  SPL  10log10     SPL  10log10 c
R Ic SPL ：補正前の音圧レベル
SPLc ：補正後の音圧レベル

X

先駆的科学計算に関するフォーラム2014

まとめ21

低マッハ数の流れ場と音響場の直接計算に対する格子ボルツマン法の有効性を示した
N-Sおよび熱流体LBMに対して，非熱流体LBMの計算速度における優位性を示すとともに，低マッハ数流れから発生する空力音の計算おいて同等の計算精度が得られることを示した
乱流音の広帯域スペクトルが高精度に予測できることを示した
今後の展開実用問題として，低速プロペラファンの空力騒音解析に適用