波動方程式の導出(音波 縦波)

# はじめに
波動方程式の導出をすっかり忘れてしまったので勉強ついでに記録します。
この波動方程式を差分化して数値計算すればこのようなシミュレーションもできます。
![4x10m_1kHz.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/a64a4492-2f62-c860-83a6-f443d3a7db06.gif)

# 参考
[波動方程式から理解する音響学]
(https://www.jstage.jst.go.jp/article/jasj/66/9/66_KJ00006579568/_pdf)

# 波動方程式の導出
音波は音圧、密度、粒子速度という物理量が密接に関わっています。
まずはこれらの関係を見てみます。続いて状態方程式、連続の式、運動方程式について考えていきます。
注意点として、ここで扱うのは微小圧力で断熱過程が仮定できる無限小振幅音波の場合で、非常に高い音圧であったり、波長が著しく短い場合は有限振幅音波として扱う必要があります。

## 0. 変数一覧
$p_0$  大気圧$\left[ \rm{Pa} \right]$
$p$  音圧(変動分)$\left[ \rm{Pa} \right]$
$\rho_0$  密度(大気圧中)$\left[ \rm{kg/m^3} \right]$
$\rho$  密度(変動分)$\left[ \rm{kg/m^3} \right]$
$u$  粒子速度$\left[ \rm{m/s} \right]$
$c$  音速$\left[ \rm{m/s} \right]$
$t$ 時間$\left[ \rm{s} \right]$
$S$  シリンダー断面積$\left[ \rm{m^2} \right]$

## 1. 圧力と密度の関係
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/78b9e099-4977-5b97-3137-e1a840e15927.png)


図でピストンが右へ動き、媒質が圧縮され、圧力は$p_{tot}=p_0+p$、密度は$\rho_{tot}=\rho_0+\rho$に変化します。
ここで、$p_{tot}=f\left( \rho_0+\rho \right)$なる関数があるとします。
$\rho$が十分に小さいとして、テイラー展開して2次以降を無視し近似します。

```math
p_{tot}=f(\rho_{tot}) \approx f(\rho_0)+\rho \frac{df(\rho_0)}{d \rho_{tot}} \\
\therefore p=\rho \frac{df(\rho_0)}{d \rho_{tot}}
```

となります。ここで微分部分は定数になるので、圧力と密度は線型関係にあります。

## 2. 圧力と粒子速度の関係
微小時間$t$秒後のピストンが気体を押す力は$pS$[N]で、作用している気体の質量は$\rho_0 ctS$[kg]なので、運動量保存則から

```math
\frac{d(\rho_0 ctS)u}{dt}=pS
```

となり、両辺を$t$について積分すれば、

```math
\begin{eqnarray} \rho_0 ctSu &=& pSt \\
p &=& \rho_0 cu \end{eqnarray}
```

またもや線型関係にあることがわかります。

## 3. 密度と粒子速度の関係
微小時間$t$秒後のピストンの変位は$ut$[m]、圧力変化の先端は$ct$[m]で、ピストンが動く前と動いた後の質量保存から、

```math
\begin{eqnarray} \rho_0 ctS &=& (\rho_0 + \rho)(ct-ut)S \\
&=& (\rho_0 ct - \rho_0 ut + \rho ct - \rho ut)S \end{eqnarray}
```

となりますが、ここで$\rho$と$u$の積は微小量同士なので無視して近似します。すると、

```math
u=\frac{\rho}{\rho_0} c
```

という線型関係が導けます。$p = \rho_0 cu$に代入すると、

```math
p=c^2 \rho
```

が導けます。

## A. 状態方程式
先程の式をそのまま使います。

```math
p=c^2 \rho
```

## B. 連続の式
![図2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/6294dd87-1ca1-16a3-48c6-e90208075f84.png)


媒質が湧き出したり消滅したりすることはないとします。
図のように微小体積$dxS$に単位時間に流入・流出する媒質の質量が保存される条件式を考えます。
左から単位時間に流入する媒質の質量は$\rho_{tot}(x)u(x)S$と表され、右側から単位時間に流出する質量は$\rho_{tot}(x+dx)u(x+dx)S$です。
微小体積中に単位時間に留まる質量は両者の差ですから、$\rho_{tot}(x)u(x)S - \rho_{tot}(x+dx)u(x+dx)S$となります。  
ここで$x+dx$に関わる項をテイラー展開して2乗項以降を無視して近似します。

```math
\rho_{tot}u(x+dx) \approx \rho_{tot}u(x) + \frac{\partial \rho_{tot}u(x)}{\partial x}dx
```

従って、単位時間に微小体積中に留まる質量は、

```math
\rho_{tot}(x)u(x)S - \rho_{tot}(x+dx)u(x+dx)S = - \frac{\partial \rho_{tot}u(x)}{\partial x}dxS
```

質量が入ってきたら密度が大きくなり、出ていったら密度は小さくなるはずです。  
今度は密度の変化を時間の関数$\rho_{tot}(t)$であると考えると、[密度の変化率]×[体積]=[単位時間に留まる質量]ですから、

```math
\frac{\partial \rho_{tot}}{\partial t}dxS = - \frac{\partial \rho_{tot}u}{\partial x}dxS \\
\therefore \frac{\partial \rho_{tot}}{\partial t} + \frac{\partial \rho_{tot}u}{\partial x} = 0
```

ここで$\rho_{tot}u=\rho_0 u + \rho u$ですから、微小量同士の積を無視して、

```math
\frac{\partial \rho}{\partial t} + \rho_0 \frac{\partial u}{\partial x} = 0
```

となります。

## C. 運動方程式
音圧によって動かされる媒質を記述する運動方程式を書きます。
図2の質量保存の時と同様に考えます。微小体積を動かす外力は左右に働く力の差です。
両側に働く力の差は$p_{tot}(x)S-p_{tot}(x+dx)S$です。$x+dx$に関係する項をテイラー展開して2乗項以降を無視して近似します。

```math
p_{tot}(x+dx) \approx p_{tot}(x) + \frac{\partial p_{tot}(x)}{\partial x}dx
```

従って、

```math
p_{tot}(x)S-p_{tot}(x+dx)S = -\frac{\partial p_{tot}(x)}{\partial x}dxS
```

となります。運動方程式$ma=F$に従って考えると、$m=\rho_0dxS$、$a=\frac{\partial u}{\partial t}$なので、

```math
\rho_0dxS \frac{\partial u}{\partial t} = -\frac{\partial p_{tot}(x)}{\partial x}dxS \\
\therefore \frac{\partial p}{\partial x} + \rho_0 \frac{\partial u}{\partial t} = 0
```

## D. 波動方程式の組み立て
ようやく準備ができました。
まずは連続の式の両辺を$t$で偏微分します。

```math
\frac{\partial^2 \rho}{\partial t^2} + \rho_0 \frac{\partial^2 u}{\partial x \partial t} = 0
```

同様にして運動方程式の両辺を$x$で偏微分します。

```math
\frac{\partial^2 p}{\partial x^2} + \rho_0 \frac{\partial^2 u}{\partial x \partial t} = 0
```

それぞれの2項目が等しいとおき、更に状態方程式$p=c^2 \rho$を代入して、

```math
\frac{\partial^2 p}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0
```

となり、一次元の波動方程式が導けました。
圧力、密度、粒子速度は線型関係にあるので、上式の圧力を密度、粒子速度に置き換えても成り立ちます。
二次元、三次元の導出も似たような手順でできます。