(57)【特許請求の範囲】
【請求項１】
地盤の変形解析を粒子法を用いて行う地盤変形解析装置であって、
立体形状を有する要素領域に関する要素情報として、少なくとも前記要素領域の識別情報、前記要素領域を形成する節点の位置を記憶し、粒子に関する粒子情報として、少なくとも粒子の識別情報、粒子の位置する要素領域の識別情報、位置、質量、密度、応力を記憶する記憶手段と、
前記粒子の位置と要素情報に基づき、前記粒子が、前記立体形状を有する要素領域の一面に対し、前記一面を形成する所定の節点を原点とし、前記一面を形成し前記節点の隣に位置する２つの節点のそれぞれへ向かう２つのベクトルの外積で示される方向に位置するかを判定し、前記要素領域の全ての面について前記粒子が前記外積で示される方向に位置する場合に、前記粒子の位置する要素領域として前記要素領域を算出し、前記要素領域内の前記粒子の数を算出する初期要素領域算出手段と、
内挿関数を用いて前記粒子の質量を基に前記節点の初期節点力を算出し、勾配関数を用いて前記節点の初期節点力、前記粒子の質量、密度を基に前記粒子の初期鉛直方向応力を算出し、前記粒子の初期鉛直方向応力を基に前記粒子の初期応力を算出する初期応力算出手段と、
前記要素領域内について、内挿関数を用いて前記粒子の質量を基に前記節点の質量を算出し、内挿関数および勾配関数を用いて前記粒子の質量、密度、応力を基に前記節点の節点力を算出し、前記節点の質量と前記節点の節点力を基に前記節点の加速度を算出し、前記節点の加速度を基に前記節点の第１の速度を算出する節点速度算出手段と、
前記要素領域内について、内挿関数を用いて前記節点の第１の速度を基に前記粒子の位置を算出し、内挿関数を用いて前記節点の加速度を基に前記粒子の速度を算出し、勾配関数を用いて前記粒子の速度を基に前記粒子のひずみ増分を算出する位置・ひずみ増分算出手段と、
前記要素領域内について、前記粒子のひずみ増分を基に前記粒子のひずみ、応力、密度を算出するひずみ・応力・密度算出手段と、
前記粒子の位置と前記要素情報を基に、前記粒子の平面位置に対応し最も下方に位置する要素領域の底面の鉛直方向座標と前記要素領域の鉛直方向の長さから、前記粒子が位置するとする要素領域を算出し、前記粒子が前記要素領域の外にある場合には、前記粒子の位置を、前記粒子の位置に応じて、前記要素領域の外面上の所定の位置へと変換する要素領域算出手段と、
を具備することを特徴とする地盤変形解析装置。
【請求項２】
前記粒子のひずみまたは前記粒子のひずみ増分を用いて、水圧を算出する水圧算出手段を更に具備することを特徴とする請求項１に記載の地盤変形解析装置。
【請求項３】
前記節点速度算出手段は、一部の節点について、ダッシュポットを模した運動方程式により前記節点の加速度を算出することを特徴とする請求項１または請求項２に記載の地盤変形解析装置。
【請求項４】
立体形状を有する要素領域に関する要素情報として、少なくとも前記要素領域の識別情報、前記要素領域を形成する節点の位置を記憶し、粒子に関する粒子情報として、少なくとも粒子の識別情報、粒子の位置する要素領域の識別情報、位置、質量、密度、応力を記憶する記憶手段を有し、地盤の変形解析を粒子法を用いて行う情報処理装置が、
前記粒子の位置と要素情報に基づき、前記粒子が、前記立体形状を有する要素領域の一面に対し、前記一面を形成する所定の節点を原点とし、前記一面を形成し前記節点の隣に位置する２つの節点のそれぞれへ向かう２つのベクトルの外積で示される方向に位置するかを判定し、前記要素領域の全ての面について前記粒子が前記外積で示される方向に位置する場合に、前記粒子の位置する要素領域として前記要素領域を算出し、前記要素領域内の前記粒子の数を算出する初期要素領域算出ステップと、
内挿関数を用いて前記粒子の質量を基に前記節点の初期節点力を算出し、勾配関数を用いて前記節点の初期節点力、前記粒子の質量、密度を基に前記粒子の初期鉛直方向応力を算出し、前記粒子の初期鉛直方向応力を基に前記粒子の初期応力を算出する初期応力算出ステップと、
前記要素領域内について、内挿関数を用いて前記粒子の質量を基に前記節点の質量を算出し、内挿関数および勾配関数を用いて前記粒子の質量、密度、応力を基に前記節点の節点力を算出し、前記節点の質量と前記節点の節点力を基に前記節点の加速度を算出し、前記節点の加速度を基に前記節点の第１の速度を算出する節点速度算出ステップと、
前記要素領域内について、内挿関数を用いて前記節点の第１の速度を基に前記粒子の位置を算出し、内挿関数を用いて前記節点の加速度を基に前記粒子の速度を算出し、勾配関数を用いて前記粒子の速度を基に前記粒子のひずみ増分を算出する位置・ひずみ増分算出ステップと、
前記要素領域内について、前記粒子のひずみ増分を基に前記粒子のひずみ、応力、密度を算出するひずみ・応力・密度算出ステップと、
前記粒子の位置と前記要素情報を基に、前記粒子の平面位置に対応し最も下方に位置する要素領域の底面の鉛直方向座標と前記要素領域の鉛直方向の長さから、前記粒子が位置するとする要素領域を算出し、前記粒子が前記要素領域の外にある場合には、前記粒子の位置を、前記粒子の位置に応じて、前記要素領域の外面上の所定の位置へと変換する要素領域算出ステップと、
を具備することを特徴とする地盤変形解析方法。
【請求項５】
前記情報処理装置が、前記粒子のひずみまたは前記粒子のひずみ増分を用いて、水圧を算出する水圧算出ステップを更に行うことを特徴とする請求項４に記載の地盤変形解析方法。
【請求項６】
前記節点速度算出ステップでは、一部の節点について、ダッシュポットを模した運動方程式により前記節点の加速度を算出することを特徴とする請求項４または請求項５に記載の地盤変形解析方法。
【請求項７】
情報処理装置を、請求項１から請求項３のいずれかに記載の地盤変形解析装置として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は地盤変形についての解析を行うための地盤変形解析装置、地盤変形解析方法、プログラムに関する。
【背景技術】
【０００２】
例えば地震時等には、地盤に対する動的な作用の下、斜面崩壊や液状化など地盤の大変形を伴う災害が多発する。その際、周辺環境の安全性を確保するためには、それらの災害による被害程度を予測することが重要になる。これらの被害予測には、地盤の変形解析が不可欠であるとともに、３次元的な変形解析や、動的解析が可能であることが望ましい。
【０００３】
このような地盤の変形解析に際しては、従来、有限要素法による変形解析が行われている。このような例が、特許文献１に記載されており、土の構成則の体積変形に抵抗する水の作用を簡便にモデル化し、砂地盤等に対する液状化の影響を考慮した地震応答解析を行うことが記載されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００６−０３８４５５公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、従来用いられているような有限要素法では、地盤の大変形解析に対応できないという問題がある。例えば有限要素法では、要素を形成する節点の移動により地盤の変形を表現するが、要素の大変形により計算の破綻を生じるなどの問題から、地盤が大きく変位する上記のような大変形時の解析に用いることが難しい。
【０００６】
一方、近年、粒子法と呼ばれる、物体の大変形解析に対応可能な手法が開発されてきている。ここでは、粒子法について概略を簡単に説明する。図１４は、粒子法であるＭＰＭ（Ｍａｔｅｒｉａｌ  Ｐｏｉｎｔ  Ｍｅｔｈｏｄ）の概要を模式的に示した図である。
【０００７】
粒子法は陽解法の時刻歴解析法である。図１４（ａ）に示すように、粒子情報としての物質情報（ラグランジュ変数）を有する個別の粒子１００（ラグランジュ粒子）は、空間の要素１０５上、あるいは要素１０５間を移動する。粒子１００がもつラグランジュ変数である物質情報は、一定時間刻み毎に粒子１００の存在する要素１０５の節点１２０に内挿関数等を通して集約される。そして、この節点１２０に対し運動方程式を解いて次ステップでの節点１２０の速度増分を求める。
【０００８】
この時点で、図１４（ｂ）で示すように、節点１２０の移動により要素１０５は粒子１００を載せながら変形し、粒子１００が移動する。また、節点１２０の物質情報より内挿関数等を通して粒子１００の物質情報が更新される。そして、図１４（ｃ）に示すように、変形した要素１０５は次ステップに備え、移動した粒子１００を残して再び元の位置に戻される。以上の過程を繰り返し、粒子１００は位置を移動するとともに、その物質情報が更新されてゆく。
【０００９】
但し、上記のような地盤の３次元的変形や動的解析の問題を考慮した場合に、粒子法をそのまま用いることはできない。即ち、３次元化への対応、加えて動的解析や地下水の挙動への対応等を行うことが必要となる。
【００１０】
本発明は、前述した問題点に鑑みてなされたもので、地盤の３次元的な大変形解析や動的解析を可能とする地盤変形解析装置等を提供することを目的とし、本発明者らが鋭意検討し考案した、粒子法であるＭＰＭについて、３次元化への対応、および動的解析や地下水の挙動への対応を実現した手法を用いることで、従来のＦＥＭ解析では不可能であった地盤の３次元的な大変形解析等を行うことを可能としたものである。
【課題を解決するための手段】
【００１１】
前述した目的を達するための第１の発明は、地盤の変形解析を粒子法を用いて行う地盤変形解析装置であって、立体形状を有する要素領域に関する要素情報として、少なくとも前記要素領域の識別情報、前記要素領域を形成する節点の位置を記憶し、粒子に関する粒子情報として、少なくとも粒子の識別情報、粒子の位置する要素領域の識別情報、位置、質量、密度、応力を記憶する記憶手段と、前記粒子の位置と要素情報に基づき、前記粒子が、前記立体形状を有する要素領域の一面に対し、前記一面を形成する所定の節点を原点とし、前記一面を形成し前記節点の隣に位置する２つの節点のそれぞれへ向かう２つのベクトルの外積で示される方向に位置するかを判定し、前記要素領域の全ての面について前記粒子が前記外積で示される方向に位置する場合に、前記粒子の位置する要素領域として前記要素領域を算出し、前記要素領域内の前記粒子の数を算出する初期要素領域算出手段と、内挿関数を用いて前記粒子の質量を基に前記節点の初期節点力を算出し、勾配関数を用いて前記節点の初期節点力、前記粒子の質量、密度を基に前記粒子の初期鉛直方向応力を算出し、前記粒子の初期鉛直方向応力を基に前記粒子の初期応力を算出する初期応力算出手段と、前記要素領域内について、内挿関数を用いて前記粒子の質量を基に前記節点の質量を算出し、内挿関数および勾配関数を用いて前記粒子の質量、密度、応力を基に前記節点の節点力を算出し、前記節点の質量と前記節点の節点力を基に前記節点の加速度を算出し、前記節点の加速度を基に前記節点の第１の速度を算出する節点速度算出手段と、前記要素領域内について、内挿関数を用いて前記節点の第１の速度を基に前記粒子の位置を算出し、内挿関数を用いて前記節点の加速度を基に前記粒子の速度を算出し、勾配関数を用いて前記粒子の速度を基に前記粒子のひずみ増分を算出する位置・ひずみ増分算出手段と、前記要素領域内について、前記粒子のひずみ増分を基に前記粒子のひずみ、応力、密度を算出するひずみ・応力・密度算出手段と、前記粒子の位置と前記要素情報を基に、前記粒子の平面位置に対応し最も下方に位置する要素領域の底面の鉛直方向座標と前記要素領域の鉛直方向の長さから、前記粒子が位置するとする要素領域を算出し、前記粒子が前記要素領域の外にある場合には、前記粒子の位置を、前記粒子の位置に応じて、前記要素領域の外面上の所定の位置へと変換する要素領域算出手段と、を具備することを特徴とする地盤変形解析装置である。
【００１２】
第１の発明の地盤変形解析装置は、前記粒子のひずみまたは前記粒子のひずみ増分を用いて、水圧を算出する水圧算出手段を更に具備してもよい。
【００１３】
前記節点速度算出手段は、一部の節点については、ダッシュポットを模した運動方程式により前記節点の加速度を算出するものであってもよい。
【００１４】
前述した目的を達するための第２の発明は、立体形状を有する要素領域に関する要素情報として、少なくとも前記要素領域の識別情報、前記要素領域を形成する節点の位置を記憶し、粒子に関する粒子情報として、少なくとも粒子の識別情報、粒子の位置する要素領域の識別情報、位置、質量、密度、応力を記憶する記憶手段を有し、地盤の変形解析を粒子法を用いて行う情報処理装置が、前記粒子の位置と要素情報に基づき、前記粒子が、前記立体形状を有する要素領域の一面に対し、前記一面を形成する所定の節点を原点とし、前記一面を形成し前記節点の隣に位置する２つの節点のそれぞれへ向かう２つのベクトルの外積で示される方向に位置するかを判定し、前記要素領域の全ての面について前記粒子が前記外積で示される方向に位置する場合に、前記粒子の位置する要素領域として前記要素領域を算出し、前記要素領域内の前記粒子の数を算出する初期要素領域算出ステップと、内挿関数を用いて前記粒子の質量を基に前記節点の初期節点力を算出し、勾配関数を用いて前記節点の初期節点力、前記粒子の質量、密度を基に前記粒子の初期鉛直方向応力を算出し、前記粒子の初期鉛直方向応力を基に前記粒子の初期応力を算出する初期応力算出ステップと、前記要素領域内について、内挿関数を用いて前記粒子の質量を基に前記節点の質量を算出し、内挿関数および勾配関数を用いて前記粒子の質量、密度、応力を基に前記節点の節点力を算出し、前記節点の質量と前記節点の節点力を基に前記節点の加速度を算出し、前記節点の加速度を基に前記節点の第１の速度を算出する節点速度算出ステップと、前記要素領域内について、内挿関数を用いて前記節点の第１の速度を基に前記粒子の位置を算出し、内挿関数を用いて前記節点の加速度を基に前記粒子の速度を算出し、勾配関数を用いて前記粒子の速度を基に前記粒子のひずみ増分を算出する位置・ひずみ増分算出ステップと、前記要素領域内について、前記粒子のひずみ増分を基に前記粒子のひずみ、応力、密度を算出するひずみ・応力・密度算出ステップと、前記粒子の位置と前記要素情報を基に、前記粒子の平面位置に対応し最も下方に位置する要素領域の底面の鉛直方向座標と前記要素領域の鉛直方向の長さから、前記粒子が位置するとする要素領域を算出し、前記粒子が前記要素領域の外にある場合には、前記粒子の位置を、前記粒子の位置に応じて、前記要素領域の外面上の所定の位置へと変換する要素領域算出ステップと、を具備することを特徴とする地盤変形解析方法である。
【００１５】
第２の発明の地盤変形解析方法は、前記情報処理装置が、前記粒子のひずみまたは前記粒子のひずみ増分を用いて、水圧を算出する水圧算出ステップを更に行ってもよい。
【００１６】
また、前記節点速度算出ステップでは、一部の節点について、ダッシュポットを模した運動方程式により前記節点の加速度を算出するものであってもよい。
【００１７】
前述した目的を達するための第３の発明は、情報処理装置を、第１の発明の地盤変形解析装置として機能させるためのプログラムである。
【００１８】
上記構成により、粒子法を用いて地盤の変形解析を行い、この際、任意の立体形状の要素領域を定め、粒子の位置する要素領域の計算や、初期応力の設定をこれに対応させて行う。粒子法を用いることにより地盤の大変形解析が可能になり、任意の立体形状の要素領域を定め、粒子の位置する要素領域の計算や、初期応力の設定をこれに対応させて行うことにより、３次元的な変形解析が可能になる。また、粒子のひずみやひずみ増分を基に水圧の算出を行うことにより、地下水を考慮した解析が可能になる。また、ダッシュポットを模した運動方程式を用い、動的解析を行うことも可能でなる。
【発明の効果】
【００１９】
本発明により、地盤の３次元的な大変形解析や動的解析を可能とする地盤変形解析装置等を提供することができる。
【図面の簡単な説明】
【００２０】
【図１】地盤変形解析装置１のハードウェア構成を示す図
【図２】粒子１００、要素領域１０５について示す図
【図３】地盤変形解析装置１による地盤変形解析方法の全体の流れを示すフローチャート
【図４】地盤変形解析装置１による地盤変形解析方法の流れを示すフローチャート
【図５】地盤変形解析装置１による地盤変形解析方法の処理について示す図
【図６】地盤変形解析装置１による地盤変形解析方法の処理について示す図
【図７】地盤変形解析装置１による地盤変形解析方法の処理について示す図
【図８】地盤変形解析装置１による地盤変形解析方法の処理について示す図
【図９】地盤変形解析装置１による地盤変形解析方法の全体の流れを示すフローチャート
【図１０】地盤変形解析装置１による地盤変形解析方法の処理について示す図
【図１１】地盤変形解析装置１による地盤変形解析の例を示す図
【図１２】地盤変形解析装置１による地盤変形解析の例を示す図
【図１３】地盤変形解析装置１による地盤変形解析の例を示す図
【図１４】粒子法の概要を示す図
【発明を実施するための形態】
【００２１】
以下、図面を参照しながら、本発明の地盤解析装置等の実施形態について説明する。
まず、図１を参照して、本実施形態の地盤変形解析装置１のハードウェア構成について説明する。地盤変形解析装置１は、例えば、制御部１１、記憶部１３、メディア入出力部１５、周辺機器Ｉ／Ｆ部１７、通信部２１、入力部２３、表示部２５等がバス１９を介して接続されて構成される。
【００２２】
制御部１１は、ＣＰＵ、ＲＯＭ、ＲＡＭ等により構成される。ＣＰＵは、記憶部１３、ＲＯＭ、記録媒体等に格納されるプログラムをＲＡＭ上のワークメモリ領域に呼び出して実行し、バス１９を介して接続された各部を駆動制御する。ＲＯＭは、コンピュータのブートプログラムやＢＩＯＳ等のプログラム、データ等を恒久的に保持する。ＲＡＭは、ロードしたプログラムやデータを一時的に保持するとともに、制御部１１が各種処理を行うために使用するワークエリアを備える。
【００２３】
記憶部１３は、ハードディスクドライブであり、制御部１１が実行するプログラムや、プログラム実行に必要なデータ、オペレーティング・システム等が格納されている。これらのプログラムコードは、制御部１１により必要に応じて読み出されてＲＡＭに移され、ＣＰＵに読み出されて実行される。記憶部１３には、後述する地盤変形解析処理に係るプログラム、データ等が格納される。
【００２４】
メディア入出力部１５（ドライブ装置）は、例えば、フロッピー（登録商標）ディスクドライブ、ＣＤドライブ、ＤＶＤドライブ、ＭＯドライブ等のメディア入出力装置であり、データの入出力を行う。
【００２５】
周辺機器Ｉ／Ｆ（インタフェース）部１７は、周辺機器を接続させるためのポートであり、周辺機器Ｉ／Ｆ部１７を介して周辺機器とのデータの送受信を行う。周辺機器Ｉ／Ｆ部１７は、ＵＳＢ等で構成されており、通常複数の周辺機器Ｉ／Ｆを有する。周辺機器との接続形態は有線、無線を問わない。
【００２６】
通信部２１は、通信制御装置、通信ポート等を有し、ネットワーク等との通信を媒介する通信インタフェースであり、通信制御を行う。
【００２７】
入力部２３は、例えば、キーボード、マウス等のポインティング・デバイス、テンキー等の入力装置であり、入力されたデータを制御部１１へ出力する。
【００２８】
表示部２５は、例えば液晶パネル、ＣＲＴモニタ等のディスプレイ装置と、ディスプレイ装置と連携して表示処理を実行するための論理回路（ビデオアダプタ等）で構成され、制御部１１の制御により入力された表示情報をディスプレイ装置上に表示させる。
【００２９】
バス１９は、各装置間の制御信号、データ信号等の授受を媒介する経路である。
【００３０】
前述したように、地震時等には、地盤に対する動的な作用の下、斜面崩壊や液状化など地盤の大変形を伴う災害が多発する。これらの被害予測には、地盤の変形解析が不可欠であり、その際３次元的な挙動を解析できることが望ましい。しかし、地盤の変形解析に従来用いられているような有限要素法では、地盤の大変形解析に対応することが困難である。
【００３１】
本発明の地盤解析装置１は、上記のような問題を解決し、粒子法であるＭＰＭを３次元的挙動を考慮した地盤の大変形解析に適合させ、地盤の３次元的挙動の解析やその際の動的解析、地下水の挙動の解析等に適用可能な処理を行うものである。以下、その処理について説明する。以下示す各処理は、予め記憶部１３等に記憶されたプログラム、データ等を制御部１１のＣＰＵが読み出すことにより、地盤変形解析装置１により実行される。
【００３２】
本実施形態の地盤変形解析装置１の行う処理においては、図２（ａ）に示すように節点１２０により形成される要素１０５の柱からなる土柱３０の平面的分布により地盤の３次元のモデル化を行う。当該モデルにおいて図３等のフローに示す計算処理を行うことにより、ＭＰＭを用いた地盤の３次元解析を可能としている。また、各粒子１００は、ラグランジュ変数である粒子情報１１０として、図２（ｂ−１）に示すように、粒子の識別情報である粒子番号と紐付いて、粒子の位置する要素番号、位置、速度、応力、ひずみ、密度、質量等の情報を有し、この粒子情報１１０がＲＡＭまたは記憶部１３に記憶されている。
図３に地盤変形解析装置１による地盤変形解析方法におけるアルゴリズムのフローを示す。
【００３３】
図３に示すように、まず、ステップＳ１００において、計算に必要な初期値の入力など、計算に係る初期条件の設定を行う。
【００３４】
即ち、はじめに時刻歴計算を始める前の処理が必要となる。ここでは、位置Ｘｐ（ｋ）＝｛ｘｐ（ｋ）、ｙｐ（ｋ）、ｚｐ（ｋ）｝（ｋは時刻ステップ数を表す）、速度ｖｐ（ｋ）＝｛ｖｘｐ（ｋ）、ｖｙｐ（ｋ）、ｖｚｐ（ｋ）｝Ｔ、ひずみεｐ（ｋ）＝εｉｊ  ｐ（ｋ）、および密度ρｐ（ｋ）、質量Ｍｐ等の初期値（ｋ＝０）を各粒子１００を示す識別情報である粒子番号等と紐付けて入力し、ＲＡＭ、あるいは記憶部１３に記憶する。
【００３５】
また、質量Ｍｐについては、
【数１】
0005221603.tif 000002
の関係がある。Ｖｐ＝Ｌ３／ｎｓ３であり、Ｌ３は要素領域の体積、ｎｓは１要素一辺あたりの粒子の個数である。ＭＰＭでは、物体を一定の質量をもった粒子としてモデル化しているので、質量Ｍｐは時刻歴計算中不変の値となる。
【００３６】
なお、後述する計算処理に必要な定数や、弾塑性モデル等履歴依存性のあるモデルを用いた場合で初期値の設定が必要な値がある場合は、それらの定数や初期値も同様にステップＳ１００で設定しＲＡＭ等に記憶する。
【００３７】
また、各要素１０５の要素領域の設定も行われている。前述したように、各要素領域は節点１２０により形成され、各要素領域は、図２（ａ）に示すように、節点１２０により構成される要素（セル）１０５の柱からなる土柱３０が平面に分布するように設定される。
このため、例えば図２（ｂ−２）に示すように、この要素領域を形成する節点１２０を識別する節点番号が、各要素を識別する要素番号、および各要素領域についての後述する局所節点番号と紐付けて記憶される。また、各節点番号は、図２（ｂ−３）に示すように、その位置と紐付けて記憶される。
以上の情報は、要素情報１１１としてＲＡＭや記憶部１３等に記憶される。本実施形態では、各節点は図２（ａ）に示すＸＹ平面においては矩形状の要素領域をもつオイラー格子を形成するものとする。なお、説明の簡便化のため図２（ａ）では要素１０５の要素領域はＸＺ平面、あるいはＹＺ平面においてもオイラー格子を形成しているが、本実施形態ではこれに限られない。
【００３８】
次に、ステップＳ１１０では、粒子の位置にかかる諸計算が行われる。ステップＳ１１０における処理を示したものが、図４である。
【００３９】
図４に示すように、ステップ１１０では、まず、時刻暦計算を始める前（ｋ＝０）の場合（ステップＳ１１１のＹｅｓ）、ステップＳ１１２において、粒子の位置と要素情報に基づき、各粒子に対応する要素領域（要素番号）を算出し、これを粒子情報１１０として記憶する。これにより各要素領域に位置する粒子が特定され、例えば要素番号と紐付けて粒子番号を記録することもできる。
【００４０】
即ち、初期設定として粒子の存在を定義するために、各要素領域内にある粒子を抽出する必要がある。要素領域は任意の立体形状をしているので、要素情報と粒子の位置を用いて、要素領域を構成する各面の法線を利用し粒子を抽出する。
具体的には、各粒子について、その粒子が、要素領域の一面に対し、当該一面を形成する所定の節点を原点とし、当該一面を形成し当該所定の節点の隣に位置する２つの節点のそれぞれへ向かう２つのベクトルの外積で示される方向に位置するかを判定する。要素領域の全ての面に対し上記の判定を行い、粒子が全ての面に対し上記の外積で示される方向に位置する場合に、その粒子の粒子情報として当該要素領域（初期要素領域）の要素番号を書き込み、その粒子が当該要素領域の内部にあることを特定する。これを以下説明する。
【００４１】
要素領域は６面体で構成するので、図５を参照しながら、６面体の内部にある粒子を抽出する方法について説明する。
【００４２】
まず、６面体の一面である節点１〜４による四辺形の局所座標を外積を用いて計算する。例えば、局所座標系を示す図５（ａ）におけるベクトルｍ１、ｍ２、ｍ３は、図５（ａ）におけるベクトルａ、ｂを用いて、それぞれ
【数２】
0005221603.tif 000003
と算出される。
【００４３】
また、図５（ｂ）に示すベクトルｘ＝（ｘ１、ｘ２、ｘ３）で表される粒子の位置を、ベクトルｍｉによる局所座標系と、ベクトルｅｉによる全体座標系の２つの座標系で表示すると、
【数３】
0005221603.tif 000004
となる。ベクトルｍiで内積を取る事で、局所座標変更後のベクトルＸ＝（Ｘ１、Ｘ２、Ｘ３）については、
【数４】
0005221603.tif 000005
となり、従って、次式
【数５】
0005221603.tif 000006
で、座標ｘ＝（ｘ１、ｘ２、ｘ３）で表される粒子位置を局所座標系における位置Ｘ＝（Ｘ１、Ｘ２、Ｘ３）で算出することができる。ここで、ｍｉｊはｍｉのｅｊ成分であることを示す。
【００４４】
ここで、図５（ｃ）で示すように、粒子ｐの位置が全体座標系でベクトルＳで表され、節点１〜４により６面体である要素領域の１面１０５ｂが構成されているとする。このとき、当該１面１０５ｂを構成する点の一つである点１の全体座標系のベクトルＴ、上式（５）で要素ｍ１１〜ｍ３３によりあらわされる変換行列である行列Ｍを用いて、次式
【数６】
0005221603.tif 000007
により表されるベクトルＬを算出し、上式（５）の成分Ｘ３に対応する成分Ｌ３がＬ３≧０なら、上記ベクトルＳで表される任意の点が、平面１０５ｂに対して図５（ａ）で示したベクトルｍ３側にあることになる。なお、図５（ｃ）で点１が原点、点１から当該点１と隣り合う点２に向かう点１→２のベクトルが図５（ａ）におけるベクトルａ、点１から当該点１と隣り合う点４に向かう点１→４がベクトルｂに対応し、これらを用いて行列Ｍの各要素ｍ１１〜ｍ３３は式（２）で計算される。即ち、上記のベクトルｍ３側とは、点１→２のベクトルと点１→４のベクトルの外積で示される方向である。
【００４５】
下記の条件全てで、ｍ３側に粒子があれば、上記６面体の内部に粒子が存在することになる。なお、（１→２、３、４）とは、点１を前述のベクトルａ、ｂの始点とすることを指す。また、ｍ１（１→２）とは、図５（ａ）に示すベクトルｍ１の向きを図５（ｄ）で示す点１から点２へ向かう方向とすることを指す。その他も同様である。
【数７】
0005221603.tif 000008
例えば２）の場合は点５に原点を移動し、前述の方法により、ｍ３側に粒子が存在するかを判定する。要素領域の全ての面に対し上記の判定を行う１）から６）の条件すべてでｍ３側に粒子があれば上記六面体の内部に粒子が存在することになる。各粒子の位置と各要素領域の要素情報に基づき上記の判定を行い、各粒子が含まれる初期要素領域の要素番号を算出し、各粒子の粒子情報としてＲＡＭ等に記憶する。
【００４６】
上記のステップＳ１１２の後、図４のステップＳ１１４において、各要素領域に含まれる粒子の数を算出する。例えば、粒子情報として当該要素領域の要素番号を有する粒子の数を集計することにより、各要素領域の粒子の数が算出される。
【００４７】
ここで、以降の処理において用いる、各粒子の粒子情報に基づき節点の物理情報を求める、あるいはその逆の計算を行うための内挿関数および勾配関数について説明する。
【００４８】
今、粒子ｐがある要素領域に存在しているとする。図６に示すように、内挿関数等の計算を行うために、全体座標系（ｘ、ｙ、ｚ軸）より局所座標系（ｒ１、ｒ２、ｒ３軸）への変換を行い、変換後の要素をアイソパラメトリック正方要素とすることを考える。局所座標系に変換した要素領域の局所節点番号は図６や図２（ｂ−２）の番号１〜８で示すように各要素領域について予め定められており、粒子ｐの座標Ｘｐ（ｋ）の関数として形状関数Ｓｉ（Ｘｐ（ｋ））（＝Ｓｉｐ（ｋ）＝｛Ｓｘｉｐ（ｋ）、Ｓｙｉｐ（ｋ）Ｓｚｉｐ（ｋ）｝Ｔ）が次式で計算される。
【数８】
0005221603.tif 000009
また、勾配関数∇Ｓｉｐ（ｋ）（＝Ｇｉｐ（ｋ）＝｛Ｇｘｉｐ（ｋ）、Ｇｙｉｐ（ｋ）Ｇｚｉｐ（ｋ）｝Ｔ）が次式で計算される。
【数９】
0005221603.tif 000010
ここで、ｉは局所節点番号である。また、ｒ１、ｒ２、ｒ３はアイソパラメトリック正方要素による局所座標系の３軸を示し、全体座標系に対し次式に示す関係がある。特にＸＹ平面では、
【数１０】
0005221603.tif 000011
の関係がある。ｘｃ、ｙｃはそれぞれ、アイソパラメトリック正方要素の局所座標系の原点の全体座標系における位置を示す。ｌｘ、ｌｙはそれぞれ、アイソパラメトリック正方要素の全体座標系におけるＸ軸方向、Ｙ軸方向の長さの１／２の値である。
以上の内挿関数ならびに勾配関数は、上記した各粒子の粒子位置ならびに各粒子が位置する要素領域の要素情報に基づき求めることができる。アイソパラメトリック正方要素により計算された形状関数は後述の計算処理の際、内挿関数として用いることができる。このため、以降形状関数を示す記号Ｓｉｐ、やＮｉｐｊを指して内挿関数と呼ぶことがある。
【００４９】
図４に戻って、時刻暦計算を始める前（ｋ＝０）では、さらに、前述のステップＳ１１４の後、ステップＳ１１５において、各粒子の応力σｐ（ｋ）の初期値（ｋ＝０、初期応力）を算出し、粒子情報１１０としてＲＡＭ等に記憶する。
【００５０】
本実施形態では、ステップＳ１１５において、粒子の応力を土柱３０の要素１０５の要素領域の応力分布より求める。土柱３０の要素領域では、内部の上下で応力が一定とし、上側の４つの節点の応力で要素領域内の粒子の応力を決定する。初期段階では、上側の節点４つに配置する応力は、粒子の重力に基づくものであり、要素領域１０５内の粒子の初期応力がこれにより決定される。
【００５１】
具体的には、内挿関数を用いて粒子の質量を基に節点の初期節点力を算出し、その後、勾配関数を用いて節点の初期節点力、粒子の質量、密度を基に粒子の初期鉛直応力を算出し、初期鉛直応力より粒子の初期応力を算出する。これを以下説明する。
【００５２】
ステップＳ１１５では、まず、粒子の質量と内挿関数を用いて、粒子の自重に伴う節点力を求める。即ち、ひとつの要素領域に注目すると、粒子の各節点１２０への影響（力）は要素領域の上方の節点１２０に対しては以下の式で、
【数１１】
0005221603.tif 000012
下方の節点１２０に対しては、以下の式で、
【数１２】
0005221603.tif 000013
それぞれ表される。ここで、Ｎｉｐｊは節点ｉ（ｉは局所節点番号）に対する粒子ｊの影響関数であり、節点ｉの変位の形状関数を使用する。前述のようにアイソパラメトリック正方要素により計算された形状関数は内挿関数に等しい。
【００５３】
また、ＭＰｊは粒子ｊの質量、ｇｚは重量加速度である。上記両式は、通常使われているＦＥＭで利用されている式と等価なものである。上式で求められる節点の節点力を鉛直方向の要素領域間で重ね合わせる（足し合わせる）ことで粒子の自重による土柱の初期節点力ｆｉ（ｉ＝１〜４）が算出される。
【００５４】
一方で、通常ＦＥＭで利用されている式
【数１３】
0005221603.tif 000014
より、ＭＰＭにより応力と節点力の関係が以下の式で記述できる。下式で微分記号∂はｚ方向の辺微分とする。なお、上式でεｖｐは粒子のひずみ、σｐｉは粒子の応力である。
【数１４】
0005221603.tif 000015
節点の応力と粒子の応力の関係については、
【数１５】
0005221603.tif 000016
であり、節点力と粒子の応力については
【数１６】
0005221603.tif 000017
の関係がある。
【００５５】
式（１５）と式（１６）を連結すると、
【数１７】
0005221603.tif 000018
となる。行列Ａは応力を節点力に変換する行列を表し、正方行列となるのでＡの逆行列を求めることができ、これを解いて節点力ｆより節点位置の応力σｉが算出できる。ステップＳ１１５では、このようにして（初期）節点力ｆｉより求めた節点位置の（初期）応力σｉを上式（１４）に代入して粒子の初期鉛直応力σｐｉを求める。
【００５６】
なお、Ａの逆行列が求まるためには、ｎｘ≧２、ｎｙ≧２である必要がある。ここで、ｎｘ、ｎｙは、それぞれ一要素内のｘ方向、ｙ方向の粒子数であり、要素領域中の粒子数および要素領域のｘ方向長さおよびｙ方向長さを用いて求めることができる。逆行列が求まらない場合はＡ行列の成分を次式のように対角化して逆行列とする。
【数１８】
0005221603.tif 000019
【００５７】
その後、初期鉛直応力σｐｉと、予め入力条件として地盤変形解析装置１に入力される値である側圧係数Ｋ０より、要素領域内の粒子の初期応力がσｚ＝σＰｉ、σＸ＝σＰｉＫ０、σＹ＝σＰｉＫ０と算出され、粒子情報としてＲＡＭ等に記憶される。ここで、せん断成分はすべてゼロとする。
【００５８】
なお、側圧係数Ｋ０は、ゼロ以外が設定されていれば優先的に利用する。入力がゼロの場合は、弾性計算の場合、ポアソン比νから決定されるＫ０＝ν／（１−ν）としておく。弾塑性解析の場合は、土の内部摩擦角φから決定される静止土圧係数Ｋ０＝１−ｓｉｎφが利用される。
【００５９】
次に、図４に戻り、時刻暦計算中の要素番号の算出について説明する。
後述するが、ＭＰＭでは、時刻ステップ毎に粒子の位置が更新される。従って、ｋ≠０である時刻暦計算中（ステップＳ１１１でＮｏ）では、ステップＳ１１３において、粒子の位置を基に、粒子の含まれる要素領域と要素領域内の局所座標を更新する。
本実施形態では、粒子の位置と要素情報を用いて、粒子のＸＹ平面位置に対応し最も下方に位置する要素領域の底面の鉛直方向座標と、要素領域のＺ方向（鉛直方向）の長さにより、粒子が位置する（とする）要素領域を算出する。そして、粒子が実際には当該要素領域の外にある場合には、粒子の位置に対応して、粒子の位置を要素領域の外面上の所定の位置に変換する。以下図７、図８を参照してこれを説明する。
【００６０】
粒子ｐのＸＹ平面での要素領域の位置は、本実施形態では要素領域がＸＹ平面においてはオイラー格子を形成するため比較的簡単に決定できる。図７（ａ）に示すように、例えば位置ｘ、ｙの粒子のＸ、Ｙ方向の要素領域の位置は、それぞれ（ｘ−Ｘｓ）／Ｌｘ＋１、
（ｙ−Ｙｓ）／Ｌｙ＋１として算定できる。ここで、Ｘｓ、Ｙｓは要素番号（１、１、Ｎｚ）の左下の節点の座標、Ｌｘ、ＬｙはＸ、Ｙそれぞれの方向の要素領域の長さである。
【００６１】
一方、前述したように、本実施形態では図７（ｂ）に示すように粒子の含まれる要素領域はＸＹ平面ではオイラー格子でＺ方向に同じ形状の要素領域が重なっているが、要素領域自体はＺ方向にオイラー格子を形成しない場合もある。そこで、Ｚ方向位置ｚの粒子の、Ｚ方向の要素領域の位置Ｎｚは以下に示す方法により算定する。
【００６２】
即ち、ＸＹ平面は直交格子なので、前述の方法でＸＹ平面上の要素領域を決定することができるが、Ｚ方向については、底面の起伏を考慮するため、ＸＹ平面と同じ手法は使用できない。Ｚ方向については、下記の手順で粒子が含まれる要素領域と要素領域内での局所座標を特定する。
【００６３】
まず、前述の（ｘ−Ｘｓ）／Ｌｘ＋１、（ｙ−Ｙｓ）／Ｌｙ＋１と要素情報を用いて、各土柱の最も底にある要素領域のうち、底面のＸＹ領域が粒子ｐの平面位置（ｘ、ｙ）を含む要素領域１０５ａ−１を算出する。
次に、当該要素領域１０５ａ−１の要素情報を用いて、図７（ｂ）の×印で示す、当該要素領域１０５ａ−１の底面上で、粒子ｐの平面位置（ｘ、ｙ）に対応するＺ方向座標Ｚｅを算出する。そして、粒子ｐの鉛直方向の位置Ｚを用いて以下の式
【数１９】
0005221603.tif 000020
で示す条件が成り立つｎを算出する。Ｌｚは要素領域１０５ａ−１（１０５ａ−２）のＺ方向の高さである。これにより、最も下方の要素領域１０５ａ−１からｎ番目の要素領域に粒子が含まれることが特定できる。このように、粒子の位置と要素情報を基に粒子の位置する要素領域を算出し更新する。例えばｎ＝２の場合、最も下方の要素領域１０５ａ−１から２番目の要素領域１０５ａ−２に粒子ｐが含まれるとする。
【００６４】
要素領域はＺ方向にオイラー格子を形成しないことがあることは前述したが、そのため、図７（ｂ）に示すように、以上のようにして求めた要素領域内に実際に粒子ｐが存在しない場合がある。
【００６５】
そこで、本実施形態では、粒子ｐが含まれる（とする）要素領域の算出後、粒子位置と要素情報を用いて、図６等で説明したようにニュートンラプソン法により要素領域がアイソパラメトリック正方要素になるよう座標変換を行い、その際、図８（ａ）に示す要素領域の底面１３０（領域ＡＡ）に対する、Ｚ軸正側から見た領域ごとに場合わけし、粒子ｐが位置する領域ごとに、以下示すような位置変換を便宜的に行う。
【００６６】
例えば図８（ａ）でＸＹ平面上ＡＡ面内に粒子ｐがある場合、図８（ｂ）のＰ１で示すように、粒子ｐが変換後の要素領域１０５より上にあるときは、粒子ｐから要素領域１０５の上面に垂線を下ろした、当該上面上の点にあるものとする。逆に要素領域１０５より下にあるときは、粒子ｐから要素領域１０５の下面（底面１３０）に向かって垂線を下ろした、当該下面上の点にあるものとする。図８（ｂ）のＰ２で示すように、粒子ｐの高さが要素領域中にあるときは、粒子ｐは要素領域１０５に含まれることになるので、その位置はそのままとする。
【００６７】
一方、図８（ａ）で粒子ｐがＡＡ面内になく、Ａ１〜Ａ４面内にあるとき、例えば図８（ｂ）のＰ３で示すように、粒子ｐが変換後のアイソパラメトリック要素より上にあるときは、粒子ｐから要素領域１０５の上面の対応する辺（上辺）に垂線を下ろした、当該上辺上の点にあるものとする。逆に要素領域１０５より下にあるときは、粒子ｐから要素領域１０５の下面の対応する辺（下辺）に向かって垂線を下ろした、当該下辺上の点にあるものとする。図８（ｂ）のＰ４で示すように、粒子ｐの高さが要素領域中にあるときは、要素領域１０５の対応する側面に向かって垂線を下ろした、当該側面上の点にあるものとする。
【００６８】
さらに、図８（ａ）でＡＡ面内になく、Ｂ１〜Ｂ４面内にあるとき、例えば図８（ｂ）のＰ５で示すように、粒子ｐが変換後のアイソパラメトリック要素より上にあるときは、粒子ｐは要素領域１０５の上面の対応する角（上角）にあるものとする。逆に要素領域１０５より下にあるときは、粒子ｐは要素領域１０５の下面の対応する角（下角）にあるものとする。図８（ｂ）のＰ６で示すように、粒子ｐの高さが要素領域中にあるときは、要素領域１０５の対応する側辺（上下辺）に向かって垂線を下ろした、当該側辺上の点にあるものとする。
【００６９】
このように、変換後の要素領域に対する粒子の位置に応じて、粒子の位置が要素領域の外面上にくるように位置の変換を行う。以上の変換は粒子位置と要素情報を用いて行うことができる。その後、ステップＳ１１４で各要素領域中の粒子数を前述したように算出、更新する。
【００７０】
地盤変形解析のフローについての説明に戻る。ステップＳ１２０では、要素領域内において、節点の質量と節点力を算出し、節点の質量と節点力より節点の加速度を算出し、節点の加速度より節点の第１の速度を算出する。
【００７１】
ステップＳ１２０では、まず、節点の質量ｍｉ（ｋ）を、粒子ｐの質量Ｍｐと内挿関数Ｓｉｐを用いて次式で算出する。
【数２０】
0005221603.tif 000021
ここで、Ｎｐは要素領域に含まれる粒子の総数である。
【００７２】
次に、節点力として、内力と外力を、粒子の質量、応力、密度、体積力、接触力を基に内挿関数、勾配関数を用いて算出する。
【００７３】
まず、内力としての節点力ｆｉｉｎｔ（ｋ）＝｛ｆｘｉｉｎｔ（ｋ）、ｆｙｉｉｎｔ（ｋ）、ｆｚｉｉｎｔ（ｋ）｝Ｔを粒子の質量、応力σｉｊ  ｐ（ｋ）、密度ρｐ（ｋ）を基に勾配関数Ｇｎｉｐを用いて次式
【数２１】
0005221603.tif 000022
で算出する。
【００７４】
一方、外力としての節点力ｆｉｅｘｔ（ｋ）＝｛ｆｘｉｅｘｔ（ｋ）、ｆｙｉｅｘｔ（ｋ）、ｆｚｉｅｘｔ（ｋ）｝ＴをＭＰＭにおける体積力ｂｉ（ｋ）＝｛ｂｘｉ（ｋ）、ｂｙｉ（ｋ）｝Ｔと接触力τｉ（ｋ）＝｛σｐ（ｋ）Ｔ・ｎ｝を用いて、
【数２２】
0005221603.tif 000023
で算出する。ここで、
【数２３】
0005221603.tif 000024
【数２４】
0005221603.tif 000025
である。ｎはｄＳの法線ベクトルを示す。Ｘｐは前述の粒子の位置ベクトルを、∂Ωは接触応力の作用する領域を、ｄＳは単位面積を表す。
【００７５】
次に、節点の質量と節点力を基に節点の加速度を求める。即ち、節点の加速度ａｉ（ｋ）＝｛ａｘｉ（ｋ）、ａｙｉ（ｋ）、ａｚｉ（ｋ）｝Ｔを、次に示す節点の運動方程式
【数２５】
0005221603.tif 000026
より算出する。
【００７６】
最後に、節点の加速度を基に次ステップｋ＋１での節点の第１の速度ｖｉ（ｋ＋１）を、
【数２６】
0005221603.tif 000027
により算出し、更新する。Δｔは時刻ステップｋの長さであり、予め適宜定めておくことができる。
【００７７】
次に、ステップＳ１３０において、粒子の位置Ｘｐ（ｋ＋１）、速度ｖｐ（ｋ＋１）、ひずみ増分Δεｉｊ  ｐ（ｋ＋１）を算出する。
【００７８】
ステップＳ１３０では、まず、粒子の位置を、節点の第１の速度を基に内挿関数を用いて次式で算出し、更新する。
【数２７】
0005221603.tif 000028
【００７９】
また、節点の加速度を基に内挿関数を用いて粒子の速度を次式で算出し、更新する。
【数２８】
0005221603.tif 000029
【００８０】
また、Ｓｕｌｓｋｙにより提案された速度ｖi（ｋ＋１）の平滑化を行い、節点の速度（第２の速度）ｖi（ｋ＋１）を次式で再び求める。即ち、節点の速度（第２の速度）を、粒子の質量、速度、節点の質量を基に内挿関数を用いて次式で求める。
【数２９】
0005221603.tif 000030
【００８１】
最後に、粒子のひずみ増分Δεｉｊ  ｐ（ｋ＋１）を、節点の速度を基に、勾配関数を用いて次式から算出する。なお、ｍ、ｎの添え字は速度の成分を表す。
【数３０】
0005221603.tif 000031
但し、工学ひずみΔγｍｎに関しては、
【数３１】
0005221603.tif 000032
である。
【００８２】
続いて、ステップＳ１４０において、粒子のひずみ増分を基に粒子のひずみεｉｊ  ｐ（ｋ＋１）、応力σｉｊ  ｐ（ｋ＋１）、密度ρｐ（ｋ＋１）を算出し、更新する。
【００８３】
ステップＳ１４０では、まず、物体の剛性マトリックスＴｉｊｋｌを用いて粒子のひずみ増分を基に粒子の応力増分Δσｉｊ  ｐ（ｋ＋１）を次式で算出する。添え字のｋ、ｌはひずみの成分を表す。
【数３２】
0005221603.tif 000033
【００８４】
また、粒子のひずみを、粒子のひずみ増分を基に次式で算出し、更新する。
【数３３】
0005221603.tif 000034
【００８５】
また、粒子の応力を、粒子の応力増分を基に算出し、更新する。
大変形解析を行うときは回転の影響も考慮することが好ましいことが知られている。ここでは、この影響を取り入れるためＪａｕｍａｎｎ率を用いる。このとき、粒子の応力増分に加え、粒子の速度、内挿関数等も用いてステップｋ＋１での粒子ｐの応力σｉｊ  ｐ（ｋ＋１）が次式で算出され、更新される。
【数３４】
0005221603.tif 000035
ここで
【数３５】
0005221603.tif 000036
である。
【００８６】
また、粒子の密度が、粒子のひずみ増分を用いて次式で算出され、更新される。
【数３６】
0005221603.tif 000037
【００８７】
以降、全時刻ステップでの計算を終えるまで（ステップＳ１５０のＮｏ）、時刻ステップ数ｋの値をインクリメントし（ステップＳ１６０）、ステップＳ１１０に戻って次の時刻暦計算を行う。この処理の過程で、各粒子のもつラグランジュ変数である粒子情報が更新される。例えば、粒子のひずみは、ステップＳ１４０において更新される。このようにして予め定めた時刻ステップ数だけ時刻暦計算を繰り返して時刻暦計算を終了すると（ステップＳ１５０のＹｅｓ）、地盤解析を終了し、表示部に位置、ひずみ、応力等、各粒子についての諸値を色分けするなどして表示する。なお、１つの時刻ステップにおいて上記の一連の処理を終了すると、上記したようにステップＳ１１０に戻り、ステップＳ１１３にて各粒子の位置する要素番号を更新し各要素領域内の粒子を特定するが、この際、各粒子の（元の）要素番号を用いて、各要素領域について、その近傍の要素領域の要素番号を有する粒子のみ当該処理を行うことができ、これにより計算量の低減につながる。
【００８８】
本実施形態では、ＭＰＭとして、上記したステップＳ１２０〜Ｓ１４０の処理にかえて、ＵＳＡＶＧと呼ばれる、１つの時刻ステップにおける時間増分を半分にし、粒子のひずみ、応力、密度を２回に分けて算出する方法をとることもできる。以下これを図９を用いて簡単に説明する。
【００８９】
図９に示すステップＳ２００における初期条件の設定、ステップＳ２１０における粒子位置に係る諸計算の処理は、上記したステップＳ１００、Ｓ１１０の処理と同様である。ＵＳＡＶＧでは、まず、ステップＳ２２０において、節点の速度を算出する。
【００９０】
ステップＳ２２０では、まず、節点の質量ｍｉ（ｋ）を前述の式（２０）により算出した後、次式により、節点の速度を算出する。
【数３７】
0005221603.tif 000038
これは、前述のステップＳ１３０で節点の速度（第２の速度）を求めたのと同様の処理である。
【００９１】
次に、ステップ２３０で、粒子のひずみ増分Δεｉｊ  ｐ（ｋ）を、次式から算出する。
【数３８】
0005221603.tif 000039
これは、前述のステップＳ１３０で粒子のひずみ増分を求めたのと同様の処理である。
【００９２】
次に、ステップＳ２４０において、粒子のひずみ増分より、時刻ステップｋ＋１／２における粒子のひずみ、応力、密度をそれぞれ算出する。
【００９３】
まず、粒子の応力増分Δσｉｊ  ｐ（ｋ）を次式で算出する。
【数３９】
0005221603.tif 000040
また、粒子のひずみを次式で算出し、更新する。
【数４０】
0005221603.tif 000041
そして、粒子の応力を次式で算出する。
【数４１】
0005221603.tif 000042
ここで、
【数４２】
0005221603.tif 000043
である。
さらに、粒子の密度が次式で算出される。
【数４３】
0005221603.tif 000044
以上の処理は、前述のステップＳ１４０における処理と同様の処理である。
【００９４】
ステップＳ２５０では、要素領域内において、節点力を算出し、節点の質量と節点力より節点の加速度を算出し、節点の加速度より節点の第１の速度を算出する。
【００９５】
まず、内力としての節点力ｆｉｉｎｔ（ｋ＋１）を粒子の質量、応力、密度を基に勾配関数Ｇｎｉｐを用いて次式
【数４４】
0005221603.tif 000045
で算出し、外力としての節点力ｆｉｅｘｔ（ｋ＋１）をＭＰＭにおける体積力ｂｉ（ｋ＋１）と接触力τｉ（ｋ＋１）を用いて次式
【数４５】
0005221603.tif 000046
で算出する。ここで、
【数４６】
0005221603.tif 000047
【数４７】
0005221603.tif 000048
である。次に、節点の質量と節点力を基に次式に示す節点の運動方程式
【数４８】
0005221603.tif 000049
で節点の加速度ａｉ（ｋ＋１）を求める。最後に、節点の加速度を基に次ステップｋ＋１での節点の第１の速度ｖｉ（ｋ＋１）を、
【数４９】
0005221603.tif 000050
により算出する。以上の処理は、前述したステップＳ１２０の処理に対応する。
【００９６】
次に、ステップＳ２６０において、粒子の位置、速度、ひずみ増分を算出する。
【００９７】
ステップＳ２６０では、まず、粒子の位置Ｘｐ（ｋ＋１）を、節点の第１の速度を基に内挿関数を用いて次式で算出し、更新する。
【数５０】
0005221603.tif 000051
また、節点の加速度を基に内挿関数を用いて粒子の速度ｖｐ（ｋ＋１）を次式で算出し、更新する。
【数５１】
0005221603.tif 000052
さらに、Ｓｕｌｓｋｙにより提案された速度ｖi（ｋ＋１）の平滑化を行い、節点の速度（第２の速度）ｖi（ｋ＋１）を次式で再び求める。即ち、節点の速度（第２の速度）を、粒子の質量、速度、節点の質量を基に内挿関数を用いて次式で算出、更新する。
【数５２】
0005221603.tif 000053
最後に、粒子のひずみ増分Δεｉｊ  ｐ（ｋ＋１）を、節点の速度を基に、勾配関数を用いて次式から算出する。
【数５３】
0005221603.tif 000054
以上の処理は、前述したステップＳ１３０の処理に対応する。
【００９８】
続いて、ステップＳ２７０において、粒子のひずみ増分を基に粒子のひずみ、応力、密度を再び算出し、更新する。
【００９９】
ステップＳ２７０では、まず、物体の剛性マトリックスＴｉｊｋｌを用いて粒子のひずみ増分を基に粒子の応力増分Δσｉｊ  ｐ（ｋ＋１）を次式で算出する。
【数５４】
0005221603.tif 000055
また、粒子のひずみεｉｊ  ｐ（ｋ＋１）を、粒子のひずみ増分を基に次式で算出し、更新する。
【数５５】
0005221603.tif 000056
また、粒子の応力σｉｊ  ｐ（ｋ＋１）を粒子の応力増分を基に算出し更新する。粒子ｐの応力σｉｊ  ｐ（ｋ＋１）は先程と同様、次式で算出され、更新される。
【数５６】
0005221603.tif 000057
ここで
【数５７】
0005221603.tif 000058
である。また、粒子の密度ρｐ（ｋ＋１）が、粒子のひずみ増分を用いて次式で算出され、更新される。
【数５８】
0005221603.tif 000059
以上の処理は、前述したステップＳ１４０の処理に対応する。
【０１００】
以降、全時刻ステップでの計算を終えるまで（ステップＳ２８０のＮｏ）、時刻ステップ数ｋの値をインクリメントし（ステップＳ２９０）、ステップＳ２１０に戻って次の時刻暦計算を行う。この処理の過程で、各粒子のもつラグランジュ変数である粒子情報が更新される。このようにして予め定めた時刻ステップ数だけ時刻暦計算を繰り返して時刻暦計算を終了すると（ステップＳ２８０のＹｅｓ）、地盤解析を終了し、表示部に位置、ひずみ、応力等、各粒子についての諸値を色分けするなどして表示する。
【０１０１】
次に、本実施形態のＭＰＭによる解析手法を地盤解析に用いる際の、水圧の考慮について説明する。ＭＰＭを地盤解析に用いるには、これを水圧を考慮したものとすることにより、例えば液状化現象や浸透流など、地下水の挙動を解析に取り込むことが可能になる。
【０１０２】
以下に示すｕ−Ｐ形式で表したダルシー則方程式においては、使用する流速の項において、加速度項はとりのぞいた方が計算が安定しているとされている。
【数５９】
0005221603.tif 000060
ここで、ｋは透水係数、ｇは重力加速度、ｎは間隙率、ｋｆは水の体積弾性係数、ρｆは水の密度を示す。非排水を仮定すると透水に関する加速度がなくなり、
【数６０】
0005221603.tif 000061
となる。さらに水は弾性体であることを考慮すると、式（６０）を時間に関して積分し、
【数６１】
0005221603.tif 000062
と表せる。ただしｔ＝０で、εν＝０、Ｐ＝０とする。ενは体積ひずみであり、粒子のひずみより求めることができる。
【０１０３】
ここで、ＭＰＭにおいて、粒子のひずみまたは粒子のひずみ増分より水圧を求めることができる。水圧は、節点力を求める際に内力として含めることができる。例えば上記のＵＳＡＶＧにおいては、上記のステップＳ２４０において時刻ステップｋ＋１／２における粒子のひずみを求めた後、次式
【数６２】
0005221603.tif 000063
により時刻ステップｋ＋１／２における水圧Ｐ（ｋ＋１／２）を算出するか、あるいは、上記のステップＳ２３０において算出した粒子のひずみ増分を用いて、次式
【数６３】
0005221603.tif 000064
で水圧Ｐ（ｋ＋１／２）を算出し更新する。これらの水圧はＲＡＭ等に記憶する。ステップＳ２５０では、発生した水圧を次式で節点の内力に含めるとよい。
【数６４】
0005221603.tif 000065
即ち、上記説明した通常のＵＳＡＶＧにＰ（ｋ＋１／２）の項を加える。ここで、δｎｍはクロネッカーのデルタを表す。
【０１０４】
水の弾性係数ｋｆが十分に大きいと土の体積ｔｒεはほぼゼロとなり、体積変化の発生しない計算になる。ＦＥＭでは要素の中心点で平均主応力と水圧を評価することで、非線形解析に使用する応力と水圧の評価点を一致させているが、粒子法においては、応力の評価点は粒子上であるが、水圧の評価点を粒子の位置によらず要素領域の中心位置にすることで、すべての粒子において式（６１）で示される状態を維持することが可能となる。具体的には、ひずみ計算と内力計算において体積ひずみεＶのみ要素領域の中心で計算する選択積分となる。
【０１０５】
次に、粒子法による解析手法を地盤解析に用いる際の、動的解析への拡張について説明する。例えば、本実施形態の地盤変形解析方法では、下方向の粘性境界の設定を行うことができる。下方向の粘性境界は、図１０（ａ）に示すようなダッシュポットを節点に付ける、即ち次式を節点の運動方程式として用いることで表現される。ダッシュポットにおいては、地震加速度ａ”、重力加速度ｇ、相対加速度ｕ”、バネ定数ｋとすると運動方程式は下記となる。
【数６５】
0005221603.tif 000066
本実施形態においては、必要な節点については節点の運動方程式として上式（６５）を用いることができる。前述のフローにおいては、上式の相対加速度を基に節点の第１の速度を求める流れとなる。
【０１０６】
また、粘性境界の節点については、図１０（ｂ）に示すように、
【数６６】
0005221603.tif 000067
により節点の内力が求められる。図１０（ｂ）に関して、Ｒｉは反力であり、節点１での反力がＲ１である。上記した初期応力時にはfｉｉｎｔとＲｉはつり合っているので作用力がゼロとなり節点は移動しない。
【０１０７】
従って、粘性境界による作用力は、図１０（ｃ）に示すように、動的な応力変化分（σ−σｉ）より求まるので、粘性境界の必要な節点では、初期応力と重力分をとりのぞいた、応力変化分を作用力として算出される粘性力が運動方程式において節点の内力に加わる。
【０１０８】
このほか、地盤の固定境界については、所定の要素領域の底面の節点において、前述のフローにおける速度ｖｉ、加速度aｉを０とすることで表現できる。
【０１０９】
さらに、側面境界条件として、要素領域を指定して、当該要素領域内で粒子の移動できない方向を指定することもできる。各要素領域内の粒子の移動について、例えばＺ方向（上下方向）への移動不可を設定し、Ｚ方向の粒子位置を固定すると粒子は上下に移動しなくなる。粒子の移動を強制的にゼロとしているので、移動している場合と同じひずみが発生し応力を発生するが、粒子がとどまることにより接触力に伴う周辺粒子への反力を発生する。
【０１１０】
また、各節点にレイリー減衰を作用させることができる。レイリー減衰の場合は、質量行列と剛性行列それぞれに係数α、βを掛けた値の減衰行列をもつダッシュポットを模した運動方程式を節点の加速度の算出に用いる。例えば上記のＵＳＡＶＧの場合では、ｋ＋１／２ステップで次式のいずれかを節点の内力として運動方程式を解くことができる。
【数６７】
0005221603.tif 000068
ここで、Ｄは剛性行列を表す。
【０１１１】
次に、本地盤解析装置により、３次元ＭＰＭを用いた解析例を示す。図１１は、本地盤解析装置により乾燥砂の流動の３次元解析を行った例である。図１１（ａ）は流動前の状態、図１１（ｂ）は流動発生５秒経過後の状態、図１１（ｃ）は流動後の状態を示す。スリットから流れ出る乾燥砂の流動の様子を３次元上で再現できていることがわかる。
【０１１２】
一方、図１２は、本地盤変形解析装置により、斜面の大変形に関する３次元動的解析を行った例で、斜面が地震動を受ける際の大変形挙動について、地盤の変形とともに最大せん断ひずみを示したものである。図１２（ａ）〜（ｃ）は斜面勾配が緩やかな場合であり、図１２（ａ）は崩壊前の状態、図１２（ｂ）は崩壊開始０．５秒経過した状態、図１２（ｃ）は崩壊後の状態を示す。また、図１３（ａ）〜（ｃ）は斜面勾配が急な場合であり、図１３（ａ）は崩壊前の状態、図１３（ｂ）は崩壊開始０．５秒経過した状態、図１３（ｃ）は崩壊後の状態を示す。斜面勾配の大きさによって、斜面の大変形挙動が異なることが示されている。斜面勾配が緩やかな場合では，斜面土塊が緩速に変形し，斜面勾配が急な場合では急速に挙動する様子が再現できていることがわかる。
【０１１３】
以上、添付図面を参照しながら、本発明に係る地盤変形解析装置等の好適な実施形態について説明したが、本発明はかかる例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
【符号の説明】
【０１１４】
１………地盤変形解析装置
１００………粒子
１０５………要素
１１０………粒子情報
１２０………節点
【図１】
0005221603.tif 000069
【図２】
0005221603.tif 000070
【図３】
0005221603.tif 000071
【図４】
0005221603.tif 000072
【図５】
0005221603.tif 000073
【図６】
0005221603.tif 000074
【図７】
0005221603.tif 000075
【図８】
0005221603.tif 000076
【図９】
0005221603.tif 000077
【図１０】
0005221603.tif 000078
【図１４】
0005221603.tif 000079
【図１１】
0005221603.tif 000080
【図１２】
0005221603.tif 000081
【図１３】
0005221603.tif 000082