(57)【要約】【課題】媒質の複雑な境界面における弾性場を簡易に解析できるようにする。【解決手段】この発明の弾性波解析方法は、各格子点がどの媒質に属するか否かの属否判定を行い各格子点に識別値を与える格子点属否判定過程と、識別値に基づき各格子点の隣接格子点の有無を判別して隣接格子点情報を与える隣接格子点有無判定過程と、隣接格子点情報の組み合わせに応じて弾性場を求める解法を判定する解法選択過程と、解法に基づいて各格子点の弾性場解析値を算出する解析値算出過程とを含み、複雑な境界面上の格子点の弾性波を簡単に解析することができる。【選択図】図７
2011112612.tif 000002
【特許請求の範囲】
【請求項１】
現時刻又は過去における弾性場を示す解析値から次時刻における各格子点の解析値を算出する解法を、
処理部において、
記憶部に記録された弾性場を伝搬する媒質が占める領域に属すか否かを格子点ごとに示す媒質属否情報から、格子点ごとに当該媒質が占める各座標軸方向ごとの隣接格子点の有無を判別し、当該媒質にかかる隣接格子点の有無情報を当該記憶部に記録する格子点属否判定過程と、
当該記憶部に記録された、各座標軸の正又は負方向ごとの隣接格子点の有無情報の組み合わせに対応する解法からなる解法選択情報テーブルを参照し、当該記憶部に記録された格子点での隣接格子点の有無情報に基づいて解法を選択し、記憶部に記録する解法選択過程を、
実行することを特徴とする弾性波解析方法。
【請求項２】
現時刻又は過去における弾性場を示す解析値から次時刻における各格子点の解析値を算出する解法を、
処理部において、
記憶部に記録された格子点ごとの解法により、上記解析値を算出する解析過程を有し、
上記解法の少なくとも一が、
一の格子点又は他の格子点が少なくとも一の座標軸方向の隣接格子点が同一媒質について存せず、その逆方向に隣接格子点が存する境界面上の格子点における次時刻の上記弾性場を、
上記逆方向の隣接格子点から入射する成分と当該境界面における反射率に応じた上記境界面上の格子点から反射する成分から合成した推定値と、
上記境界面上の格子点における現時刻の上記解析値とから、
現時刻と次時刻の間の上記弾性場の変化に対する上記一の座標軸方向についての寄与として、上記境界面上の格子点における次時刻の上記解析値を算出することを特徴とし、当該解法は当該座標軸方向の隣接格子点が存しないことに対応付けられていることを特徴とする弾性波解析方法。
【請求項３】
請求項１又は２に記載した弾性波解析方法において、
上記境界面上の格子点における弾性場が上記反射率が正値の場合は粒子速度であり、上記反射率が負値の場合は圧力である、
ことを特徴とする弾性波解析方法。
【請求項４】
現時刻又は過去における弾性場を示す解析値から次時刻における各格子点の解析値を算出する解法を実行する弾性波解析装置であって、
記憶部に記録された弾性場を伝搬する媒質が占める領域に属すか否かを格子点ごとに示す媒質属否情報から、格子点ごとに当該媒質が占める各座標軸方向ごとの隣接格子点の有無を判別し、当該媒質にかかる隣接格子点の有無情報を当該記憶部に記録する格子点属否判定部と、
当該記憶部に記録された、各座標軸の正又は負方向ごとの隣接格子点の有無情報の組み合わせに対応する解法からなる解法選択情報テーブルを参照し、当該記憶部に記録された格子点での隣接格子点の有無情報に基づいて解法を選択し、記憶部に記録する解法選択部を、
具備する弾性波解析装置。
【請求項５】
現時刻又は過去における弾性場を示す解析値から次時刻における各格子点の解析値を算出する解法を実行する弾性波解析装置であって、
処理部において、
記憶部に記録された格子点ごとの解法により、上記解析値を算出する解析過程を有し、
上記解法の少なくとも一が、
一の格子点又は他の格子点が少なくとも一の座標軸方向の隣接格子点が同一媒質について存せず、その逆方向に隣接格子点が存する境界面上の格子点における次時刻の上記弾性場を、
上記逆方向の隣接格子点から入射する成分と当該境界面における反射率に応じた上記境界面上の格子点から反射する成分から合成した推定値と、
上記境界面上の格子点における現時刻の上記解析値とから、
現時刻と次時刻の間の上記弾性場の変化に対する上記一の座標軸方向についての寄与として、上記境界面上の格子点における次時刻の上記解析値を算出する解析値算出部を備え、当該解法は当該座標軸方向の隣接格子点が存しないことに対応付けられていることを特徴とする弾性波解析装置。
【請求項６】
請求項４又は５に記載した弾性波解析装置において、
上記境界面上の格子点における弾性場が上記反射率が正値の場合は粒子速度であり、上記反射率が負値の場合は圧力である、
ことを特徴とする弾性波解析装置。
【請求項７】
請求項１乃至３の何れか１項に記載した弾性波解析処理方法をコンピュータに実行させるための弾性波解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、空間を伝播する音波その他の弾性波の波動現象を解析する弾性波解析方法とその装置と、プログラムに関する。
【背景技術】
【０００２】
計算機の計算速度の向上に伴い、音波その他の弾性波のシミュレーションが試みられている。その中でも時間領域差分（ＦＤＴＤ：Finite Difference Time Domain）法は、弾性場を直接的に求めるものとして注目されている。
【０００３】
ＦＤＴＤ法は、例えば非特許文献１に開示されており、ある座標軸点上における次時刻の解析値を現時刻の解析値と隣接する格子点上の解析値に基づき算出するものである。ＦＤＴＤ法では、図１７に示すように弾性波が伝搬する解析対象となる領域を複数の格子点の集合で指定する。ここで、その位置で圧力を計算する格子点（以下、「圧力格子点」ｐと示す）と、その位置で粒子速度を計算する格子点（以下、「速度格子点」、ｕ，ｖ，ｗと示す）が交互に配置される。但し、簡単のため図１７は２次元で格子点の配置を示しており、ｚ軸方向の粒子速度ｗをとる格子点を省略している。
【０００４】
例えば、３次元空間において空気中を伝搬する音波の場合、ある座標（ｉｘ，ｉｙ，ｉｚ）の圧力格子点ｐにおける次時刻ｎ＋１/２の圧力（音波を扱う場合は特に音圧ともいう）ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）は、現時刻ｎ−１/２の圧力ｐｎ−１/２（ｉｘ，ｉｙ，ｉｚ）と、その圧力格子点と隣接する速度格子点の現時刻の粒子速度・ｕｎ（ｉｘ＋１/２，ｉｙ，ｉｚ）等の各座標軸方向の差分に基づいて式（１）で求められる。
【０００５】
【数１】
2011112612.tif 000003
【０００６】
ここで、κは体積弾性率、Δｔは離散時間間隔、Δｘ，Δｙ，Δｚはｘ，ｙ，ｚ方向の
格子点間隔、・ｕｎ（ｉｘ＋１/２,ｉｙ,ｉｚ），・ｖｎ（ｉｘ,ｉｙ＋１/２,ｉｚ），・ｗｎ（ｉｘ,ｉｙ,ｉｚ＋１/２）はｘ，ｙ，ｚ方向の現時刻の粒子速度を示す。本文でｕｎ（ｉｘ＋１/２,ｉｙ,ｉｚ）等に前置している「・」の表記は、本来式中又は図中に示すようにｕの上に置く表記が正しい。
【０００７】
また、時刻ｎ＋１/２とは時刻ｎからΔｔ/２経過した時刻ｎとｎ＋１の中間の時刻を表す。ｎは現実の時刻ｔをｔ＝ｎΔｔと離散化するものである。同様にｉｘ＋１/２とは、ｘ軸方向のｉｘ番目の格子点から正方向にΔｘ/２ずれた位置であり、ｉｘ番目とｉｘ＋１番目の中間の格子点を示し、現実のｘ軸方向の座標ｘをｘ＝（ｎ＋１/２）Δｘと離散化している。つまり、圧力ｐｎ＋１／２（ｉｘ,ｉｙ,ｉｚ）、粒子速度・ｕｎ（ｉｘ＋１/２,ｉｙ,ｉｚ）, ・ｖｎ（ｉｘ,ｉｙ＋１/２,ｉｚ）, ・ｗｎ（ｉｘ,ｉｙ,ｉｚ＋１/２）の格子点のｘ,ｙ,ｚ座標は各々（ｉｘΔｘ,ｉｙΔｙ,ｉｚΔｚ）,（（ｉｘ＋１/２）Δｘ, ｉｙΔｙ,ｉｚΔｚ）, （ｉｘΔｘ, （ｉｙ＋１/２）Δｙ,ｉｚΔｚ）, （ｉｘΔｘ,ｉｙΔｙ, （ｉｚ＋１/２）Δｚ）に位置することになる。
【０００８】
次時刻の粒子速度・ｕｎ+1（ｉｘ＋１/２,ｉｙ,ｉｚ），・ｖｎ+1（ｉｘ,ｉｙ＋１/２,ｉｚ），・ｗｎ+1（ｉｘ,ｉｙ,ｉｚ＋１/２）は、現時刻ｎにおける粒子速度・ｕｎ（ｉｘ＋１/２,ｉｙ,ｉｚ），・ｖｎ（ｉｘ,ｉｙ＋１/２,ｉｚ），・ｗｎ（ｉｘ,ｉｙ,ｉｚ＋１/２）と、その隣接する圧力格子点の現時刻ｎ＋１/２の圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）と圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）の差分に基づいて式（２）、式（３）、式（４）で求められる。
【０００９】
【数２】
2011112612.tif 000004
【００１０】
ここでρは媒質の密度を示す。
【００１１】
図１７に格子点間の関係を矢印で模式的に示す。実線の矢印で、座標（ｉｘ，ｉｙ）の圧力格子点ｐの次時刻の圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）が、その座標と各座標軸方向で隣接する現時刻の速度格子点間の粒子速度の差によって求められる関係を表している。求められた次時刻の圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）によって、ｘ軸方向で隣接する速度格子点ｕの粒子速度が求められる関係を、破線の矢印で示している。
【００１２】
このようにＦＤＴＤ法では、圧力と粒子速度とを交互に求める。この求め方を蛙飛び差分スキームと称することもある（非特許文献１）。図１８にＦＤＴＤ法の計算フローを示す。ステップＳ８０で初期値を設定する。初期値は、現在時刻の各圧力格子点の圧力ｐｎ−１/２（ｉｘ，ｉｙ，ｉｚ）と速度格子点の粒子速度・ｕｎ（ｉｘ＋１/２,ｉｙ,ｉｚ），・ｖｎ（ｉｘ,ｉｙ＋１/２,ｉｚ），・ｗｎ（ｉｘ,ｉｙ,ｉｚ＋１/２）である。  ステップＳ８１で、各圧力格子点の次時刻の圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）を式（１）で算出する。式（１）の計算は対象となる全ての圧力格子点について行われる（ステップＳ８２）。そして時刻をΔｔ/２増加させる（ステップＳ８３）。ステップＳ８４で、 ステップ８１で求めた圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）を用いて式（２）〜（４）で速度粒子点の次時刻ｎ＋１の粒子速度を対象となるすべての速度格子点について算出する（ステップＳ８５）。その後、再び時刻をΔｔ/２増加させる（ステップＳ８６）。ステップ８１と８５の計算は、時刻が予め定めた解析時間に達するまで繰り返される（ステップＳ８７のＮＯ）。
【００１３】
ＦＤＴＤ法では、隣接する格子点の媒質の体積弾性率や密度が一定であることを前提にしている。そのため、隣接格子点が解析対象となる領域外となる場合、つまり隣接格子点が存在しない境界面や、体積弾性率や密度が異なる媒質の境界面においてはそのまま適用できない。
【００１４】
複雑な形状の領域の弾性波を解析するには、媒質の境界条件、例えば、空気中を伝搬する音波に対する壁面による反射の条件を与える必要がある。境界面の形状や反射率などは波面の反射などの波動現象に大きな影響を与える反面、とりわけ３次元空間内では境界面の形状を設定することは困難である。
【００１５】
従来、複雑な形状の領域において解析を行うには、直方体などの単純な形状の領域内に含まれる格子点同士で解析を行い、部分領域が重複する部分において圧力や粒子速度などの解析値を共通に利用する方法がある。この方法は、例えば非特許文献２に開示されている。
【００１６】
当該方法について図１９を用いて説明する。例えば、解析領域が部分領域Ａと、部分領域Ｂとから構成され、両者が重複領域Ｓで重複するものとする。
【００１７】
部分領域Ａで次時刻の圧力を求める際には、原則として部分領域Ａの現時刻の圧力と粒子速度を用いるが、重複領域Ｓにおける粒子速度として部分領域Ｂで算出した値を用いる。また、部分領域Ｂで次時刻の粒子速度を求める際には、原則として部分領域Ｂ内の現時刻の圧力と粒子速度を用いるが、重複領域Ｓにおける圧力として部分領域Ａで算出した値を用いる。
【００１８】
このように、部分領域が重複する領域については一の部分領域で算出した現時刻の圧力や粒子速度を他の部分領域で次時刻の粒子速度や圧力を求めるために共通に用いるため、部分領域間で重複領域における圧力や粒子速度の受け渡しを要していた。
【先行技術文献】
【非特許文献】
【００１９】
【非特許文献１】佐藤雅弘「ＦＤＴＤ法による弾性振動・波動の解析入門」森北出版
【非特許文献２】橋本  修「実践ＦＤＴＤ時間領域差分法」森北出版
【発明の概要】
【発明が解決しようとする課題】
【００２０】
しかし、従来のＦＤＴＤ法では、解析対象の媒質の形状が複雑になると、多くの領域を接続して、それらの重複領域における解析値を共通化しなければならない。その解析値の受け渡しが煩雑となる。また、境界面の形状が複雑になることで、境界条件の設定が煩雑になり、シミュレーションにおける利便性を低下させる要因になっていた。
【００２１】
この発明は、このような問題点に鑑みてなされたものであり、複雑な形状の解析領域でも、シミュレーションの利便性を低下させない弾性波解析方法とその装置と、プログラムを提供することを目的とする。
【課題を解決するための手段】
【００２２】
この発明の弾性波解析法は、現時刻又は過去における弾性場を示す解析値から次時刻における各格子点の解析値を算出する解法であって、処理部において格子点属否判定過程と解法選択過程を実行する。格子点属否判定過程は、記憶部に記録された弾性場を伝搬する媒質が占める領域に属すか否かを格子点ごとに示す媒質属否情報から、格子点ごとに当該媒質が占める各座標軸方向ごとの隣接格子点の有無を判別し、当該媒質にかかる隣接格子点の有無情報を当該記憶部に記録する。解法選択過程は、当該記憶部に記録された格子点での隣接格子点の有無情報に基づいて解法を選択し、記憶部に記録する。
【発明の効果】
【００２３】
この発明の弾性波解析方法によれば、各格子点がどの媒質に属するか否かを表す識別値に基づいて、同一媒質に属する各格子点の隣接する格子点を持つか否かを判定し、各格子点の隣接格子点有無情報を生成する。そして、その隣接格子点有無情報に対応した弾性波の解析方法を選択する。これにより境界面上の格子点であることを識別し、解析方法の一環として当該格子点における境界条件が設定される。したがって、従来技術のように境界面において解析値の受け渡しを行う必要が無いので、複雑な形状の解析空間における弾性波解析をより容易に行うことができる。
【図面の簡単な説明】
【００２４】
【図１−１】圧力格子点と速度格子点を模式的に示す図であり、（ａ）はｘ，ｙ，ｚ軸上の正負方向ともに隣接格子点を有する場合（Ｓｐ＝１）、（ｂ）はｘ，ｙ，ｚの一つの座標軸上の正負いずれか一方について隣接格子点を有しない場合（Ｓｐ＝２）である。
【図１−２】圧力格子点と速度格子点を模式的に示す図であり、（ｃ）はｘ，ｙ，ｚのうち２つの座標軸上の正負両方向ともいずれも隣接格子点を有する場合（Ｓｐ＝３）、（ｄ）はｘ，ｙ，ｚのうち２つの座標軸上の正負いずれか一方について隣接格子点を有しない場合（Ｓｐ＝４）である。
【図１−３】圧力格子点と速度格子点を模式的に示す図であり、（ｅ）はｘ，ｙ，ｚのうち１つの座標軸上の正負両方向とも隣接格子点を有するが、１つの座標軸上の正負いずれか一方のみ隣接格子点を有する場合（Ｓｐ＝５）、（ｆ）はｘ，ｙ，ｚ各座標軸ともに正負いずれか一方について隣接格子点を有しない場合（Ｓｐ＝６）である。
【図１−４】圧力格子点と速度格子点を模式的に示す図であり、（ｇ）はｘ，ｙ，ｚのうち１つの座標軸上の正負両方向とも隣接格子点を有する場合（Ｓｐ＝７）、（ｈ）はｘ，ｙ，ｚのうち２つの座標軸上の正負いずれか一方のみ隣接格子点を有する場合（Ｓｐ＝８）である。
【図１−５】圧力格子点と速度格子点を模式的に示す図であり、（ｉ）はｘ，ｙ，ｚのうち１つの座標軸上の正負いずれか一方のみ隣接格子点を有する場合（Ｓｐ＝９）、（ｊ）は隣接格子点を有しない場合（Ｓｐ＝１０）である。
【図２】境界面上の圧力格子点における入射波と反射波の関係を示す図。
【図３】この発明の弾性波解析装置１００の機能構成例を示す図。
【図４】弾性波解析装置１００の動作フローを示す図。
【図５】属否判定の方法の一例を示す図。
【図６】圧力格子点上に境界面をとるモデルの例を示す図。
【図７】隣接格子点有無判定部２３、解法選択部２４、解析値算出部２５の機能構成例を示す図。
【図８】圧力格子点の隣接格子点有無判定部２３１の動作フローを示す図。
【図９】速度格子点（粒子速度・ｕ）の隣接格子点有無判定部２３２の動作フローを示す図。
【図１０】速度格子点（粒子速度・ｖ）の隣接格子点有無判定部２３２の動作フローを示す図。
【図１１】速度格子点（粒子速度・ｗ）の隣接格子点有無判定部２３２の動作フローを示す図。
【図１２−１】解法選択テーブルを示す図であり、（ａ）は項番１〜１３、（ｂ）は項番１４〜２６を示す図である。
【図１２−２】解法選択テーブルを示す図であり、（ｃ）は項番２７〜３９、（ｄ）は項番４０〜５２を示す図である。
【図１２−３】解法選択テーブルを示す図であり、（ｅ）は項番５３〜６４を示す図である。
【図１３】速度格子点の解法選択テーブルを示す図であり、（ａ）はｘ軸方向、（ｂ）はｙ軸方向、ｚはｚ軸方向を示す図である。
【図１４】速度格子点と圧力格子点を模式的に示す図であり、（ａ）は正負両方向とも隣接格子点を有する場合（Ｓｕ＝１）、（ｂ）は正負いずれか一方に隣接格子点を有する場合（Ｓｕ＝２）、（ｃ）は隣接格子点を有しない場合（Ｓｕ＝３）である。
【図１５】解析値算出部２５の動作フローを示す図。
【図１６】境界面上の速度格子点における入射波と反射波の関係を示す図。
【図１７】従来のＦＤＴＤ法で弾性波を解析するのに解析対象を格子点でモデル化した例を示す図。
【図１８】従来のＦＤＴＤ法の動作フローを示す図。
【図１９】部分領域Ａと部分領域Ｂとから構成される解析領域の例を示す図。
【発明を実施するための形態】
【００２５】
以下、この発明の実施の形態を図面を参照して説明する。複数の図面中同一のものには同じ参照符号を付し、説明は繰り返さない。実施例の説明の前に、この発明における圧力、粒子速度の取り扱いについて説明する。
【００２６】
媒質内部からその物質の境界面に入射する弾性波（圧力、粒子速度で示される）は境界面を反射する。裏返せば、境界面への入射成分と反射成分を求め、それらを反射の度合いに応じて合成する条件を境界面上で与えることができれば、境界面上の格子点における弾性波の解析値（圧力、粒子速度）を求めることが可能になる。
【００２７】
図１−１〜１−５に、一の圧力格子点ｐに隣接するｘ,ｙ,ｚ軸方向の速度格子点ｕ,ｖ,ｗを模式的に示す。図１−１（ａ）は圧力格子点ｐの周りの全ての速度格子点が同一媒質に存在する場合を示す。図１−１（ａ）の場合は、従来のＦＤＴＤ法で次時刻の弾性波解析値を求めることが可能である。
【００２８】
図１−１（ｂ）に示すように、ｘ軸上の負方向に、解析領域内における隣接格子点（ｉｘ，ｉｙ−１/２，ｉｚ）が存在せず、ｙ，ｚ軸方向に、解析領域内における隣接格子点（ｉｘ，ｉｙ−１/２，ｉｚ），（ｉｘ，ｉｙ＋１/２，ｉｚ），（ｉｘ，ｉｙ，ｉｚ−１/２），（ｉｘ，ｉｙ，ｉｚ＋１/２）が、いずれも存在する場合、次時刻の圧力を式（５）で算出できる。
【００２９】
【数３】
2011112612.tif 000005
【００３０】
κは体積弾性率である。式（５）のΔｐｘｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）は、推定値ｐｅｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）と現時刻の圧力との差分であって、次時刻と現時刻の圧力の差分におけるｘ軸方向についての寄与を示すにすぎない（式（６））。
【００３１】
【数４】
2011112612.tif 000006
【００３２】
推定値ｐｅｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）は、この発明は最も簡単な一次形式の例として式（７）に基づいて、その圧力格子点における次時刻ｎ＋１/２の圧力を算出する。
【００３３】
【数５】
2011112612.tif 000007
【００３４】
ここで、Ｒ１は反射率（但しＲ１≦０）、ｂｘはｘ方向の伝搬係数（ｂｘ＝（ｃΔｔ−Δｘ）/（ｃΔｔ＋Δｘ）、ｃは波面進行速度である。ρは密度。ｄ１は解の発散を抑制するための安全係数であり、零以上の微小な正値である。解析の条件にもよるが好ましくはｄ１＜０．０１である。また、ｄ１は境界面で波面が反射する際に生じる減衰を解析する際の減衰率としても利用できる。
【００３５】
式（７）の内、式（８）を満たす成分は、式（９）が示すｘ方向へ波面進行速度ｃで反射波が生じずに伝搬する条件、即ち入射波だけを与え反射波を生じさせない条件を時刻ｎにおいて格子点（ｉｘ,ｉｙ,ｉｚ）,（ｉｘ−１,ｉｙ,ｉｚ）の中間点の座標（ｉｘ−１/２,ｉｙ,ｉｚ）において離散化することにより与えられる。
【００３６】
【数６】
2011112612.tif 000008
【００３７】
また、上記した式（７）の内、格子点（ｉｘ,ｉｙ,ｉｚ）上の圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）等に着目すると式（１０）が与えられる。
【００３８】
【数７】
2011112612.tif 000009
【００３９】
式（１０）を満たす成分は、時間経過と共に零に近似する。これは、ｃΔｔ、Δｘ共に正値であるから、伝搬係数ｂｘの絶対値が常に１未満となるからである。もっとも、圧力の初期値ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）が零であれば時間を経過しても零のままであるため、圧力の初期値ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）を当初から零としてもよい。
【００４０】
次時刻の圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）は零になると、圧力格子点（ｉｘ，ｉｙ，ｉｚ）ではｘの正方向に波面進行速度ｃで伝搬する入射波成分ｐ（ｃｔ−ｋ（ｘ−ｉｘΔｘ））と、その逆の負方向に波面進行速度ｃで伝搬する反射波成分ｐ（ｃｔ＋ｋ（ｘ−ｉｘΔｘ））が逆位相で生じる。ここでｋは波数（ｋ＝２π/λ、λは波長）である。図２に境界面上の圧力格子点における入射波と反射波の関係を示す。図２の横方向は波面進行波の進行方向、縦方向は進行波の振幅を表す。入射波を実線、反射波を破線で示す。
【００４１】
即ち、式（８）の部分により境界面への入射波に相当する成分が与えられる。また式（１０）により、又は境界面上の格子点（ｉｘ,ｉｙ,ｉｚ）で圧力ｐ（ｉｘ,ｉｙ,ｉｚ）＝０とおくことで入射波を受けて境界面から媒質内部に位相が反転した反射波に相当する成分が与えられる。その成分に反射率Ｒ１を乗じ、入射波に相当する成分と合成するという条件を与える。よって、入射波ｐ（ｃｔ−ｋ（ｘ−ｉｘΔｘ））に応じ反射率Ｒ１で逆位相の反射波成分−ｐ（ｃｔ＋ｋ（ｘ−ｉｘΔｘ））が合成されているとの条件のもとで、境界面上の格子点における圧力ｐ（ｉｘ,ｉｙ,ｉｚ）を求めることができる。このようにして得た図２中に太線で示す次時刻の圧力ｐ（ｉｘ,ｉｙ,ｉｚ）は、波動方程式を満たすので境界面上の格子点の圧力を安定して推定することができる。
【００４２】
即ち、式（７）に示したように少なくとも一の方向の隣接格子点が同一媒質について存
せず、その逆方向に隣接格子点が存する境界面上の格子点における次時刻の弾性波変数を、逆方向の隣接格子点から入射する成分と当該境界面における反射の度合いに応じた当該境界面上の格子点から反射する成分と合成することと等価となる。
【００４３】
上記と異なり図１−１（ｂ）に示すように、ｘ軸の負方向に隣接する速度格子点（座標（ｉｘ−１/２,ｉｙ,ｉｚ））が同一媒質について存せず、逆方向に隣接する速度格子点（座標（ｉｘ＋１/２,ｉｙ,ｉｚ））が存する場合は式（７）における格子点（ｉｘ−１,ｉｙ,ｉｚ）を（ｉｘ＋１,ｉｙ,ｉｚ）に置き換えて適用すればよい。この場合の解法を、例えばＳｐ＝２とする。
【００４４】
図１−２（ｄ）に示すように、ｙ,ｚ軸の正方向に隣接する速度格子点（座標（ｉｘ,ｉｙ＋１/２,ｉｚ）, （ｉｘ,ｉｙ,ｉｚ＋１/２））が同一媒質について存しない場合も、式（５）における格子点（ｉｘ−１,ｉｙ,ｉｚ）を（ｉｘ,ｉｙ−１,ｉｚ）, （ｉｘ,ｉｙ,ｉｚ−１）と各々置き換えて適用すればよい。負方向の隣接する速度格子点がない場合も式（５）における格子点（ｉｘ−１,ｉｙ,ｉｚ）を（ｉｘ,ｉｙ＋１,ｉｚ）, （ｉｘ,ｉｙ,ｉｚ＋１）と各々置き換えて適用すればよい。この場合の解法を、例えばＳｐ＝４とする。
【００４５】
ｚ軸方向のみ隣接格子点が正負ともに存在し、ｘ，ｙ軸各々について正負何れか一方の隣接格子点が存在しないときは、各方向軸に対する寄与を合成して次時刻の圧力ｐｎ+１/２（ｉｘ，ｉｙ，ｉｚ）を求めることができる。
【００４６】
【数８】
2011112612.tif 000010
【００４７】
また、図１−３（ｆ）に示すようにｘ,ｙ,ｚ軸各々について正負の何れか一方の隣接格子点が各々存在しないときは、次のように求めることができる。
【００４８】
【数９】
2011112612.tif 000011
【００４９】
このように方向成分ごとの寄与を合成して次時刻の圧力を求めることにより、境界面への波動の入射角θiが任意の場合でも、圧力の推定をおこなうことができる。隣接格子点が存在しない座標軸が異なる場合でもこの方法は適用できる。
【００５０】
なお、図１−２（ｃ）及び図１−３〜１−５については、解析値算出部２５の説明の際に参照する。また、速度格子点の粒子速度についての解法については後述する。
【００５１】
なお、式(７)は、反射率Ｒ１が−１≦Ｒ１≦０の場合でなければ適用できない。何故ならば反射率Ｒ１が正値の場合、式（７）により求められる圧力ｐｎ＋１/２（ｉｘ,ｉｙ,ｉｚ）は不安定になるからである。
【００５２】
その理由を説明する。ここで、空間上の座標ｒにおける圧力ｐ（ｒ，ｔ）が時刻ｔにおける各格子点上の値ｐｒ＝（ｉｘΔｘ,ｉｙΔｙ，ｉｚΔｚ）（ｔ）を与える平面波ｅ−ｊｋｒで表されるとする。つまりｐ（ｒ,ｔ）＝ｐｒ（ｔ）ｅ−ｊｋｒとおく。このとき、式（７）は、式（７ｂ）のように変形される。
【００５３】
【数１０】
2011112612.tif 000012
【００５４】
両辺を整理すると式（７ｃ）となる。
【００５５】
【数１１】
2011112612.tif 000013
【００５６】
ここで、ｐｎ＋１/２を安定に求めるためには式（７ｄ）で表す関係を満たす必要がある。
【００５７】
【数１２】
2011112612.tif 000014
【００５８】
ここで、０＜ｄ１≪１,０＜ｂｘ＜１であることを考慮すると式（７ｅ）を満たせばよい。
【００５９】
【数１３】
2011112612.tif 000015
【００６０】
式（７ｅ）によれば、そのようなＲ１は、−２≦Ｒ１≦０となるが、一般に反射率は−１≦Ｒ１≦１であるため、両者が重複する−１≦Ｒ１≦０の場合が式（７）により安定して解を求める条件となる。
【００６１】
そこで、反射率Ｒ１＞０の場合の解析は、粒子速度・ｕ，・ｖ，・ｗにおいて対処する。その場合については後述する（粒子速度の計算方法の説明、式（２３）参照）。
【００６２】
なお、反射率Ｒ１は一の媒質１における密度ρ１と波面進行速度ｃ１、他の媒質２における密度ρ２と波面進行速度ｃ２により式（１３）に示すように定められる。
【００６３】
【数１４】
2011112612.tif 000016
【００６４】
ここでφは媒質２における回折角を表す。また、スネルの法則により回折角φは式（１４）で与えられる。但し、ｃ２＞ｃ１であって、式（１５）に示す関係がある場合は媒質２への回折は起こらず境界面で全反射するためＲ１＝１となる。
【００６５】
【数１５】
2011112612.tif 000017
【実施例１】
【００６６】
図３にこの発明の弾性波解析装置１００の機能構成例を示す。弾性波解析装置１００は、データ入出力部１０と、処理部２０と、記憶部３０と、解析値表示部４０と、を具備する。弾性波解析装置１００は、例えば記憶部３０としてＲＡＭ、処理部２０としてＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、処理部２０がそのプログラムを実行することで実現することができる。
【００６７】
データ入出力部１０は、弾性波解析に必要な各種データを入力すると共に解析結果を出力する。各種データとしては、解析する空間領域、ｘ，ｙ，ｚの各方向の格子点距離、媒質の体積弾性率κ、密度ρ（または比重）、境界面における反射率Ｒがある。各種データと解析結果は、記憶部３０に記憶される。解析値表示部４０は、各時刻において解析された各圧力格子点の圧力と速度格子点の粒子速度を、圧力記憶領域３１、ｘ軸粒子速度記憶領域３２〜ｚ軸粒子速度記憶領域３４から読み出して表示する。
【００６８】
処理部２０は、解析領域設定部２１と、格子点属否判定部２２と、隣接格子点有無判定部２３と、解法選択部２４と、解析値算出部２５と、を備える。図４に弾性波解析装置１００の動作フローを示す。
【００６９】
データ入出力部１０より解析に必要なデータが入力される（ステップＳ１０）。必要なデータとしては、弾性波の振る舞いを解析する領域Ｓ、ｘ，ｙ，ｚ各方向の格子点距離、または格子点、波動を伝搬する媒質の体積弾性率κ、媒質の密度ρ（または比重）、境界面における反射率Ｒ１等である。これらのデータは記憶部３０の各種係数記憶領域３５に記憶される。
【００７０】
解析領域設定部２１は、解析対象とする格子点を定める（ステップＳ２１）。圧力ｐｎ＋１/２（ｉｘ,ｉｙ，ｉｚ）、粒子速度・ｕｎ（ｉｘ＋１/２,ｉｙ，ｉｚ）、・ｖｎ（ｉｘ,ｉｙ＋１/２，ｉｚ）、・ｗｎ（ｉｘ,ｉｙ，ｉｚ＋１/２）の格子点のｘ,ｙ,ｚ座標は各々（ｉｘΔｘ,ｉｙΔｙ,  ｉｚΔｚ）, （（ｉｘ＋１/２）Δｘ,ｉｙΔｙ,  ｉｚΔｚ）, （ｉｘΔｘ,（ｉｙ＋１/２）Δｙ,  ｉｚΔｚ）, （ｉｘΔｘ,ｉｙΔｙ,  （ｉｚ＋１/２）Δｚ）に位置することとなる。ここで、１≦ｉｘ≦Ｎｘ、１≦ｉｙ≦Ｎｙ、１≦ｉｚ≦Ｎｚ、Ｎｘ,Ｎｙ,Ｎｚは各々ｘ,ｙ,ｚ軸方向への格子点数であり、全解析領域はＮｘ×Ｎｙ×Ｎｚ点からなる３次元直方体とする。
【００７１】
格子点属否判定部２２は、波動が伝播する媒質ごとに空間領域Ｓが与えられている場合にどの領域に属するか否かを判定する（ステップS２２）。その判定は、格子点ごとにどの媒質に属すか、あるいはいずれの媒質にも属さないかを示す属否判定情報ｃ（ｉｘ，ｉｙ，ｉｚ）を与えて行う。ｃは任意の整数である。
【００７２】
図５を参照して属否判定の方法を説明する。領域Ｓに属するか否かは、ｘ，ｙ，ｚのうち、ある平面座標（例えばｉｘ,ｉｙ）において、対応する座標軸（例えばｚ軸方向）を解析領域全体の一端から（例えばｉｚ＝１）から他端（ｉｚ＝Ｎｚ）までの間を走査し、境界面に交わった回数を数えることで判断できる。即ち、交わった回数が偶数回である格子点は領域外（例えば属否判定情報＝０）、奇数回の格子点は領域内（例えば属否判定情報＝１）と判断できる。これを他のｘ，ｙ座標に亘って実行することで、全格子点について、ある媒質の解析領域への属否判定が可能になる。
【００７３】
媒質によって、密度、体積弾性率が異なるため、２以上の媒質について解析を行う場合は、媒質ごとに属否判定を行い各格子点ごとに、例えば、何れの媒質にも属さない（０）,媒質１（１）,媒質２（２）を示す属否判定情報を与える。
【００７４】
内外判定の方法として、ある座標面（例えばｘ,ｙ座標平面）内である媒質（媒質１）が占める領域の外縁を与える多角形の頂点座標群が対応する座標（この場合はｚ座標）ごとに領域空間情報Ｓとして与えられている場合には、既知の解法によりその多角形内に属する格子点を領域内（属否判定情報＝１）、属さない格子点を領域外（属否判定情報＝０）と判断することができる。これを他のｚ座標に亘って実行することで、全格子点について媒質１の解析領域への属否判定を成しうる。
【００７５】
なお、領域空間情報Ｓと座標軸が接する点、領域を画する多角形の頂点又は辺上にある点は特異点となる。特異点で交わった場合は、交わった回数として加えないことで正しく属否判定を行うことができる。求められた格子点毎の領域への属否判定情報ｃ（ｉｘ,ｉｙ,ｉｚ）は解析領域記憶領域３６に記憶される。
【００７６】
この発明の弾性波解析方法においては、反射率Ｒ１が零以上の場合は境界面として粒子速度・ｕ，・ｖ，・ｗをとる速度格子点を、反射率Ｒ１が零以下の場合は圧力を取る圧力格子点を境界面とする必要がある。もっとも、反射率Ｒ１が零の場合はどちらの格子点でもよい。
【００７７】
図６に圧力格子点上に境界面（ｘ軸座標がｉｘ、ｘ軸方向に垂直）をとる場合の例を示す。この例では、境界面上にｙ軸方向（ｘ軸方向に垂直）の粒子速度・ｖをとる速度格子点を通る（太い実線）が、何れも各々隣接する圧力格子点を有するため粒子速度・ｖ，・ｗを従来のＦＤＴＤ法で算出しても特に問題にならない。このようなことは、圧力格子点上の境界面が一つの座標軸方向に垂直な場合であって、その座標軸に垂直な方向の粒子速度の速度格子点を含む場合についていえる。
【００７８】
しかし、一の座標軸方向の速度格子点上の境界面が、その座標軸に垂直であっても圧力格子点を含む場合（例えば、ｙ座標がｉｙの粒子速度・ｕをとる速度格子点上に境界面（ｙ軸方向に垂直）が存する場合）は、その圧力格子点は同一の媒質にかかる隣接格子点の何れかが存しない。従って、そのような境界面をとらないように、境界面に垂直な方向と同一の座標軸方向の速度格子点をとるように領域設定を行う。例えば、図６の例（太い破線）のようにｙ軸方向に垂直な境界面を取る場合は図６の例（太い破線）のように最も近接したｙ座標ｉｙ＋１/２の格子点を通るようにｙ軸粒子速度・ｖをとる速度格子点上にとればよい。
【００７９】
隣接格子点有無判定部２３は、格子点毎の領域への属否判定情報ｃ（ｉｘ,ｉｙ,ｉｚ）に基づき、隣接格子点の有無を判別する（ステップＳ２３）。解法選択部２４は、各格子点における隣接格子点の有無情報の組み合わせに応じて圧力ｐｎ＋１/２（ｉｘ,ｉｙ,ｉｚ）や粒子速度・ｕｎ＋１（ｉｘ＋１/２,ｉｙ,ｉｚ）, ・ｖｎ＋１（ｉｘ,ｉｙ＋１/２,ｉｚ）, ・ｗｎ＋１（ｉｘ,ｉｙ,ｉｚ＋１/２）を求める際の解法を、解法選択テーブル記憶領域３８に記憶されたテーブルを参照して選択する（ステップＳ２４）。選択された解法は、解法情報記憶部３７に記憶される。
【００８０】
解法値算出部２５は、格子点毎に選択された解法に基づいて次時刻の圧力と粒子速度を算出する（ステップＳ２５）。以上述べたように、この発明の弾性波解析方法は、各格子点が、どの媒質に属するか否かを表す属否情報に基づいて、同一の媒質に属す各格子点の隣接する格子点を持つか否かを判定し、その判定結果に対応した解析方法を適用する。よって、従来技術のように境界面において解析値の受け渡しを行う必要が無いので、複雑な形状の解析空間における弾性波解析をより容易に行うことができる。
【００８１】
図７に、隣接格子点有無判定部２３、解法選択部２４、解析値算出部２５の機能構成例を示して更に詳しく説明する。
【００８２】
〔隣接格子点有無判定部〕
隣接格子点有無判定部２３は、圧力格子点の隣接格子点有無判定部２３１と、速度格子点の隣接格子点有無判定部２３２と、を備える。図８に圧力格子点の隣接格子点有無判定部２３１の動作フローを示す。
【００８３】
圧力格子点の隣接格子点有無判定部２３１は、属否情報ｃ（ｉｘ,ｉｙ,ｉｚ）を入力として、圧力格子点（ｉｘ,ｉｙ,ｉｚ）のｘ,ｙ,ｚ各座標軸の正負６方向について、隣接格子点の有無を判定する。最初に任意の圧力格子点（ｉｘ,ｉｙ,ｉｚ）を選択する（ステップＳ２３１ａ）。
【００８４】
圧力格子点（ｉｘ,ｉｙ,ｉｚ）のｘ軸の正方向で隣接する速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）が、同一の媒質が占める領域内であるか否かをその速度格子点の属否情報ｃ（ｉｘ＋１/２,ｉｙ,ｉｚ）から判定（ステップＳ２３１ｂ）し、当該圧力格子点（ｉｘ,ｉｙ,ｉｚ）と同一の媒質が占める領域内であれば隣接格子点の有無情報Ｘ正＝１（ステップＳ２３１ｄ）とし、同一の媒質が占める領域内で無ければＸ正＝０とする（ステップＳ２３１ｃ）。
【００８５】
次に、圧力格子点（ｉｘ,ｉｙ,ｉｚ）のｘ軸の正方向で隣接する速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）が、同一の媒質が占める領域内であるか否かを判定する（ステップＳ２３１ｅ）。そして、当該圧力格子点（ｉｘ,ｉｙ,ｉｚ）と同一の媒質が占める領域内であれば隣接格子点の有無情報Ｘ負＝１（ステップＳ２３１ｇ）とし、同一の媒質が占める領域内で無ければＸ負＝０とする（ステップＳ２３１ｆ）。
【００８６】
同様の処理をｙ軸、ｚ軸方向の計６方向について行う。その結果、各圧力格子点（ｉｘ,ｉｙ,ｉｚ）について、６個の隣接格子点有無情報（Ｘ正, Ｘ負, Ｙ正, Ｙ負, Ｚ正, Ｚ負）が出力される。
【００８７】
ステップＳ２３１ｂ〜ステップＳ２３１ｕまでの動作は、未処理の圧力格子点が存在しなくなるまで圧力格子点（ｉｘ,ｉｙ,ｉｚ）を変更して繰り返される（ステップＳ２３１ｖ）。
【００８８】
速度格子点の隣接格子点有無判定部２３２は、属否情報ｃ（ｉｘ,ｉｙ,ｉｚ）を入力として、粒子速度・ｕを取る格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）の正負２方向について、隣接格子点の有無を判定する。図９に、粒子速度・ｕを取る速度格子点の隣接格子点有無判定動作の動作フローを示す。
【００８９】
最初に任意の粒子速度・ｕにかかる速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）を選択する（ステップＳ２３２ａ）。速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）のｘ軸の正方向で隣接する圧力格子点（ｉｘ＋１,ｉｙ,ｉｚ）が、同一の媒質が占める領域内であるか否かをその圧力格子点の属否情報ｃ（ｉｘ＋１,ｉｙ,ｉｚ）から判定（ステップＳ２３２ｂ）し、当該速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）と同一の媒質が占める領域内であれば隣接格子点の有無情報ＵＸ正＝１（ステップＳ２３２ｄ）とし、同一の媒質が占める領域内で無ければＵＸ正＝０とする（ステップＳ２３２ｃ）。
【００９０】
次に、速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）のｘ軸の負方向で隣接する圧力格子点（ｉｘ,ｉｙ,ｉｚ）が、同一の媒質が占める領域内であるか否かを判定する（ステップＳ２３２ｅ）。そして、当該速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）と同一の媒質が占める領域内であれば隣接格子点の有無情報ＵＸ負＝１（ステップＳ２３２ｇ）とし、同一の媒質が占める領域内で無ければＵＸ負＝０とする（ステップＳ２３２ｆ）。
【００９１】
このように、粒子速度・ｕにかかる格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）について、２個の隣接格子点有無情報（ＵＸ正, ＵＸ負）が出力される。他のｙ軸、ｚ軸についても、図１０と図１１に示すように同様な処理が行われ、それぞれの軸に対して２個の隣接格子点有無情報（ＶＹ正, ＶＹ負）, （ＷＺ正, ＷＺ負）が出力される。図１０と図１１の動作フローの説明は、図９と同様なので省略する。
【００９２】
〔解法選択部〕
解法選択部２４は、圧力格子点解法選択部２４１と、速度格子点解法選択部２４２と、を備える。圧力格子点解法選択部２４１は、各圧力格子点（ｉｘ,ｉｙ,ｉｚ）について、６個の隣接格子点有無情報（Ｘ正, Ｘ負, Ｙ正, Ｙ負, Ｚ正, Ｚ負）を入力として、圧力の解析方法である解法Ｓｐを選択する。
【００９３】
速度格子点解法選択部２４２は、隣接格子点有無情報（ＵＸ正, ＵＸ負），（ＶＹ正, ＶＹ負）,（ ＷＺ正, ＷＺ負）を入力として粒子速度の解析方法である解法Ｓｕ，Ｓｖ，Ｓｗを選択する。
【００９４】
〔圧力格子点解法選択部〕
図１２−１〜１２−３に、圧力格子点解法選択部２４１で行われる処理において用いる解法テーブルの内容を示す。項番は単なる識別用の番号である。隣接格子点有無情報（Ｘ正, Ｘ負, Ｙ正, Ｙ負, Ｚ正, Ｚ負）に解法Ｓｐが対応している。即ち、圧力格子点解法選択部２４１では、この解法の選択は、隣接格子点有無情報（Ｘ正, Ｘ負, Ｙ正, Ｙ負, Ｚ正, Ｚ負）に基づき、記憶部３０に予め記憶した解法テーブルを参照して行われる。この解法テーブルは、解法選択テーブル記憶領域３８に記憶されている。
【００９５】
ここで、解法Ｓｐ＝１は隣接格子点有無情報（Ｘ正, Ｘ負, Ｙ正, Ｙ負, Ｚ正, Ｚ負）が全て「１」の場合に選択される（図１−１（ａ）参照）。解法Ｓｐ＝２は、例えば項番７のようにｘ軸の負方向の隣接格子点（ｉｘ−１/２，ｉｙ，ｉｚ）が存在せず、ｙ，ｚ軸方向における隣接格子点が存在する場合に選択される（図１−１（ｂ））。選択した解法Ｓｐを解法情報記憶領域３７に記録する。
【００９６】
このように隣接格子点有無情報（Ｘ正, Ｘ負, Ｙ正, Ｙ負, Ｚ正, Ｚ負）と解法Ｓｐとの関係は、例えば図１２−１〜１２−３に示す通りなので、個々の解法と隣接格子点有無情報との関係の説明は省略する。
【００９７】
〔速度格子点解法選択部〕
図１３に、速度格子点解法選択部２４２における処理で用いられる解法テーブルの内容を示す。この解法テーブルも解法選択テーブル記憶領域３８に記憶されている。項番は単なる識別用の番号である。図１３（ａ）に隣接格子点有無情報（ＵＸ正, ＵＸ負）に対応して選択される解法Ｓｕ、図１３（ｂ）に隣接格子点有無情報（ＶＹ正, ＶＹ負）に対応して選択される解法Ｓｖ、図１３（ｃ）に隣接格子点有無情報（ＷＺ正, ＷＺ負）に対応して選択される解法Ｓｗが示されている。
【００９８】
解法Ｓｕ＝１は、隣接格子点有無情報（ＵＸ正, ＵＸ負）が全て「１」の場合に選択される。解法Ｓｕ＝２は、ｘ軸上の正負何れか一方にのみ隣接格子点が存する場合（例えば図（図１４（ｂ））に選択される。解法Ｓｕ＝３は、隣接格子点が存しない場合（例えば図１４（ｃ））に選択される。選択した解法Ｓｕを解法情報記憶領域３７に記録する。
【００９９】
他のｙ，ｚ軸方向の解法Ｓｖ，Ｓｗも同様に選択され記憶される。各解法の具体的な計算方法については後述する。
【０１００】
〔解析値算出部〕
解析値算出部２５は、現在の圧力ｐｎ−１/２（ｉｘ，ｉｙ，ｉｚ）と粒子速度・ｕｎ（ｉｘ＋１/２,ｉｙ,ｉｚ），・ｖｎ（ｉｘ,ｉｙ＋１/２,ｉｚ），・ｗｎ（ｉｘ，ｉｙ，ｉｚ＋１/２）とを用いて、解法選択部２４で選択された解法（Ｓｐ，Ｓｕ，Ｓｖ，Ｓｗ）を適用して次時刻の弾性波解析値を離散時間毎に求め、記憶部３０に記憶する。図１５に解析値算出部２５の動作フローを示す。
【０１０１】
この解析値算出過程の前の段階で、圧力格子点解法選択部２４１と速度格子点解法選択部２４２による解法選択過程（ステップＳ２４、図４参照）を実行しておき、各格子点（ｉｘ，ｉｙ，ｉｚ）の初期値又は現時刻の圧力ｐｎ−１/２（ｉｘ，ｉｙ，ｉｚ）、初期値又は現時刻の粒子速度・ｕｎ（ｉｘ＋１/２,ｉｙ,ｉｚ），・ｖｎ（ｉｘ,ｉｙ＋１/２,ｉｚ），・ｗｎ（ｉｘ，ｉｙ，ｉｚ＋１/２）と、各格子点の弾性波解析値を求める解法Ｓｐ，Ｓｕ，Ｓｖ，Ｓｗとを設定し、解法情報記憶領域３７に記録させておく（ここでｎ＝１）。
【０１０２】
圧力算出過程（ステップＳ２５０）は、解法Ｓｐ読出ステップＳ２５０ａと圧力算出ステップＳ２５０ｂを含む。解法Ｓｐ読出ステップＳ２５０ａは、解法情報記憶領域３７に記録されている解法Ｓｐを読み出す。圧力算出ステップＳ２５０ｂは、その読み出した解法にしたがって次時刻の圧力を算出する。各解法については後述する。圧力算出の後、時刻がΔｔ/２増加させる（ステップＳ２５１）。
【０１０３】
ｘ軸粒子速度算出過程は、解法Ｓｕ読出ステップＳ２５２ａとｘ軸粒子速度算出ステップＳ２５２ｂを含む。解法Ｓｕ読出ステップＳ２５２ａは、解法情報記憶領域３７に記録されている解法Ｓｕを読み出す。ｘ軸粒子速度算出ステップＳ２５２ｂは、時刻ｎ＋１における粒子速度ｕにかかる速度格子点の粒子速度ｕを読み出した解法（Ｓｕ）に従って算出する（ステップＳ２５２）。  
【０１０４】
同様に、ｙ軸粒子速度算出過程は、時刻ｎ＋１における粒子速度ｖにかかる速度格子点の粒子速度を読み出した解法（Ｓｖ）に従って算出する（ステップＳ２５３）。ｚ軸粒子速度算出過程は、時刻ｎ＋１における粒子速度ｗにかかる速度格子点の粒子速度を読み出した解法（Ｓｗ）に従って算出する（ステップＳ２５４）。
【０１０５】
その後、時刻がΔｔ/２増加（ステップＳ２５５）させ、予め定めた解析時間に達するまで一連（ステップＳ２５０〜Ｓ２５５）の処理を繰り返す（ステップＳ２５６のＮＯ）。各過程で算出された弾性波解析値は、記憶部３０に記憶される。
【０１０６】
このように、格子点ごとに選択した解法を解法情報記憶領域３７に記憶した解法に従って次時刻の圧力や粒子速度を、逐一解法を選択せずに効率的に弾性波解析値を求めることができる。
【０１０７】
〔圧力算出ステップ〕
解法Ｓｐ＝２は、例えば項番７（図１２−１（ａ））のようにｘ軸の負方向の隣接格子点（ｉｘ,ｉｙ−１/２,ｉｚ）が存在せず、ｙ，ｚ軸方向における隣接格子点が存在する場合である（図１−１（ｂ）参照）。この場合は、式（１６）で次時刻の圧力を算出する。
【０１０８】
【数１６】
2011112612.tif 000018
【０１０９】
即ち、ｘ軸方向の寄与については式（７）、ｙ，ｚ軸方向の寄与については従来のＦＤ
ＴＤ法を用いることで算出できる。
【０１１０】
解法Ｓｐ＝２が選択された他の隣接格子点有無情報（Ｘ正, Ｘ負, Ｙ正, Ｙ負, Ｚ正, Ｚ負）の組み合わせの場合（項番２〜６）も、式（１６）と同様の考え方を他の座標軸、正負方向を違えて適用することで次時刻の圧力を求めることができる。
【０１１１】
解法Ｓｐ＝３は、例えば、項番２２のようにｘ軸の正負方向の隣接格子点が無い場合（図１−２（ｃ）参照）は、隣接格子点の存する座標軸方向についてのみＦＤＴＤ法で解析する（式（１７））。
【０１１２】
【数１７】
2011112612.tif 000019
【０１１３】
解法Ｓｐ＝３が選択された他の組み合わせの場合（項番８、１３）も、式（１７）と同様の考え方で座標軸方向を違えて次時刻の圧力を求めることができる。
【０１１４】
解法Ｓｐ＝４は、２座標軸方向について、各々正負方向何れかの隣接格子点がない場合（例えば図１−２（ｄ））に選択される。解法Ｓｐ＝４は、正負何れか隣接格子点の座標軸方向への成分について圧力を推定してその方向への時間差分の寄与分を各々求める。その他の座標軸方向について１次元ＦＤＴＤ法で求めた値に基づき寄与分を求め、各座標軸に亘り合成する。
【０１１５】
項番２１のようにｘ,ｙ軸上の負方向に各々隣接格子点がなく、ｚ軸上の正負両方向に隣接格子点がある場合（例えば図１−２（ｄ））には、上記した式（１４）に基づき次時刻の圧力を求めることができる。解法Ｓｐ＝４が選択された他の組み合わせの場合（項番９〜１２、１４〜２０）も、式（１４）と同様の考え方で次時刻の圧力を求めることができる。
【０１１６】
解法Ｓｐ＝５は、１座標軸方向のみ正負両方向とも隣接格子点がなく、１座標軸上の正負何れかの方向に隣接格子点があり、さらに１座標軸上の正負両方向とも隣接格子点がある場合（例えば図１−３（ｅ））に選択される。ここでは、正負何れかの方向に隣接格子点がある座標軸方向への音圧を推定して、その方向への時間差分の寄与分を求める。正負両方向とも隣接格子点のある座標軸方向については、その方向への寄与分をＦＤＴＤ法で求め、これらを合成する。例えば、項番４２のように、ｙ軸上の負方向の隣接格子点がなく、ｚ軸上の両方向の隣接格子点が存する場合には、式（１８）に基づき次時刻の圧力を求めることができる。
【０１１７】
【数１８】
2011112612.tif 000020
【０１１８】
解法Ｓｐ＝５が選択された他の組み合わせの場合（項番２３〜２７、３２、３３、３８〜４１）も、式（１８）と同様の考え方で次時刻の圧力を求めることができる。
【０１１９】
解法Ｓｐ＝６は、ｘ,ｙ,ｚ各座標軸について正負方向何れか一方の隣接格子点が存する場合（例えば図１−３（ｆ））に選択される。ここでは、各座標軸方向への圧力を推定してその方向への時間差分を求め、これらを合成する。例えば、項番３７のようにｘ,ｙ,ｚ軸上各々について正方向の隣接格子点が存する場合には、上記した式（１５）に基づき次時刻の圧力を求めることができる。
【０１２０】
解法Ｓｐ＝６が選択された他の組み合わせの場合（項番２８〜３１、３４〜３６）も、式（１５）と同様の考え方で次時刻の圧力を求めることができる。
【０１２１】
解法Ｓｐ＝７は、１座標軸方向のみ正負方向両者の隣接格子点が存する場合（例えば図１−４（ｇ））に選択される。ここでは、隣接格子点の存する座標軸方向についてのみＦＤＴＤ法で解析する。例えば、項番５７のようにｚ軸方向についてのみ隣接格子点が存する場合には、式（１９）に基づき次時刻の圧力を求めることができる。
【０１２２】
【数１９】
2011112612.tif 000021
【０１２３】
解法Ｓｐ＝７が選択された他の組み合わせの場合（項番４３、５２）も、式（１９）と同様の考え方で次時刻の圧力を求めることができる。
【０１２４】
解法Ｓｐ＝８は、ｘ,ｙ,ｚ座標軸の内２つの座標軸上の正負方向何れか一方の隣接格子点が存する場合（例えば図１−４（ｈ））に選択される。ここでは、各座標軸方向への圧力を推定してその方向への時間差分の寄与分を求め、これらを合成する。例えば、項番５３のようにｙ,ｚ軸上各々について正方向の隣接格子点が存する場合には、式（２０）に基づき次時刻の圧力を求めることができる。
【０１２５】
【数２０】
2011112612.tif 000022
【０１２６】
解法Ｓｐ＝８が選択された他の組み合わせの場合（項番４４〜５１、５４〜５６）も、式（２０）と同様の考え方で次時刻の圧力を求めることができる。
【０１２７】
解法Ｓｐ＝９は、ｘ,ｙ,ｚ座標軸の内、１の座標軸上の正負方向何れか一方の隣接格子
点が存する場合（例えば図１−５（ｉ））に選択される。ここでは、その座標軸方向への
圧力を推定する。例えば、項番６３のようにｚ軸上の正方向の隣接格子点のみが存する場
合には、式（２１）に基づき次時刻の圧力を求めることができる。解法Ｓｐ＝９が選択さ
れた他の組み合わせの場合（項番５８〜６２）も、式（２１）と同様の考え方で次時刻の
圧力を求めることができる。
【０１２８】
【数２１】
2011112612.tif 000023
【０１２９】
解法Ｓｐ＝１０は、ｘ,ｙ,ｚ何れの座標軸上にも隣接格子点が存在しない場合に選択される。この場合（図１−５（ｊ））、格子点が孤立しているため波動の伝搬を解析することができない。したがって、式（２２）のように現時刻の圧力を次時刻の圧力としても良い。
【０１３０】
【数２２】
2011112612.tif 000024
【０１３１】
この場合（項番６４）は、不適切な領域設定であるとして解析値表示部４０にその旨を表示したうえで、動作を停止するようにしても良い。以上述べた１０種の解法に基づいて、全ての圧力格子点の次時刻の圧力を求めることが可能である。
【０１３２】
粒子速度・ｕｎ＋１（ｉｘ＋１/２,ｉｙ,ｉｚ），・ｖｎ＋１（ｉｘ,ｉｙ＋１/２,ｉｚ），・ｗｎ＋１（ｉｘ，ｉｙ，ｉｚ＋１/２）についても、同様に隣接格子点の有無の組み合わせに応じた解法の選択により求めることができる。
【０１３３】
〔粒子速度の計算方法〕
ここで、ｘ軸方向で隣接する圧力格子点が解析領域外である境界面上の速度格子点の次時刻の粒子速度の計算方法について説明する。この実施例では、境界面上の速度格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）の正方向に隣接する圧力格子点（ｉｘ＋１,ｉｙ,ｉｚ）が存在せず、その逆方向に隣接する圧力格子点（ｉｘ,ｉｙ,ｉｚ）が、さらに隣接して速度格子点（ｉｘ−１/２,ｉｙ,ｉｚ）が存する場合に、その境界面上の速度格子点（ｉｘ＋１/２，ｉｙ，ｉｚ）の粒子速度を、最も簡単な一次形式として式（２３）に基づき算出する。
【０１３４】
【数２３】
2011112612.tif 000025
【０１３５】
ここで、式（２３）の内、式（２４）を満たす成分は、時刻ｎ＋１/２においてｘ軸の正方向へ波面進行速度ｃで反射波が生じずに伝搬する条件である式（２５）を離散化して得られる。
【０１３６】
【数２４】
2011112612.tif 000026
【０１３７】
また、式（２３）の内、式（２６）を満たす成分は、式（８）と同様に時間経過と共に粒子速度・ｕｎ＋１（ｉｘ＋１/２，ｉｙ，ｉｚ）は零に近似する。もっとも初期値が零の場合には常に零となり、初期値の設定において・ｕｎ＋１（ｉｘ＋１/２,ｉｙ,ｉｚ）＝０としてもよい。
【０１３８】
【数２５】
2011112612.tif 000027
【０１３９】
この値が零となると境界面上の速度格子点（ｉｘ＋１/２，ｉｙ，ｉｚ）では、ｘの正方向に進行速度ｃで伝搬する入射波成分ｐ（ｃｔ−ｋ（ｘ−（ｉｘ＋１/２）Δｘ））が生じる場合に、その逆の負方向に進行速度ｃで伝搬する反射波成分ｐ（ｃｔ＋ｋ（ｘ−（ｉｘ＋１/２）Δｘ））が同位相で生じる条件を離散化したものである。
【０１４０】
即ち、ｘ軸方向の粒子速度・ｕ（ｘ，ｔ）の時間微分（つまり単位体積当りの力）は、一般にその方向についての圧力ｐ（ｘ，ｔ）の傾き（式（２７））で表される。
【０１４１】
【数２６】
2011112612.tif 000028
【０１４２】
これを境界面において常に零、つまり・ｕ（（ｉｘ＋１/２）Δｘ，ｔ）＝０とすると、式（２８）に示す関係が成り立つ。つまり境界面への入射波ｐ（ｃｔ−ｋ（ｘ−（ｉｘ＋１/２）Δｘ）））があれば、同相の反射波成分ｐ（ｃｔ＋ｋ（ｘ−（ｉｘ＋１/２）Δｘ）））が生じる。
【０１４３】
【数２７】
2011112612.tif 000029
【０１４４】
ここでｋは波数（ｋ＝２π/λ、λは波長）である。
【０１４５】
以上まとめると、式（２３）のうち式（２４）の部分により境界面への入射波に相当する成分が与えられる。また、式（２６）により、又は境界面上の格子点（ｉｘ＋１/２,ｉｙ,ｉｚ）で粒子速度・ｕ（ｉｘ＋１/２,ｉｙ,ｉｚ）＝０と置くことで入射波を受けて境界面から媒質内部に同位相の反射波に相当する成分が与えられる。
【０１４６】
図１６にこの入射波と反射波を模式的に示す。横方向は波面の進行方向、縦方向は波面の振幅を示す。入射波を実線、反射波を破線で示す。式（２３）は、同位相の反射波に相当する成分に反射率Ｒ１を乗じ、入射波に相当する成分と合成するという条件を与える。よって、入射波ｐ（ｃｔ−ｋ（ｘ−（ｘ＋１/２）Δｘ））に応じた反射率Ｒ１で同位相の反射波成分ｐ（ｃｔ＋ｋ（ｘ−（ｘ＋１/２）Δｘ））が合成されているとの条件のもとで、境界面上の格子点における粒子速度・ｕ（ｉｘ＋１/２,ｉｙ,ｉｚ）を求めることができる。
【０１４７】
式（２３）は、反射率Ｒ１が０≦Ｒ１≦１の場合でなければ適用できない。何故ならば反射率Ｒ１が負値の場合、式（２３）により求められる粒子速度・ｕｎ＋１（ｉｘ＋１/２,ｉｙ,ｉｚ）は不安定になるからである。
【０１４８】
その理由を説明する。ここで、空間上の任意の点ｒにおける粒子速度・ｕ（ｒ,ｔ）が時刻ｔにおける格子点上の値・ｕｒ＝（ｉｘ＋１/２）Δｘ,ｉｙΔｙ,ｉｚΔｚ）（ｔ）に応じた平面波ｅ−ｊｋｒで表わされるものとする。つまり、・ｕ（ｒ,ｔ）＝・ｕｒ＝（ｉｘ＋１/２）Δｘ,ｉｙΔｙ,ｉｚΔｚ）（ｔ）ｅ−ｊｋｒとおく。このとき、式（２７）は式（２７ａ）と変形される。
【０１４９】
【数２８】
2011112612.tif 000030
【０１５０】
両辺を整理すると式（２７ｂ）となる。
【０１５１】
【数２９】
2011112612.tif 000031
【０１５２】
即ち、・ｕｎ＋１を安定に求めるには式（２７ｃ）を満たすことが必要となる。
【０１５３】
【数３０】
2011112612.tif 000032
【０１５４】
ここで、０＜ｄ１≪１,０＜ｂｘ＜１であることを考慮すると式（２７ｄ）を満たせばよい。
【０１５５】
【数３１】
2011112612.tif 000033
【０１５６】
式（２７ｄ）によれば、そのようなＲ１は、０≦Ｒ１≦２となるが、一般に反射率は−１≦Ｒ１≦１であるため、両者が重複する範囲をとると０≦Ｒ１≦１となる。このようにして得た境界面上の速度格子点の次時刻の粒子速度は、波動方程式を満たすので境界面上の格子点の解析値を安定して推定することができる。
【０１５７】
〔粒子速度算出ステップ〕
解法Ｓｕ＝１が選択された場合は、通常のＦＤＴＤ法により次時刻の粒子速度を求める（式（２））。
【０１５８】
解法Ｓｕ＝２が選択された場合は、ｘ軸の負方向に反射率Ｒ１で同位相の反射波が生じる状況のもとで、式（２５）を格子点（ｉｘ＋１/２，ｉｙ，ｉｚ）を基準に正負方向を逆に適用することで、次時刻の粒子速度を算出する。
【０１５９】
解法Ｓｕ＝３が選択された場合は、速度格子点が孤立しているため、波動の伝播を解析することができない。便宜上、現時刻の粒子速度と等しい一定値としても良い（式（２９））。
【０１６０】
【数３２】
2011112612.tif 000034
【０１６１】
または、このような場合は、不適切な領域設定であるとして解析値表示部４０にその旨を表示したうえで、動作を停止するようにしても良い。ｙ，ｚ軸方向の粒子速度も同様に計算できる。よって、説明は省略する。
【０１６２】
なお、境界面と波面の入射方向により境界面に接する両媒質に応じた反射率Ｒ１（式（１１））を一の媒質への格子点の属否判定又は隣接格子点の有無判定、つまり境界面が判別されるまでに設定しておく必要がある。
【０１６３】
境界面における解析値の推定は式（７）による圧力推定、式（２３）による粒子速度推定は各々、反射率Ｒ１が負値の場合、正値の場合とで使い分ける必要がある。反射率Ｒ１を与える式（１１）によれば一方の媒質１から他方の媒質２への入射波の反射率Ｒ１と、媒質２から媒質１への入射波の反射率Ｒ２とは、Ｒ１＝−Ｒ２となる関係がある。
【０１６４】
そのため、何れか一方の媒質側から境界面上における解析値を式(７)、式（２３）では求められない。例えば、境界面が圧力を取る格子点で媒質１から２への反射率が正値をとる場合が該当する。その場合でも、他方の媒質２から１の反射率は負値をとるため、その解析値（圧力）を式（７）で求めることができるため、その値を用いれば特に問題は生じない。式（２３）の場合も同様である。むしろ、両方向から相矛盾する解析値を推定することが避けられ、一つの境界面における安定した解析値を一意に求められる。
【０１６５】
そのために、解法を判定する際、一方の媒質１と他方の媒質２の間の反射率Ｒ１が正値で一方の媒質１において式（７）を用いて圧力ｐを求めることが不適切な場合、その旨を示す上記の１〜６４の値とは別途の識別値を解法Ｓｐ（ｉｘ，ｉｙ，ｉｚ）として設定しておく。又は反射率Ｒ１が負値で式（２３）を用いて粒子速度・ｕ，・ｖ，・ｗを求めることが不適切な場合、その旨を示す上記の１〜３とは別の識別値を解法Ｓｕ（ｉｘ＋１/２，ｉｙ，ｉｚ）等として設定しておく。そして、一方の媒質１において弾性波解析値を算出する際、その解法Ｓｐ（ｉｘ，ｉｙ，ｉｚ）、解法Ｓｕ（ｉｘ＋１/２，ｉｙ，ｉｚ）等で設定した識別値が式（７）、式（２３）を用いて算出することが不適切であることを示す場合、他方の媒質２において算出された圧力ｐｎ＋１/２（ｉｘ，ｉｙ，ｉｚ）、粒子速度・ｕｎ＋１（ｉｘ＋１/２，ｉｙ，ｉｚ）等を利用する。
【０１６６】
解法Ｓｐ＝１が選択された圧力格子点（ｉｘ,ｉｙ,ｉｚ）の次時刻の圧力ｐｎ＋１/２（ｉｘ,ｉｙ,ｉｚ）は、通常の３次元ＦＤＴＤ法に基づいて算出される。
【０１６７】
以上述べたように、この発明の弾性波解析方法によれば、各格子点がどの媒質に属するか否かを表す識別値に基づいて、同一媒質に属す各格子点の隣接する格子点を持つか否かを判定し、各格子点の隣接格子点有無情報を生成する。そして、その隣接格子点有無情報に対応した弾性波の解析方法を選択する。したがって、従来技術のような解析値の受け渡しを行う必要が無いので、複雑な形状の解析空間における弾性波解析をより容易に行うことができる。
【０１６８】
また、上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。
【０１６９】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）/ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１７０】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１７１】
また、各手段は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図２】
2011112612.tif 000035
【図３】
2011112612.tif 000036
【図４】
2011112612.tif 000037
【図７】
2011112612.tif 000038
【図８】
2011112612.tif 000039
【図９】
2011112612.tif 000040
【図１０】
2011112612.tif 000041
【図１１】
2011112612.tif 000042
【図１２−１】
2011112612.tif 000043
【図１２−２】
2011112612.tif 000044
【図１２−３】
2011112612.tif 000045
【図１３】
2011112612.tif 000046
【図１６】
2011112612.tif 000047
【図１８】
2011112612.tif 000048
【図１−１】
2011112612.tif 000049
【図１−２】
2011112612.tif 000050
【図１−３】
2011112612.tif 000051
【図１−４】
2011112612.tif 000052
【図１−５】
2011112612.tif 000053
【図５】
2011112612.tif 000054
【図６】
2011112612.tif 000055
【図１４】
2011112612.tif 000056
【図１５】
2011112612.tif 000057
【図１７】
2011112612.tif 000058
【図１９】
2011112612.tif 000059