音波伝播のシミュレーションと風車騒音への応用をめざして理学専攻 情報科学コース

三木 麻梨子（指導教員：河村 哲也）

１．はじめに近年、問題になっている、エネルギー源の枯渇や地球温暖化の対策のひとつとして注目されているのが、風力・太陽・光等の地球にやさしいエネルギーである。その中でも、風車による風力発電は、自然の風を利用したクリーンな発電方法として、注目されている。
その一方で、騒音問題・電波障害・景観障害などの欠点が挙げられている。風車に限らず、騒音問題や振動問題がますます深刻になってきている昨今、音の発生や伝播の予測が必要とされている。その方法として、数値シミュレーションが活用されている。
そこで、本研究では、音波の基本現象と共に、回転する音源周りの解析を行った。モデルとして、Fig.1 に示す垂直軸揚力型のジャイロミル型風車に近いものを考える。波動性の考慮を重視するため、波動方程式を差分法によって解いた。

(Fig.1) ジャイロミル型風車

２．モデル化
本研究では、回転する音源周りの音波を観察するため、円柱座標系の３次元空間領域を考えた。円柱内の一部を音源と考え、それを回転させる。（Fig.2）

z

∂2 u

∂2u1 ∂u∂t∂rr ∂r= c22+21 ∂2u∂2 ur ∂θ∂z 2+ 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=TT=t（ω：回転角速度）
上記関係式によって式（3.1.1）を座標変換し、以下の式を得た。
∂2 u∂2 u ∂2 u ∂ 2 u∂2 u ∂u= c2+ 2 + 2 − ω2 Y Y 2 −22∂T∂X∂Y∂Z∂X∂Y22∂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.1.3）を解くために、同上の極座標変換を用いて、次式を導出した。

∂2 u∂2 uc2∂2 u c 2 ∂u
∂2 u
∂2 u222=c+−ω++c−2ω∂T 2∂r 2r2∂θ2 r ∂r
∂Z 2∂θ ∂T（3.1.4）

３－２．差分方法

音源

波動方程式の数値解法には、中心差分法を用いる。時間方向・空間方向共に、２次精度の中心差分（陽解法）を用いて近似した。ただし、（3.1.4）の最後の項

r(Fig.2)

３－１．基礎方程式
基礎方程式として、下記の３次元波動方程式を用いた。
∂t 2= c2∂u
∂Tについては、後退差分を用いて近似した。

３－３．格子生成３．計算方法∂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：時間）
以下の２通りについて検証した。
(1)固定座標系
式（3.1.1）を以下の極座標変換を用いると次式が得られる。
x = r cos θ , y = r sin θ(Fig.3)３－４．境界条件(1)r 方向
円周（r = rmax ）においては、反射要素のない空間を考えるため、完全吸収条件（3.4.1）を用いた。また、中心軸（r = r1）に関しては、回転中心軸条件（3.4.2）を用いた。
∂1

∂r

∂rr

∂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 の１つ外側の格子点 ）（u1 = u r1 u2 = u(r2 )）(2)𝛉方向周方向には、周期境界条件を与えた。
(3)z 方向上面・下面（z=0，zmax ）共に、r 方向と同様に完全吸収条件（（3.4.1）において r を z に置き換えた式）を用いた。

４．計算結果音 源 を k=17 ～ 25 の 長 さ に 固 定 し 、 ２ か 所 対 称(i,j,k)=(14,21,k),(14,61,k)に配置し、10000step 計算した。
角速度としてω = 0.05 × π(反時計回り)を設定した。

４－１．固定座標系分かる。近接点(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)

４－２．回転座標系
固定座標系では、音源が 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

５．まとめと今後の課題
本研究では、円柱座標の３次元空間に音源を２つ設置し、差分法による数値計算を用いて、回転する音源周りを解析した。音源を回転させる固定座標系による解析、また、回転座標系による解析共に、良い結果を得ることができた。
また、ドップラー効果もとらえることができた。
今後の課題として、風速や音の反射など他の要因を考慮し、音波の伝わりをシミュレーションしていきたい。また、風車の構造を厳密に再現し、ブレード数による違いも観察していきたい。

謝辞
本研究を進めるにあたり、ご尽力くださいました指導教員の河村哲也先生・諸先輩方に深く感謝いたします。
(Fig.6) 7880step(Fig.7)7880step（音源静止）（音源回転）
Fig.5 ,7 では、２つの音源からの波が時間的に干渉し、周辺に広がっていく様子が現われている。
音源が静止している Fig.4,6 と比較すると、Fig.5,7 共に、回転の前方向に波が圧縮され間隔が密になるのに対して、後方向は間隔が疎になり、ドップラー効果が起きていることが分かる。
次の Fig.8 は地面に近い水平面(k=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〔森北出版株式会社〕