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

音響理論の基礎ということで，初最も重要な音響学の方程式である波動方程式の導出についてお話します。

まず波動方程式が何かということから解説して簡単な導出をご紹介したいと思います。

できるだけわかり易く書くつもりですが，偏微分などの簡単な数学の基礎は知らないと難しいかもしれません。

また波動方程式は量子力学や電磁気学など様々な分野で登場しますが，今回は音響分野の視点でのお話になります。

未熟ゆえ何かと間違いがあるかもしれません。誤りに気づきましたらご指摘頂ければ幸いです。

それでは本題です。

波動方程式とは？

波動方程式とは”波動現象を表す偏微分方程式”なのですが，初めて触れる方はよくわかりませんよね？

まず「波動」というのは”動く波”だとイメージするのが一番わかり易いんじゃないかと思っています。

「振動」と言えば皆さんイメージが湧きますよね？振動といえばギターの弦の振動や振り子の振動など色々ありますが，どれも振動する場所は固定されているかと思います。

これに対して「波動」は振動する位置も時間とともに動いていく(伝搬していく)ようなイメージです。

では実際に音圧についての3次元波動方程式を見てみましょう。以下のような形です。

$$\Delta p(\boldsymbol{r},t) – \frac{1}{c^{2}} \frac{\partial^{2} p(\boldsymbol{r},t)}{\partial t^{2}}= 0$$

ここで$\Delta$はラプラシアンと呼ばれる偏微分微分演算子で，$\Delta := \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}}$のように定義されています。

$\boldsymbol{r}$は3次元空間の位置で，$t$は勿論時刻です。

この方程式は3次元空間において時間とともに音(音圧)が従う振る舞いを表しています。

音を支配する方程式ですので音響分野にとって非常に重要なわけです。

では早速導出に移りましょう。

波動方程式の導出

まず音は空気の振動なので空気の運動を考える必要があります。

空気の運動と言うと「流体力学」を思い浮かべる方もいるかと思いますがその通りです。

音の波動方程式は流体力学の方程式から出てくるもので，流体力学の視点で考えると音は非常に微小な変動です。

実際大気圧が1013hPaなのに対して，音(として認識できる圧力)は20μPa〜2Pa程度しかありません。

そういうわけで流体力学に則って方程式を構築していきます。

空気を物理的に扱う時，原子や分子1つ1つに着目すると非常に乱雑な運動をするため上手く解析できません。

そこで流体力学ではある微小領域の分子の集まりを”粒子”として，その粒子達の振る舞いを解析します。

流体を表す変数は，圧力$P$，密度$\rho$，粒子速度$\boldsymbol{v}$，温度$T$，エントロピー$S$，内部エネルギー$U$等があるのですが，独立している変数は粒子速度(3変数)と他の2つだけです。

つまり合計で5つの変数を考えれば良いので5本の方程式を立てることで解析していきます。
具体的にこの5本の方程式は

連続方程式
運動方程式(3本)
状態方程式

です。これから1つずつ見ていきましょう。

※ここで流体力学の基礎方程式の導出には今井先生の「流体力学」を参考にさせて頂きました。流体力学の入門書としても非常に良い本だと思いますのでチェックしてみてください。

流体力学 (物理テキストシリーズ 9)

created by Rinker

連続方程式

連続方程式はある(有界閉)領域の流体の流入・流出の関係から導かれる式です。

単位時間あたりに領域$V$から出ていく(減少する)流体の質量は境界$\partial V$における粒子速度の外側法線方向成分を考えればいいので

$$\iint_{\partial V} \rho \boldsymbol{v} \cdot \boldsymbol{n} dS$$

のように書けます。ここで$\boldsymbol{n}$は境界における外側単位法線ベクトルです。

一方で領域全体の質量の単位時間の変化(増加分)は以下のように表せます。

$$\frac{\partial}{\partial t} \iiint_{V} \rho dV$$

2つの式から

$$\iint_{\partial V} \rho \boldsymbol{v} \cdot \boldsymbol{n} dS = – \frac{\partial}{\partial t} \iiint_{V} \rho dV$$

となります(符号に注意)。更にガウスの定理から

$$\iint_{\partial V} \rho \boldsymbol{v} \cdot \boldsymbol{n} dS = \iiint_{V} \boldsymbol{\nabla} (\rho \boldsymbol{v}) dV$$

なので整理すると

$$\iiint_{V} \left \{ \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} (\rho \boldsymbol{v}) \right \} dV = 0 $$

ここで上式について領域形状$V$は任意なので常に成り立つためには被積分関数が0になるしかありません。つまり

$$\frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} (\rho \boldsymbol{v}) = 0 \tag{1}$$

これが(オイラーの)連続方程式です。

運動方程式

運動方程式($ma=F$)は馴染み深いと思いますが，流体力学においても当然同様の意味です。

連続方程式の時のように領域$V$について考えて，力の釣り合いから導きます。

慣性力($ma$部分)については次のように書けます。

$$\iiint_{V} \rho \frac{D \boldsymbol{v}}{D t} dV$$

ここで$D/Dt$はラグランジュ微分と呼ばれる演算子で，$t$以外の変数も$t$に依存する場合に現れます。
あまり難しく考えず，位置を時刻に依存する変数として普通に全微分すれば良いと思います。
外力については，まず圧力について考えると，$\boldsymbol{n}$と逆向きに働くので

$$\iint_{\partial V} -P \boldsymbol{n} dS$$

と書けます。その他に単位質量あたりに働く外力を$\boldsymbol{F}$とすると

$$\iiint_{V} \rho \boldsymbol{F} dV$$

となります。これらをまとめて

$$\iiint_{V} \rho \frac{D \boldsymbol{v}}{D t} dV = \iiint_{V} \rho \boldsymbol{F} dV -\iint_{\partial V} P \boldsymbol{n} dS$$

となります。連続方程式と同様に体積積分=0の形に持っていきたいので，ガウスの定理から

$$\iint_{\partial V} P \boldsymbol{n} dS = \iiint_{V} \boldsymbol{\nabla} P dV$$

として代入すれば

$$\iiint_{V} \rho ( \boldsymbol{F} – \frac{D \boldsymbol{v}}{D t} – \frac{1}{\rho} \boldsymbol{\nabla} P) dV = 0$$

$V$は任意より

$$\frac{D \boldsymbol{v}}{D t} = \boldsymbol{F} – \frac{1}{\rho} \boldsymbol{\nabla} P \tag{2}$$

これが(オイラーの)運動方程式です。

状態方程式

状態方程式は熱力学の理論から導かれる方程式です。

厳密な議論には熱力学や統計力学の知識が必要となりますので本記事では深入りしません。

空気など普通の気体では理想気体の仮定が成り立ち，以下の方程式を満たします(ここまでは高校で習うレベル)。

$$PV = nRT$$

ここで$V$は体積，$n$は物質量，$R$は気体定数，$T$は絶対温度でした。

空気中の音の伝搬は熱の伝搬に比べて十分速いので，断熱過程とみなすことができます。よって断熱過程における理想気体の状態方程式

$$P = C\rho^{\gamma} \tag{3}$$

が成り立ちます。ここで$\gamma$は比熱比，$C$は定数です。

これが状態方程式です。ちなみにこのように圧力が密度の関数になっている流体をバロトロピー流体と言います。

3つの式の線形化

ここまでで得られた方程式は流体力学分野の方程式です。

音の方程式は当然この方程式に基づくのですが，先に述べたように音は大気圧に比べて非常に小さな変化ですので線形化することで単純な式に帰着することができます。

この節では波動方程式導出のため，3つの式を線形化します。

定常状態の圧力を$P_{0}$とすると，音圧$p$はそこからの変化分ですので

$$p = P – P_{0}$$

と書けます。同様に定常状態の密度を$\rho_{0}$として，密度変化を$\delta \rho = \rho – \rho_{0}$としておきます。

今風が吹かないような場を仮定すると音による粒子速度$\boldsymbol{u}$は$\boldsymbol{u} \approx \boldsymbol{v}$となります。

更に空気に対する外力(重力など)が十分に小さいことを仮定すると3つの式は

\begin{align} &\frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} (\rho \boldsymbol{u}) = 0 \tag{$1^{\prime}$} \\ &\frac{D \boldsymbol{u}}{D t} = – \frac{1}{\rho} \boldsymbol{\nabla} P \tag{$2^{\prime}$} \\\ &\frac{P}{P_{0}} = \left(\frac{\rho}{\rho_{0}}\right)^{\gamma} = (1 + s)^{\gamma} \tag{$3^{\prime}$} \end{align}

ここで

$$s := \frac{\delta \rho}{\rho_{0}}$$

は圧縮度と呼ばれる微小量です。

ここから線形化していきます。まず式($3^{\prime}$)から音圧と密度変化の関係を導出します。
式($3^{\prime}$)をテイラー展開すると

$$\frac{P}{P_{0}} = 1 + \gamma s + \frac{1}{2} \gamma(\gamma – 1)s^{2} + \cdots$$

これから

$$P = P_{0} + Ks + O(s^{2})$$

となります。ここで$K:=\rho_{0}\gamma$は体積弾性率と呼ばれる定数で，$O(s^{2})$はランダウの記号で$s^{2}$以降の高次項をまとめたものです。

$s$は微小量なので$O(s^{2})$は非常に小さいことから無視すると

$$p = P – P_{0} = Ks = \frac{K}{\rho_{0}} \delta \rho \tag{4}$$

ここで以下の定数を定義します。

$$c := \sqrt{\frac{dP}{d\rho}} = \sqrt{\frac{d}{d\rho}(P_{0} + K \frac{\rho-\rho_{0}}{\rho_{0}})} = \sqrt{\frac{K}{\rho_{0}}}$$

実はこの定数は音速に対応します。両辺2乗して

$$c^{2} = \frac{K}{\rho_{0}}$$

式(4)から

$$p = c^{2} \delta \rho \tag{$3^{\prime \prime}$}$$

次に式($1^{\prime}$)について

\begin{align} &\frac{\partial (\rho_{0} + \delta \rho)}{\partial t} + \boldsymbol{\nabla} {(\rho_{0} + \delta \rho) \boldsymbol{u}} \\ &= \frac{\partial \delta \rho}{\partial t} + \rho_{0} \boldsymbol{\nabla} \cdot \boldsymbol{u} + \boldsymbol{\nabla} (\delta \rho) \cdot \boldsymbol{u} + \delta \rho \boldsymbol{\nabla} \cdot \boldsymbol{u} = 0 \end{align}

2次の微小量(微小量×微小量)の項を無視して

$$\frac{\partial \delta \rho}{\partial t} + \rho_{0} \boldsymbol{\nabla} \cdot \boldsymbol{u} = 0 \tag{$1^{\prime \prime}$}$$

となります。式($2^{\prime}$)について

\begin{align} &(\rho_{0} + \delta \rho) \frac{D \boldsymbol{u}}{D t} + \boldsymbol{\nabla} (P_{0} + p) \\ &= (\rho_{0} + \delta \rho) \left \{ \frac{\partial \boldsymbol{u}}{\partial t} + u_{x} \frac{\partial \boldsymbol{u}}{\partial x} + u_{y} \frac{\partial \boldsymbol{u}}{\partial y} + u_{z} \frac{\partial \boldsymbol{u}}{\partial z} \right \} + \boldsymbol{\nabla} p = 0 \end{align}

先程同様2次の微小量を無視すると

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

となります。これで波動方程式を導出する準備は完了です。
では次はついに波動方程式を導出します。

波動方程式の組み立て

\begin{align} &\frac{\partial \delta \rho}{\partial t} + \rho_{0} \boldsymbol{\nabla} \cdot \boldsymbol{u} = 0 \tag{$1^{\prime \prime}$} \\ &\rho_{0} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} p = 0 \tag{$2^{\prime \prime}$} \\ &p = c^{2} \delta \rho \tag{$3^{\prime \prime}$} \end{align}

から波動方程式を導出します。
式($1^{\prime \prime}$)，($3^{\prime \prime}$)から

$$\frac{1}{c^{2}} \frac{\partial p}{\partial t} + \rho_{0} \boldsymbol{\nabla} \cdot \boldsymbol{u} = 0$$

$t$で偏微分して

$$\frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}} + \boldsymbol{\nabla} \cdot \left(\rho_{0} \frac{\partial \boldsymbol{u}}{\partial t} \right) = 0$$

式($2^{\prime \prime}$)より

$$\frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}} – \boldsymbol{\nabla} \cdot \boldsymbol{\nabla} p = 0$$

よって波動方程式

$$\Delta p – \frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}} = 0$$

が導出できました。長い！笑