OpenAcoustics CookBook内のコードについて

[sakakibara]

お世話になります。

FEM、BEMで閉管路内で定常状態の音圧・粒子速度の分布を比較してみようと考えています。
OpenAcoustics　CookBook内のコードを利用して計算を試みていますが、判らない点をご教授頂けたらと思います。

FEM、BEMともに、速度ポテンシャルの計算結果までしか出されていないように見えますが、
何か意図があってのことなのでしょうか？
速度ポテンシャルを音圧に換算するには、時間微分（2πｆを掛ける）するので良いと思うのですが、
粒子速度に換算する場合は、2点間のレベル差と距離から算出する事になると思います。
3方向のベクトル算出になると思いますが、具体的にはどの様に算出するべきなのでしょうか？
メッシュ間隔などもバラバラですので、どの様にすればいいのか判りません。

また、ＢＥＭのコードがうまく動いていない気がします。
使用するPythonのバージョンなどで、不具合などが発生するようなものなのでしょうか？
OpenAcoustics　CookBookで使用されていました、Pythonのバージョンを教えていただけませんでしょうか？
当方、Ver2.6.6 にて実施しています。
よろしくお願いいたします。

[Hisaharu SUZUKI]

榊原さん

はじめまして。日本エヴィクサーの鈴木と申します。
主にFEM, FDTDの部分を担当しております。

> FEM、BEMで閉管路内で定常状態の音圧・粒子速度の分布を比較してみようと考えています。
> OpenAcoustics　CookBook内のコードを利用して計算を試みていますが、判らない点をご教授頂けたらと思います。

とりあえず、FEMの部分は了解です！

> FEM、BEMともに、速度ポテンシャルの計算結果までしか出されていないように見えますが、
> 何か意図があってのことなのでしょうか？

説明が簡便なので、速度ポテンシャルにしていました。

> 速度ポテンシャルを音圧に換算するには、時間微分（2πｆを掛ける）するので良いと思うのですが、

おっしゃるとおり、時間微分すれば計算できます。

> 粒子速度に換算する場合は、2点間のレベル差と距離から算出する事になると思います。
> 3方向のベクトル算出になると思いますが、具体的にはどの様に算出するべきなのでしょうか？
> メッシュ間隔などもバラバラですので、どの様にすればいいのか判りません。

粒子速度は、FEMの場合は要素自体を空間微分して
計算することもできると思いますし、
近傍の2点で空間差分をとって粒子速度値とすることもできます。

要素内の音圧の値は、補間して計算しますが
その辺の理論説明もコードも省いてしまっているので
ちょっと大変かも知れませんね。
資料とコードを検索してみますので、見つかったらお送りいたします。

BEMは得意でないのであれですが、
解析的に微分した式を離散化もできると思いますし
上記の近傍2点の計算でもできると思います。
近傍2点の計算はBEMのほうがやりやすそうですね。

補足があれば、お願いいたします→みなさん

先日著者で集まったのですが、粒子速度が必要かどうかは
けっこう議論が分かれていました。
個人的には音響インテンシティの計算とかしたいひともいると思うので
あってもいいかなと思っています。

> また、ＢＥＭのコードがうまく動いていない気がします。
> 使用するPythonのバージョンなどで、不具合などが発生するようなものなのでしょうか？
> OpenAcoustics　CookBookで使用されていました、Pythonのバージョンを教えていただけませんでしょうか？
> 当方、Ver2.6.6 にて実施しています。
> よろしくお願いいたします。

エラーメッセージ等お送りいただくことはできますか？
私のほうでも、後で検証してみます。

[Takuya OSHIMA]

大嶋です。

> 補足があれば、お願いいたします→みなさん

FEMであれば、勾配の計算は汎用のポストプロセッサでも出来てしまう、という
のはありますね。

大嶋

[sakakibara]

鈴木様、大嶋先生

　早速のご返答ありがとうございます。
この掲示板にふさわしい内容かどうかわからなかったのですが、
いきなりややこしい質問をしてしまったのかもしれません。
申し訳ありません。


あまり解析のことがわかっていませんで、ポスト処理で勾配計算ができるというなども知りませんでした。
ポストプログラムに持っていくファイルフォーマットも、まだよく判っていないのですが・・・・（汗

著者の方々の意見が分かれている粒子速度の有無に関しては、「可能であれば機能として欲しい」と考えます。
やはり、そのうち音響インテンシティも計算して、ポスト処理でインテンシティの表示などができたら良いな等と思ったりしてます。

BEMのプログラムは、もう少し見直してみます。
当方のミスで書き間違いなどがある場合は、お手数を煩わせるのは申し訳ありませんので。

また幼稚な質問を含め、色々とご質問させていただくこともあるかと思いますが、よろしくお願いいたします。

[Hisaharu SUZUKI]

榊原様

鈴木です。

> 　早速のご返答ありがとうございます。
> この掲示板にふさわしい内容かどうかわからなかったのですが、
> いきなりややこしい質問をしてしまったのかもしれません。
> 申し訳ありません。

とんでもないです。

> あまり解析のことがわかっていませんで、ポスト処理で勾配計算ができるというなども知りませんでした。
> ポストプログラムに持っていくファイルフォーマットも、まだよく判っていないのですが・・・・（汗

僕も知りませんでした(汗
ParaViewでできると思っていいのでしょうか？→大嶋さん

> 著者の方々の意見が分かれている粒子速度の有無に関しては、「可能であれば機能として欲しい」と考えます。
> やはり、そのうち音響インテンシティも計算して、ポスト処理でインテンシティの表示などができたら良いな等と思ったりしてます。

FEMだと何個か関数を追加すればできると思うので
検討してみます。できたらアップしますね。

> BEMのプログラムは、もう少し見直してみます。
> 当方のミスで書き間違いなどがある場合は、お手数を煩わせるのは申し訳ありませんので。

さきほど試してみたのですが

WindowsXP
Python 2.5.4
matplotlib 1.0.1
numpy 1.5.1
scipy 0.9.0

の組み合わせで、私の環境だと動作していました。

もしエラーが出たら、下記のようなコードで
バージョンを確認していただくと何かわかるかも知れません。

import numpy,scipy,matplotlib
print numpy.__version__
print scipy.__version__
print matplotlib.__version__

> また幼稚な質問を含め、色々とご質問させていただくこともあるかと思いますが、よろしくお願いいたします。

お待ちしています。

[ISHIZUKA Takashi]

榊原様

BEM担当の石塚と申します。
宜しくお願い致します。

我々のプロジェクトに興味をお持ち頂きありがとうございます。

【BEMコードの動作】
BEMコードの動作確認済み環境ですが、
鈴木さんのWindowsに加え、
当方では、

Linux（Ubuntu10.04）
python 2.6.5
numpy 1.3.0
scipy 0.7.0
matplotlib 0.99.1.2

で動作しております。
また、MacOS上でも音線法担当の方が動作を確認していますが、
私の方では、バージョン情報までは把握しておりません。

鈴木さんのコメントの通り、エラーコードを頂ければ
何かお役に立てるかも知れません。

【粒子速度の計算】
サンプルコードはわかり易さを最優先していますので、
速度ポテンシャルのみを計算しています。

粒子速度は、
近接２点間の音圧差から求めても良いと思います。
（私自身はやったことが無いので、どの程度近接させるか等のノウハウはありませんが。）

また、連立方程式の解（境界表面上の速度ポテンシャル）を
境界積分方程式を方向微分した形式（NDF型境界積分方程式）
に代入すれば、１点の計算で粒子速度が求まります。

サンプルコードには実装していませんが、
実用コード段階にあたるSimpleAPIには、
粒子速度や音響インテンシティの計算機能も実装したいと考えています。
公開にはまだ少し時間がかかりますが、ご期待頂ければ幸いです。

その他、ご不明な点やご要望がありましたら、
当フォーラムにあげて頂ければ、可能な限り対応させて頂きます。
（タイムラグはご容赦下さい。）

[Takuya OSHIMA]

大嶋です。

> 僕も知りませんでした(汗
> ParaViewでできると思っていいのでしょうか？→大嶋さん

ParaView、Gmshともできます。
ただ、複素数の扱いは工夫が必要かもしれません。

大嶋

[sakakibara]

大嶋先生殿、鈴木様、石塚様

お世話になります。
質問しておいて、ご返答が遅くなってすみません。
サンプルコードは、あくまでサンプルコードで、実用コードはSimpleAPIを作成することになったのですね。
あとで、Code　ＡｓｔｅｒとElmerの計算結果と比較していた
「各種オープンソースFEMコードによる周波数領域音響伝搬解析精度のベンチマークテスト:Code_Aster・Elmer・
OpenAcoustics」
を見させていただいて1st　STEPだったんですね。

先は長そうですが、とりあえずサンプルコードの件に着いてご返答させていただきます。
現在使用しているPythonのバージョンをご連絡いたします。
Python　2.6.6
numpy　1.6.0
scipy　　0.9.0
matplotlib　1.0.1
以上を使用しております。

BEMがうまく動作していない点に関しては、少し原因について進展（推測）がありました。
計算コードのタイプミスを疑い、CookBook内のList3.17にあるモデルを計算して見ました。

lc = 0.49;（←簡単化のために0.1より変更しています）
// Definition of points
Point(1) = {0.5, 0.5, 0, lc};
Point(2) = {-0.5, 0.5, 0, lc};
Point(3) = {-0.5, -0.5, 0, lc};
Point(4) = {0.5, -0.5, 0, lc};

Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};

Physical Line("rigid") = {1,2,4};
Physical Line("absorbing") = {3};
Line Loop(1) = {1,2,3,4};

上記のモデルで計算すると、例題と同じようになるのは確認させていただきました。
簡単化のために各Lineを0.49mピッチで2分割にして、8要素になります。

これを勝手なのですが、例えばLine(3)をrigidに変更し、周囲を完全剛壁状態にしたいと思い
Line(3)が吸音要素になっておりますが、剛壁に変更するために
//Physical Line("rigid") = {1,2,4};
//Physical Line("absorbing") = {3};
Physical Line("rigid") = {1,2,3,4};
というコードに変更して計算していました。

すると、なぜか要素数が16になり、以下のようなエラーメッセージが出てきています。
Warning (from warnings module):
File "C:\Documents and Settings\B000\BEM.py", line 33
dr_dn = - (r_x * element['normalVector'][0] + r_y *
element['normalVector'][1]) / abs_r
RuntimeWarning: invalid value encountered in divide

Warning (from warnings module):
File "C:\Documents and Settings\B000\BEM.py", line 17
return special.j1(x) + 1j * special.y1(x)
RuntimeWarning: invalid value encountered in multiply

Warning (from warnings module):
File "C:\Documents and Settings\B000\BEM.py", line 35
return h11(abs_r * waveNumber) * dr_dn
RuntimeWarning: invalid value encountered in multiply

Warning (from warnings module):
File "C:\Documents and Settings\B000\BEM.py", line 13
return special.j0(x) + 1j * special.y0(x)
RuntimeWarning: invalid value encountered in multiply

Warning (from warnings module):
File "C:\Python26\lib\site-packages\scipy\integrate\quadrature.py",
line 58
return (b-a)/2.0*sum(w*func(y,*args),0), None
RuntimeWarning: invalid value encountered in multiply

Warning (from warnings module):
File "C:\Python26\lib\site-packages\scipy\integrate\quadrature.py",
line 168
AccuracyWarning)
AccuracyWarning: maxiter (50) exceeded. Latest difference = -1.#IND00e
+00

現状では、モデルの中に減衰が無いのでエネルギーが無限大になり発散するのでは無いかと思っています。
BEMでは完全反射の閉空間問題は、計算できないのでしょうか？

以上、よろしくお願いいたします。

追記
色々なコードに触れてみようと考えており、Code　Asterも使用すべく、
とりあえずUbuntu＋Salome　MECAをPCにインストールしてみました。
こちらの件に着いても、質問などさせていただいてもよろしいのでしょうか？？？
プロジェクトの見当違いでしたら、申し訳ありません。

[Takuya OSHIMA]

大嶋です。

遅くなりましたが、全周rigidにした場合の不具合報告ありがとうございます。
私どもでも再現し、原因は判明しました。
目下修正中ですので、今しばらくお待ちください。

> 現状では、モデルの中に減衰が無いのでエネルギーが無限大になり発散するのでは無いかと思っています。
> BEMでは完全反射の閉空間問題は、計算できないのでしょうか？

私どものコードではバグのためお恥ずかしながら計算できていませんが、一般
的には全境界が完全反射の場合でも、発散することなく計算できます。物理的
には、音源から音響エネルギーが音場内に流れていかない、という状態になり
ます(そういう理解でいいですよね？＞石塚さん)。

> 色々なコードに触れてみようと考えており、Code　Asterも使用すべく、
> とりあえずUbuntu＋Salome　MECAをPCにインストールしてみました。
> こちらの件に着いても、質問などさせていただいてもよろしいのでしょうか？？？
> プロジェクトの見当違いでしたら、申し訳ありません。

私どももSalome-Mecaは注目していますが、現状はインストール位までしかサー
ベイしていません。一方で手前味噌で恐縮ですが、小生が役員をしております
OpenCAE学会にSalome-Mecaを熱心にスタディしているグループがおりますので、
そちらのGoogleグループに投稿されることをお勧めします。
http://groups.google.com/group/openfoam
「OpenFOAM」というグループ名がついていますが、Salome-Mecaでも遠慮なく質
問して頂いて結構です。

大嶋

[sakakibara]

お世話になります。
原因が判明してよかったです。
全周Rigidで発生するBEMのバグだったのですね・・・・・
今はFDTDを利用して、教科書の内容と比較したりして遊んで（？）いますので、
BEMでも解ける様になりましたら、その内やってみたいと思います。

> OpenCAE学会にSalome-Mecaを熱心にスタディしているグループがおりますので、
> そちらのGoogleグループに投稿されることをお勧めします。http://groups.google.com/group/openfoam
> 「OpenFOAM」というグループ名がついていますが、Salome-Mecaでも遠慮なく質
> 問して頂いて結構です。

情報を頂き有難うございます。
色々と触ってみたのですが、なかなか思うように進まないのが実情で断念しかけていました。
（短期間でやろうとするのがいけないのですが・・・・）
時間をかけて、少しずつ進めていこうと思います。

有難うございました。

[Hisaharu SUZUKI]

鈴木です。

> お世話になります。
> 原因が判明してよかったです。
> 全周Rigidで発生するBEMのバグだったのですね・・・・・
> 今はFDTDを利用して、教科書の内容と比較したりして遊んで（？）いますので、
> BEMでも解ける様になりましたら、その内やってみたいと思います。

BEMというより、私の書いたメッシャクラスのバグでした。とほほ。。。
修正したらまたご連絡いたしますので、少しお時間ください。

[ISHIZUKA Takashi]

石塚です。

> 遅くなりましたが、全周rigidにした場合の不具合報告ありがとうございます。
> 私どもでも再現し、原因は判明しました。
> 目下修正中ですので、今しばらくお待ちください。

メッシャクラスの修正は鈴木さんにお願いするとしまして、

とりあえず全周rigidでBEMコードを動くようにするには、
サンプルコード「sampleBem2.py」の119行目で、
境界条件（名称と値）を設定するときに、
以下のように使われない境界条件（今回はabsorbing）を定義しないようにすれば
計算が最後まで走ります。
（結果の図まで表示されることを確認しております。）

[sampleBem2.py, l.119]
# boundaryCondition = {'rigid':0., 'absorbing':1.}
boundaryCondition = {'rigid':0.}

試して頂ければ幸いです。

様々な境界条件への柔軟な対応は、検討課題とさせて下さい。

>> 現状では、モデルの中に減衰が無いのでエネルギーが無限大になり発散するのでは無いかと思っています。
>> BEMでは完全反射の閉空間問題は、計算できないのでしょうか？
>
> 私どものコードではバグのためお恥ずかしながら計算できていませんが、一般
> 的には全境界が完全反射の場合でも、発散することなく計算できます。物理的
> には、音源から音響エネルギーが音場内に流れていかない、という状態になり
> ます(そういう理解でいいですよね？＞石塚さん)。

そのように理解して頂ければ良いと思います。
以下は補足です。
このコードは、周波数領域の境界要素法ですので、
音源からある周波数の純音が放射されている「定常音場」を解析しております。
私の理解では、全周rigidの場合の「定常音場」とは、
・場の音圧が0の時点から音波が放射される。
・全周rigidでエネルギーを吸収する要素が無い為、
　音波の放射に従い、場に音響エネルギーがどんどん蓄積される。
・場の音響エネルギーが十分に蓄積された状態になると、
　音源からのエネルギーが場に供給されなくなる。
　直感的な言い方をすると、
　音源からエネルギーを流し込もうとしても、場のエネルギーに押し返される。
・この段階で、場の音響エネルギー（音圧も）は平衡状態（定常音場）になる。
と言うことであり、
全周rigidの条件では、境界要素法はこの平衡状態を求めていることになります。

と、言うわけで計算上は解が求まります。

--
ISHIZUKA

[sakakibara]

鈴木様、石塚様

　ご連絡有難うございます。

> とりあえず全周rigidでBEMコードを動くようにするには、

> # boundaryCondition = {'rigid':0., 'absorbing':1.}

そうですね。言われるとおり
#boundaryCondition = {'rigid':0., 'absorbing':0.}
と言う形にしても、全周rigidで計算できるかもしれませんね。
また後で試して見ます。

> 様々な境界条件への柔軟な対応は、検討課題とさせて下さい。
やはりサンプルコードと言う事ですので、どこまで対応するのかが問題になるのでしょうね。
サンプルコードを修正する事でsimpleAPIに反映できるのであれば、あまり問題ないのかもしれませんが、
反映できないのであれば、時間のロスになりますから・・・・お手数をおかけしてすみません。
ふと思ったのですが、サンプルコードのFEMで開空間の計算って可能なのでしょうか？？？

> このコードは、周波数領域の境界要素法ですので、
> 音源からある周波数の純音が放射されている「定常音場」を解析しております。
> 私の理解では、全周rigidの場合の「定常音場」とは、
> ・場の音圧が0の時点から音波が放射される。
> ・全周rigidでエネルギーを吸収する要素が無い為、

＞・・・・・・・・省略

> ・この段階で、場の音響エネルギー（音圧も）は平衡状態（定常音場）になる。
> と言うことであり、
> 全周rigidの条件では、境界要素法はこの平衡状態を求めていることになります。

そうですよね。
平衝状態を求めているのですから、発散するのはおかしいですよね。
とんちんかんな質問をして、お恥ずかしい限りです。
また詳しくご説明していただき、非常に勉強になります。
有難うございました。

今後とも、よろしくお願いいたします。
sakakibara