FEM for Sound Field Analysis

有限要素法による音場解析

Last Update:2009年05月28日

目次
概要
波動方程式とHelmholtz方程式
境界条件
有限要素法離散化
- Helmholtz方程式の弱形式化と有限要素法離散化
  - 弱形式化
  - 有限要素法離散化
- 吸収境界条件の弱形式化と有限要素法離散化
ヘルムホルツ方程式における連立一次方程式の反復解法
- COCG法のアルゴリズム
  - アルゴリズム(COCG)
- 前処理について
軸対称問題における定式化
参考にしたもの

概要

FEMによる音場解析の基礎について解説する．音の解析には波動方程式を非定常問題として直接解く方法と，Helmholtz方程式を用いて特定の周波数に着目して空間的な振幅の分布について静的に解く方法の２つがある．過答応答に興味がない限り，Helmholtz方程式を用いた方法の方が効率が良いので，ここではHelmhotlz方程式を用いた方法について解説する．

Helmhotlz方程式は，音場を空間分布する振幅の項と注目する周波数の振動成分の積で表現し，波動方程式に変数分離法を適応することで得られる．

波動方程式とHelmholtz方程式

静止した空気を伝わる波

密度の２階時間微分と圧力のラプラシアンの関係

質量保存則
$\frac{\partial\rho}{\partial t} = -\nabla\cdot(\rho\b{v})$

圧縮性流体の運動方程式は

圧縮性流体の運動方程式
$\frac{\partial (\rho \b{v})}{\partial t} + \b{v}\cdot\{\nabla \otimes (\rho \b{v})\} = -\nabla p + \nabla\cdot\b{T}$

のように書くことができる．但し， $\b{T}$ は粘性の効果によるCauchy応力テンソルであるとする．ここで移流と粘性の効果が無視できるとする．つまり左辺第２項と右辺第二項を無視する．

$\frac{\partial (\rho \b{v})}{\partial t} = -\nabla p$

両辺の発散を計算して，微分の順序を入れ替え，質量保存の式を適応する．

$\;\;\;\nabla\cdot\frac{\partial (\rho \b{v})}{\partial t} = -\nabla\cdot(\nabla p)\\\Leftright\frac{\partial}{\partial t}\{\nabla\cdot(\rho\b{v})\} = -\nabla\cdot(\nabla p)\\\Leftright \frac{\partial^2 \rho}{\partial t^2} = -\nabla\cdot(\nabla p)$

圧力と密度の線形近似

状態方程式から，

$p_0=\rho_0 RT_0$

が成り立っている．は空気に対する気体定数とする(簡単のため一般的なモル数を用いた定義ではない)． p_0 $\rho_0$ T_0 はそれぞれ標準圧力，標準密度，標準温度であるとする．断熱圧縮の元では，

$p = \rho^\gamma RT_0$

が成り立つ．但し， $\gamma$ は比熱比である．標準状態における密度と圧力に対する応答を調べてみよう．

$\left.\frac{\partial p}{\partial \rho}\right|_{\rho_0} = \gamma \rho^{\gamma-1}_0 RT_0 =\frac{\gamma\rho^{\gamma}_0 RT_0}{\rho_0} = \frac{\gamma p_0}{\rho_0} = \frac{\gamma\rho_0 RT_0}{\rho_0} = \gamma RT_0$

圧力の波動方程式

上の式から密度の微小な変動 $d\rho$ と圧力の微小な変動に次のような関係があることが分かる．

$d\rho = \frac{1}{\gamma RT_0} dp$

この微小な変化を時間に対する変化であるとすると，次の式が成り立つ．

$\frac{\partial^2 \rho}{\partial t^2} = \frac{1}{\gamma RT_0} \frac{\partial^2 p}{\partial t^2}$

$\frac{\partial^2 p}{\partial t^2} = -\gamma RT_0 \nabla\cdot(\nabla p)$

これは圧力に対する波動方程式であり，音速は次のように表される

音速
$c = \sqrt{\gamma RT_0 }$

波動方程式とHelmholtz方程式

次の波動方程式を解くことを考える．

波動方程式
$\nabla\cdot(\nabla U) - \frac{1}{c^2} \frac{\partial^2 U}{\partial t^2} = 0$

変数分離法によるHelmholtz方程式の導出

解が次のように振動数 $\omega$ での単振動の部分と空間的な振幅の分布の積によって分けられると仮定する．

$U(\b{x},t) = u(\b{x}) e^{i\omega t}$

これを以下の式に代入すると

$\nabla\cdot(\nabla u) e^{i\omega t} - \frac{1}{c^2} (-\omega^2) u e^{i\omega t}= 0$

となる．結局

ヘルムホルツ方程式
$\nabla\cdot(\nabla u) + k^2 u = 0$

を満たす．但し，

$k=\frac{\omega}{c}$ は波数(wave number)と呼ばれる.

カーネル関数

次のようなデルタ関数を右辺に持つHelmholtz方程式を考えよう．

$\nabla\cdot(\nabla u) + k^2 u = \delta(r)$

r=|x|

右辺が点対称なので，解も点対称となる．カーネル数が存在するので，境界要素法(BEM)を用いたシミュレーションも可能である．

解は２次元と３次元の場合で異なり，次のようになる．

２次元の時

２次元の場合，もしくは軸対称問題の場合，解析的な解は次のように第一種Hankel関数で書くことができる．

２次元の場合の解析解
$g(r)=\frac{iH^{(1)}_0(kr)}{4}$

ここで出てくるHankel関数とは次のように第一種ベッセル関数 J_0 と複素単位をかけた第２種ベッセル関数 Y_0 を足したものである．

$H^{(1)}_0 = J_0 + iY_0$

解析解のベッセル関数については以下でもう少し詳しく説明した．

ベッセル関数

３次元の時

３次元の場合，解析的な解は次のように第一種Hankel関数で書くことができる．

３次元の場合の解析解
$g(r)=\frac{e^{ikr}}{4\pi r}$

境界条件

全反射境界条件(Hard wall)

波が境界に完全に反射されることを考えよう．つまり，境界がx=0平面であるとし，境界条件の法線ベクトルをｘ軸に，波数ベクトルの境界に接する成分をｙ軸に取る．

図１．境界と波数ベクトルの位置関係

この時音圧の空間分布は次のような式で表される．

$p(x,y,z)=\bar{p}e^{i(\alpha x+\beta y)} + \bar{p} e^{i(-\alpha x+\beta y)}$

音場の法線方向における微分を考える．

$\left.\frac{\partial p}{\partial x}\right|_0 = \bar{p}i\alpha e^{i(\alpha x+\beta y)} - \bar{p}i\alpha e^{i(-\alpha x+\beta y)}\\\;\;\; = \bar{p}i\alpha e^{i\beta y}$\;e^{i\alpha x} - e^{i\alpha x}$\\\;\;\; = \bar{p}i\alpha e^{i\beta y}$1-1$ = 0$

座標系を一般的な物に戻すと次のようになる．

全反射の境界条件
$\frac{\partial p}{\partial n}= 0$

吸収境界条件(Absorbing Boundary Condition)

波が境界を反射せずに通過することを考えよう．境界がx=0平面であるとし，境界条件の法線ベクトルをｘ軸に，波数ベクトルの境界に接する成分をｙ軸に取る．

図２．境界と波数ベクトルの位置関係

このとき，音圧は次のような式で表される．

$p(x,y,z)=\bar{p}e^{i(\alpha x+\beta y)}$

$k = \sqrt{\alpha^2+\beta^2},\;\;\; \theta = \tan^{-1}$\frac{\beta}{\alpha}$$

$\sigma=\sin\theta = \frac{\beta}{k},\;\;\;\alpha = k\cos\theta = k\sqrt{1-\sigma^2}$

$\frac{\partial p}{\partial y} = i\beta p\;\;\;\Leftright\;\;\; \frac{\partial^2 p}{\partial y^2} = i\beta\frac{\partial p}{\partial y}=-\beta^2p \;\;\;\Leftright\;\;\; \beta^2 = -\frac{1}{p}\frac{\partial^2 p}{\partial y^2}$

$\frac{\partial p}{\partial x} = i\alpha p = i k $\sqrt{1-\sigma^2}$p \\\;\;\; \simeq ik$1-\frac{\sigma^2}{2}$p = ik$1+\frac{\beta^2}{2k^2}$p = $ik+\frac{i}{2k}\frac{\partial^2 }{\partial y^2}$p$

$\frac{\partial p}{\partial z}=0$ であるから，

$\frac{\partial p}{\partial x} \simeq \{ik+\frac{i}{2k}$\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}$\}p$

となる．上の式では法線ベクトルをｘ軸に，波数ベクトルの境界に接する成分をｙ軸に取った．これを法線とそれに接する方向微分に直すと次のようになる．

$\frac{\partial p}{\partial n} \simeq $ik+\frac{i}{2k}\Delta_t$p$

但し， $\Delta_t$ は接平面におけるラプラシアンである．数値解析では上の近似式を用いて吸収境界条件を表す．

吸収境界条件
$\frac{\partial p}{\partial n} = $ik+\frac{i}{2k}\Delta_t$p$

波数ベクトルが境界の法線に対して大きな角度を持つ場合はこの近似の誤差は大きくなる．より正確な近似には，より複雑で高次の導関数が必要であり，ここでは触れない．

固体からの音の放射

固体が角速度で周期的に振動しているとすると固体の速度 $\b{v}_s$ は

$\b{v}_s = \bar{\b{v}}e^{iwt}$

となる．この時粘性を無視するとすると，流体の速度の法線方向成分は固体の法線方向の速度と一致するために，次のようになる．

$\b{v}_f\cdot\b{n} = \bar{\b{v}}\cdot\b{n}e^{iwt}$

圧力が流体の速度変化に応じて周期的に変化するとすると次のようになる．

$p = \bar{p} e^{iwt}$

以上を移流と粘性を無視したNavier-Stokes方程式

$\frac{\partial (\rho \b{v})}{\partial t} = -\nabla p$

へ代入して， $\b{n}$ の内積を計算すると次の通りになる

$\rho iw\bar{\b{v}}\cdot\b{n} = - (\nabla \bar{p}\cdot\b{n})$

が成り立つ．

固体表面からの音の放射
$\frac{\partial p}{\partial n} = -\rho iw(\bar{\b{v}}\cdot\b{n})$

有限要素法離散化

Helmholtz方程式の弱形式化と有限要素法離散化

弱形式化

次のような式を弱形式化してみよう．

$-\{\nabla\cdot(\nabla p)+k^2 p\} = f$

この両辺に $\delta p$ を掛けて全領域で積分する

$\int_V -\delta p \{\nabla\cdot(\nabla p) + k^2 p \} dV = \int_V \delta p f dV\;\;\;\forall \delta p$

上式第一項にGreen-Gaussの定理を適応すると，

$\;\;\; \int_S -n\cdot (\delta p\nabla p) dS + \int_V \nabla\delta p \cdot\nabla p - \delta p k^2 p dV = \int_V \delta p f dV\;\;\;\forall \delta p\\\Leftright \int_S -\delta p\frac{\partial p}{\partial n} dS + \int_V \nabla\delta p \cdot\nabla p - \delta p k^2 p dV = \int_V \delta p f dV\;\;\;\forall \delta p$

弱形式のHelmholtz方程式
$\int_V \nabla\delta p \cdot\nabla p - \delta p k^2 p dV = \int_V \delta p f dV + \int_S \delta p\frac{\partial p}{\partial n} dS \;\;\;\forall \delta p$

有限要素法離散化

体積積分に寄与する項のみ離散化を行う（表面積分に関する項は後から足し合わせればよい）．メッシュ分割されているとすると体積分は各要素の和で表すことができる．

$\sum_e\int_{V_e} \nabla\delta p \cdot\nabla p - \delta p k^2 p dv = \sum_e\int_{V_e} \delta p f dv\;\;\;\forall \delta p$

要素内で内挿関数を用いて次のように補間されているとする．

$\delta p = N^p \delta p^p\;\;\;\; p = N^q p^q$

上式に内挿関数の補間式を代入して，式変形すると連立一次方程式が得られる．

$\;\;\;\sum_e\int_{V_e} \nabla(N^p\delta p^p) \cdot\nabla (N^q p^q) - k^2(N^p\delta p^p) (N^q p^q) dV = \sum_e\int_{V_e} (N^p\delta p^p) f dV\;\;\;\forall \delta p\\\Leftright \sum_e \delta p^p \[\int_{V_e} \nabla N^p \cdot\nabla N^q - k^2N^pN^qdV\] p^q = \sum_e\delta p^p\int_{V_e} N^p f dV\;\;\;\forall \delta p\\\Leftright \sum_e $\delta p^p A_e^{pq} p^q $ = \sum_e$\delta p^p f^p$\;\;\;\forall \delta p\\\Leftright \delta p^a$\sum_e A_e^{pq}$ p^b = \delta p^a $\sum_ef^p$\;\;\;\forall \delta p^a\\\Leftright A^{ab}\ p^b = F^a$

但し，要素内での節点の番号づけ $\{p,q\}$ が全体の節点の番号付け $\{a,b\}$ に対応しているとした．

要素剛性行列は以下のようになることがわかる．

要素剛性行列
$A_e^{pq} = \int_{V_e} \nabla N^p \cdot\nabla N^q - k^2N^pN^qdV$

吸収境界条件の弱形式化と有限要素法離散化

弱形式化

境界の積分式

$\int_S \delta p \frac{\partial p}{\partial n} dS$

に境界条件の式

$\frac{\partial p}{\partial n} = $ik+\frac{i}{2k}\Delta_t$p$

を代入して

$\int_S \delta p $ik+\frac{i}{2k}\Delta_t$p ds = \int_S \delta p (ik p)-\frac{i}{2k}(\nabla_t\delta p)\cdot(\nabla_t p) ds + \int_{\partial S}\frac{i}{2k}\delta p \frac{\partial p}{\partial n_t} dt$

以下の議論では，簡単のために $\partial S = \emptyset$ として，右辺第３項を無視する．

有限要素法離散化

表面がメッシュ分割されているとすると，積分は各要素の和で表すことができる．

$\int_S \delta p (ik p)-\frac{i}{2k}(\nabla_t\delta p)\cdot(\nabla_t p) dS = \sum_e \int_{S_e} \delta p (ik p)-\frac{i}{2k}(\nabla_t\delta p)\cdot(\nabla_t p) dS$

表面の要素内で内挿関数を用いて次のように補間されているとする．

$\delta p = N^p \delta p^p\;\;\;\; p = N^q p^q$

これを上に代入すると境界での要素剛性行列が得られる．

$A^{pq}_e = \int_{S_e} ik N^p N^q -\frac{i}{2k}(\nabla N^p)\cdot(\nabla N^q) dS$

吸収境界条件の次数による精度の比較

吸収境界条件の精度によって，どの程度解が変化するのかを以下の２式で確認した．式(1)は１次の吸収境界条件，式(2)は２次の吸収境界条件である．

$\begin{array}{ll}\frac{\partial p}{\partial n} = ikp & \;\;\;\cdots\;\;\; (1)\\ \frac{\partial p}{\partial n} = $ik+\frac{i}{2k}\frac{\partial^2 }{\partial y^2}$p & \;\;\;\cdots\;\;\; (2)\end{array}$

図３の左の図は式(1)を，右の図は式(2)を境界条件として用いた場合の，解の実部をコンターで表示している．

図３．吸収境界条件の精度の違いによる解への影響

１次２次ともに、波源から点対称になっていないのがわかる．これはまだまだ放射境界条件が完全ではなく壁の影響があるからである．低次の放射境界条件は境界と波の進行方向が垂直に近いという仮定に基いているので、画面下隅のように波の進行方向と境界が平行に近い場合は誤差が出ている．とはいえ２次は１次に比べれば改善されていることがわかる．１次では画面の上の方に壁と反射した波との干渉模様が見えるが、２次ではほとんど見られない．

ヘルムホルツ方程式における連立一次方程式の反復解法

有限要素法のような疎行列を解くのに、一番効率が良いのは反復解法だが、係数行列が非エルミートで不定値であるので、難しい

吸収境界条件や粘性を考えると有限要素法離散化によって得られる係数行列は複素数行列となる。この行列は対称であるがエルミートではないので、CG法のような方法を用いることができない。実非対称行列で使うようなBiCG法(複素数の場合はCOCG[5]),CGNR法,BiCGStab法,GMRes法を複素数に拡張した方法を用いて行列を解く。

COCG法のアルゴリズム

COCG(Conjugate orthogonal-conjugate gradient)法はBiCG法の複素数の場合の特殊な場合で，BiCG法のシャドーベクトルを共役なベクトルとしたものである．COCG法のアルゴリズムは次の通りになる．

アルゴリズム(COCG)

1.   Compute $\b{r}_0=\b{b}-\b{A}\b{x}_0$
2.   Set , $\b{p}_0=\b{r}_0$
3.   For $k=0,1,\ldots,m\quad$ Do:
4.      $\alpha_k=\frac{(\bar{\b{r}}_k,\b{r})}{(\bar{\b{A}\b{p}}_k,\b{p}_k)}$
5.      $\b{x}_{k+1}=\b{x}_k+\alpha_k\b{p}_k$
6.      $\b{r}_{k+1}=\b{r}_k-\alpha_k\b{A}\b{p}_k$
7.      $\beta_k=\frac{\bar{(\b{r}}_{k+1},\b{r}_{k+1})}{(\bar{\b{r}}_,\b{r}_k)}$
8.      $\b{p}_{k+1}=\b{r}_{k+1}+\beta_k\b{p}_k$
9. $End \quad Do$

********************************

ここでの記号 $(\;,\;)$ は複素数の内積を表す．つまり， $(a,b)=\bar{\b{a}}^T\b{b}$ である．エルミート行列に対してはCG法を適応したときのアルゴリズムではこの内積を使うが，COCG法では内積の代わりに，成分どうしの掛け算の足し合わせ $(\bar{a},b)=\b{a}^T\b{b}$ を用いていることに注意されたい．

前処理について

また、波数の大きなHelmholtzを離散化した時に得られる係数行列は不定値なので、COCG法やBiCGStab法で安定して解くために、特殊な前処理手法(複素シフトラプラス前処理)を使うこともある。詳しくは[2][3][4]を参考にして頂きたい。

軸対称問題における定式化

以下の図のように軸対称な３次元問題について考える，つまり，形状，境界条件，物性値全てが軸対称な場合であり，解も軸対称となる．この場合は断面における２次元的な問題に帰着でき，３次元解析のコストを大幅に削減することができる．

弱形式化

軸対称な場合について，断面について２次元的にHelmholtz方程式を有限要素法離散化する．Helmholtz方程式の弱形式は以下の通りであった．

$\int_V \nabla\delta p \cdot\nabla p - \delta p k^2 p dV = \int_V \delta p f dV + \int_S \delta p\frac{\partial p}{\partial n} dS \;\;\;\forall \delta p$

積分の変数を以下のように， $\{x,y,z\}\rightarrow \{r,\theta,z\}$ に変換しよう．

$\{\begin{array}{l}x\rightarrow rsin\theta\\ y\rightarrow rcos\theta\\ z\rightarrow z\end{array}$

この場合のヤコビアンの絶対値は

$\det J(r,\theta,z)=\det \[\begin{array} \sin\theta & r \cos\theta & 0 \\ \cos\theta & -r\sin\theta & 0 \\ 0 & 0 & 1 \end{array}\] = r$

となるので， $dxdydz = r drd\theta dz$ となる．また，この場合下の図のように，

$dS = rdsd \theta$

が成り立つので，これを上に代入すると，

$\int_\theta\int_r\int_z r\nabla\delta p \cdot\nabla p - r\delta p k^2 p dzdrd\theta = \int_\theta\int_r\int_z \delta p f rdzdrd\theta + \int_\theta \int_s \delta p\frac{\partial p}{\partial n} rdsd\theta\;\;\;\forall \delta p$

$\theta$ 方向に被積分関数は変化しないので， $\theta$ の積分を取り除くことができる．また，断面の２次元的な積分領域を,断面の境界をのように書くと，下のような軸対称問題についての弱形式化が得られる．

軸対称問題についてのHelmholtz方程式の弱形式
$\int_r\int_z r\nabla\delta p \cdot\nabla p - r\delta p k^2 p dzdr = \int_r\int_z \delta p f rdzdr + \int_s \delta p\frac{\partial p}{\partial n} rds\;\;\;\forall \delta p$

有限要素法離散化

断面の２次元領域がメッシュ分割されているとすると，上式の積分は各要素内での積分を全て足し合わせて得ることができる．つまり，

$\int_{v_e} r\nabla\delta p \cdot\nabla p - r\delta p k^2 p dv = \int_{v_e} \delta p f rdv + \int_{s_e} \delta p\frac{\partial p}{\partial n} rds\;\;\;\forall \delta p$

ここで，体積積分に関係する１項と２項について計算する．，ガラーキン法で変分を行うとして，要素内で次のように補間関数を用いてと $\delta T$ が離散化されているとする．

p = N^q p^q ， $\delta p = N^p \delta p^p$

これらを弱形式に代入すると以下のように，要素剛性行列，全体剛性行列を導出することができる．但し，式変形では要素内での番号付け $\{p,q\}\rightarrow \{a,b\}$ が成り立つとした．

左辺第一項

$\sum_e\int_{v_e} r\nabla N^p\delta p^p\nabla N^qp^q - rk^2N^p\delta p^pN^q p^qdv\\ \Leftright \sum_e$\delta p^p\int_{v_e} r\nabla N^p\nabla N^q - rk^2 N^pN^qdv p^q$\\ \Leftright\delta p^a$\sum_e A_e^{pq} $ p^b\\ \Leftright\delta p^a A^{ab} p^b$

係数行列
$A_e^{pq}=\int_{v_e}r\nabla N^p\nabla N^q - rk^2 N^pN^q dV$

右辺第一項

$\sum_e\int_{v_e} \delta p f rdv\\ \Leftright \sum_e$\delta p^p\int_{v_e}r N^p f dv$\\ \Leftright\delta p^a$\sum_e F^p_e $\\ \Leftright\delta p^a F^a$

これらは要素単位での次のような行列の足し合わせで得られることが分かる．

外力ベクトル
$F_e^p = \int_{V_e} rN^p f dv$

開口端補正を求める例題

軸対称問題のHelmhltz方程式の解析精度の検証のために，開口端補正を求める実験を行った．開口端補正とはパイプの共鳴を考えるときに，パイプの長さの４分の１とかの整数倍のように計算される共鳴する波長よりも少し波長の方が長いところで共鳴が起るということ．パイプのように一端が開いていると，その周りの空気も一緒に振動するので，その影響を加味しなければならないのだ．