超音波流体屋のプログラム備忘録
分布点音源法(DPSM)
最終更新：
usapfrog
-
view
分布点音源法(Distributed Point Source Method, DPSM)とは、線形(音)場の解析手法の一つである。
手法としては、球面波の重ね合わせで音場を再現する点音源法に、境界要素法的な音場散乱の勘定が加わったものである。
複雑な定式化を有する境界要素法に比較して、原理・背景の理解が平易な計算手法であるので、
境界要素法の入門代わりに使うのも良いかもしれない。
手法としては、球面波の重ね合わせで音場を再現する点音源法に、境界要素法的な音場散乱の勘定が加わったものである。
複雑な定式化を有する境界要素法に比較して、原理・背景の理解が平易な計算手法であるので、
境界要素法の入門代わりに使うのも良いかもしれない。
多分2015年の現時点で日本語の教科書は無いので、興味があればページ末の英語の教科書を参照されたい。
このページでは流体中の線形音場について記述するが、
固体内の音場伝播やら線形空間微分方程式なら解けることになっているので必要に応じてページ末の教科書を参照すること。
このページでは流体中の線形音場について記述するが、
固体内の音場伝播やら線形空間微分方程式なら解けることになっているので必要に応じてページ末の教科書を参照すること。
本ページではMatlabまたはOctaveで実行可能なサンプルプログラムを付す。
https://github.com/usairpump/dpsmからダウンロードされたい。
https://github.com/usairpump/dpsmからダウンロードされたい。
ホイヘンスの原理・レイリー積分
速度$v$(または加速度$\alpha$)で振動する振動子から放射される音場(調和振動音場$e^{j\omega t}$)は、
散乱体がない場合、波動方程式を解く必要なく、
振動板$S$上からの数値積分(球面波の重ね合わせ)で算出することができる。
散乱体がない場合、波動方程式を解く必要なく、
振動板$S$上からの数値積分(球面波の重ね合わせ)で算出することができる。
振動板上の点の位置ベクトルを${\bf a}(x,y,z)$、求めたい空間上の点を${\bf b}(x,y,z)$とすると、音圧$p$は下式のとおり。
$\L p({\bf b})= \rho \int \alpha({\bf a}) \frac{\exp(-jkr)}{2\pi r} dS
=\frac{j\omega\rho}{2\pi} \int v({\bf a}) \frac{\exp(-jkr)}{r} dS$
$\L r = |{\bf b} - {\bf a} | $
$\L p({\bf b})= \rho \int \alpha({\bf a}) \frac{\exp(-jkr)}{2\pi r} dS
=\frac{j\omega\rho}{2\pi} \int v({\bf a}) \frac{\exp(-jkr)}{r} dS$
$\L r = |{\bf b} - {\bf a} | $
二次元の場合は下記のようにハンケル関数を使用する。
$\L p({\bf b})= \rho \int \alpha({\bf a}) \frac{-j}{2} H_0(kr) dS
=\frac{\omega\rho}{2} \int v({\bf a}) H_0(kr) dS$
$\L G_{ij} = H_0(kr) = J_0(kr) - jY_0(kr)$
$\L p({\bf b})= \rho \int \alpha({\bf a}) \frac{-j}{2} H_0(kr) dS
=\frac{\omega\rho}{2} \int v({\bf a}) H_0(kr) dS$
$\L G_{ij} = H_0(kr) = J_0(kr) - jY_0(kr)$
数値的にも振動板表面を適当な微小領域に分割して、結果を合計すれば良い。
例題01
空気中(音速 340 m/s, 密度1.29 kg/m3)で
20kHz, 1 m/sで振動する半径100 mmの放射源からの音圧分布を算出する。
ソースファイル:q01_rayleigh_integral.m
20kHz, 1 m/sで振動する半径100 mmの放射源からの音圧分布を算出する。
ソースファイル:q01_rayleigh_integral.m
最大振幅: 943 Pa
点音源法
微小面積$dS$を点音源とみなしても計算する式は同じである。
ただ、円の分割などで担当する面積ができるだけ同じになるようにコード上で配慮する必要がある。
$\L p({\bf b}_j)= \sum_i A_i G_{ij}$
$\L A_i = \frac{j\omega\rho}{2\pi} v({\bf a}_i) dS,\ G_{ij} = \frac{\exp(-jkr_{ij})}{r_{ij}}$
$G_{ij}$はグリーン関数と呼ばれる。
ただ、円の分割などで担当する面積ができるだけ同じになるようにコード上で配慮する必要がある。
$\L p({\bf b}_j)= \sum_i A_i G_{ij}$
$\L A_i = \frac{j\omega\rho}{2\pi} v({\bf a}_i) dS,\ G_{ij} = \frac{\exp(-jkr_{ij})}{r_{ij}}$
$G_{ij}$はグリーン関数と呼ばれる。
例題02
最大振幅: 761 Pa
分布自体は例題01と遜色ないが、振幅が0.6倍程度に減じる欠点がある【1】。
分布自体は例題01と遜色ないが、振幅が0.6倍程度に減じる欠点がある【1】。
未知音源
レイリー積分では振動板の前に散乱体、いわば二次音源がある場合この影響を考慮できない。
そこで、DPSMやBEMは元々の一次音源も二次音源を一度振幅未知の音源とみなす。
その後、一次音源の振動や二次音源の性状から決定される境界条件を満たすように、
連立方程式を解いて音源振幅を決定する。
音源振幅が分かれば、点音源法により音場が計算できるというわけである。
そこで、DPSMやBEMは元々の一次音源も二次音源を一度振幅未知の音源とみなす。
その後、一次音源の振動や二次音源の性状から決定される境界条件を満たすように、
連立方程式を解いて音源振幅を決定する。
音源振幅が分かれば、点音源法により音場が計算できるというわけである。
点音源法ではあまり利用されないが、点音源法式に粒子速度${\bf u}$が算出できることが以降で重要になってくる。
$\L u({\bf b})= \frac{-\nabla p}{j\omega\rho} = -\frac{1}{j\omega\rho} \sum_i \left[ \frac{\partial r}{\partial {\bf x}} \frac{\partial}{\partial r} A_i G_i \right] = \sum_i A_i {\bf M}_{ij} $
$\L {\bf M}_{ij} = \frac{({\bf r}/r)}{j\omega\rho} \left(jk +\frac{1}{r} \right) \frac{\exp(-jkr)}{r}$
$\L u({\bf b})= \frac{-\nabla p}{j\omega\rho} = -\frac{1}{j\omega\rho} \sum_i \left[ \frac{\partial r}{\partial {\bf x}} \frac{\partial}{\partial r} A_i G_i \right] = \sum_i A_i {\bf M}_{ij} $
$\L {\bf M}_{ij} = \frac{({\bf r}/r)}{j\omega\rho} \left(jk +\frac{1}{r} \right) \frac{\exp(-jkr)}{r}$
また、二次元空間で解析する場合はグリーン関数はハンケル関数になる。
$\L G_{ij} = H_0(kr) = J_0(kr) - jY_0(kr)$
$\L {\bf M}_{ij} = \frac{({\bf r}/r)}{j\omega\rho} kH_1(kr)$
$\L G_{ij} = H_0(kr) = J_0(kr) - jY_0(kr)$
$\L {\bf M}_{ij} = \frac{({\bf r}/r)}{j\omega\rho} kH_1(kr)$
分布点音源
DPSMとBEMで異なるのはここから先の取り扱いである。
グリーン関数を扱う上で問題となるのは、単純に扱うと自音源が自身の界面に形成する音場が無限大になってしまう点である。
BEMは界面をメッシュ分割し数値積分することによりこれを処理しているが、
DPSMは音源と界面をずらすことで解決する。
DPSMは点音源を半径$r_s$の小球と近似し、界面から$r_s$だけ後退させた位置に配置する。
すると球の表面が界面に接する形で並ぶこととなる。
グリーン関数を扱う上で問題となるのは、単純に扱うと自音源が自身の界面に形成する音場が無限大になってしまう点である。
BEMは界面をメッシュ分割し数値積分することによりこれを処理しているが、
DPSMは音源と界面をずらすことで解決する。
DPSMは点音源を半径$r_s$の小球と近似し、界面から$r_s$だけ後退させた位置に配置する。
すると球の表面が界面に接する形で並ぶこととなる。
※ところで、教科書には隣接する球は、上図のように接して並べろと書いてあるが、
先の低振幅問題を鑑みて、下図のように両者が貫通するくらいの配置で並べるほうが個人的には良いと思っている。
(dSの定義が曖昧になるが、未知音源の振幅に換算されるので問題ないはず)
先の低振幅問題を鑑みて、下図のように両者が貫通するくらいの配置で並べるほうが個人的には良いと思っている。
(dSの定義が曖昧になるが、未知音源の振幅に換算されるので問題ないはず)
音源振幅の決定
音源の数(N個)だけある未知音源振幅$A_i$を放置したまま、
音源群が音源$j$の近傍界面に形成する法線方向粒子速度は下式の通り。${\bf n}_j$が放射方向向き法線ベクトル。
$\L u_j =\sum_i A_i M_{ij} = v_0, \ \ (M_{ij} = {\bf n}_j\cdot {\bf M}_{ij}) $
音源$j$が主音源である場合、この法線方向粒子速度$v_0$は既知である。つまり、方程式がひとつできる。
音源群が音源$j$の近傍界面に形成する法線方向粒子速度は下式の通り。${\bf n}_j$が放射方向向き法線ベクトル。
$\L u_j =\sum_i A_i M_{ij} = v_0, \ \ (M_{ij} = {\bf n}_j\cdot {\bf M}_{ij}) $
音源$j$が主音源である場合、この法線方向粒子速度$v_0$は既知である。つまり、方程式がひとつできる。
音源$j$が二次音源の場合、下記のように境界条件から方程式ができる。
- 音源$j$が音響的に剛壁とみなせる場合なら、$\L u_j =\sum_i A_i M_{ij} = 0 $
- 媒質が水などで音源$j$が気泡のような圧力が0になる場合なら、 $\L p_j =\sum_i A_i G_{ij} = 0 $
あとは、音源の数だけ方程式を用意すればよい。
音源の数だけ近傍界面があるため何ら問題なく数を揃えることができるであろう。
音源の数だけ近傍界面があるため何ら問題なく数を揃えることができるであろう。
境界条件の集合を$\L \{V_j\}=\{ v_0, \cdots, v_0, 0, \cdots, 0 \}$ (音源面で振動速度, 反射面では0)とすると、
境界条件に関して下記の行列表現が可能である。
$\L \left[ {\cal M}_{ij} \right] \{ A_i \} = \{ V_j \} $
この逆行列演算を実施すれば、音源振幅が既知となるのであとは適宜場の計算に使用すれば良い。
境界条件に関して下記の行列表現が可能である。
$\L \left[ {\cal M}_{ij} \right] \{ A_i \} = \{ V_j \} $
この逆行列演算を実施すれば、音源振幅が既知となるのであとは適宜場の計算に使用すれば良い。
行列の各要素は境界条件によって適切なグリーン関数を選択する。
$\L {\cal_M}_{ij} = M_{ij} \ \ or \ \ G_{ij} $
$\L {\cal_M}_{ij} = M_{ij} \ \ or \ \ G_{ij} $
例題03:分布点音源法による円盤からの放射
例題02と同内容の音場をDPSMで計算する。
$\L \left[ M_{ij} \right] \{ A_i \} = v_0 \{\bf 1 \} $
$\L \left[ M_{ij} \right] \{ A_i \} = v_0 \{\bf 1 \} $
ソースファイル:q03_dpsm_disc.m
前述の円盤音源配置関数:def_circ3.mとともに使用すること。
前述の円盤音源配置関数:def_circ3.mとともに使用すること。
青い点群が放射表面(例題02の音源と同一)、緑の点群が未知音源群。
最大振幅: 830 Pa
例題04:二円盤間の定在波
反射体のある計算例として、半径50 mmの円盤を
距離2波長の位置に二枚付きあわせて配置し、下面側を1 m/sで駆動、
上面は単純な剛壁である場合の円盤間の定在波音場を算出する。
下面のの主音源(l), 上面の二次音源(u)を用いて音源決定の行列は下記のようになる。
$\L \left[ \begin{array}{cc} M_{l\rightarrow l} & M_{u\rightarrow l} \\ M_{l\rightarrow u} & M_{u\rightarrow u} \end{array} \right] \left\{ \begin{array}{c} {\bf A}_l \\ {\bf A}_u \end{array} \right\} = \left\{ \begin{array}{c} {\bf V_0} \\ {\bf 0} \end{array} \right\} $
ソースコード:q04_dpsm_2disc.m
距離2波長の位置に二枚付きあわせて配置し、下面側を1 m/sで駆動、
上面は単純な剛壁である場合の円盤間の定在波音場を算出する。
下面のの主音源(l), 上面の二次音源(u)を用いて音源決定の行列は下記のようになる。
$\L \left[ \begin{array}{cc} M_{l\rightarrow l} & M_{u\rightarrow l} \\ M_{l\rightarrow u} & M_{u\rightarrow u} \end{array} \right] \left\{ \begin{array}{c} {\bf A}_l \\ {\bf A}_u \end{array} \right\} = \left\{ \begin{array}{c} {\bf V_0} \\ {\bf 0} \end{array} \right\} $
ソースコード:q04_dpsm_2disc.m
青い点群が放射表面、緑の点群が未知音源群。
最大振幅: 10.3 kPa
FEMとの比較検証結果を準備中
例題05:反射と屈折
音波が境界に気液境界に入射する場合、気相側用の音源と液相側の音源が必要となる。
そして異なる流体面での音の透過には圧力共通と法線速度共通の境界条件が2本必要になる。
そして異なる流体面での音の透過には圧力共通と法線速度共通の境界条件が2本必要になる。
気相で放射された音場が気液界面で反射・屈折し液相に入射する二次元音場を計算する。
気相の主音源(s), 界面から気相に放射する二次音源(a), 界面から液相に放射する二次音源(w)を用いて音源決定の行列は下記のようになる。
$\L \left[ \begin{array}{cc} M_{s\rightarrow s} & M_{a\rightarrow s} & O \\ M_{s\rightarrow a} & M_{a\rightarrow a} & -M_{w\rightarrow a} \\ G_{s\rightarrow w} & G_{a\rightarrow w} & -G_{w\rightarrow w} \end{array} \right] \left\{ \begin{array}{c} {\bf A}_s \\ {\bf A}_a \\ {\bf A}_w \end{array} \right\} = \left\{ \begin{array}{c} {\bf V_0} \\ {\bf 0} \\ {\bf 0} \end{array} \right\} $
気相の主音源(s), 界面から気相に放射する二次音源(a), 界面から液相に放射する二次音源(w)を用いて音源決定の行列は下記のようになる。
$\L \left[ \begin{array}{cc} M_{s\rightarrow s} & M_{a\rightarrow s} & O \\ M_{s\rightarrow a} & M_{a\rightarrow a} & -M_{w\rightarrow a} \\ G_{s\rightarrow w} & G_{a\rightarrow w} & -G_{w\rightarrow w} \end{array} \right] \left\{ \begin{array}{c} {\bf A}_s \\ {\bf A}_a \\ {\bf A}_w \end{array} \right\} = \left\{ \begin{array}{c} {\bf V_0} \\ {\bf 0} \\ {\bf 0} \end{array} \right\} $
二次元空間、20 kHz, 幅60 mm, 1m/sの振動子から4波長程度先の気液界面に10度傾けて音波を入射する。
気相は音速340 m/s,密度1.3 kg/m3, 液相は音速1500 m/s, 密度 998 kg/m3。
気液界面は16波長程度確保する。
気相は音速340 m/s,密度1.3 kg/m3, 液相は音速1500 m/s, 密度 998 kg/m3。
気液界面は16波長程度確保する。
音源配置。界面の両側に気液音源が並んでいるのが特徴的。界面の青点はダブって配置されている。
メモ: w側界面の法線ベクトルはw音源の定義にのみ使用し、Gの算出では使用されない 法線粒子速度共通の境界条件はw音源からの音波もa側の法線ベクトルを用いて算出される
M行列の絶対値の対数をコンター表示した結果。上の部分行列の分割具合がよくわかる。
音圧分布。y=0の界面で連続・微分可能が保持されている。
利点
- 有限要素法や有限差分法に比べて空間を離散化する必要がない。単純に計算節点が三乗 (空間) から二乗 (表面) に減少する上、音源が小さければ格安の計算負荷で広い空間音場を知ることができる。
- 開領域をそのまま扱える。吸収境界条件の適用やPMLなど減衰で吸収させる必要がない。
- 遠距離場で精度が低下しない
- 離散化されているのは (音) 源のみであるため解析に必要ない点の (音) 場を計算する必要がない。また、必要とあらばどんな細かさででも空間の (音) 場を計算することができる。
- (音)場の空間微分が何回でもできる。物理量の微分がメッシュによる精度落ち無しに何階でも手に入る。(空間格子など数値差分で計算する手法は、数値微分演算自体がノイズ発生源となる。)
欠点
以下の場合は有限差分法・有限要素法を使うべき。
- 密行列演算によるメモリ枯渇。疎行列 (対角成分とその周囲のみ非零の行列) を扱う有限要素法と比較して、使用メモリが多い傾向にある。このため、計算の可否は計算時間ではなく使用メモリ制限により依存すると言っていい。例としては 有限要素法なら現行50万節点程度では非力なPCでも計算できるのに対して、分布点音源法では数Gメモリかつ流体音場の解析では6000音源あたりが限度である。
- 閉領域では使用メモリが激増する。境界には波長に対して最低3音源は配置する必要がある。
- 波動方程式が一様でない不均質媒質を扱えない。
- 非線形方程式は解けない。
境界要素法との相違
- 境界要素法は境界を要素分割し、定式化が複雑な境界積分方程式のもと要素内で数値積分を要する。
メッシュ生成の手間があるが、計算精度に関する評価は確立されている。
- 分布点音源法は点音源を境界(のちょっと後方)に並べる。
背景定式化は幾分か容易であるが、精度の方は点音源法相応なので多少の誤差がある。
教科書
リンクはAmazon
【1】 D. Placko, T. Kundu 編 "DPSM for Modeling Engineering Problems" (Wiley, New York, 2007).
【2】 T. Kundu 編 "Ultrasonic Nondestructive Evaluation: Engineering and Biological Material Characterization" (CRC Press, Boca Ration, FL, 2004).
【1】 D. Placko, T. Kundu 編 "DPSM for Modeling Engineering Problems" (Wiley, New York, 2007).
【2】 T. Kundu 編 "Ultrasonic Nondestructive Evaluation: Engineering and Biological Material Characterization" (CRC Press, Boca Ration, FL, 2004).
