BCM およびメッシュレス法を用いたハイブリッド LBM による空力騒音の直接解析
山田和豊（九大），文 吉周（九大院），岡 明宏（九大院），黒川 靖（九大院）

1.

緒

言

2.2

コンピュータに代表される電子機器や空調システムは居住環境で作動することから，それらに組み込まれた冷却用ファンや送風用ファンなどから放射される空力騒音が大きな問題となっている．このように，居住環境下で作動するファンには低騒音化が求められ，空力騒音の低減や予測に関する研究が精力的に行われている．ファンのように空力騒音が問題となる流れは低マッハ数である
場合が多く，空力騒音の予測手法として音響アナロジー
を用いた分離解法がよく用いられる．しかし，高レイノ
ルズ数になると乱流渦構造が十分に解像されないことに加え，高周波数の乱流騒音に対してコンパクトの仮定を用いる分離解法では，空力騒音の予測精度は十分でない．
流体音の直接計算が望ましいが，計算規模が大きく，また低マッハ数流れでは難しい．
格子ボルツマン法 (LBM: Lattice Boltzmann Method)
は，Navier-Stokes 方程式を離散化して解く従来の計算手法とは異なる，新しい流れの計算手法である (1) ．LBM
は，計算アルゴリズムが単純なことから，従来の計算手法に比べて計算負荷が小さく高速に計算が行え，また極めて低いマッハ数の流れを高精度に計算できる．流体音の直接計算において，Navier-Stokes 解析では空間 6 次精度・時間 4 次精度のスキームが必要といわれるが，LBM
では空間 3 次精度・時間 2 次精度でも音圧に関して同程度の結果が得られている (2) ．したがって，LBM はター
ボ機械で問題となるような低マッハ数流れの空力騒音の予測に適した計算手法であるといえ，その応用が期待される．そこで本研究では，LBM をターボ機械における空力音の直接解析に応用することを目的として，NACA 翼から発生する広帯域騒音の予測精度を実験結果と比較することにより検証した．

2.

数値計算手法2.1

格子ボルツマン法

格子 BGK モデル

本研究では，衝突則にすべての粒子が同じ割合で平衡状態へ向かうと仮定した格子 BGK(Bhatnagar-GrossKrook) モデルを使用した．このとき，格子ボルツマン方程式は次のように表される．

fi (x + ci ∆t, t + ∆t)= fi (x, t) −]1[
fi (x, t) − fieq (x, t) (1)τ

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

i

局所平衡分布関数 f eq は，Boltzmann 方程式の平衡解である Maxwell 分布を基にして，Taylor 展開により次式で近似される．
[
ci ·u(ci ·u)2
u·u ]fieq = ρwi 1 +
+
−
(3)
2
RT
2(RT )
2RT
ここで，R は気体定数であり，T は温度である．粒子速
度数および Taylor 展開の打ち切り次数が限られるため，
D3Q15 モデルでは RT = 1/3 の制約が課される．このよ
うに，本モデルでは温度は固定されるため，圧力は状態方√
程式から p = ρ/3 の関係で与えられ，音速は cs = c/ 3
で一定となる．最終的に，局所平衡分布関数は次式で与
えられる．
]
[
3
9
(4)
fieq = ρwi 1 + 3ci ·u + (ci ·u)2 − u·u
2
2

LBM は，流体を仮想粒子の集合体と考え，仮想粒子の衝突および並進から巨視的な流体運動を再現する数値計算手法である．LBM では，流体の等方性を満足するように仮想粒子の速度は定義され，格子を併進する仮想粒子にあわせて規則格子が使用される．仮想粒子の衝突則を適切に与えることで，LBM の基礎方程式は低マッハ数の条件下で Navier-Stokes 方程式に漸近する．ただし，仮想粒子速度の数が少ない場合，限られた粒子速度数の制約からエネルギー式が満足されない．このため，そのようなモデルは非熱流体モデルといわれる．しかしながら，
通常の非圧縮性解法とは異なり，LBM の基礎式は密度に対応する分布関数に関する方程式であるため，導出されるマクロな保存則としては圧縮性の連続の式となる．し
たがって，非熱流体モデルであっても，それを陽的に解く LBM は，極低マッハ数流れを取り扱う計算手法でありながら，流体音の計算に応用できる可能性を有する．

ここで，wi は重み係数であり，Navier-Stokes 方程式導出の条件から，D3Q15 モデルの場合，wi = 2/9 (i = 0)，
wi = 1/9 (i = 1 ∼ 6)，wi = 1/72 (i = 7 ∼ 14) と決定される．また，動粘性係数 ν と緩和係数 τ は，次式のように関係づけられる．

ν=11(τ − )c∆x
32

ここで，∆x は格子幅を表す．

1

(5)

q

q

Coarse gridx2

x1

2q 1 − 2q Wall(a) q < 1/2

x0

x2

x1

x0

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

Fine gridFig. 2 Illustration of interpolated Bounce-Back scheme
Interpolation

Scaling

Streaming

Fig. 1 Computational procedure for interface between
two blocks with diﬀerent lattice spacing in multiscale
model

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

2.3 マルチスケールモデル
流体音を直接計算する場合，流れ計算において最小渦
スケールまで解像する必要がある．そのため高レイノル
ズ数流れにおいて，等間隔規則格子を用いる格子ボルツ
マン法では，格子数が計算不可能なほどに膨大となる．そこで，計算領域を複数のブロックに分割し，各ブロックで独立に格子幅を設定して等間隔規則格子を生成する．これにより，壁面に向かって段階的にブロックの格子幅を小さくすれば，計算領域全体の格子数を削減できる．この場合，各ブロックで粒子速度の大きさが異なるため，ブロック境界で以下のような処理を行う（図 1）(4) ．
細かい格子ブロックにおける緩和係数 τf は，粗い格子ブロックにおける緩和係数 τc を用いて，ブロック境界で動粘性係数が一致するという条件から，次式のように決
定される．
τf =

1
1
+ n(τc − )
2
2

1
2
1
1 (10)
2q − 1fi (x0 , t) =fi (x1 , t) +fi (x1 , t), q ≥2q2q2fi (x0 , t) = 2qfi (x1 , t) + (1 − 2q)fi (x2 , t), q <ここで，i は i と反対向きの速度を持つ粒子を指す．
後述する単独翼の計算においては，スパン方向の境界
条件には周期条件を用いた．また，計算領域の境界は翼から十分に離れていると仮定し，遠方境界条件として一様な密度および速度を設定してそれらのマクロな物理量から計算される局所分布関数を固定して与えた．このとき，計算領域外縁での音波の反射を防ぐため，計算領域の外側には Damping region を設け，音波を減衰させた．
Damping region では，式 (1) に変えて次式を用いた (7) ．

(6)

ここで，n はブロック間での格子幅の比 ∆xc /∆xf を表し，本研究では n = 2 とした．また，ブロック境界で巨視的な物理量およびその導関数が連続であるという条件から，以下の関係が導出される．
′= F̃ieq + (F̃i − F̃ieq )/Ω′= fieq + (fi − fieq )Ω

fi
Fi

′

′

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) (11)
−α
Ld

(7)
(8)

ここで，大文字 F は粗い格子ブロック，小文字 f は細か′
い格子ブロックでの分布関数を表し，上付き添え字 は衝突後の分布関数であることを表す．また，Ω は次式で定義されるパラメータである．

Ω=nτc − 1
τf − 1

ここで，f target は分布関数の目標値であり，主流の密度および速度から計算される平衡分布関数を与えた．Ld はDamping region 入口の位置から計算領域外縁までの距離であり，xd は Damping region 入口からの距離である．
係数 α は 0.05 とした．

(9)

2.5 Building-Cube Method
マルチスケールモデルにおいて，各ブロックのサイズ
および各ブロック間での n を任意に設定できるが，並列計算を行う上では各ブロックの格子点数をそろえるのが望ましい．そこで，本研究では Building-Cube Method（BCM）(6) を用いて格子の細分化を行った．BCM では，計算領域が多数の”Cube”と呼ばれる同じ格子点数を持つ
直交等間隔格子のブロックに分割される．物体近傍など
格子解像度が必要な個所では，この Cube のサイズを局
所的かつ段階的に小さくすることで格子解像度を変更する．つまり，格子細分化によって Cube に含まれる格子
点数自体は変化しない．そのため，各 Cube の計算負荷は等しく，BCM は並列化に適した格子細分化を行える．

図 1 にマルチスケールモデルにおける計算手順について示す．粗い格子と細かい格子の境界においては同図に示すように格子点を重複させる．まず，衝突演算後，粗い格子において分布関数の空間補間を行う．次に，式 (7),(8) を用いて，分布関数について格子スケールの変換を行う．その後，移流演算を行い，1 時間ステップが完了する．また，粗い格子において時刻 t + 1/2∆tc の分布関数
は，時刻 t および t + ∆tc の分布関数から補間して得られる．

2.4 境界条件
物体表面の境界条件として，Interpolated BounceBack(5) を適用した．Interpolated Bounce-Back は比較
2

(×10−7 )Tab. 1 Computational time ratio of each scheme to theisothermal LBM

10.0
N-S (2nd-oder central)

8.0

p′

N-S (4th-order central)

LBM

6.0

N-S (6th-order compact)

4.0

LBM (thermal)
LBM (isothermal)

2.0

Navier-Stokes

Isothermal

ThermalC2-RK4C4-RK4

CP6-RK4

1.0

2.5

6.0

12.3

16.3

0.0
-2.0

図 3 に時刻 t = 0.0176s のときの圧力の変動量分布をそれぞれの格子について示す．まず，六次精度のコンパクト
スキームの結果を両格子間で比較すると，上下のピークの位置や大きさはほとんど変わらないことがわかり，比較的粗い G201 の格子でも十分な精度で計算できている．
なお，同スキームの結果においてシンボルは計算格子点
における値を示しており，G201 では上下のピーク間は約
7 格子点で計算されている．次に二次精度中心差分の結果に着目すると，他の計算結果と比べて，G201 では振幅および位相のずれが大きいことがわかる．G401 では
G201 に比べて改善されているが，六次精度のコンパクト
スキームの結果との差異が見られる．四次精度中心差分
および熱流体 LBM，非熱流体 LBM の結果を比べると，ほぼ等しい結果であることがわかる．さらに，G201 を用いた四次精度中心差分および LBM の結果は，G401 の二次精度中心差分の結果と同等であることがわかる．G401の格子では，四次精度中心差分および LBM の結果は六次精度のコンパクトスキームの結果とほぼ一致している．
以上より，LBM の計算結果は時間および空間二次精度のNavier-Stokes 解よりも低散逸な特性を示し，四次精度のNavier-Stokes 計算と同程度の結果を与えることがわかった．また，本計算では非熱流体モデルと熱流体モデルの計算結果に違いは認められない．
次に，計算時間の比較を示す．表 1 に各計算スキーム
で要した計算時間を非熱流体 LBM に対する比率により示す．本計算は 1CPU による計算であり，並列化は行っていない．熱流体 LBM は非熱流体 LBM よりも 2.5 倍の計算時間を要する．LBM と同精度の離散化スキーム
である Navier-Stokes 計算の空間二次精度の中心差分および時間二次精度の Runge-Kutta 法は 6.0 倍の計算時間を要する．一方，本計算結果において，LBM と同等の結果が得られた四次精度の中心差分および時間四次精度のRunge-Kutta 法では約 12 倍の計算時間がかかる．さらに，空力音の直接計算に広く用いられる六次精度のコン
パクトスキームおよび四次精度の Runge-Kutta 法では約16 倍の時間を要する．以上により，計算速度において，Navier-Stokes 計算 S および熱流体 LBM に対する非熱流体 LBM の優位性が示された．

-4.0
-6.0
-8.0
4.0

5.0

6.0

x [m]

7.0

8.0

(a) G201 (coarse mesh)(×10 −7 )10.0
N-S (2nd-oder central)

8.0

p′

N-S (4th-order central)

6.0

N-S (6th-order compact)

4.0

LBM (thermal)
LBM (isothermal)

2.0
0.0
-2.0
-4.0
-6.0
-8.0
4.0

5.0

6.0

x [m]

7.0

8.0

(b) G401 (ﬁne mesh)Fig. 3 Instantaneous distributions of pressure ﬂuctuation on the x-axis

3.

パルス伝播問題による検証

LBM と差分法による圧縮性 Navier-Stokes 解の予測精度および計算速度の比較のために，一様流中を伝播するパルスの計算を行った．LBM では，非熱流体モデルに加えて，熱流体モデル (8) による計算についても実施した．
一方，Navier-Stokes 計算では，空間の離散化には二次精度および四次精度の中心差分スキーム，六次精度のコン
パクトスキーム (9) の 3 種類のスキームについて比較調査した．時間積分には，空間二次精度の中心差分に対しては二次精度の Runge-Kutta 法，空間四次精度の中心差
分および六次精度のコンパクトスキームに対しては四次精度の Runge-Kutta 法をそれぞれ用いた．
一様流のマッハ数は 0.1 とし，初期条件として密度の変動量を次式で与えた．
[]
ρ′ = ϵ exp −(ln 2)(x2 + y 2 )/b
(12)

4.

単独翼から発生する広帯域騒音の予測

4.1 解析対象および計算格子
空力音の検証として，NACA0012 翼から発生する広帯域空力騒音について解析を実施した．流れ場および空力音の評価においては，宮澤らの実験結果 (10) と比較した．
実験が行われた条件は，レイノルズ数 Re = 2.0 × 105 ，
迎え角 9.5◦ ，主流流速 20m/s，コード長 C = 0.15m である．この条件では，流れは翼負圧面側前縁付近で一端はく離し，すぐに乱流へと遷移して再付着する．計算では，領域全体で温度が 20◦ C で一様であると仮定し，音速一定として主流マッハ数 0.058 で計算を行った．
計算領域は，流れ方向，垂直方向，スパン方向の長さ

ここで，ϵ = 1.0 × 10−5 および b = 0.1 とした．また，温度は T = 20◦ C で一様に与えた．
計算領域は −10m≤ x ≤ 10m，−10m≤ y ≤ 10m の正方形領域であり，計算格子は各方向等間隔に 201 点および 401 点を配置した二種類の格子を用いた．以下では，それぞれ G201 および G401 と称する．LBM ではクーラン数は常に 1 であり，格子幅によって時間刻み幅は固定される．Navier-Stokes 計算の時間刻み幅は熱流体 LBMの値に一致させ，∆t = 8.8 × 10−5 s とした．

3

-4.0

EXP (Miyazawa)
LBM (Fine grid)
LBM (Coarse grid)

-3.0

-2.0

Cp
-1.0

y

0.0

x
1.0
0.0

Fig. 4 Computational mesh for NACA0012 (block lines
are only shown)

0.2

0.4

0.6

0.8

1.0

x /C
(a) MeanϬ͘ϱTab. 2 Details of computational mesh
Points in Cube

Coarse

Fine

173

173

Number of Cubes (×104 )7

24

Minimum Spacing (y + )

4.6

2.3

Total Points

340 million

1.2 billion

yW;DŝǇĂǌĂǁĂͿ
Ϭ͘ϰ

>D

Ϭ͘ϯ

C p′
Ϭ͘Ϯ

Ϭ͘ϭ

がそれぞれ 24C ，24C ，0.125C の直方体領域とした．図4 に翼周りの計算格子を示す．同図には，計算ブロックに当たる BCM の Cube が示されている．各 Cube には，
17 × 17 × 17 の格子点を設定した．翼表面近傍の最小格子
幅は壁座標で 2.3 である．主流部の最大格子幅は 0.0078C
であり，これは 15kHz の音波についてその波長を 20 格
子点で解像することができる．総格子点数は約 12 億点である．また本研究では，表 2 に示すように，格子幅が約 2倍の粗い格子についても計算を実施し，計算で解像される最小渦スケールの影響について調査した．なお，本計算
では，計算領域中心から直径 16C の円の外側を Damping
region とした．

Ϭ
Ϭ͘Ϭ

Ϭ͘ϭ

Ϭ͘Ϯ

Ϭ͘ϯ

Ϭ͘ϰ

Ϭ͘ϱ

Ϭ͘ϲ

Ϭ͘ϳ

Ϭ͘ϴ

Ϭ͘ϵ

ϭ͘Ϭ

x/C
(b) RMS
Fig. 5 Pressure distributions around the airfoil

′
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

4.2 結果および考察
はじめに，流れ場に関して実験結果 (10) との比較を示す．図 5 に翼面静圧の時間平均および分散分布を示す．上
述のとおり，負圧面側の前縁近傍で境界層の剥離が生じているため，時間平均圧力がコード方向に変化していな
い領域を確認できる．しかし，その剥離は小規模であり，
x/C = 0.1 までに再付着している．x/C = 0.1 のすぐ上流では，剥離遷移が起こっており，これに伴って圧力のRMS が高くなっている．LBM 結果は，実験結果と比較して概ね良く一致しているが，再付着がわずかに遅れている．粗い格子の結果をみると，層流剥離泡が過大に評価されているのがわかる．剥離せん断層における K-H 不安定性をうまく予測できていない．
図 6 は，翼下流における時間平均速度分布および RMS
速度分布の比較を示している．時間平均速度分布をみると，LBM 結果は後流の幅および後流の深さともに実験結果と非常によく一致している．RMS 速度をみると，実験
結果には主流中でもわずかながら速度変動がみられ，主流乱れの存在を確認できる．LBM では，主流乱れは導入していないため，主流の RMS 速度はゼロとなっている．
後流中について実験と比較すると，速度変動は過大評価されている．翼負圧面側の境界層は，後縁でかなり厚くなっており，境界層外縁付近の格子解像度が十分でない

LBM

y/C

0.1
0.0
-0.1
-0.2
-0.2

0.0

0.2

0.4

umean

Fig. 6 Time-averaged velocity and RMS velocity distributions at the downstream of airfoil
ことが考えられる．
図 7 は瞬時の流れ場を示す．同図 (a) は等スパン面上の渦度分布を，同図 (b) は翼前縁近傍の渦流れ構造を示している．前縁近傍で剥離したせん断層が巻き上がり，それが崩壊して下流へと移流されている様子が示されている．
同図 (c) には翼負圧面上における圧力変動量の分布，同図(d) には同時刻における翼負圧面近くの流れ方向速度の変動量分布を示している．同図 (c) をみると，x/C = 0.07
付近で，スパン方向に連なった低圧領域が観察され，は4

Separated boundary layer

Observation point
0

ω 300

6.7C

(a) Vorticity magnitude contour

3.5 ×10−4ps′

0.5 ρ∞U ∞ 2−3.5 ×10−4

1.0

Hn

Fig. 8 Instantaneous pressure ﬂuctuation ﬁeld

−1.0(b) Vortex cores (∇2 p = 150)

では，高い圧力変動が確認できる．また，同図 (d) の翼負圧面近傍（y + = 5 付近）の流れ方向速度の変動分布をみると，x/C > 0.10 の領域で流れ方向に筋状に連なった低速域（低速ストリーク）がスパン方向に複数観察される．
そのスパン方向の間隔は壁座標で 150 程度である．本計算では，十分な解像度で壁乱流を支配する縦渦構造が捉えられており，壁まで流れ場を直接計算できていると考
えられる．

0.1

p′
0.5ρ ∞U ∞ 2
−0.1

0.00

x/C

0.05

0.10

0.15

0.20

図 8 は，等スパン断面における瞬時の圧力変動分布を示している．前縁および後縁付近から不規則な圧力波が伝播している様子が見られる．前縁付近から発生している圧力波は，前縁近傍で生じた剥離の非定常的な変動によると考えられる．剥離せん断層は，上述のとおり周期
的に巻き上がり三次元的に崩壊しながら下流へと放出される過程で，前縁近傍で圧力変動が生じている．また後縁では，翼負圧面上の発達した乱流境界層内の渦が干渉して騒音を引き起こしている．図 9 には，x/C = 0.075，
x/C = 0.25 および x/C = 0.9 における翼面圧力変動のスペクトル密度の分布を示す．同図 (a) の x/C = 0.075
についてみると，計算結果の圧力変動は全体的に過小評
価されている．また，計算では St = 50 付近に実験結果にはみられないピークを確認できる．このピークは剥離
せん断層の巻き上がりに対応する．x/C = 0.075 は，実際の再付着点付近であると考えられ，計算上の再付着点である x/C = 0.1 をみると実験結果によく一致することがわかる．上述したように，計算では剥離が過大評価されているとみられ，これはスパン長の違いだけではなく主流乱れの有無とも関連すると推察される．同図 (b) に示す x/C = 0.25 における実験結果の翼面圧力変動スペ
クトルについてみると，St = 60 付近にピークがみられる．文献 (10) によると，このピークは静圧孔のヘルムホルツ共鳴に起因しており，流れの変動とは直接関係がな
い．したがって，このピークを除けば，本計算結果は広い周波数域に渡って実験結果と良い一致が得られている．
同図 (c) についても同様のことがいえ，音源となる翼面近くの乱流境界層内の乱流渦まで精度良く再現できていると考えられる．

0.25(c) Pressure ﬂuctuation contour

0.2
u′

U∞
−0.2

0.00

x/C

0.05

0.10

0.15

0.20

0.25

(d) Streamwise velocity ﬂuctuation contour
Fig. 7 Instantaneous ﬂow ﬁeld

く離せん断層が二次元的に巻き上がっていると考えられる．剥離せん断層の巻き上がりはじめはある程度二次元性を有しており，スパン方向に軸を持つ組織的な大規模渦構造が形成される．このことは，同図 (b) の渦コアから確認できる．しかし，同図 (b) に示されるように，剥離せん断層と翼表面の間には比較的小規模な流れ方向に軸を持つ縦渦構造が多数発生しており，剥離せん断層の巻
き上がりはすぐに崩壊していることがわかる．この剥離
せん断層の巻き上がりが崩壊して三次元化する結果，乱流に遷移して再付着する．再付着点の x/C = 0.10 付近

5ϭ͘ϬнϬϬ50ϭ͘ϬͲϬϮ
ϭ͘ϬͲϬϯ
ϭ͘ϬͲϬϰ
ϭ͘ϬͲϬϱ
ϭ͘ϬͲϬϲ
ϭ͘ϬͲϬϳ

yW;DŝǇĂǌĂǁĂͿϭ͘ϬͲϬϴ

>D;ǆͬĐсϬ͘ϬϳϱͿϭ͘ϬͲϬϵ

>D;ǆͬĐсϬ͘ϭϬϬͿ

EXP (Kato)

]z
40
H
/
B
[d30
yt
is
n20
e
D
l 10
a
rt
c
e0
p
S

ϭ͘ϬͲϬϭLBM (Fine grid)
LBM (Coarse grid)-10ϭ͘ϬͲϭϬ

100

ϭ

ϭϬ

St

1000

10000

f [Hz]

ϭϬϬ(a) x/C = 0.075

Fig. 10 Frequency spectra of far-ﬁeld soundϭ͘ϬнϬϬ

yW;DŝǇĂǌĂǁĂͿϭ͘ϬͲϬϭ
ϭ͘ϬͲϬϮ
ϭ͘ϬͲϬϯ

的に大きく評価されている．補正後の結果について見ると，Fine grid の計算結果は，低周波数域で実験結果（黒い線）よりやや大きく評価されているものの，700Hz より高周波数域で良い一致が得られており，乱流音の広帯
域スペクトルの定量的に予測できている．一方，Coarse
grid の計算結果は，600Hz 付近より高周波数域で急激に音圧レベルが下がっており，600Hz から 1100Hz の周波
数域で音圧レベルを過小評価している．乱流音の広帯域
スペクトルを定量的に予測するためには，Fine grid 程度の格子解像度が必要であることがわかる．
700Hz 以下の低周波数域では，実験結果と比較して計算結果に過大評価がみられる．これは波長の長い低周波の音波には計算領域が十分でないことが考えられる．しかしながら，700Hz から 2kHz の高周波数の領域については，LBM 結果は実験結果と良い一致が得られており，
空力騒音を定量的に予測できていることがわかる．

>D

Slope -5/3

ϭ͘ϬͲϬϰ
ϭ͘ϬͲϬϱ
ϭ͘ϬͲϬϲ
ϭ͘ϬͲϬϳ
ϭ͘ϬͲϬϴ
ϭ͘ϬͲϬϵ
ϭ͘ϬͲϭϬ
ϭ

ϭϬ

St

ϭϬϬ

(b) x/C = 0.25ϭ͘ϬнϬϬ

yW;DŝǇĂǌĂǁĂͿϭ͘ϬͲϬϭ
ϭ͘ϬͲϬϮ

>D

Slope -5/3

ϭ͘ϬͲϬϯ

結 言
本研究では，LBM をターボ機械における空力音の直接解析に応用することを目的として，NACA 翼から発生する広帯域騒音の予測精度を実験結果と比較することにより検証した．
格子解像度による比較では，粗い格子では解像可能な乱流境界層内の最小渦スケールが大きくなるだけでなく，微小擾乱の成長を再現する必要のある剥離遷移現象をうまく予測できない．その結果，遠方場の音圧スペクトルが過小評価された．
十分密な計算格子では，翼面上の圧力分布および後流速度分布に関する流れ場の比較において，LBM 結果は実
験結果と良い一致を示した．遠方場の空力音スペクトルにおいても，700Hz から 2kHz の高周波数帯域で実験結
果と良い一致が得られた．以上のように，非熱流体の等温モデルを使用した LBM であっても，温度変化の小さ
い低速流れであれば空力騒音の予測が可能であることが
示された．

5.

ϭ͘ϬͲϬϰ
ϭ͘ϬͲϬϱ
ϭ͘ϬͲϬϲ
ϭ͘ϬͲϬϳ
ϭ͘ϬͲϬϴ
ϭ͘ϬͲϬϵ
ϭ͘ϬͲϭϬ
ϭ

ϭϬ

St

ϭϬϬ

(c) x/C = 0.90
Fig. 9 Power spectral densities of blade surface pressurecoeﬃcient

図 10 は，翼後縁から流れに対して上方に 6.7C だけ離れた位置における空力音のスペクトルを示している．同図には，比較のために計算結果に加えて実験結果 (11) が示されている．本研究では，FFT による周波数分析には，
サンプリング点数 8192 個のデータを用いた．サンプリン
グ周波数は 250kHz である．さらに，20 組の周波数分析
結果に対しアンサンブル平均を施した．

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

赤線は Fine grid を用いて計算した結果であり，青線はFine grid に対して各方向に格子解像度を 1/2 に粗くした
Coarse grid を用いて計算した結果である．計算結果にはスパン方向長さに関する補正を施している．補正前を点線で，補正後を実線で示している．スパン方向に周期境界条件を用いて計算しているため，補正前の結果は全体(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,6

(2006), pp. 2999–3006

Boundary Condition,” AIAA Journal, Vol.45, No.7
(2007), pp.1703-1712.

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

(8) Li, Q., Luo, K. H. He, Y. L., Gao, Y. J. and Tao, W.
Q., “Coupling lattice Boltzmann model for simulation of thermal ﬂows on standard lattices,” PhysicalReview E, Vol.85, 016719 (2012).

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

(9) Lele, S., “Compact ﬁnite diﬀerence schemes with
spectral-like resolution,” Journal of Computational
Physics, Vol.103, No.1 (1992), pp.16-42.

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

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

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

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

(7) Kam, E. W. S., et al., “Lattice Boltzmann Method
Simulation of Aeroacoustics and Nonreﬂecting

7