音響インテンシティとpythonでの可視化【音響理論基礎】

本記事では音響インテンシティとは何かを解説して，その後平面波と球面波の音響インテンシティ(ベクトル場)を可視化してみます。

平面波及び球面波については以下の記事をご覧ください。

【音響理論基礎】平面波と球面波〜pythonによる可視化〜

まずは簡単な理論を述べてからpython言語での可視化を試みます。

音響インテンシティとは？

音響インテンシティ(acoustic intensity)とは流体中における音のエネルギー移動の割合を表したベクトル量で，音響パワー束密度とも呼ばれます。

イメージとしては音のエネルギーの大きさと流れ(方向)をベクトルで表現したものだと思ってもらって良いと思います。

音響インテンシティには大きく分けて瞬時音響インテンシティ(instantaneous acoustic intensity)と平均音響インテンシティ(average acoustic intensity)の2つがあります。

それぞれ簡単にご紹介します。

瞬時音響インテンシティ

瞬時音響インテンシティはある時刻の瞬時的な音響インテンシティ(エネルギーの移動)を表すベクトル量です。

瞬時音響インテンシティ$\boldsymbol{I}(\boldsymbol{r},t)$は音圧$p(\boldsymbol{r},t)$と粒子速度$\boldsymbol{u}(\boldsymbol{r},t)$を用いて

$$\boldsymbol{I}(\boldsymbol{r},t) = p(\boldsymbol{r},t) \boldsymbol{u}(\boldsymbol{r},t)$$

のように表されます。

また，音響インテンシティ$\boldsymbol{I}(\boldsymbol{r},t)$とエネルギー密度$e(\boldsymbol{r},t)$の関係は

$$\frac{\partial e}{\partial t}=-\boldsymbol{\nabla} \cdot \boldsymbol{I}$$

となっています。ここでエネルギー密度は

$$e=\frac{\rho_{0}}{2} |\boldsymbol{u}|^{2}+\frac{1}{2 \rho_{0} c^{2}} p^{2}$$

で与えられています。

平均音響インテンシティ

定常状態(単一周波数)の場では1周期に渡る全体のインテンシティを指標として用いる場合が多く，以下のように書かれます。

$$\boldsymbol{I}(\boldsymbol{r},\omega) = \frac{1}{2} {\rm Re} \{ p(\boldsymbol{r},\omega) \boldsymbol{u}(\boldsymbol{r},\omega)^{\ast} \}$$

また，オイラーの運動方程式より音響と粒子速度の関係は

$$\rho_{0} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} p = 0 \tag{$2^{\prime \prime}$}$$

で与えられています。 ※必要に応じて以下の記事をご参照ください

【音響理論基礎】波動方程式の導出

この式を両辺フーリエ変換して(周波数領域で考えて)みると

$$i \omega \rho_{0} \boldsymbol{u}=\boldsymbol{\nabla} p$$

と整理できます。これより平均音響インテンシティは

$$\boldsymbol{I} = \frac{1}{2} {\rm Re} \left\{ p \left(\frac{1}{i \omega \rho_{0}} \boldsymbol{\nabla} p \right)^{\ast} \right\}$$

のように簡単に書き換えることが出来ます。

今回扱うのは単一周波数における平面波及び球面波ですので，こちらの式を用いて定式化，可視化していきます。

平面波・球面波の平均音響インテンシティ

上で述べたインテンシティの式から平面波と球面波の平均音響インテンシティを導出します。
平面波及び球面波の式は

\begin{align}
p_{p} &= e^{i \boldsymbol{k} \cdot \boldsymbol{r}} \notag \\
p_{s} &= \frac{e^{ikr}}{r} \notag
\end{align}

とします。
では早速導出していきましょう！単純な計算ですので皆さんもやってみて下さい。

平面波の平均音響インテンシティ

$$\boldsymbol{\nabla} p_{p} = i \boldsymbol{k} e^{i \boldsymbol{k} \cdot \boldsymbol{r}}$$

となることに注意すると

\begin{aligned}
\boldsymbol{I}_{p} &=\frac{1}{2} {\rm Re}\left\{e^{i \boldsymbol{k} \cdot \boldsymbol{r}}\left(\frac{1}{i \omega \rho_{0}} \boldsymbol{\nabla} e^{i \boldsymbol{k} \cdot \boldsymbol{r}}\right)^{\ast}\right\} \notag \notag \\
&=\frac{1}{2} {\rm Re}\left\{ e^{i \boldsymbol{k} \cdot \boldsymbol{r}} \left(\frac{i \boldsymbol{k}}{i \omega \rho_{0}} e^{i \boldsymbol{k} \cdot \boldsymbol{r}}\right)^{\ast}\right\} \notag \notag \\
&=\frac{1}{2 \omega \rho_{0}} \boldsymbol{k} \notag
\end{aligned}

のようになります。
この式から平面波の平均音響インテンシティは位置に依存しないことが分かりますね。

球面波の平均音響インテンシティ

極座標での微分演算子$\boldsymbol{\nabla}$は

$$\boldsymbol{\nabla} =\boldsymbol{e}_{r} \frac{\partial}{\partial r} + \boldsymbol{e}_{\theta} \frac{1}{r} \frac{\partial}{\partial \theta} + \boldsymbol{e}_{\phi} \frac{1}{r \sin \theta} \frac{\partial}{\partial \phi}$$

ですので

$$\boldsymbol{\nabla} p_{s} = \frac{(ik -1)e^{ikr}}{r} \boldsymbol{e}_{r}$$

となって

\begin{aligned}
\boldsymbol{I}_{p} &=\frac{1}{2} {\rm Re}\left\{\frac{e^{i k r}}{r}\left(\frac{1}{i \omega \rho_{0}} \boldsymbol{\nabla} \frac{e^{i k r}}{r}\right)^{\ast}\right\} \notag \\
&=\frac{1}{2} {\rm Re}\left\{\frac{e^{i k r}}{r}\left(\frac{k r-i}{\omega \rho_{0} r} \frac{e^{i k r}}{r} \boldsymbol{e}_{r}\right)^{\ast}\right\} \notag \\
&=\frac{1}{2 c \rho_{0} r^{2}} \boldsymbol{e}_{r} \notag
\end{aligned}

のように導くことが出来ます。ここで$\omega = ck$を使いました。

音源からの距離が離れるにつれてインテンシティが小さくなっていくことが見て取れますね。

平面波・球面波のインテンシティ可視化

ここではpythonによる2次元平面のインテンシティの可視化を試みます。

例として上で導出した平面波・球面波のインテンシティをmatplotlibのquiver関数で可視化しますが，同じように他の場にも適用出来ますので是非使ってみて下さい。

まず音響インテンシティと音圧を引数にして可視化する関数を以下に示します。

def plot_intensity(I, P, x, y, disp):
fig, ax = plt.subplots()
ax.quiver(x, y, I[0], I[1], color='black', angles='xy', scale_units='xy')
p = P.real
ax.imshow(p, cmap=plt.cm.bwr, vmin=-1, vmax=1, interpolation='bicubic', origin='lower', extent=disp)
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
fig.tight_layout()
fig.show()

Iがインテンシティ，Pが複素音圧です。音圧については可視化しなくても良いのですが，理解を助けるために一緒に可視化することにしました。

平面波と球面波を算出するにあたって使用したパラメータもお示ししておきます。

freq = 1000
omega = 2 * np.pi * freq
rho_0 = 1.293
kappa = 142.0e3
c = np.sqrt(kappa/rho_0)
k = omega / c
row = 21
col = 21
width = 0.05
y, x = np.mgrid[0:row:1+width, 0:col:1+width] * width
y = y - width(row/2) x = x - width(col/2)
disp = np.array([np.min(x), np.max(x), np.min(y), np.max(y)])

平面波の可視化

平面波の平均音響インテンシティ(と音圧)を算出して可視化するコードは以下のように書けます。

# plane wave
theta = np.pi/2
k_vec = k*np.array([np.cos(theta), np.sin(theta)])
I_p = k_vec / (2*omega*rho_0)
P_p = np.exp(1j * (x*k_vec[0] + y*k_vec[1]))
# drawing
plot_intensity(I_p, P_p, x, y, disp)

結果として出力される画像は以下です。

同じ方向に同じエネルギーが伝搬している様子が分かりますね。

球面波の可視化

同様にして球面波も可視化してみましょう。

# spherical wave
source = np.array([0, -1])
tmp = np.empty([2,x.shape[0],x.shape[1]], dtype=np.float64)
tmp[0] = x - source[0]
tmp[1] = y - source[1]
r = np.linalg.norm(tmp, axis=0)
e_r = tmp / r
I_s = e_r / (2*c*rho_0*np.power(np.pi*r,2))
P_s = np.exp(1j*k*r) / r
P_s[np.isinf(P_s)] = 0
P_s[np.abs(P_s)&amp;gt;1] /= np.abs(P_s[np.abs(P_s)&amp;gt;1])
# drawing
plot_intensity(I_s, P_s, x, y, disp)

結果は以下の通りです。