種々の風車から発生する音波伝播の数値シミュレーション
理学専攻 情報科学コース 寺町幸希子 (指導教員：河村哲也)
はじめに

1

風力発電は代替エネルギー源として近年注目されているが，日本ではあまり実用化されておらず，その理由の一つとして風車騒音が挙げられている．
そのため本研究では，種々の大型風車において、風車ブレードの先端から発生する音を回転移動する音源とみなし，発生する音波がどのように伝播するかを数値シミュレーションすることで，風車騒音への応用を考えることを目的とした．

空間のモデル化・格子生成

2プロペラ型 (Fig.2.1) と直線翼型 (Fig.2.2) の 2 種類の風車について考えた. プロペラ型では, 縦横それぞれ 1，奥行き 0.5 の直方体領域を考え, 直線翼型では半径 1, 高さ 1 の円柱領域 (円柱座標系) を考えた. 本研
究では, 以下の 3 つの場合において音波伝播のシミュレーションを行った.

Fig.2.1: Propeller type2.1

境界条件

3.2.1

Case 1

PML 法（Perfectly Matched Layer）(2)(3) を利用し⃗ に対する補正項を V⃗1 , P に対て吸収層を設定した. Vする補正項を P1 とすると，∫∫ ∫⃗⃗⃗ dt)dt
V1 = m4 u + n4 V dt + n5 ( V∫∫ ∫∂∂P dt + n1( P dt)dt
∂x
∂x
∫
∫ ∫
P1 = m4 P + n4 P dt + n5 ( P dt)dt∫∫ ∫∂∂+ m1u dt + n1( u dt)dt
∂x
∂x
∫
∫ ∫∂∂
v dt + n2( v dt)dt+ m2∂y∂y∫∫ ∫∂∂+ m3w dt + n3( w dt)dt∂z∂z+ m1

Case 1Case 2

吸収層は設定せず，少ない格子数で遠方境界を扱えるようにするため，音源の回転半径を 1 とした時，x軸方向が 16，y 軸方向が 16，z 軸方向が 8 となるよう，音源から遠ざかる程格子の間隔が広くなる不等間隔格子を用いる. 格子数は 200 × 200 × 100 とした.

2.3

3.2

Fig.2.2:
Straight wing type

格子数は 200 × 200 × 100 とし, 各方向の比率がx : y : z = 1 : 1 : 1 である直交等間隔格子を用いる．
また地面を想定した面以外の面に吸収層を設置することで遠方境界を扱った.

2.2⃗ : 速度ベクトル
P : 圧力，VQ は音源項で，音源から発せられる波から数値誤差により余分な波が発生しないように，空間的に数点を使ったガウス分布の振幅を時間的に振動させることに⃗ 1 , P1 はよってなめらかな音源を再現した. (1) また, V吸収層を用いた際に生じる項である. (吸収層を用いない場合は 0 である)

Case 3

音源の回転半径を 1 とした時，円柱の半径が 72，高さが 72 となるよう, 音源から遠ざかる程格子の間隔が広くなる不等間隔格子を用いて遠方境界を扱った. 格
子数は 200 × 200 × 100 とした．

3

計算方法

3.1

基礎方程式

(3)となる．ここで補正項に現れるパラメータは，
m1 = σy + σz , m2 = σz + σx , m3 = σx + σy ,m1 + m2 + m3
m4 =, n1 = σy σz , n2 = σz σx , (4)2n3 = σx σy , n4 = n1 + n2 + n3 , n5 = σx σy σzである． σx , σy , σz はそれぞれ x 方向，y 方向，z 方向の吸収率を表しており，本研究ではそれぞれ 0.4 とし，吸収層が領域全体の 15 ％ となるように設定した.
式 (1)(2)(3)(4) より Case 1 では下記の方程式を解いた.
⃗∂V
⃗1(5)
= −∇P + V∂t∂P⃗ + Q + P1(6)= −∇·V∂t

3.2.2

Case 2 および Case 3

壁面での流速および圧力は 0 とした。壁面での微分値は２次精度の片側差分で近似した. また, Case 3 においては周方向に周期条件を課した.

以下の線形化オイラー方程式を使用する．

3.3
3.3.1

⃗
∂V= −∇P∂t(1)∂P⃗ +Q= −∇·V
∂t(2)

解法

Case 1

式 (5)(6) を以下の式を用いた陰的な差分法であるコンパクト差分法 (4)fi+ 23 − fi− 23fi+ 12 − fi− 12′′αfi+1+ fi′ + αfi−1=b+a3hh(7)Case 3 では, 音波は放射状に伝播した. 音源付近では, Case 1, 2 とほぼ同様な分布になった.

で解く. この式の係数を

a=322α − 11(3 − 2α), b =,α =8822とすることで，1 階微分を 4 次精度で近似した．

3.3.2

Case 2式 (1)(2) を以下の陽的な差分法で解く．5 点 (x −s, x − r, x, x + q, x + p) を用いた点 x における f の 1 階微分は，
fx = af (x + p) + bf (x + q) + ef (x)(8)+ cf (x − r) + df (r − s)Fig.4.3: Pressure distribution obtained by case 3であり，この式の係数を

prs
qrs,b = −p(p − q)(s + p)(r + p)q(p − q)(s + q)(r + q)pqs
pqr
c=−,d =r(r − s)(r + p)(r + q)s(r − s)(s + p)(s + q)a=

また, 風車の中心からの距離を等しくとったある地点での圧力 P の時間変化をグラフ表示した．青色の曲線は Case 1, 赤色の曲線は Case 2 であり, 縦軸は圧力P の大きさを, 横軸は時間を表している．

e = −(a + b + c + d)とすることで，1 階微分を 4 次精度で近似した．

3.3.3

Fig.4.4: Zoom up ﬁgure of Fig.4.3

Case 3式 (1)(2) を円柱座標系に変換する.

∂Vr∂P=−∂t∂r∂Vθ1 ∂P=−∂tr ∂θ∂P∂Vz=−∂t∂z(9)∂P∂Vr1 ∂Vθ
∂Vz= −(++)+Q∂t∂rr ∂θ∂z(10)Fig.4.5: Time history of pressure P

Case 2 と同様に, 式 (9)(10) を 4 次精度の陽的な差分法で解く．

全体を見ると曲線の推移は概ね一致していることから, 不等間隔格子を用いた方法で, 風車から遠い場所の結果を求めることは妥当であると考えられる.

54

結果と考察

音源を一定速度で 7/4 回転させた時の圧力 P を等値線で表した結果を示す.

今後の課題

今後は，風車を実際に設置する環境により近づくような不等間隔格子の格子間隔を検討し, また風車に対して様々な方向からの音波伝播を調べ, 風車騒音への応用を考えていきたい．

参考文献

Fig.4.1:
Pressure distribution
obtained by case 1

Fig.4.2:
Pressure distribution
obtained by case 2

Case 1 では, 音源が吸収層の近くを回転移動しているため，ほとんどの音波がすぐに吸収されている. 一方で, 1 つの音源が通過してから次の音源が通過するまでの時間が短いため, 地面からの反射波以外に, 吸収層の表面からのわずかな反射波の影響を受けている.
Case 2 では, 波は音源から遠く離れた場所へは, 時間が経つにつれて半円状を描いて音が伝播していくことが分かった.

[1] Tomonori Nihei, Katsuya Ishii, ”A fast solver ofthe shallow water equations on a sphere using a
combined compact diﬀerence scheme”, Journal ofComputational Physics 187, 639-659 (2003)[2] Berenger, J-P. ”A Perfectly Matched Layer for
the Absorption of Electomagnetic Waves” Journal of Computational Physics 114, 185-200 (1994)[3] Tam, C.K.W. T. ”Computational Acoustics : A
Wave Number Approach (Cambridge AerospaceSeries)” , Cambridge University Press (2012)[4] 鶴秀生, 岩津玲磨, 音響設計のための時間領域差分法の高精度化, 京都大学 数理解析研究所講究録, 第1645 巻, 177-188 (2009)