音波伝播のシミュレーションと風車騒音への応用をめざして
理学専攻 情報科学コース
三木 麻梨子(指導教員:河村 哲也)
1.はじめに
近年、問題になっている、エネルギー源の枯渇や地球温
暖化の対策のひとつとして注目されているのが、風力・太
陽・光等の地球にやさしいエネルギーである。その中でも、
風車による風力発電は、自然の風を利用したクリーンな発
電方法として、注目されている。
その一方で、騒音問題・電波障害・景観障害などの欠点
が挙げられている。風車に限らず、騒音問題や振動問題が
ますます深刻になってきている昨今、音の発生や伝播の予
測が必要とされている。その方法として、数値シミュレー
ションが活用されている。
そこで、本研究では、音波の基本現象と共に、回転する
音源周りの解析を行った。モデルとして、Fig.1 に示す垂
直軸揚力型のジャイロミル型風車に近いものを考える。波
動性の考慮を重視するため、波動方程式を差分法によって
解いた。
(Fig.1) ジャイロミル型風車
2.モデル化
本研究では、回転する音源周りの音波を観察するため、
円柱座標系の3次元空間領域を考えた。円柱内の一部を音
源と考え、それを回転させる。(Fig.2)
z
∂2 u
∂t 2
= c2
∂2u
∂r 2
+
1 ∂u
r ∂r
+
1 ∂2u
∂2 u
r 2 ∂θ
∂z 2
+
2
音源が、時間とともに角速度ωで回転するとして計算する。
(2)回転座標系
音源の回転とともに回転する座標系を考え、音源を格子
に固定した状態で計算する。固定座標系(x,y,z)から回転座
標系(X,Y,Z)への変換には以下の関係式を用いた。
x = X cos(ωT) + Y sin( ωT)
X = x cos(ωt) − y sin(ωt)
y = −X sin(ωT) + Y cos(ωT)
Y = x sin(ωt) + y cos(ωt)
t=T
T=t
(ω:回転角速度)
上記関係式によって式(3.1.1)を座標変換し、以下の式
を得た。
∂2 u
∂2 u ∂2 u ∂ 2 u
∂2 u ∂u
= c2
+ 2 + 2 − ω2 Y Y 2 −
2
2
∂T
∂X
∂Y
∂Z
∂X
∂Y
2
2
∂
u
∂u
∂
u
− ω2 X X 2 −
+ 2ω2 XY
∂Y
∂X
∂X ∂Y
∂2 u
∂2 u
− 2ωX
+ 2ωY
∂T ∂Y
∂T ∂X
(3.1.3)
次に、3次元円柱座標系にて式(3.1.3)を解くために、
同上の極座標変換を用いて、次式を導出した。
∂2 u
∂2 u
c2
∂2 u c 2 ∂u
∂2 u
∂2 u
2
2
2
=
c
+
−
ω
+
+
c
−
2ω
∂T 2
∂r 2
r2
∂θ2 r ∂r
∂Z 2
∂θ ∂T
(3.1.4)
3-2.差分方法
音源
波動方程式の数値解法には、中心差分法を用いる。時間
方向・空間方向共に、2次精度の中心差分(陽解法)を用
いて近似した。ただし、
(3.1.4)の最後の項
r
(Fig.2)
3-1.基礎方程式
基礎方程式として、下記の3次元波動方程式を用いた。
∂t 2
= c2
∂u
∂T
について
は、後退差分を用いて近似した。
3-3.格子生成
3.計算方法
∂2u
(3.1.2)
∂2u
∂x 2
+
∂2u
∂y 2
+
∂2 u
∂z 2
格子数は、半径方向(r)に 50、周方向(θ)に 81、軸方向
(z)に 40 とした。円柱の中心は回転軸としている。計算に
おいては、それぞれを i,j,k としている。 音源
(3.1.1)
(u:空気密度または圧力、c:音波の伝播速度、t:時間)
以下の2通りについて検証した。
(1)固定座標系
式(3.1.1)を以下の極座標変換を用いると次式が得られ
る。
x = r cos θ , y = r sin θ
(Fig.3)
3-4.境界条件
(1)r 方向
円周(r = rmax )においては、反射要素のない空間を考
えるため、完全吸収条件(3.4.1)を用いた。また、中心軸
(r = r1)に関しては、回転中心軸条件(3.4.2)を用いた。
∂
1
∂
r
∂r
r
∂r
∂u
∂r
−
∂
∂T
r=r 1
u
=
=0
r=r max
−u 2 +u 1
−r 22 4+r 21 4
(3.4.1)
(3.4.2)
( r1 :最も中心に近い格子点, r2 :r1 の1つ外側の格子点 )
(u1 = u r1 u2 = u(r2 ))
(2)方向
周方向には、周期境界条件を与えた。
(3)z 方向
上面・下面(z=0,zmax )共に、r 方向と同様に完全吸収
条件(
(3.4.1)において r を z に置き換えた式)を用いた。
4.計算結果
音 源 を k=17 ~ 25 の 長 さ に 固 定 し 、 2 か 所 対 称
(i,j,k)=(14,21,k),(14,61,k)に配置し、10000step 計算した。
角速度としてω = 0.05 × π(反時計回り)を設定した。
4-1.固定座標系
分かる。近接点(i=25)に比べて、遠方の点(i=35,45)の方が、
波の干渉で密度変化が大きく現われている。
2.50E-01
①I=25,J=41,K=3
2.00E-01
②I=35,J=41,K=3
空1.50E-01
気
密
1.00E-01
度
③I=45,J=41,K=3
5.00E-02
0.00E+00
100 1100 2100 3100 4100 5100 6100 7100 8100 9100
(Fig.8)
4-2.回転座標系
固定座標系では、音源が Fig.3 に示した同心円(i=14)上
の格子点を時間的に移動するとした。音源周りの特性を精
度良く取り扱うため、回転座標系で音源近傍を計算した。
Fig.9 , 10 は音源が静止した場合と回転した場合を、同一
の時間ステップで比較したものである。Fig.9 では、音波が
左右対称になっているのに対して、Fig.10 では、回転によ
って音波が非対称に歪んでいる様子が現われている。
Fig.4,5 は、軸方向で音源中央(k=21)の水平断面での密度
の時間変化を示したものである。
(Fig.9)3720step
(Fig.4) 7880step
(Fig.5)7880step
(音源静止)
(音源回転)
Fig.6,7 は、半径方向音源位置(i=14)での円筒面での密度
の時間変化を示したものである。
(Fig.10)3720step
5.まとめと今後の課題
本研究では、円柱座標の3次元空間に音源を2つ設置し、
差分法による数値計算を用いて、回転する音源周りを解析
した。音源を回転させる固定座標系による解析、また、回
転座標系による解析共に、良い結果を得ることができた。
また、ドップラー効果もとらえることができた。
今後の課題として、風速や音の反射など他の要因を考慮
し、音波の伝わりをシミュレーションしていきたい。また、
風車の構造を厳密に再現し、ブレード数による違いも観察
していきたい。
謝辞
本研究を進めるにあたり、ご尽力くださいました指導教
員の河村哲也先生・諸先輩方に深く感謝いたします。
(Fig.6) 7880step
(Fig.7)7880step
(音源静止)
(音源回転)
Fig.5 ,7 では、2つの音源からの波が時間的に干渉し、周
辺に広がっていく様子が現われている。
音源が静止している Fig.4,6 と比較すると、Fig.5,7 共に、
回転の前方向に波が圧縮され間隔が密になるのに対して、
後方向は間隔が疎になり、ドップラー効果が起きているこ
とが分かる。
次の Fig.8 は地面に近い水平面(k=3)上で、半径方向の異
なる3点(i=25,35,45)における密度変化を比較したもので
ある。音源に近い点から順次、音波が伝わっていることが
参考文献
[1]桑原邦郎・河村哲也:流体計算と差分法 2005〔朝倉書店〕
[2]高見穎郎・河村哲也:偏微分方程式の差分解法 1994〔東京大
学出版会〕
[3]割田真弓:「閉空間内における音波の伝播の数値シミュレーシ
ョン」お茶の水女子大学大学院数理・情報科学専攻修士論文平成
16 年度
[4]河村哲也:エネルギーと風車 2003〔山海堂〕
[5]Kunio Kuwahara and Yuko Oshima:Thermal Convection Caused
by Ring-Type Heat Source 〔JPS Vol.51,No.11〕
[6]牛山泉:風車工学入門 2002〔森北出版株式会社〕
ダウンロード

音波伝播のシミュレーションと風車騒音への応用をめざして