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