(57)【要約】      （修正有）【課題】解析対象領域の場を正確にモデル化することが可能な物理量シミュレーション方法を提供する。【解決手段】シミュレーション対象格子点における第１物理量を、前記シミュレーション対象格子点以外の複数の格子点の第２物理量と第１定数とに基づいて演算し（Ｓ１０４）、シミュレーション対象格子点における第２物理量を、前記シミュレーション対象格子点以外の複数の格子点の第１物理量と第２定数とに基づいて演算し（Ｓ１０８）、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第１媒質の第１定数と第２媒質の第１定数との加重調和平均を用い（Ｓ１０３）、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第１媒質の第２定数と第２媒質の第２定数との加重平均を用いる（Ｓ１０７）。【選択図】図３
2013235506.tif 000002
【特許請求の範囲】
【請求項１】
それぞれ異なる第１定数と第２定数とを有する第１媒質と第２媒質とが混在する解析対象領域に設定された複数の格子点における第１物理量と第２物理量をシミュレーションする物理量シミュレーション方法であって、
シミュレーション対象格子点における第１物理量を、前記シミュレーション対象格子点以外の複数の格子点の第２物理量と第１定数とに基づいて演算する第１物理量演算ステップと、
シミュレーション対象格子点における第２物理量を、前記シミュレーション対象格子点以外の複数の格子点の第１物理量と第２定数とに基づいて演算する第２物理量演算ステップと、を有し、
前記第１物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点のうち、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第１定数として、媒質境界と前記第１媒質に属する格子点との間の距離と、媒質境界と前記第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第１定数と第２媒質の第１定数との加重調和平均を用い、
前記第２物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点のうち、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第２定数として、媒質境界と前記第１媒質に属する格子点との間の距離と、媒質境界と前記第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第２定数と第２媒質の第２定数との加重平均を用いることを特徴とする物理量シミュレーション方法。
【請求項２】
前記第１物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点が、前記シミュレーション対象格子点から等しい距離離れた２つの格子点であることを特徴とする請求項１に記載の物理量シミュレーション方法。
【請求項３】
前記第２物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点が、前記シミュレーション対象格子点から等しい距離離れた２つの格子点であることを特徴とする請求項１又は請求項２に記載の物理量シミュレーション方法。
【請求項４】
シミュレーションがＦＤＴＤ法に基づいて行われることを特徴とする請求項１乃至請求項３のいずれか１項に記載の物理量シミュレーション方法。
【請求項５】
シミュレーションがＷＥ−ＦＤＴＤ法に基づいて行われることを特徴とする請求項１乃至請求項３のいずれか１項に記載の物理量シミュレーション方法。
【請求項６】
前記第１物理量が音圧であり、前記第２物理量が粒子速度であることを特徴とする請求項１乃至請求項５のいずれか１項に記載の物理量シミュレーション方法。
【請求項７】
前記第１係数が体積弾性率であり、前記前記第２係数が密度であることを特徴とする請求項１乃至請求項６のいずれか１項に記載の物理量シミュレーション方法。
【請求項８】
請求項１乃至請求項７のいずれか１項に記載の方法を実行することを特徴とする物理量シミュレーションシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピューター上で音波や電磁波の時間に伴う伝搬をシミュレーションするために好適な物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムに関する。
【背景技術】
【０００２】
従来より、コンピューター上で音波や電磁波の時間に伴う伝搬を波動理論に基づきシミュレーションを行う手法として時間領域差分法（Ｆｉｎｉｔｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｔｉｍｅ  Ｄｏｍａｉｎ  Ｍｅｔｈｏｄ：  以下、ＦＤＴＤ法）が提案されている（例えば、非特許文献１）。
【０００３】
このＦＤＴＤ法では、解析対象領域である場を時間および空間について離散化し格子点を設定し、波動の支配方程式の微分を差分表現して伝搬解析を行う。一般的に用いられるＹｅｅのＦＤＴＤ法では、解析対象領域である場の物理量である音圧と粒子速度を空間的に１／２グリッドずらして配置し、時間更新も音圧と粒子速度を１／２ステップずつ互い違いに計算する（Ｓｔａｇｇｅｒｅｄ−ｇｒｉｄ、非特許文献１参照）。以下、１／２グリッドずらして設定された空間の各点を格子点として表現する。
【０００４】
図５はＦＤＴＤ法に基づくシミュレーションを説明する図である。図５において、ｘ軸の下方に記載された・・・，ｉ−１／２，ｉ，ｉ＋１／２，ｉ＋１，・・・がそれぞれ位置を示しており、それぞれが本明細書では格子点として定義されるものである。また、格子点のうち、○は音圧をシミュレーションするために定義される格子点であり、△は粒子速度をシミュレーションするために定義される格子点である。ｐは音圧を、ｕは粒子速度を示している。
【０００５】
図５では、解析対象領域である場は均質な媒質が想定されており、この媒質の体積弾性率をＫ、密度をρとしている。ＦＤＴＤ法に基づくシミュレーションでは、音圧ｐ及び粒子速度ｕはそれぞれ下式（１）及び（２）によって演算を行う。
【０００６】
【数１】
2013235506.tif 000003
【０００７】
【数２】
2013235506.tif 000004
【０００８】
ここで、Δｘは空間離散化幅であり、Δｔは時間離散化幅である。Δｘ／２が相隣る格子点の間隔である。また、実時間ｔとΔｔとの間にはｔ＝ｎ・Δｔの関係がある。ここで、ｎは離散更新時間として定義されるものである。音圧ｐ及び粒子速度ｕの上付き文字は、この離散更新時間を示している。
【０００９】
式（１）及び図５を参照すると、音圧ｐは、離散更新時間が１つ前の音圧ｐ（第１項）と、着目している音圧定義点（ｘ＝ｉ）を挟む２つの格子点における粒子速度ｕの値と体
積弾性率Ｋなど（第２項）によって決定されることがわかる。
【００１０】
同じように、式（２）及び図５を参照すると、粒子速度ｕは、離散更新時間が１つ前の粒子速度ｕ（第１項）と、着目している粒子速度定義点（ｘ＝ｉ＋１／２）を挟む２つの格子点における音圧ｐの値と密度ρなど（第２項）によって決定されることがわかる。
【００１１】
ところで、解析対象領域に第１の媒質と、これと異なる第２の媒質とが混在する場合、すなわち、体積弾性率がＫ1、密度がρ1である第１媒質と、体積弾性率がＫ2、密度がρ2である第２媒質と、が存在する解析対象領域に対し、ＦＤＴＤ法に基づくシミュレーションを適用する場合には、従来技術では、格子点上に、第１媒質と第２媒質と間の媒質境界を設定し、媒質に関する値を以下のように設定した上でシミュレーションが実行される（非特許文献１参照）。
【００１２】
図６に示すように音圧定義点（ｘ＝ｉ）である格子点に媒質境界を設定した場合においては、その両側の格子点における体積弾性率の調和平均である
【００１３】
【数３】
2013235506.tif 000005
【００１４】
をｘ＝ｉの位置の音圧ｐをシミュレーションする際に利用する体積弾性率とする。
また、図７に示すように粒子速度定義点（ｘ＝ｉ＋１／２）である格子点に媒質境界を設定した場合においては、その両側の格子点における密度の相加平均である
【００１５】
【数４】
2013235506.tif 000006
【００１６】
をｘ＝ｉ＋１／２の位置の粒子速度ｕをシミュレーションする際に利用する密度とする。
【非特許文献１】「ＦＤＴＤ音響解析では媒質間の境界における媒質定数はどう与えればよいのか―ディリクレ・ノイマン混合境界条件に基づく音響数値計算のためのインピーダンス境界の取り扱い―」，日本音響学会建築音響研究会，建築音響研究資料、資料番号ＡＡ２０１１−３６，２０１１年１０月１１日
【発明の概要】
【発明が解決しようとする課題】
【００１７】
従来技術では、媒質境界を設定する位置が格子点上に制限され、解析対象領域の場を正確にモデル化することができない、という問題があった。例えば、壁面が対向しているような解析対象領域をシミュレーションする場合にその間の距離を厳密に与えることができず、モデル化誤差により解析精度が低下する原因となる。
【００１８】
上記のような問題に対処するために、空間離散化幅をより小さくすれば境界位置の誤差をある程度小さくすることができるが、シミュレーションの際、媒質境界を設定する位置については、やはり空間離散化幅の制限を受けてしまうことには変わらない。また、空間離散化幅をより小さくすることは計算容量・時間の増大を招く。さらに、空間離散化幅に合わせて時間離散化幅も小さくする必要があり、これも計算時間の増加につながる。
【００１９】
媒質境界の近傍を小さい空間離散化幅で計算し、より大きい空間離散化幅で計算する均質媒質中と連結するサブグリッド手法（例えば、特開２０１２−１４２４６号）により、計算コストを抑えつつ境界位置の誤差をある程度小さくすることはできる。しかしながら、この場合でも境界位置がサブグリッドの空間離散化幅に制限を受けることには変わらない。
【００２０】
また、媒質境界が設定された格子点をシミュレーションする際、両媒質の体積弾性率の調和平均や、密度の相加平均を用いる従来の手法は、媒質境界が設定された格子点の両側の２つの格子点の値を用いて、２次中心差分で表現するＦＤＴＤ法に対してのみ適用可能である。
【００２１】
すなわち、計算精度を向上するために、それ以上の点の値を用いて空間微分をより高次の差分で表現する手法（例えば、後述する両側４点を用いるＦＤＴＤ（２，４）法）があるが、従来技術では、そのような手法には適用できない、という問題があった。
【００２２】
なお、ＦＤＴＤ（２，４）法を適用するにあたっては、媒質境界の近傍については、両側２つの格子点の値を用いて２次中心差分で計算し、一様媒質部分ではより高次の差分で計算し、２つの計算結果を連結する方法が考えられるが、媒質境界部分での計算精度が低下する。また、この場合も媒質境界の位置はグリッドの離散化幅の制限を受ける、という問題が解消されるわけではない。
【課題を解決するための手段】
【００２３】
この発明は、上記課題を解決するものであって、請求項１に係る発明は、それぞれ異なる第１定数と第２定数とを有する第１媒質と第２媒質とが混在する解析対象領域に設定された複数の格子点における第１物理量と第２物理量をシミュレーションする物理量シミュレーション方法であって、シミュレーション対象格子点における第１物理量を、前記シミュレーション対象格子点以外の複数の格子点の第２物理量と第１定数とに基づいて演算する第１物理量演算ステップと、シミュレーション対象格子点における第２物理量を、前記シミュレーション対象格子点以外の複数の格子点の第１物理量と第２定数とに基づいて演算する第２物理量演算ステップと、を有し、前記第１物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点のうち、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第１定数として、媒質境界と前記第１媒質に属する格子点との間の距離と、媒質境界と前記第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第１定数と第２媒質の第１定数との加重調和平均を用い、前記第２物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点のうち、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第２定数として、媒質境界と前記第１媒質に属する格子点との間の距離と、媒質境界と前記第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第２定数と第２媒質の第２定数との加重平均を用いることを特徴とする。
【００２４】
また、請求項２に係る発明は、請求項１に記載の物理量シミュレーション方法において、前記第１物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点が、前記シミュレーション対象格子点から等しい距離離れた２つの格子点であることを特徴とする。
【００２５】
また、請求項３に係る発明は、請求項１又は請求項２に記載の物理量シミュレーション方法において、前記第２物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点が、前記シミュレーション対象格子点から等しい距離離れた２つの格子点であることを特徴とする。
【００２６】
また、請求項４に係る発明は、請求項１乃至請求項３のいずれか１項に記載の物理量シミュレーション方法において、シミュレーションがＦＤＴＤ法に基づいて行われることを特徴とする。
【００２７】
また、請求項５に係る発明は、請求項１乃至請求項３のいずれか１項に記載の物理量シミュレーション方法において、シミュレーションがＷＥ−ＦＤＴＤ法に基づいて行われることを特徴とする。
【００２８】
また、請求項６に係る発明は、請求項１乃至請求項５のいずれか１項に記載の物理量シミュレーション方法において、前記第１物理量が音圧であり、前記第２物理量が粒子速度であることを特徴とする。
【００２９】
また、請求項７に係る発明は、請求項１乃至請求項６のいずれか１項に記載の物理量シミュレーション方法において、前記第１係数が体積弾性率であり、前記前記第２係数が密度であることを特徴とする。
【００３０】
また、請求項８に係る発明は、請求項１乃至請求項７のいずれか１項に記載の方法を実行することを特徴とする物理量シミュレーションシステムである。
【発明の効果】
【００３１】
本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、演算ステップにおいて、定数として、媒質境界と第１媒質に属する格子点との間の距離と、媒質境界と第２媒質に属する格子点との間の距離と、に応じた適切な平均値を用いることで、複数の格子点のいずれも第１媒質と第２媒質と間の媒質境界上に設けられていなくてもよく、解析対象領域の場を正確にモデル化することが可能となる。
【００３２】
すなわち、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、解析対象領域の場の媒質が不均質な場合においても、空間離散化幅に関わらず媒質境界を正確に設定でき、媒質境界の位置に関するモデル化誤差を無くし計算精度を向上できるようになる。
【００３３】
また、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムを適用するにあたっては、空間離散化幅及び時間離散化幅は特に変更する必要は無く、計算コストは変わらない。
【００３４】
また、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、空間微分をより多くの点の値を用いて高次差分で表現する高精度な手法用いて、不均質媒質中の波動伝搬を解析することが可能になる。
【００３５】
本発明の提供する物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、音圧の波動方程式に差分法を適用したＷＥ−ＦＤＴＤ法（参考文献１）にも適用することが可能である。
【００３６】
本発明の提供する物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムは、音波の伝搬解析を例のみならず、電磁波や振動等の波動一般の解析に適用することが可能である。
【図面の簡単な説明】
【００３７】
【図１】本発明の実施形態に係る物理量シミュレーションシステムを構成するコンピューターの一例を示す図である。
【図２】本発明の実施形態に係る物理量シミュレーションを説明する図である。
【図３】本発明の実施形態に係る物理量シミュレーションのフローチャートの一例を示す図である。
【図４】空間微分を両側の４点の格子点の値を用いて、４次中心差分で表すＦＤＴＤ（２，４）法に本発明を適用した場合を説明する図である。
【図５】ＦＤＴＤ法に基づくシミュレーションを説明する図である。
【図６】従来技術に係るシミュレーションを説明する図である。
【図７】従来技術に係るシミュレーションを説明する図である。
【図８】実施例を説明する図である。
【発明を実施するための形態】
【００３８】
以下、本発明の実施の形態を図面を参照しつつ説明する。図１は本発明の実施形態に係る物理量シミュレーションシステムを構成するコンピューターの一例を示す図である。図１において、１０はシステムバス、１１はＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）、１２はＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）、１３はＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、１４は外部情報機器との通信を司る通信制御部、１５はキーボードコントローラなどの入力制御部、１６は出力制御部、１７は外部記憶装置制御部、１８はキーボード、ポインティングデバイス、マウスなどの入力機器からなる入力部、１９は印刷装置などの出力部、２０はＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）等の外部記憶装置、２１はグラフィック制御部、２２はディスプレイ装置をそれぞれ示している。
【００３９】
図１において、ＣＰＵ１１は、ＲＯＭ１３内のプログラム用ＲＯＭ、或いは、大容量の外部記憶装置２０に記憶されたプログラム等に応じて、外部機器と通信することでデータを検索・取得したり、また、図形、イメージ、文字、表等が混在した出力データの処理を実行したり、更に、外部記憶装置２０に格納されているデータベースの管理を実行したり、などといった演算処理を行うものである。
【００４０】
また、ＣＰＵ１１は、システムバス１０に接続される各デバイスを統括的に制御する。ＲＯＭ１３内のプログラム用ＲＯＭあるいは外部記憶装置２０には、ＣＰＵ１１の制御用の基本プログラムであるオペレーティングシステムプログラム（以下ＯＳ）等が記憶されている。また、ＲＯＭ１３あるいは外部記憶装置２０には出力データ処理等を行う際に使用される各種データが記憶されている。メインメモリーであるＲＡＭ１２は、ＣＰＵ１１の主メモリ、ワークエリア等として機能する。
【００４１】
入力制御部１５は、キーボードや不図示のポインティングデバイスからの入力部１８を制御する。また、出力制御部１６は、プリンタなどの出力部１９の出力制御を行う。
【００４２】
外部記憶装置制御部１７は、ブートプログラム、各種のアプリケーション、フォントデータ、ユーザーファイル、編集ファイル、プリンタドライバ等を記憶するＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）や、或いはフロッピーディスク（ＦＤ）等の外部記憶装置２０へのアクセスを制御する。本発明の物理量シミュレーション方法を実現するシステムプログラムは、上記のような外部記憶装置２０に記憶されている。また、グラフィック制御部２１は、ディスプレイ装置２２に表示する情報を描画処理するための構成である。
【００４３】
また、通信制御部１４は、ネットワークを介して、外部機器と通信を制御するものであり、これによりシステムが必要とするデータを、インターネットやイントラネット上の外部機器が保有するデータベースから取得したり、外部機器に情報を送信したりすることができるように構成される。
【００４４】
外部記憶装置２０には、ＣＰＵ１１の制御プログラムであるオペレーティングシステムプログラム（以下ＯＳ）以外に、本発明の物理量シミュレーションシステムをＣＰＵ１１上で動作させるシステムプログラム、及びこのシステムプログラムで用いるデータなどがインストールされ保存・記憶されている。
【００４５】
本発明の物理量シミュレーション方法を実現するシステムプログラムで利用されるデータとしては、基本的には外部記憶装置２０に保存されていることが想定されているが、場合によっては、これらのデータを、通信制御部１４を介してインターネットやイントラネット上の外部機器から取得するように構成することも可能である。また、本発明の物理量シミュレーション方法を実現するシステムプログラムで利用されるデータを、ＵＳＢメモリやＣＤ、ＤＶＤなどの各種メディアから取得するように構成することもできる。
【００４６】
図２は本発明の実施形態に係る物理量シミュレーションを説明する図である。以下、本明細書においては、ｘ軸方向の１次元を例として説明するが、本発明は３次元のいずれの方向にも拡張することが可能であり、本発明は２次元及び３次元空間の波動解析に対しても有効である。
【００４７】
図２において、ｘ軸の下方に記載された・・・，ｉ−１／２，ｉ，ｉ＋１／２，ｉ＋１，・・・がそれぞれ位置を示しており、それぞれが本明細書では格子点として定義されるものである。また、格子点のうち、○は音圧をシミュレーションするために定義される格子点であり、△は粒子速度をシミュレーションするために定義される格子点である。ｐは音圧を、ｕは粒子速度を示している。
【００４８】
図２に示す解析対象領域は、体積弾性率がＫ1、密度がρ1である第１媒質と、体積弾性率がＫ2、密度がρ2である第２媒質と、が存在するものであり、点線で示すラインが第１媒質と第２媒質との媒質境界である。なお、特許請求の範囲では、体積弾性率を第１定数、密度を第２定数として表現している。すなわち、特許請求の範囲では、
Ｋ1；第１媒質の第１定数
ρ1；第１媒質の第２定数
Ｋ2；第２媒質の第１定数
ρ2；第２媒質の第２定数
である。
【００４９】
また、解析対象領域において定義されている複数の格子点のいずれも第１媒質と第２媒質と間の媒質境界上に設けないような条件が設定されている。本発明では、後述するように、演算ステップにおいて、定数として、媒質境界と第１媒質に属する格子点との間の距離と、媒質境界と第２媒質に属する格子点との間の距離と、に応じた適切な平均値を用いることで、このような条件設定を可能としている。これにより、本発明では、解析対象領域の場を正確にモデル化することが可能となる。
【００５０】
本発明においても、ＦＤＴＤ法に基づくシミュレーションでは、音圧ｐ及び粒子速度ｕはそれぞれ下式（１）及び（２）によって演算を行う。
【００５１】
【数１】
2013235506.tif 000007
【００５２】
【数２】
2013235506.tif 000008
【００５３】
ここで、Δｘは空間離散化幅であり、Δｔは時間離散化幅である。Δｘ／２が相隣る格子点の間隔である。また、実時間ｔとΔｔとの間にはｔ＝ｎ・Δｔの関係がある。ここで、ｎは離散更新時間として定義されるものである。音圧ｐ及び粒子速度ｕの上付き文字は、この離散更新時間を示している。
【００５４】
なお、特許請求の範囲においては、式（１）に基づく演算を「第１物理量演算ステップ」、式（２）に基づく演算を「第２物理量演算ステップ」と表現している。また、特許請求の範囲においては、「第１物理量」を音圧、「第２物理量」を粒子速度として表現している。
【００５５】
式（１）及び図２を参照すると、音圧ｐは、離散更新時間が１つ前の音圧ｐ（第１項）と、着目している音圧定義点（ｘ＝ｉ）を挟む２つの格子点における粒子速度ｕの値と体積弾性率Ｋなど（第２項）によって決定されることがわかる。
【００５６】
式（１）に基づく演算を第１物理量演算ステップで、格子点ｘ＝ｉにおける音圧ｐを演算する際には、上記の第２項に用いる体積弾性率Ｋが必要となるが、格子点ｘ＝ｉ−１／２は第１媒質に属し、格子点ｘ＝ｉ+１／２は第２媒質に属しているために、このときに
用いる体積弾性率をどのように定義するかが問題となる。
【００５７】
本発明では、このような場合、体積弾性率（第１定数）として、媒質境界と第１媒質に属する格子点との間の距離と、媒質境界と前記第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第１定数と第２媒質の第１定数との加重調和平均を用いる。
【００５８】
より具体的には、媒質境界と第１媒質に属する格子点との間の距離をｌu1とし、媒質境界と第２媒質に属する格子点との間の距離をｌu2とした場合に、格子点ｘ＝ｉにおける音圧ｐを演算する際に用いる体積弾性率（第１定数）としては、
【００５９】
【数５】
2013235506.tif 000009
【００６０】
を用いるようにする。ただし、数式を簡略にするために、Δｘ＝ｌu1＋ｌu2としている。
【００６１】
また、式（２）及び図２を参照すると、粒子速度ｕは、離散更新時間が１つ前の粒子速度ｕ（第１項）と、着目している粒子速度定義点（ｘ＝ｉ＋１／２）を挟む２つの格子点における音圧ｐの値と密度ρなど（第２項）によって決定されることがわかる。
【００６２】
式（２）に基づく演算を第２物理量演算ステップで、格子点ｘ＝ｉ+１／２における粒
子速度ｕを演算する際には、上記の第２項に用いる密度ρが必要となるが、格子点ｘ＝ｉは第１媒質に属し、格子点ｘ＝ｉ+１は第２媒質に属しているために、このときに用いる
密度をどのように定義するかが問題となる。
【００６３】
本発明では、このような場合、密度（第２定数）として、媒質境界と第１媒質に属する
格子点との間の距離と、媒質境界と第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第２定数と第２媒質の第２定数との加重平均を用いる。
【００６４】
より具体的には、媒質境界と第１媒質に属する格子点との間の距離をｌp1とし、媒質境界と第２媒質に属する格子点との間の距離をｌp2とした場合に、格子点ｘ＝ｉ＋１／２における粒子速度ｕを演算する際に用いる密度（第２定数）としては、
【００６５】
【数６】
2013235506.tif 000010
【００６６】
を用いるようにする。ただし、数式を簡略にするために、Δｘ＝ｌp1＋ｌp2としている。
【００６７】
媒質境界近傍の格子点における演算において、上記のように体積弾性率（第１定数）と密度（第２定数）とを決定すると、媒質境界におけるディリクレ条件及びノイマン条件を満たし、境界条件を理論的に正確に記述することが可能となる。
【００６８】
本発明においては、媒質境界近傍の格子点については、以上の要領で演算するが、次に、解析対象領域に属する全ての格子点を演算するための処理手順について説明する。図３は本発明の実施形態に係る物理量シミュレーションのフローチャートの一例を示す図である。
【００６９】
図３において、ステップＳ１００で処理が開始されると、続くステップＳ１０１では、ｕ（ｉ−１／２）、ｕ（ｉ＋１／２）が同じ媒質に属するか否かが判定される。
【００７０】
ｕ（ｉ−１／２）、ｕ（ｉ＋１／２）が同じ媒質に属し、ステップＳ１０１における判定がＹＥＳであるときには、ｕ（ｉ−１／２）、ｕ（ｉ＋１／２）が属する媒質に応じて、体積弾性率Ｋ1又はＫ2を選択（ステップＳ１０２）し、ステップＳ１０４で音圧の更新計算（第１物理量演算ステップ）を実行する。
【００７１】
一方、ｕ（ｉ−１／２）、ｕ（ｉ＋１／２）が異なる媒質に属し、ステップＳ１０１における判定がＮＯであるときには、体積弾性率Ｋを式（５）に基づいて演算（ステップＳ１０３）し、ステップＳ１０４で音圧の更新計算（第１物理量演算ステップ）を実行する。
【００７２】
ステップＳ１０５では、ｐ（ｉ）、ｐ（ｉ＋１）が同じ媒質に属するか否かが判定される。
【００７３】
ｐ（ｉ）、ｐ（ｉ＋１）が同じ媒質に属し、ステップＳ１０５における判定がＹＥＳであるときには、ｐ（ｉ）、ｐ（ｉ＋１）が属する媒質に応じて、密度ρ1又はρ2を選択（ステップＳ１０６）し、ステップＳ１０８で粒子速度の更新計算（第２物理量演算ステップ）を実行する。
【００７４】
一方、ｐ（ｉ）、ｐ（ｉ＋１）が異なる媒質に属し、ステップＳ１０５における判定がＮＯであるときには、密度ρを式（６）に基づいて演算（ステップＳ１０７）し、ステップＳ１０８で粒子速度の更新計算（第２物理量演算ステップ）を実行する。
【００７５】
ステップＳ１０９では、離散更新時間ｎを１インクリメントし、ステップＳ１１０では最後の離散更新時間に到達したか否かが判定される。ステップＳ１１０の判定がＮＯであ
る場合にはステップＳ１０１に戻り、判定がＹＥＳである場合には、ステップＳ１１１に進み、処理を終了する。
【００７６】
ステップＳ１０３及びステップＳ１０７の計算は、計算開始前にあらかじめ処理して結果を保持しておけば、本発明の実施による計算量の増加は生じない。
【００７７】
以上をまとめると、本発明は、それぞれ異なる第１定数（体積弾性率）と第２定数（密度）とを有する第１媒質と第２媒質とが混在する解析対象領域に設定された複数の格子点における第１物理量と第２物理量をシミュレーションする物理量シミュレーション方法、物理量シミュレーションシステムであって、シミュレーション対象格子点における第１物理量を、前記シミュレーション対象格子点以外の複数の格子点の第２物理量と第１定数（体積弾性率）とに基づいて演算する第１物理量演算ステップと、シミュレーション対象格子点における第２物理量を、前記シミュレーション対象格子点以外の複数の格子点の第１物理量と第２定数（密度）とに基づいて演算する第２物理量演算ステップと、を有し、前記第１物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点のうち、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第１定数（体積弾性率）として、媒質境界と前記第１媒質に属する格子点との間の距離と、媒質境界と前記第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第１定数（体積弾性率）と第２媒質の第１定数（体積弾性率）との加重調和平均を用い、前記第２物理量演算ステップで用いる前記シミュレーション対象格子点以外の複数の格子点のうち、第１媒質に属する格子点と、第２媒質に属する格子点とが存在する場合には、第２定数（密度）として、媒質境界と前記第１媒質に属する格子点との間の距離と、媒質境界と前記第２媒質に属する格子点との間の距離と、を重み係数とした、第１媒質の第２定数（密度）と第２媒質の第２定数（密度）との加重平均を用いることを特徴としている。
【００７８】
媒質境界近傍の格子点における演算において、上記のように体積弾性率（第１定数）と密度（第２定数）とを決定すると、媒質境界におけるディリクレ条件及びノイマン条件を満たし、境界条件を理論的に正確に記述することが可能となる。
【００７９】
したがって、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、演算ステップにおいて、定数として、媒質境界と第１媒質に属する格子点との間の距離と、媒質境界と第２媒質に属する格子点との間の距離と、に応じた適切な平均値を用いることで、複数の格子点のいずれも第１媒質と第２媒質と間の媒質境界上に設けられていなくてもよく、解析対象領域の場を正確にモデル化することが可能となる。
【００８０】
すなわち、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、解析対象領域の場の媒質が不均質な場合においても、空間離散化幅に関わらず媒質境界を正確に設定でき、媒質境界の位置に関するモデル化誤差を無くし計算精度を向上できるようになる。
【００８１】
また、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムを適用するにあたっては、空間離散化幅及び時間離散化幅は特に変更する必要は無く、計算コストは変わらない。
【００８２】
次に、より多くの格子点の値を用いて空間微分を表現する手法には以下のように適用する。
【００８３】
例えば、空間微分を両側の４点の格子点の値を用いて、４次中心差分で表すＦＤＴＤ（
２，４）法においては、音圧ｐ及び粒子速度ｕは一般にそれぞれ式（７）（第１物理量演算ステップ）及び式（８）（第２物理量演算ステップ）で更新する。
【００８４】
【数７】
2013235506.tif 000011
【００８５】
【数８】
2013235506.tif 000012
【００８６】
ただし、Ａ＝９／８である。
【００８７】
本発明においては、式（７）及び式（８）を、以下のようにΔｘ離れた２点の２次中心差分と３Δｘ離れた２点の２次中心差分に分けて、式（９）及び式（１０）の形で考える。
【００８８】
【数９】
2013235506.tif 000013
【００８９】
【数１０】
2013235506.tif 000014
【００９０】
図４は空間微分を両側の４点の格子点の値を用いて、４次中心差分で表すＦＤＴＤ（２，４）法に本発明を適用した場合を説明する図である。以下、この図を参照して、ＦＤＴＤ（２，４）法における媒質境界近傍における演算で、体積弾性率（第１定数）と密度（第２定数）を決定する際の考え方を説明する。
【００９１】
式（９）に基づく演算を第１物理量演算ステップで、格子点ｘ＝ｉにおける音圧ｐを演算する際には、式（９）の第２項に用いる体積弾性率Ｋ'と、式（９）の第３項に用いる
体積弾性率Ｋ' 'が必要となる。本発明では、Ｋ'とＫ' 'とをそれぞれ下式（１１）及び
（１２）によって決定する。（図４（ａ）参照）
【００９２】
【数１１】
2013235506.tif 000015
【００９３】
【数１２】
2013235506.tif 000016
【００９４】
ただし、体積弾性率Ｋ'を求める際の、媒質境界と第１媒質に属する格子点との間の距
離をｌ'u1とし、媒質境界と第２媒質に属する格子点との間の距離をｌ'u2とする。また、数式を簡略にするために、Δｘ＝ｌ'u1＋ｌ'u2である。また、体積弾性率Ｋ' 'を求める
際の、媒質境界と第１媒質に属する格子点との間の距離をｌ' 'u1とし、媒質境界と第２
媒質に属する格子点との間の距離をｌ' 'u2とする。３Δｘ＝ｌ' 'u1＋ｌ' 'u2である。
【００９５】
また、式（１０）に基づく演算を第２物理量演算ステップで、格子点ｘ＝ｉ＋１／２における粒子速度ｕを演算する際には、式（１０）の第２項に用いる密度ρ'と、式（１０
）の第３項に用いる密度ρ' 'が必要となる。本発明では、ρ'とρ' 'とをそれぞれ下式
（１３）及び（１４）によって決定する。（図４（ａ）参照）
【００９６】
【数１３】
2013235506.tif 000017
【００９７】
【数１４】
2013235506.tif 000018
【００９８】
ただし、密度ρ' を求める際の、媒質境界と第１媒質に属する格子点との間の距離をｌ'p1とし、媒質境界と第２媒質に属する格子点との間の距離をｌ' p2とする。ただし、数
式を簡略にするために、Δｘ＝ｌ'p1＋ｌ' p2である。また、密度ρ' 'を求める際の、媒質境界と第１媒質に属する格子点との間の距離をｌ' ' p1とし、媒質境界と第２媒質に属する格子点との間の距離をｌ' ' p2とする。３Δｘ＝ｌ' 'p1＋ｌ' ' p2である。
【００９９】
また、図４（ｂ）に示す場合においては、式（９）や式（１０）の第２項に関連する格子点は同一媒質に属し、また第３項に関連する格子点は異なる媒質に属している状態となっている。この場合、式（９）に基づく演算を第１物理量演算ステップで、格子点ｘ＝ｉにおける音圧ｐを演算する際の、式（９）の第２項に用いる体積弾性率Ｋ'と、式（９）
の第３項に用いる体積弾性率Ｋ' 'とはそれぞれ下式（１５）及び（１６）によって決定
する。
【０１００】
【数１５】
2013235506.tif 000019
【０１０１】
【数１６】
2013235506.tif 000020
【０１０２】
また、式（１０）に基づく演算を第２物理量演算ステップで、格子点ｘ＝ｉ＋１／２における粒子速度ｕを演算する際の、式（１０）の第２項に用いる密度ρ'と、式（９）の
第３項に用いる密度ρ' 'とはそれぞれ下式（１７）及び（１８）によって決定する。
【０１０３】
【数１７】
2013235506.tif 000021
【０１０４】
【数１８】
2013235506.tif 000022
【０１０５】
以上のように、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、空間微分をより多くの点の値を用いて高次差分で表現する高精度な手法用いて、不均質媒質中の波動伝搬を解析することが可能になる。
【０１０６】
また、本発明の提供する物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムによれば、音圧の波動方程式に差分法を適用したＷＥ−ＦＤＴＤ法（参考文献１）にも適用することが可能である。
【０１０７】
また、本発明の提供する物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムは、音波の伝搬解析を例のみならず、電磁波や振動等の波動一般の解析に適用することが可能である。
（実施例）
図８は、第１媒質と第２媒質との間の境界条件を組み込んだ音波伝搬の物理量シミュレーションを、従来技術に基づいて実施したものと、本発明に基づいて実施したものを示す。さらに、図８においては、第１媒質と第２媒質との間の境界条件を組み込んだ音波伝搬の理論解も示している。
【０１０８】
図８を参照すると、本発明に係る物理量シミュレーションは、従来技術に係る物理量シミュレーションに比べて、より理論解に近く、計算精度が良いことがわかる。
【符号の説明】
【０１０９】
１０・・・システムバス
１１・・・ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）
１２・・・ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）
１３・・・ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）
１４・・・通信制御部
１５・・・入力制御部
１６・・・出力制御部
１７・・・外部記憶装置制御部
１８・・・入力部
１９・・・出力部
２０・・・外部記憶装置
２１・・・インターフェイス部
２１・・・グラフィック制御部
２２・・・ディスプレイ装置
【図１】
2013235506.tif 000023
【図２】
2013235506.tif 000024
【図３】
2013235506.tif 000025
【図４】
2013235506.tif 000026
【図５】
2013235506.tif 000027
【図６】
2013235506.tif 000028
【図７】
2013235506.tif 000029
【図８】
2013235506.tif 000030