(57)【要約】  地盤の各土要素に対する粘土・中間土・砂などの土の設定や地盤立体構造モデルの設定、計算条件の設定に基づいて土骨格の時間体積変化率及び接線剛性並びに土の透水性、透気性、水分特性及び質量に関する各種マトリクスを計算し（Ｓ１４０〜Ｓ２２０）、これらを用いた連立方程式を作成し、変形または応力速度等に関する境界条件、間隙水の水理条件及び間隙空気の空気境界条件を与えて、未知の「加加速度（躍度）場」と間隙水圧場及び間隙空気圧場の求解を行なう（Ｓ２３０）。これにより、あらゆる状態・種類の土からなる地盤についての解析を精度よく行なうことができる。
2014002977.tif 000002
【特許請求の範囲】
【請求項１】
土骨格の状態に関するデータと外力に関するデータとを入力するデータ入力手段と、
前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を構築し、該連立速度型運動方程式に変位、変位速度又は加速度に関する幾何的境界条件と応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう解析手段と、
該解析結果を出力する結果出力手段と、
を備える空気と水と土骨格の連成計算装置。
【請求項２】
請求項１記載の空気と水と土骨格の連成計算装置であって、
前記解析手段は、前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水・間隙空気の時間的体積変化率に関する第５項目と、間隙水圧の荷重変化率寄与に関する第６項目と、間隙空気圧の荷重変化率寄与に関する第７項目と、間隙水の流量変換に関する第８項目と、間隙空気の流量変換に関する第９項目と、質量に関する第１０項目と、土骨格の接線剛性に関する第１１項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する手段である、
空気と水と土骨格の連成計算装置。
【請求項３】
請求項１記載の空気と水と土骨格の連成計算装置であって、
前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、により構築される方程式である、
空気と水と土骨格の連成計算装置。
【請求項４】
請求項３記載の空気と水と土骨格の連成計算装置であって、
前記空気と水と土骨格からなる三相混合体の運動方程式は、式(Ａ１)により表わされ、
前記水と土骨格間での質量保存式は、式(Ｂ１)により表わされ、
前記空気と土骨格間での質量保存式は、式(Ｃ１)により表わされる、
空気と水と土骨格の連成計算装置。
【数１】
2014002977.tif 000129
【請求項５】
請求項１記載の空気と水と土骨格の連成計算装置であって、
前記連立速度型運動方程式は、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧と飽和度に関するものであり、
前記土骨格の変形解析は、前記幾何的境界条件と前記力学的境界条件と前記水理境界条件と前記空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場と飽和度場の時刻歴応答を求めることにより行なうものである、
空気と水と土骨格の連成計算装置。
【請求項６】
請求項５記載の空気と水と土骨格の連成計算装置であって、
前記解析手段は、前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水の時間的体積変化率に関する第５項目と、間隙空気の時間的体積変化率に関する第６項目と、飽和度の時間的変化率に関する第７項目、間隙水圧の荷重変化率寄与に関する第８項目と、間隙空気圧の荷重変化率寄与に関する第９項目と、飽和度の荷重変化率寄与に関する第１０項目と、間隙水の流量変換に関する第１１項目と、間隙空気の流量変換に関する第１２項目と、質量に関する第１３項目と、土骨格の接線剛性に関する第１４項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する手段である、
空気と水と土骨格の連成計算装置。
【請求項７】
請求項５記載の空気と水と土骨格の連成計算装置であって、
前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、土の水分特性に関する式と、により構築される方程式である、
空気と水と土骨格の連成計算装置。
【請求項８】
請求項７記載の空気と水と土骨格の連成計算装置であって、
前記空気と水と土骨格からなる三相混合体の運動方程式は、式(Ａ２)により表わされ、
前記水と土骨格間での質量保存式は、式(Ｂ２)により表わされ、
前記空気と土骨格間での質量保存式は、式(Ｃ２)により表わされ、
前記土の水分特性に関する式は、式(Ｄ２）により表わされる、
空気と水と土骨格の連成計算装置。
【数２】
2014002977.tif 000130
【請求項９】
前記解析手段は、有限要素法、差分法等の数値解析法を用いて前記連立速度型運動方程式を構築する手段である請求項１ないし８のうちのいずれか１つの請求項に記載の空気と水と土骨格の連成計算装置。
【請求項１０】
前記土骨格の状態に関するデータは、土骨格を構成する複数の要素における各要素の性状に関するデータおよび隣接する要素との関係に関するデータとが含まれるデータである請求項１ないし９のうちのいずれか１つの請求項に記載の空気と水と土骨格の連成計算装置。
【請求項１１】
前記各要素の性状に関するデータは、土の性状に関するデータが含まれるデータである請求項１０記載の空気と水と土骨格の連成計算装置。
【請求項１２】
前記外力に関するデータは、空気と水と土骨格の連成系に与える荷重及び変位に関するデータが含まれるデータである請求項１ないし１１のうちのいずれか１つの請求項に記載の空気と水と土骨格の連成計算装置。
【請求項１３】
前記外力に関するデータは、空気と水と土骨格の連成系に与える振動に関するデータが含まれるデータである請求項１ないし１２のうちのいずれか１つの請求項に記載の空気と水と土骨格の連成計算装置。
【請求項１４】
前記解析手段は、前記外力に関するデータを複数回変更して解析する手段である請求項１ないし１３のうちのいずれか１つの請求項に記載の空気と水と土骨格の連成計算装置。
【請求項１５】
前記解析手段は、前記外力に関するデータを変更するときには変更前の解析結果を初期状態の土骨格の状態に関するデータとして用いて解析する手段である請求項１４記載の空気と水と土骨格の連成計算装置。
【請求項１６】
土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を構築し、該連立方程式に変位、変位速度又は加速度に関する幾何的境界条件と応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場、間隙水圧場及び間隙空気圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう、
ことを特徴とする空気と水と土骨格の連成計算方法。
【請求項１７】
請求項１６記載の空気と水と土骨格の連成計算方法であって、
前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水・間隙空気の時間的体積変化率に関する第５項目と、間隙水圧の荷重変化率寄与に関する第６項目と、間隙空気圧の荷重変化率寄与に関する第７項目と、間隙水の流量変換に関する第８項目と、間隙空気の流量変換に関する第９項目と、質量に関する第１０項目と、土骨格の接線剛性に関する第１１項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する、
ことを特徴とする空気と水と土骨格の連成計算方法。
【請求項１８】
請求項１６記載の空気と水と土骨格の連成計算方法であって、
前記連立速度型運動方程式は、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧と飽和度に関するものであり、
前記土骨格の変形解析は、前記幾何的境界条件と前記力学的境界条件と前記水理境界条件と前記空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場と飽和度場の時刻歴応答を求めることにより行なうものである、
空気と水と土骨格の連成計算方法。
【請求項１９】
請求項１８記載の空気と水と土骨格の連成計算方法であって、
前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水の時間的体積変化率に関する第５項目と、間隙空気の時間的体積変化率に関する第６項目と、飽和度の時間的変化率に関する第７項目、間隙水圧の荷重変化率寄与に関する第８項目と、間隙空気圧の荷重変化率寄与に関する第９項目と、飽和度の荷重変化率寄与に関する第１０項目と、間隙水の流量変換に関する第１１項目と、間隙空気の流量変換に関する第１２項目と、質量に関する第１３項目と、土骨格の接線剛性に関する第１４項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する、
ことを特徴とする空気と水と土骨格の連成計算方法。
【請求項２０】
土骨格の変形解析を行なう空気と水と土骨格の連成計算プログラムであって、
土骨格の状態に関するデータと外力に関するデータとを入力するデータ入力モジュールと、
前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を構築し、該連立速度型運動方程式に変位、変位速度又は加速度に関する幾何的境界条件と応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう解析モジュールと、
を備える空気と水と土骨格の連成計算プログラム。
【請求項２１】
請求項２０記載の空気と水と土骨格の連成計算プログラムであって、
前記解析モジュールは、前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水・間隙空気の時間的体積変化率に関する第５項目と、間隙水圧の荷重変化率寄与に関する第６項目と、間隙空気圧の荷重変化率寄与に関する第７項目と、間隙水の流量変換に関する第８項目と、間隙空気の流量変換に関する第９項目と、質量に関する第１０項目と、土骨格の接線剛性に関する第１１項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築するモジュールである、
空気と水と土骨格の連成計算プログラム。
【請求項２２】
請求項２０記載の空気と水と土骨格の連成計算プログラムであって、
前記連立速度型運動方程式は、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧と飽和度に関するものであり、
前記土骨格の変形解析は、前記幾何的境界条件と前記力学的境界条件と前記水理境界条件と前記空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場と飽和度場の時刻歴応答を求めることにより行なうものである、
空気と水と土骨格の連成計算プログラム。
【請求項２３】
請求項２２記載の空気と水と土骨格の連成計算プログラムであって、
前記解析モジュールは、前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水の時間的体積変化率に関する第５項目と、間隙空気の時間的体積変化率に関する第６項目と、飽和度の時間的変化率に関する第７項目、間隙水圧の荷重変化率寄与に関する第８項目と、間隙空気圧の荷重変化率寄与に関する第９項目と、飽和度の荷重変化率寄与に関する第１０項目と、間隙水の流量変換に関する第１１項目と、間隙空気の流量変換に関する第１２項目と、質量に関する第１３項目と、土骨格の接線剛性に関する第１４項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築するモジュールである、
空気と水と土骨格の連成計算プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラムに関し、詳しくは、土骨格の変形解析を行なう空気と水と土骨格の連成計算を行なう装置および土骨格の変形解析を行なう空気と水と土骨格の連成計算を行なう方法並びに土骨格の変形解析を行なう空気と水と土骨格の連成計算を行なうためのプログラムに関する。
【背景技術】
【０００２】
従来、この種の土骨格の連成計算装置としては、土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧に関する連立速度型運動方程式を構築し、変位等の幾何的境界条件と応力等の力学的境界条件と間隙水圧等の水理境界条件とを与えて随時積分し、速度場および間隙水圧場の時刻歴応答を求めることにより土骨格の変形解析を行なうものが提案されている（例えば、特許文献１参照）。この土骨格の連成計算装置では、上述した解析により、加速度の影響が大きいときには動的解析による土骨格の変形解析を行ない、加速度の影響が小さいときには静的解析による土骨格の変形解析を行なうことができ、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことを可能としている。なお、「加加速度項」は、説明変数である速度ベクトルの時間に関する２階微分の項を意味しており、「加加速度」は「躍度」とも称される。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】国際公開ＷＯ  ２００７／０４６１７８  Ａ１
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上述の土骨格の連成計算装置では、土粒子と水のみで構成されている飽和土を対象としているため、土粒子の間隙が水により飽和されている砂から粘土に至るまでのあらゆる中間土からなる地盤について、静的解析か動的解析かに拘わらず土骨格の変形解析を行なうことが可能であるものの、土粒子の間隙が水で飽和されていない土、即ち、土粒子の間隙に水の他に空気が存在する不飽和土に対しては十分な解析を行なうことができない。盛土や斜面は不飽和土であることが多く、上述の土骨格の連成計算装置では、降雨浸透に伴う盛土や斜面の変形・破壊解析や、豪雨と地震の複合時の盛土や斜面の変形・破壊解析を行なうことができない。【０００５】
本発明の空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラムは、土粒子の間隙に水の他に空気が存在する不飽和土に対して静的解析か動的解析かに拘わらずに土骨格の変形解析を行なうことを主目的とする。即ち、本発明の空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラムは、土が飽和土であるか不飽和土であるかに拘わらずに静的または動的に土骨格の変形解析を行なうことを主目的としており、例えば、不飽和土から飽和土までを連続的に扱う解析に適用したり、降雨浸透に伴う盛土や斜面の変形・破壊解析や豪雨と地震の複合時の盛土や斜面の変形・破壊解析に適用することを可能とする。
【課題を解決するための手段】
【０００６】
本発明の空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラムは、上述の主目的を達成するために以下の手段を採った。【０００７】
本発明の空気と水と土骨格の連成計算装置は、
土骨格の状態に関するデータと外力に関するデータとを入力するデータ入力手段と、
前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を構築し、該連立速度型運動方程式に変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう解析手段と、
該解析結果を出力する結果出力手段と、
を備えることを要旨とする。【０００８】
この本発明の空気と水と土骨格の連成計算装置では、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を変位、変位速度又は加速度に関する幾何的境界条件と応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場の時刻歴応答を求めるから、土骨格の間隙に水の他に空気が存在する不飽和土に対して、加速度の影響が大きいときには動的解析による土骨格の変形解析を行なうことができ、加速度の影響が小さいときには静的解析による土骨格の変形解析を行なうことができる。したがって、土骨格の間隙に水が飽和している飽和土だけでなく、土骨格の間隙に水が飽和していない土、即ち、土骨格の間隙に水の他に空気が存在する不飽和土に対して、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。ここで、「加加速度項」は、説明変数である速度ベクトルの時間に関する２階微分の項を意味しており、「加加速度」を「躍度」とも称する。【０００９】
こうした本発明の空気と水と土骨格の連成計算装置において、前記解析手段は、前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水・間隙空気の時間的体積変化率に関する第５項目と、間隙水圧の荷重変化率寄与に関する第６項目と、間隙空気圧の荷重変化率寄与に関する第７項目と、間隙水の流量変換に関する第８項目と、間隙空気の流量変換に関する第９項目と、質量に関する第１０項目と、土骨格の接線剛性に関する第１１項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する手段である、ものとすることもできる。こうすれば、不飽和土に対して、より適正に土骨格の変形解析を行なうことができる。【００１０】
また、本発明の空気と水と土骨格の連成計算装置において、前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、により構築される方程式である、ものとすることもできる。この場合、前記空気と水と土骨格からなる三相混合体の運動方程式は式(A1)により表わされ、前記水と土骨格間での質量保存式は式(B1)により表わされ、前記空気と土骨格間での質量保存式は式(C1)により表わされる、ものとすることもできる。【００１１】
【数１】
2014002977.tif 000003
【００１２】
ここで、式中の上付きの「＼」（以下「＼（上付き）」と記載する。）は土骨格から見た物質時間微分を意味し、例えば「＼（上付き）」が２つ付いた「v」は土骨格の加加速度ベクトルを表している。なお、以下では、土骨格から見た物質時間微分を意味する「＼（上付き）」は上付きの「・」（以下「・（上付き）」と表記し、また、土骨格を意味するs（下付き）は省略して記述する。例えば「＼（上付き）」が２つ付いた「v」は「・（上付き）」が２つ付いた「v」となる。さらに簡単のため、間隙水は非圧縮性流体（Kwが無限大）として扱う。その結果、式(A1)、式(B1)および式(C1)の基礎式は、以下の式(A1)’、式(B1)’および式(C1)’になる。【００１３】
【数２】
2014002977.tif 000004
【００１４】
次に、上述の式(A1)’を弱形式化すると、すなわち下記の式（ａ１１）、式（ａ１２）を満たす仮想変位速度を用いて弱形式化を施すと、式（ａ１３）のようになる。ここでは、有効応力式として平均化骨格応力の式（Bishopの有効応力式においてχ=swとした式）、客観性のある有効応力速度としてGreen Naghdiの有効応力速度を用いると、有効応力速度テンソルは式（ａ１４）のようになる。以上をまとめると、空気と水と土骨格からなる三相混合体の運動方程式の弱形式は式（ａ１５）のようになる。【００１５】
【数３】
2014002977.tif 000005
【００１６】
【数４】
2014002977.tif 000006
【００１７】
【数５】
2014002977.tif 000007
【００１８】
ここでは、空間的離散化手法として、有限要素法を用い、さらに間隙水圧および間隙空気圧を要素中心に割り当てる田村流・Christian流に基づくと、この弱形式は次式(A1)”になる。各マトリクスの具体形は上述の実施の形態にて示す通りであり、Mマトリクスは式(48)、Kマトリクスは式(50)、LWTマトリクスは式(75)、LATマトリクスは式(77)に相当する。また右辺第１項は式(29)に相当する。【００１９】
【数６】
2014002977.tif 000008
【００２０】
次に、水と土骨格の連成式の式(B1)’と空気と土骨格の連成式の式(C1)’を空間的に離散化する。空気と水と土骨格からなる三相混合体の運動方程式の場合と同様に、田村流・Christian流に基づいて離散化を施すと、式(B1)”と式(C1)”を得る。各マトリクスの具体形は後述する実施の形態にて示す通りであり、LW1マトリクスは式(81)、LW2マトリクスは式(83)、LA1マトリクスは式(82)、LA2マトリクスは式(84)、LBWマトリクスは式(20)、LBAマトリクスは式(25)、SWAマトリクスは式(79)、SACマトリクスは式(80)、HWマトリクスは式(14)、HAマトリクスは式(15)である。また式(B1)”の右辺第１項は式(44)、式(C1)”の右辺第１項は式(46)に相当する。【００２１】
【数７】
2014002977.tif 000009
【００２２】
以上、式(A1)”、式(B1)”、式(C1)”を常微分方程式の形に整理すると、下記の式(E1)を得る。ここで、LW1、LW2、LA1およびLA2マトリクスが土骨格の時間的体積変化率に関する第１項目、HWマトリクスが土の透水性に関する第２項目、HAマトリクスが土の透気性に関する第３項目、SWAマトリクスが土の水分特性に関する第４項目、SACマトリクスが（間隙水・）間隙空気の時間的体積変化率に関する第５項目、LWTマトリクスが間隙水圧の荷重変化率寄与に関する第６項目、LATマトリクスが間隙空気圧の荷重変化率寄与に関する第７項目、LBWマトリクスが間隙水の流量変換に関する第８項目、LBAマトリクスが間隙空気の流量変換に関する第９項目、Mマトリクスが質量に関する第１０項目、Kマトリクスが土骨格の接線剛性に関する第１１項目に当たる。【００２３】
【数８】
2014002977.tif 000010
【００２４】
さらに、連立一次方程式の形に帰着させるために、Wilsonのθ法やNewmarkのβ法に準拠した時間積分方法を用いることができる。以下に、Wilsonのθ法に準拠した場合を示す。【００２５】
【数９】
2014002977.tif 000011
【００２６】
また、間隙水圧および間隙空気圧に関しては、台形公式を用いて、時間積分方法を用いることができる。以下に、台形公式を用いた時間積分方法を示す。【００２７】
【数１０】
2014002977.tif 000012
【００２８】
以上より、常微分方程式として与えられる上述の式(E1)は、下記の式(F1)の連立一次方程式に帰着させることができる。【００２９】
【数１１】
2014002977.tif 000013
【００３０】
本発明の空気と水と土骨格の連成計算装置において、前記連立速度型運動方程式は、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧と飽和度に関するものであり、前記土骨格の変形解析は、前記幾何的境界条件と前記力学的境界条件と前記水理境界条件と前記空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場と飽和度場の時刻歴応答を求めることにより行なうものである、ものとすることもできる。このように間隙水の圧縮性を考慮すると共に飽和度についてより詳細に考慮することにより、土骨格の変形解析をより詳細に行なうことができる。【００３１】
こうした間隙水の圧縮性を考慮した態様の本発明の空気と水と土骨格の連成計算装置において、前記解析手段は、前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水の時間的体積変化率に関する第５項目と、間隙空気の時間的体積変化率に関する第６項目と、飽和度の時間的変化率に関する第７項目、間隙水圧の荷重変化率寄与に関する第８項目と、間隙空気圧の荷重変化率寄与に関する第９項目と、飽和度の荷重変化率寄与に関する第１０項目と、間隙水の流量変換に関する第１１項目と、間隙空気の流量変換に関する第１２項目と、質量に関する第１３項目と、土骨格の接線剛性に関する第１４項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する手段である、ものとすることもできる。こうすれば、不飽和土に対して、より適正に土骨格の変形解析を行なうことができる。【００３２】
また、間隙水の圧縮性を考慮した態様の本発明の空気と水と土骨格の連成計算装置において、前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、土の水分特性に関する式と、により構築される方程式である、ものとすることもできる。この場合、前記空気と水と土骨格からなる三相混合体の運動方程式は、式(A2)により表わされ、前記水と土骨格間での質量保存式は、式(B2)により表わされ、前記空気と土骨格間での質量保存式は、式(C2)により表わされ、前記土の水分特性に関する式は、式(D2）により表わされる、ものとすることもできる。【００３３】
【数１２】
2014002977.tif 000014
【００３４】
ここで、式中の上付きの「＼」は、上述したとおり、土骨格から見た物質時間微分を意味し、例えば「＼（上付き）」が２つ付いた「v」は土骨格の加加速度ベクトルを表している。また、式(D2)に示すように、土の水分特性に関する式は用いるモデルに依存するため、以下では簡単のために、土の水分特性に関する式は、サクションとサクションの土骨格から見た物質時間微分と飽和度と飽和度の土骨格から見た物質時間微分の関係で表す形を用いる。サクションは間隙空気圧と間隙水圧の差を意味する。その結果、式(A2)、式(B2)、式(C2)および式(D2)の基礎式は、以下の式(A2)’、式(B2)’、式(C2)’および式(D2)’になる。【００３５】
【数１３】
2014002977.tif 000015
【００３６】
次に、上述の式(A2)’を弱形式化すると、すなわち下記の式（ａ２１）、式（ａ２２）を満たす仮想変位速度を用いて弱形式化を施すと、式（ａ２３）のようになる。ここでは、有効応力式として平均化骨格応力の式（Bishopの有効応力式においてχ=swとした式）、客観性のある有効応力速度としてGreen Naghdiの有効応力速度を用いると、有効応力速度テンソルは式（ａ２４）のようになる。以上をまとめると、空気と水と土骨格からなる三相混合体の運動方程式の弱形式は式（ａ２５）のようになる。【００３７】
【数１４】
2014002977.tif 000016
【００３８】
【数１５】
2014002977.tif 000017
【００３９】
【数１６】
2014002977.tif 000018
【００４０】
ここでは、空間的離散化手法として、有限要素法を用い、さらに田村流・Christian流に倣って、間隙水圧、間隙空気圧および飽和度を要素中心に割り当てる方法を用いると、この弱形式は次式(A2)”になる。各マトリクスの具体形は上述の実施の形態にて示す通りであり、Mマトリクスは式(249)、Kマトリクスは式(251)、LWTマトリクスは式(276)、LATマトリクスは式(278)、LSTマトリクスは式(280)に相当する。また右辺第１項は式(229)に相当する。【００４１】
【数１７】
2014002977.tif 000019
【００４２】
次に、水と土骨格の連成式の式(B2)’、空気と土骨格の連成式の式(C2)’および土の水分特性に関する式の式(D2)’を空間的に離散化する。空気と水と土骨格からなる三相混合体の運動方程式の場合と同様に、田村流・Christian流に倣って離散化を施すと、式(B2)”、式(C2)”および式(D2)”を得る。各マトリクスの具体形は後述する実施の形態にて示す通りであり、LW1マトリクスは式(287)、LW2マトリクスは式(289)、LA1マトリクスは式(288)、LA2マトリクスは式(290)、LBWマトリクスは式(220)、LBAマトリクスは式(225)、SWCマトリクスは式(282)、SACマトリクスは式(283)、SSWマトリクスは式(284)、HWマトリクスは式(214)、HAマトリクスは式(215)、SCマトリクスは式(285)、SWマトリクスは式(286)である。また式(B2)”の右辺第１項は式(244)、式(C2)”の右辺第１項は式(247)に相当する。【００４３】
【数１８】
2014002977.tif 000020
【００４４】
以上、式(A2)”、式(B2)”、式(C2)”、式(D2)”を常微分方程式の形に整理すると、下記の式(E2)を得る。ここで、LW1、LW2、LA1およびLA2マトリクスが土骨格の時間的体積変化率に関する第１項目、HWマトリクスが土の透水性に関する第２項目、HAマトリクスが土の透気性に関する第３項目、SCおよびSWマトリクスが土の水分特性に関する第４項目、SWCマトリクスが間隙水の時間的体積変化率に関する第５項目、SACマトリクスが間隙空気の時間的体積変化率に関する第６項目、SSWマトリクスが飽和度の時間的変化率に関する第７項目、LWTマトリクスが間隙水圧の荷重変化率寄与に関する第８項目、LATマトリクスが間隙空気圧の荷重変化率寄与に関する第９項目、LSTマトリクスが飽和度の荷重変化率寄与に関する第１０項目、LBWマトリクスが間隙水の流量変換に関する第１１項目、LBAマトリクスが間隙空気の流量変換に関する第１２項目、Mマトリクスが質量に関する第１３項目、Kマトリクスが土骨格の接線剛性に関する第１４項目に当たる。【００４５】
【数１９】
2014002977.tif 000021
【００４６】
さらに、連立一次方程式の形に帰着させるために、Wilsonのθ法やNewmarkのβ法に準拠した時間積分方法を用いることができる。以下に、Wilsonのθ法に準拠した場合を示す。【００４７】
【数２０】
2014002977.tif 000022
【００４８】
また、間隙水圧、間隙空気圧および飽和度に関しては、例えば台形公式を用いて、時間積分方法を用いることができる。以下に、台形公式を用いた時間積分方法を示す。【００４９】
【数２１】
2014002977.tif 000023
【００５０】
以上より、常微分方程式として与えられる上述の式(E2)は、下記の式(F2)の連立一次方程式に帰着させることができる。【００５１】
【数２２】
2014002977.tif 000024
【００５２】
また、本発明の空気と水と土骨格の連成計算装置において、前記解析手段は、有限要素法、差分法等の数値解析法を用いて前記連立速度型運動方程式を構築する手段であるものとすることもできる。こうすれば、計算を容易なものとすることができる。【００５３】
さらに、本発明の空気と水と土骨格の連成計算装置において、前記土骨格の状態に関するデータは、土骨格を構成する複数の要素における各要素の性状に関するデータおよび隣接する要素との関係に関するデータとが含まれるデータであるものとすることもできる。この場合、前記各要素の性状に関するデータは、土の性状に関するデータが含まれるデータであるものとすることもできる。こうすれば、土骨格を構成する複数の要素毎に要素の性状を設定することができると共に要素間の関係を要素間毎に設定することができる。したがって、砂から粘土に至るまでのあらゆる不飽和の中間土の地盤についての土骨格における変形解析を精度よく行なうことができる。また、粘土や砂，中間土が層をなす不飽和土としての地盤についての土骨格における変形解析を精度よく解析することができる。【００５４】
あるいは、本発明の空気と水と土骨格の連成計算装置において、前記外力に関するデータは、空気と水と土骨格の連成系に与える荷重及び変位に関するデータが含まれるデータであるものとすることもできるし、前記外力に関するデータは、空気と水と土骨格の連成系に与える振動に関するデータが含まれるデータであるものとすることもできる。こうすれば、不飽和土に対して、荷重及び変位に対する土骨格における変形解析を行なうことができると共に振動に対する土骨格における変形解析を行なうこともできる。【００５５】
また、本発明の空気と水と土骨格の連成計算装置において、前記解析手段は、前記外力に関するデータを複数回変更して解析する手段であるものとすることもできる。この場合、前記解析手段は、前記外力に関するデータを変更するときには変更前の解析結果を初期状態の土骨格の状態に関するデータとして用いて解析する手段であるものとすることもできる。こうすれば、不飽和土に対して、土の状態と載荷の状態に応じた土骨格における変形解析を繰り返し行うことができる。また、不飽和土に対して、動的な解析に連続して静的な解析を行なったり、逆に静的な解析に連続して動的な解析を行なったりすることができる。【００５６】
本発明の空気と水と土骨格の連成計算方法は、
土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を構築し、該連立方程式に変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場、間隙水圧場及び間隙空気圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう、
ことを特徴とする。【００５７】
本発明の空気と水と土骨格の連成計算方法では、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を変位、変位速度又は加速度に関する幾何的境界条件と応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場の時刻歴応答を求めるから、土骨格の間隙に水の他に空気が存在する不飽和土に対して、加速度の影響が大きいときには動的解析による土骨格の変形解析を行なうことができ、加速度の影響が小さいときには静的解析による土骨格の変形解析を行なうことができる。したがって、土骨格の間隙に水が飽和している飽和土だけでなく、土骨格の間隙に水が飽和していない土、即ち、土骨格の間隙に水の他に空気が存在する不飽和土に対して、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。【００５８】
こうした本発明の空気と水と土骨格の連成計算方法において、前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水・間隙空気の時間的体積変化率に関する第５項目と、間隙水圧の荷重変化率寄与に関する第６項目と、間隙空気圧の荷重変化率寄与に関する第７項目と、間隙水の流量変換に関する第８項目と、間隙空気の流量変換に関する第９項目と、質量に関する第１０項目と、土骨格の接線剛性に関する第１１項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する、ことを特徴とするものとすることもできる。こうすれば、不飽和土に対して、より適正に土骨格の変形解析を行なうことができる。【００５９】
また、本発明の空気と水と土骨格の連成計算方法おいて、前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、により構築される方程式である、ものとすることもできる。この場合、前記空気と水と土骨格からなる三相混合体の運動方程式は上述の式(A1)により表わされ、前記水と土骨格間での質量保存式は上述の式(B1)により表わされ、前記空気と土骨格間での質量保存式は上述の式(C1)により表わされる、ものとすることもできる。この場合、式(A1)，式(B1)，及び式(C1)は、本発明の空気と水と土骨格の連成計算装置と同様に、間隙水を非圧縮性流体（Kwが無限大）として扱うことにより上述の式(A1)’、式(B1)’および式(C1)’と示される。式(A1)’は、空間的離散化手法として有限要素法を用い、さらに間隙水圧および間隙空気圧を要素中心に割り当てる田村流・Christian流に基づくと、上述の次式(A1)”と示され、式(B1)’と式(C1)’についても空間的に離散化し、同様に田村流・Christian流に基づいて離散化を施すと、上述の式(B1)”と式(C1)”を得ることができる。式(A1)”、式(B1)”、式(C1)”を常微分方程式の形に整理すると上述の式(E1)を得ることができ、連立一次方程式の形に帰着させるためにWilsonのθ法に準拠すれば、最終的には、上述の式(F1)の連立一次方程式に帰着させることができる。【００６０】
本発明の空気と水と土骨格の連成計算方法において、前記連立速度型運動方程式は、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧と飽和度に関するものであり、前記土骨格の変形解析は、前記幾何的境界条件と前記力学的境界条件と前記水理境界条件と前記空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場と飽和度場の時刻歴応答を求めることにより行なうものである、ものとすることもできる。このように間隙水の圧縮性を考慮すると共に飽和度についてより詳細に考慮することにより、土骨格の変形解析をより詳細に行なうことができる。【００６１】
この間隙水の圧縮性を考慮する態様の本発明の空気と水と土骨格の連成計算方法において、前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水の時間的体積変化率に関する第５項目と、間隙空気の時間的体積変化率に関する第６項目と、飽和度の時間的変化率に関する第７項目、間隙水圧の荷重変化率寄与に関する第８項目と、間隙空気圧の荷重変化率寄与に関する第９項目と、飽和度の荷重変化率寄与に関する第１０項目と、間隙水の流量変換に関する第１１項目と、間隙空気の流量変換に関する第１２項目と、質量に関する第１３項目と、土骨格の接線剛性に関する第１４項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築する、ことを特徴とするものとすることもできる。こうすれば、不飽和土に対して、より適正に土骨格の変形解析を行なうことができる。【００６２】
また、間隙水の圧縮性を考慮する態様の本発明の空気と水と土骨格の連成計算方法おいて、前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、土の水分特性に関する式と、により構築される方程式である、ものとすることもできる。この場合、前記空気と水と土骨格からなる三相混合体の運動方程式は上述の式(A2)により表わされ、前記水と土骨格間での質量保存式は上述の式(B2)により表わされ、前記空気と土骨格間での質量保存式は上述の式(C2)により表わされ、前記土の水分特性に関する式は上述の式(D2)により表わされる、ものとすることもできる。この場合、式(A2)，式(B2)，式(C2)，及び式(D2)は、本発明の空気と水と土骨格の連成計算装置と同様に、土の水分特性に関する式は、サクションとサクションの土骨格から見た物質時間微分と飽和度と飽和度の土骨格から見た物質時間微分の関係で表す形を用いることにより上述の式(A2)’、式(B2)’、式(C2)’および式(D2)’と示される。式(A2)’は、空間的離散化手法として有限要素法を用い、さらに田村流・Christian流に倣って、間隙水圧、間隙空気圧および飽和度を要素中心に割り当てると、上述の次式(A2)”と示され、式(B2)’と式(C2)’と式(D2)’についても空間的に離散化し、同様に田村流・Christian流に倣って離散化を施すと、上述の式(B2)”と式(C2)”と式(D2)”を得ることができる。式(A2)”、式(B2)”、式(C2)”、式(D2)”を常微分方程式の形に整理すると上述の式(E2)を得ることができ、連立一次方程式の形に帰着させるためにWilsonのθ法に準拠すれば、最終的には、上述の式(F2)の連立一次方程式に帰着させることができる。【００６３】
本発明の空気と水と土骨格の連成計算プログラムは、
土骨格の変形解析を行なう空気と水と土骨格の連成計算プログラムであって、
土骨格の状態に関するデータと外力に関するデータとを入力するデータ入力モジュールと、
前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を構築し、該連立速度型運動方程式に変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう解析モジュールと、
を備えることを要旨とする。【００６４】
この本発明の空気と水と土骨格の連成計算プログラムでは、このプログラムをコンピュータにインストールして起動することにより、コンピュータを上述した本発明の空気と水と土骨格の連成計算装置として機能させることができる。したがって、コンピュータを、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧に関する連立速度型運動方程式を変位、変位速度又は加速度に関する幾何的境界条件と応力又は応力速度に関する力学的境界条件と間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件と間隙空気圧又は間隙空気の流量に関する空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場の時刻歴応答を求める空気と水と土骨格の連成計算装置として機能させるから、コンピュータを、土骨格の間隙に水の他に空気が存在する不飽和土に対して、加速度の影響が大きいときには動的解析による土骨格の変形解析を行なうことができると共に加速度の影響が小さいときには静的解析による土骨格の変形解析を行なうことができる装置として機能させることができる。この結果、コンピュータを、土骨格の間隙に水が飽和している飽和土だけでなく、土骨格の間隙に水が飽和していない土、即ち、土骨格の間隙に水の他に空気が存在する不飽和土に対して、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる装置として機能させることができる。【００６５】
こうした本発明の空気と水と土骨格の連成計算プログラムにおいて、前記解析モジュールは、前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水・間隙空気の時間的体積変化率に関する第５項目と、間隙水圧の荷重変化率寄与に関する第６項目と、間隙空気圧の荷重変化率寄与に関する第７項目と、間隙水の流量変換に関する第８項目と、間隙空気の流量変換に関する第９項目と、質量に関する第１０項目と、土骨格の接線剛性に関する第１１項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築するモジュールである、ものとすることもできる。こうすれば、コンピュータを、不飽和土に対して、より適正に土骨格の変形解析を行なうことができる装置として機能させることができる。【００６６】
また、本発明の空気と水と土骨格の連成計算プログラムおいて、前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、により構築される方程式である、ものとすることもできる。この場合、前記空気と水と土骨格からなる三相混合体の運動方程式は上述の式(A1)により表わされ、前記水と土骨格間での質量保存式は上述の式(B1)により表わされ、前記空気と土骨格間での質量保存式は上述の式(C1)により表わされる、ものとすることもできる。この場合、式(A1)，式(B1)，及び式(C1)は、本発明の空気と水と土骨格の連成計算装置と同様に、間隙水を非圧縮性流体（Kwが無限大）として扱うことにより上述の式(A1)’、式(B1)’および式(C1)’と示される。式(A1)’は、空間的離散化手法として有限要素法を用い、さらに間隙水圧および間隙空気圧を要素中心に割り当てる田村流・Christian流に基づくと、上述の次式(A1)”と示され、式(B1)’と式(C1)’についても空間的に離散化し、同様に田村流・Christian流に基づいて離散化を施すと、上述の式(B1)”と式(C1)”を得ることができる。式(A1)”、式(B1)”、式(C1)”を常微分方程式の形に整理すると上述の式(E1)を得ることができ、連立一次方程式の形に帰着させるためにWilsonのθ法に準拠すれば、最終的には、上述の式(F1)の連立一次方程式に帰着させることができる。【００６７】
本発明の空気と水と土骨格の連成計算プログラムにおいて、前記連立速度型運動方程式は、土骨格の加加速度項を含む土骨格の速度と間隙水圧と間隙空気圧と飽和度に関するものであり、前記土骨格の変形解析は、前記幾何的境界条件と前記力学的境界条件と前記水理境界条件と前記空気境界条件とを与えて随時積分し、速度場と間隙水圧場と間隙空気圧場と飽和度場の時刻歴応答を求めることにより行なうものである、ものとすることもできる。このように間隙水の圧縮性を考慮すると共に飽和度についてより詳細に考慮することにより、コンピュータを、土骨格の変形解析をより詳細に行なうことができる装置として機能させることができる。【００６８】
この間隙水の圧縮性を考慮する態様の本発明の空気と水と土骨格の連成計算プログラムにおいて、前記解析モジュールは、前記土骨格の状態に関するデータと外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と、土の透水性に関する第２項目と、土の透気性に関する第３項目と、土の水分特性に関する第４項目と、間隙水の時間的体積変化率に関する第５項目と、間隙空気の時間的体積変化率に関する第６項目と、飽和度の時間的変化率に関する第７項目、間隙水圧の荷重変化率寄与に関する第８項目と、間隙空気圧の荷重変化率寄与に関する第９項目と、飽和度の荷重変化率寄与に関する第１０項目と、間隙水の流量変換に関する第１１項目と、間隙空気の流量変換に関する第１２項目と、質量に関する第１３項目と、土骨格の接線剛性に関する第１４項目とを計算し、該計算した項目を用いて前記連立速度型運動方程式を構築するモジュールである、ものとすることもできる。こうすれば、コンピュータを、不飽和土に対して、より適正に土骨格の変形解析を行なうことができる装置として機能させることができる。【００６９】
また、本発明の空気と水と土骨格の連成計算プログラムおいて、前記連立速度型運動方程式は、空気と水と土骨格からなる三相混合体の運動方程式と、水と土骨格間での質量保存式と、空気と土骨格間での質量保存式と、土の水分特性に関する式と、により構築される方程式である、ものとすることもできる。この場合、前記空気と水と土骨格からなる三相混合体の運動方程式は上述の式(A2)により表わされ、前記水と土骨格間での質量保存式は上述の式(B2)により表わされ、前記空気と土骨格間での質量保存式は上述の式(C2)により表わされ、前記土の水分特性に関する式は式(D2)により表わされるものとすることもできる。この場合、式(A2)，式(B2)，式(C2)，及び式(D2)は、本発明の空気と水と土骨格の連成計算装置と同様に、サクションとサクションの土骨格から見た物質時間微分と飽和度と飽和度の土骨格から見た物質時間微分の関係で表す形を用いることにより上述の式(A2)’、式(B2)’、式(C2)’および式(D2)’と示される。式(A2)’は、空間的離散化手法として有限要素法を用い、さらに田村流・Christian流に倣って、間隙水圧、間隙空気圧および飽和度を要素中心に割り当てると、上述の次式(A2)”と示され、式(B2)’と式(C2)’と式(D2)’についても空間的に離散化し、同様に田村流・Christian流に倣って離散化を施すと、上述の式(B2)”と式(C2)”と式(D2)”を得ることができる。式(A2)”、式(B2)”、式(C2)”、式(D2)”を常微分方程式の形に整理すると上述の式(E2)を得ることができ、連立一次方程式の形に帰着させるためにWilsonのθ法に準拠すれば、最終的には、上述の式(F2)の連立一次方程式に帰着させることができる。
【図面の簡単な説明】
【００７０】
【図１】本発明の第１実施例の空気と水と土骨格の連成計算装置２０の構成の概略を示す構成図である。
【図２】第１実施例の空気と水と土骨格の連成計算プログラム３０における処理の一例を示すフローチャートである。
【図３】解析対象としての地盤を２行２列のマトリクスの４要素からなるときの具体例を示す説明図である。
【図４】各要素における局所節点を説明する説明図である。
【図５】計算条件の設定を入力する具体例としてのファイルの一例を示す説明図である。
【図６】計算条件の設定を入力する具体例としてのファイルの一例を示す説明図である。
【図７】土の設定を入力する具体例としてのファイルの一例を示す説明図である。
【図８】具体例としての地盤立体構造モデルの設定を入力するファイルの一例を示す説明図である。
【図９】２次元平面ひずみ条件での間隙水および間隙空気の流れをモデル化した説明図である。
【図１０】軸対称条件での間隙水および間隙空気の流れをモデル化した説明図である。
【図１１】３次元条件での間隙水および間隙空気の流れをモデル化した説明図である。
【図１２】出力ファイルの一例を示す説明図である。
【図１３】空気と水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて空気と水と土骨格の連成計算プログラムを複数回に亘って実行する際の処理の一例を示すフローチャートである。
【図１４】初期飽和度が最大飽和度の場合に対して行なった非排気・非排水三軸圧縮試験のシミュレーション結果の一例を示す説明図である。
【図１５】初期飽和度が最大飽和度の場合に対して行なった排気・排水三軸圧縮試験のシミュレーション結果の一例を示す説明図である。
【図１６】初期飽和度をそれぞれ変えた５つの場合に対して行なった非排気・非排水三軸圧縮試験のシミュレーション結果の一例を示す説明図である。
【図１７】第２実施例の空気と水と土骨格の連成計算プログラム２３０における処理の一例を示すフローチャートである。
【図１８】第２実施例の方法を用いて行なった、初期飽和度をそれぞれ変えた５つの場合に対して行なった非排気・非排水三軸圧縮試験のシミュレーション結果の一例を示す説明図である。
【図１９】ほぼ飽和状態にある粘土地盤上に、不飽和状態にある盛土を載荷したときの飽和度分布を示す説明図である。
【図２０】水位上昇による盛土内の浸潤面生成過程を示す説明図である。
【発明を実施するための形態】
【００７１】
次に、本発明を実施するための形態を実施例を用いて説明する。
【実施例１】
【００７２】
図１は、本発明の第１実施例の空気と水と土骨格の連成計算装置２０の構成の概略を示す構成図である。第１実施例の空気と水と土骨格の連成計算装置２０は、図示するように、一般的な汎用コンピュータ２２にアプリケーションソフトウエアとしての空気と水と土骨格の連成計算プログラム３０がインストールされたものとして構成されている。空気と水と土骨格の連成計算プログラム３０は、計算条件の設定や土の性状の設定、地盤立体構造モデルの設定などを行なうためにデータを入力するデータ入力モジュール３２と、入力されたデータを用いて不飽和土の力学的挙動を弾塑性構成式を用いて計算する計算モジュール３４と、計算した結果を出力する結果出力モジュール３６と、から構成されている。【００７３】
図２は、第１実施例の空気と水と土骨格の連成計算プログラム３０における処理の一例を示すフローチャートである。以下、このフローチャートを用いて空気と水と土骨格の連成計算プログラム３０の処理について説明する。ここで、空気と水と土骨格の連成計算プログラム３０の処理を説明するに当たり、説明を簡易にするために、解析対象としての地盤を２行２列のマトリクスの４要素とするときを具体例として用いる。具体例の４要素を図３に示す。図示するように、４要素については下行左から要素１，要素２、上行左から要素３，要素４と要素番号を割り振り、節点については下段左から節点１，節点２，節点３、中段左から節点４，節点５，節点６、上段左から節点７，節点８，節点９と節点番号を割り振るものとした。また、各要素における局所節点については、図４に示すように、左下から左回りに局所節点１，局所節点２，局所節点３，局所節点４と局所節点番号を割り振り、各要素の節点間については図４に示すように下の節点間から左回りに節点間１，節点間２，節点間３，節点間４と節点間番号を割り振るものとした。【００７４】
空気と水と土骨格の連成計算プログラム３０が実行されると、まず、計算条件の設定の入力処理が行なわれる（ステップＳ１００）。ここで、計算条件の設定の入力としては、例えば、差分法に基づいた時間数値積分方法の選択、時刻の更新回数（計算ステップ数），２次元平面ひずみや軸対称，３次元などの計算条件を入力する。計算条件の設定を入力する具体例としてのファイルの一例を図５および図６に示す。図５のファイルでは、上から、計算条件として「平面ひずみ」であるか「軸対称」であるか「３次元」であるかの設定、変形理論として「微小変形理論」を選択するか「有限変形理論」を選択するかの設定、有限変形理論時に使用する客観性のある有効応力速度として「有効Cauchy応力のJaumann速度」であるか「Green-Nagdhi速度」であるかの設定、収束の有無の設定、構成式の選択として「下負荷面・上負荷面・回転硬化付きオリジナルカムクレイモデル」を選択するか「下負荷面・上負荷面・回転硬化付き修正カムクレイモデル」を選択するかの設定、状態量出力データの方法として「Gauss点の要素平均」であるか「ある特定のGauss点の値」であるかの設定、出力データの形式として「Text」であるか「Binary」であるかの設定、自重考慮の有無の設定、を入力する。図６のファイルでは、上から、現在のプログラム使用番号、計算ステップ数、１計算ステップ当たりの計算時間DT（以下Δｔと記載する場合がある。）(sec/step)、盛土・掘削時の荷重分割数、分布荷重作用時の分布荷重分割数、土の材料定数と初期値入力ファイル、メッシュデータ入力ファイル、次のプログラム実行時に用いる計算初期条件用出力ファイル指定を入力する。【００７５】
続いて、地盤の各要素に対する粘土・中間土・砂などの土の設定を入力する（ステップＳ１１０）。この土の設定では、材料定数や初期状態が入力される。土の設定を入力する具体例としてのファイルの一例を図７に示す。図７のファイルでは、弾塑性体である土の種類の数や弾性体の種類の数、材料番号、過圧密の発展則に用いる塑性変形の進展尺度の指定、構造の発展則に用いる塑性変形の進展尺度の指定、練返し土の等方正規圧密線の切片N、限界状態定数M、圧縮指数λ、膨潤指数κ、ポアソン比ν、飽和透水係数ksw、土粒子の比重Gs、初期側圧係数K0、初期過圧密比1/R0、初期構造の程度1/R*0、正規圧密土化指数m、U*の形状、構造劣化指数a、構造劣化指数b、構造劣化指数c、構造劣化指数cs、回転硬化指数br、回転硬化限界面mb、初期異方性の程度、初期間隙比、乾燥透気係数kda、水分特性曲線のパラメータswmax，swmin，αv，mv，nv、透水係数・透気係数モデルのパラメータξ，γ、空気の温度、地盤の初期上載圧、セル圧(引張：正)、空気圧、サクション、空気の気体定数、単位系の設定（時間、力、長さ）、粘性境界の密度、圧縮波速度、せん断波速度などを入力する。【００７６】
そして、地盤立体構造モデルの設定を入力する（ステップＳ１２０）。この地盤立体構造モデルの設定では、地盤モデルの要素に対する条件などを入力する。図３および図４に示す具体例としての地盤立体構造モデルの設定を入力するファイルの一例を図８に示す。図８のファイルの第１行目の１行６列のデータは、地盤立体構造モデルの要素数、節点数、分布荷重速度作用要素数、１要素中の節点数、計算ステップ数、弾性体の要素数である。第２行目から第１０行目までの９行５列のデータは、左から順に、節点番号、ｘ方向の境界条件の種別、ｘ方向の境界条件の値、ｙ方向の境界条件の種別、ｙ方向の境界条件の値である。ここで、境界条件の種別としては、座標制御（０）、荷重速度制御（１）、変位速度制御（−１）、変位加速度制御（−２）、粘性を考慮した加速度境界（２）である。図８のファイルの第１１行目から第１４行目までの４行５列のデータは、左端が要素番号、要素番号の右から順に各要素における局所節点番号に対応する全体における節点番号である（局所節点番号と全体節点番号の対応の設定）。ファイルの第１５行目から第２３行目までの９行３列のデータは、節点番号、ｘ方向の座標、ｙ方向の座標である。ファイルの第２４行目および第２５行目については、第２４行目は反力の出力を必要とする節点数であり、第２５行目はその節点番号である。この場合、反力の出力は節点番号に対してｘ方向とｙ方向とがあるため、例えば節点番号７の場合のｘ方向については７×２−１＝１３として示し、ｙ方向については７×２＝１４として示している。ファイルの第２６行目および第２７行目の２行６列のデータは、左から順に、分布荷重速度作用要素、作用辺（作用節点間）、その左側節点の分布荷重速度のｘ方向の値、その左側節点の分布荷重速度のｙ方向の値、その右側節点の分布荷重速度のｘ方向の値、その右側節点の分布荷重速度のｙ方向の値である。ファイルの第２８行目から第３１行目までの４行６列のデータは、左から順に、要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号である（節点間番号と全体要素番号の対応の設定）。ただし、隣接する要素へは非排水になっている場合は、その要素自身の要素番号を入力する。ファイルの第３２行目から第３５行目までの４行６列のデータは、左から順に、要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号である（節点間番号と全体要素番号の対応の設定）。ただし、隣接する要素へは非排気になっている場合は、その要素番号を入力する。ファイルの第３６行目のデータは、２節点間の長さ不変の条件を与えて計算を行う場合の条件数である。２節点間の長さ不変の条件数が存在する場合には、この第３６行目のデータの下に新たな行を設けて条件番号、当該条件に関する節点間の指定がなされる。第３７行目のデータは３節点間の角度不変の条件数であり、第３８行目のデータは、この３節点間の角度不変の条件に対して、左から順に、条件番号、当該条件に関する節点の��定である。第３９行目のデータは、前２節点の相対位置ベクトルと後２節点で指定された２節点の相対速度の方向が不変の条件数であり、第４０行目のデータはこの２節点の相対速度の方向が不変の条件に対して、左から順に、条件番号、当該条件に関する節点の指定である。第４１行目のデータは、等変位（速度）の条件数である。等変位（速度）の条件が存在する場合には、この第４１行目のデータの下に新たな行を設けて条件番号、当該条件に関する事項の指定がなされる。第４２行目のデータは、土の種類である。第４３行目のデータは、左から順に、モデルに入力される規則波の入力数、モデルに入力される不規則波の入力数であり、第４４行目のデータは、モデルに入力される規則波に対して、左から順に、規則波振動の種類、作用する節点数、振幅、周期である。第４５行目のデータは、規則波振動が作用する節点番号である。この場合、規則波振動が作用する方向はx方向とｙ方向とがあるため、例えば節点番号７の場合のｘ方向については７×２−１＝１３として示し、ｙ方向については７×２＝１４として示している。ここで、モデルに入力する不規則波が存在する場合は、第４５行目の下に新たな行を設けて不規則波のファイルを指定する。【００７７】
こうしてデータの入力が終了すると、入力したデータを用いて計算処理を行なう（ステップＳ１３０〜Ｓ２７０）。計算処理は、第１実施例では、時間積分の方法としてWilsonのθ法を用いて行なったが、Newmarkのβ法などを用いて行なうものとしてもよい。また、第１実施例では、空間的離散化に対して、有限要素法を用いて間隙水圧および間隙空気圧を要素中心に割り当てる田村流・Christian流に基づいて行ったが、節点に間隙水圧を割り当てるSandhu流に基づいて行う方法やあるいはメッシュフリー法などの方法を用いて行ってもよい。以下、Wilsonのθ法と田村流に基づく有限要素法を用いた処理について説明する。まず、時刻ｔ＝ｔ＋θΔｔにおける「全体位置ベクトル」、「全体速度関連増分ベクトル」ならびに「全体加速度関連増分ベクトル」と、各種状態（有効応力と間隙水圧・間隙空気圧、構造・過圧密・異方性、飽和度等）を次式（１）〜（５）を用いて予測する（ステップＳ１３０）。この予測計算処理は、短時間に収束させるために行なうものであり、必ずしも必要な処理ではない。予測計算処理を行わない場合は、時刻ｔ＝ｔで最終的に決められた値、ただし、初期時刻の場合は初期値、を時刻ｔ＝ｔ＋θΔｔにおける値として用いる。ここで、式（１）の左辺や右辺の各項のように全節点のあるベクトル量をまとめて列ベクトル表現するときには、式（６）に示すように成分表示される。また、式（４）の左辺や右辺の各項のように全要素のあるスカラー量をまとめて列ベクトル表現するときには、式（７）に示すように成分表示される。【００７８】
【数２３】
2014002977.tif 000025
【００７９】
【数２４】
2014002977.tif 000026
【００８０】
【数２５】
2014002977.tif 000027
【００８１】
次に、節点の変位速度を土骨格の時間的体積変化率に変換するマトリクスとして、次式（８）により全体Lマトリクス（全体体積変化率マトリクス）を計算する（ステップＳ１４０）。全体Lマトリクスの各要素は式（９）により示される。全体LマトリクスはＮＥ行（ＮＰ×ｒ）列のマトリクス（ｒ：一節点あたりの成分の数）で、地盤立体構造モデルの設定時に指定した全体節点番号と局所節点番号の対応に従うとき、式（９）の要素ｉにおける要素Lマトリクスにおいて、ｍ番目の節点が全体の節点番号と一致しない列成分に対しては、数値ゼロを代入する。要素Lマトリクスを構成するBvマトリクスは、２次元平面ひずみ条件の場合は式（１０）により、軸対称条件の場合は式（１１）により、３次元条件の場合には式（１２）により表わされる。【００８２】
【数２６】
2014002977.tif 000028
【００８３】
【数２７】
2014002977.tif 000029
【００８４】
式（９）などの積分はいずれも、式（１３）のように、全体座標系における解析対象領域を局所座標系に変換した後、Gaussの数値積分法を用いて行われる。ここに、式（１３）は３次元条件の場合における体積積分を示す。なお、２次元平面ひずみ条件の場合は解析対象面に垂直な方向に単位厚さを考慮して、軸対称条件の場合は軸の周りの回転対称性を考慮して積分が行われる。【００８５】
【数２８】
2014002977.tif 000030
【００８６】
続いて、土の透水性を表す全体HWマトリクス（全体透水性マトリクス）を式（１４）により、土の透気性を表す全体HAマトリクス（全体透気性マトリクス）を式（１５）により計算する（ステップＳ１５０）。ここで、排水・非排水等の水理境界条件と排気・非排気等の空気境界条件は、それぞれ式（１６）の要素HWマトリクスの成分と式（１７）の要素HAマトリクスの成分を変化させることで与える。例えば、非排水条件の場合は、この条件が設定されている節点間に対応する要素HWマトリクスの成分の値をゼロとする。非排気条件の場合は、この条件が設定されている節点間に対応する要素HAマトリクスの成分の値をゼロとする。全体HWマトリクスと全体HAマトリクスはＮＥ行ＮＥ列のマトリクスで、地盤立体構造モデルの設定時に指定した当該要素の節点間番号と全体要素番号の対応に従うとき、次式（１４）と次式（１５）で示す要素ｉの要素HWマトリクスと要素HAマトリクスにおいてｍ番目の節点間番号が全体の要素番号と一致しない全体HWマトリクスと全体HAマトリクスの成分に対しては、数値ゼロを代入する。なお、式（１６）と式（１７）の要素ｉにおける要素HWマトリクスと要素HAマトリクスにおける最後の成分は、要素ｉ番目（すなわち、ｉ行目）のｉ番目の列成分に代入する。式（１６）の右辺の各成分は、式（１８）により示される。式（１７）の右辺の各成分は、式（１９）により示される。各記号は、２次元平面ひずみ条件の場合は図９を、軸対称条件の場合は図１０を、３次元条件の場合は図１１を参照することになる。【００８７】
【数２９】
2014002977.tif 000031
【００８８】
【数３０】
2014002977.tif 000032
【００８９】
次に、全体LBWマトリクス（全体水流量変換マトリクス）を式（２０）により、全体LBAマトリクス（全体空気流量変換マトリクス）を式（２５）により計算する（ステップS１６０）。全体LBWマトリクスの各要素は式（２１）により、全体LBAマトリクスの各要素は式（２６）により示される。全体LBWマトリクスまたは全体LBAマトリクスはＮＥ行（ＮＰ×ｒ）列のマトリクスで、地盤立体構造モデルの設定時に指定した全体節点番号と局所節点番号の対応に従うとき、式（２１）または式（２６）の要素ｉにおける要素LBWマトリクスまたは要素LBAマトリクスにおいて、ｎ番目の節点が全体の節点番号と一致しない列成分に対しては、数値ゼロを代入する。要素LBWマトリクスまたは要素ＬBAマトリクスを構成するLBマトリクスは式（２２）により表される。要素LBWマトリクスの各成分は式（２３）により表され、隣接する要素間の水理条件が排水条件か非排水条件かによって足し合わされる値が異なり、排水条件の場合に用いられるβjwiの値を式（２４）に示す。要素LBAマトリクスの各成分は式（２７）により表され、隣接する要素間の空気条件が排気条件か非排気条件かによって足し合わされる値が異なり、排気条件の場合に用いられるβjaiの値を式（２８）に示す。各記号は、２次元平面ひずみ条件の場合は図９を、軸対称条件の場合は図１０を、３次元条件の場合は図１１を参照することになる。【００９０】
【数３１】
2014002977.tif 000033
【００９１】
【数３２】
2014002977.tif 000034
【００９２】
そして、次式（２９）により時刻ｔ＝ｔ＋θΔｔにおける全体荷重増分ベクトルを、式（４４）より間隙水の全体流量増分ベクトルを、さらに式（４６）より間隙空気の全体流量増分ベクトルの計算を行なう（ステップＳ１７０）。式（２９）の右辺第１項は公称応力速度ベクトルに関する項で式（３０）から構成され、右辺第２項はGreen-Naghdi速度を用いる場合に現れる項で式（３１）から構成される。ここで、式（３０）の右辺のＮマトリクスと式（３１）の右辺のＢマトリクスは、２次元平面ひずみ条件については式（３２）および式（３３）で表わされ、軸対称条件については式（３４）および式（３５）で表わされ、３次元条件について式（３６）および式（３７）で表わされる。また、式（３０）の公称応力速度ベクトルに関する項は、式（３８）に示すように、右辺第１項に示す境界を通じて外部から与えられる表面力速度ベクトルと右辺第２項に示すこの境界の幾何形状変化に伴って現れる表面力ベクトルの増分の和からなる。このとき、右辺第１項の時刻ｔ＝ｔ＋θΔｔにおける表面力速度ベクトルは、Wilsonのθ法に従う式（３９）により、時刻ｔ＝ｔ＋Δｔにおいて境界条件として設定された表面力速度ベクトルを換算して与える。式（３８）右辺第２項に含まれる応力又は応力速度に関する力学的境界に立てた単位外向き法線ベクトルからなるnマトリクスは、２次元平面ひずみ条件または軸対称条件の場合は次式（４０）で与えられ、３次元条件の場合は式（４１）で与えられる。また、式（３１）の右辺第１項の有効応力成分を含む列ベクトルは、２次元平面ひずみ条件または軸対称条件の場合は式（４２）で与えられ、３次元条件の場合は式（４３）で与えられる。【００９３】
【数３３】
2014002977.tif 000035
【００９４】
【数３４】
2014002977.tif 000036
【００９５】
【数３５】
2014002977.tif 000037
【００９６】
【数３６】
2014002977.tif 000038
【００９７】
【数３７】
2014002977.tif 000039
【００９８】
【数３８】
2014002977.tif 000040
【００９９】
【数３９】
2014002977.tif 000041
【０１００】
【数４０】
2014002977.tif 000042
【０１０１】
また、式（４４）に示す時刻ｔ＝ｔ＋θΔｔにおける間隙水の全体流量増分ベクトルは、式（４５）で与える全体流量速度ベクトルにθΔｔを乗じたものとして与える。さらに、式（４６）に示す時刻ｔ＝ｔ＋θΔｔにおける間隙空気の全体流量増分ベクトルは、式（４７）で与える間隙空気の全体流量速度ベクトルにθΔｔを乗じたものとして与える。【０１０２】
【数４１】
2014002977.tif 000043
【０１０３】
【数４２】
2014002977.tif 000044
【０１０４】
次に、全体Mマトリクス（全体質量マトリクス）を式（４８）により計算する（ステップＳ１８０）。ここで、式（４８）右辺の要素Ｍマトリクスは、式（４９）で表される。なお、全体Mマトリクスは（ＮＰ×ｒ）行（ＮＰ×ｒ）列のマトリクスとなる。【０１０５】
【数４３】
2014002977.tif 000045
【０１０６】
そして、全体Kマトリクス（土骨格の全体接線剛性マトリクス）を式（５０）により計算する（ステップＳ１９０）。ここで、全体Kマトリクスは（ＮＰ×ｒ）行（ＮＰ×ｒ）列のマトリクスとなる。全体Kマトリクスを構成する要素Kマトリクスは式（５１）で示される。この式（５１）中、右辺第１項は特定された土材料（砂や粘土など）の力学挙動を与える接線剛性を示し、右辺第２項は刻々の土の形状変化による接線剛性への寄与を表わし、右辺第３項は衝撃荷重など土骨格に加速度が顕著に働く条件下での接線剛性への寄与を示し、右辺第４項は物体力作用条件下で生じる接線剛性への寄与を示す。材料定数の中で、特にこれらの「構造低位化指数」、「正規圧密土化指数」および「回転硬化限界定数」に対して相対的な差異を与えることによって、砂から中間土、粘土までを連続的に取り扱うことができる。式（５１）の右辺第１項に含まれるDepマトリクス（弾塑性マトリクス）は、弾塑性状態（負荷状態）では式（５２）で与えられ、さらに式（５２）は２次元平面ひずみ条件で式（５３）と式（５４）より与えられ、軸対称条件で式（５５）と式（５６）より与えられ、３次元条件で式（５７）と式（５８）により与えられる。また、弾性状態では式（５９）で与えられる。式（５３）〜式（５８）中に現れる記号は式（６０）より与えられ、式（６０）中に現れる記号は式（６１）を一個ないしは複数個利用して与えられる。【０１０７】
【数４４】
2014002977.tif 000046
【０１０８】
【数４５】
2014002977.tif 000047
【０１０９】
【数４６】
2014002977.tif 000048
【０１１０】
【数４７】
2014002977.tif 000049
【０１１１】
【数４８】
2014002977.tif 000050
【０１１２】
【数４９】
2014002977.tif 000051
【０１１３】
【数５０】
2014002977.tif 000052
【０１１４】
【数５１】
2014002977.tif 000053
【０１１５】
p'〜q応力空間での塑性圧縮と塑性膨張の敷居線の傾きと硬化と軟化の敷居線の傾きに対する式（６２）と式（６３）の表記は、土の骨格構造のうちで、構造R*（０＜R*≦１、R*が０に近いほど構造は高位であることを示す）の低位化の仕方と過圧密の程度Ｒ（０＜R≦１、Rが０に近いほど過圧密であることを示す）の解消の仕方を与える発展則において、次式（６４）および（６５）のように、塑性測度として塑性ストレッチングのノルムを用いた場合を示している。土の種類に応じて塑性ストレッチングのせん断成分や有効応力テンソルと塑性ストレッチングの内積で表す塑性仕事率などで与えることもできる。なお、微小変形解析を用いる場合、塑性ストレッチングは塑性ひずみ速度で表される。また、式（６５）に現れる「Ｕ*」および式（６６）に現れる「Ｕ」はここでは非負の関数とし、例えば、次式（６６）と（６７）で与える。また、異方性を表す回転硬化変数テンソルの発展則については、次式（６８）を用いている。【０１１６】
【数５２】
2014002977.tif 000054
【０１１７】
【数５３】
2014002977.tif 000055
【０１１８】
式（５２）の右辺第２項を構成するN’マトリクスとT1マトリクスは、２次元平面ひずみ条件では式（６９）と式（７０）より与えられ、軸対称条件では式（７１）と式（７２）より与えられ、３次元条件では式（７３）と式（７４）により与えられる。【０１１９】
【数５４】
2014002977.tif 000056
【０１２０】
【数５５】
2014002977.tif 000057
【０１２１】
【数５６】
2014002977.tif 000058
【０１２２】
次に、全体LWTマトリクス（全体間隙水圧分荷重変化率変換マトリクス）と全体LATマトリクス（全体間隙空気圧分荷重変化率変換マトリクス）を式（７５）と式（７７）により計算する（ステップＳ２００）。ここで、式（７５）右辺の要素LWTマトリクスは、式（７６）で表され、式（７７）右辺の要素LATマトリクスは、式（７８）で表される。なお、全体LWTマトリクスと全体LATマトリクスはどちらも（ＮＰ×ｒ）行ＮＥ列マトリクスとなる。【０１２３】
【数５７】
2014002977.tif 000059
【０１２４】
【数５８】
2014002977.tif 000060
【０１２５】
全体SWAマトリクス（全体水分特性マトリクス）と全体SACマトリクス（全体空気圧縮性マトリクス）をそれぞれ式（７９）と式（８０）より求める（ステップＳ２１０）。どちらもＮＥ行ＮＥ列マトリクスとなる。【０１２６】
【数５９】
2014002977.tif 000061
【０１２７】
全体透水修正原Lマトリクス、全体透気修正原Lマトリクス、全体飽和度修正原Lマトリクス、全体空気間隙比修正原Ｌマトリクス、全体間隙水修正原Lマトリクスおよび全体間隙空気修正原Lマトリクスをそれぞれ式（８１）、式（８２）、式（８３）、式（８４）、式（８５）および式（８６）より求める（ステップＳ２２０）。全体Lマトリクスと同様にして作られるＮＥ行（ＮＰ×ｒ）列のマトリクスである。ただし、全体透水修正原Lマトリクスと全体透気修正原Lマトリクスは、全体接線剛性方程式（連立１次方程式）を構築する上では、それぞれ全体間隙水修正原Lマトリクスと全体間隙空気修正原Lマトリクスに包含されるため、全体透水修正原Lマトリクスと全体透気修正原Lマトリクスは実質計算する必要はない。【０１２８】
【数６０】
2014002977.tif 000062
【０１２９】
【数６１】
2014002977.tif 000063
【０１３０】
次式（８７）に示す全体接線剛性方程式（連立１次方程式）より、変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件並びに間隙空気圧又は間隙空気の流量に関する空気境界条件を与えて未知の「加加速度（躍度）場」、間隙水圧場と間隙空気圧場の求解を行なう（ステップＳ２３０）。この式（９１）は、速度型運動方程式の弱形式の有限要素離散化と慣性項の影響を考慮した空気と水と土骨格連成式のモデル化によって最終的に得られる連立２階常微分方程式（式（８８））において、時間微分項を陰解法（差分法）に基づいて処理したもので、変形場に対してはWilsonのθ法を用いて、間隙水圧と間隙空位圧に対しては台形公式を利用した場合の（最終的に）解くべき連立１次方程式である。全体荷重増分ベクトル、間隙水の全体流量増分ベクトル及び間隙空気の全体流量増分ベクトル並びに全体Mマトリクス及び全体Kマトリクス並びに全体HWマトリクス及び全体HAマトリクス並びに全体LBWマトリクス及び全体LBAマトリクス並びに全体LWTマトリクス及び全体LATマトリクス並びに全体SWAマトリクス及び全体SACマトリクス並びに全体LW2マトリクス及び全体LA2マトリクス並びに全体LθnWマトリクス及び全体LθnAマトリクスは時刻ｔ＝ｔと時刻ｔ＝ｔ＋θΔｔとの間の時間ステップ内で繰り返し計算を行ないながら更新する。なお、式（８８）では、運動方程式に対しさらに土骨格から見た物質時間微分を施しているため、説明変数である速度ベクトルの時間に関する２階微分の項、すなわち、「加加速度（躍度とも称する）」項を有している点がもっとも特徴的である。これは、土骨格の構成式を速度型で与えるためだけでなく、地盤等の土の幾何形状変化を考慮するためで、これによって、変形から破壊まで、あるいは破壊後の挙動を対象にした土の体積変化の正確な計量を必要とする大変形解析を可能にするためである。またこの加加速度（躍度）導入のため、Wilsonのθ法の基礎となる線形加速度法の考えに基づいて、変形に対して式（８９）〜式（９１）を誘導している。なお、加速度の影響が小さいような応答が得られるときは、式（８８）は、静的解析で得られる連立１階常微分方程式と一致する。このことは、動的および静的解析の使い分けを必要としないことを意味する。また、例えば、二節点の間に距離が不変あるいは三節点がなす角度が不変の条件などの束縛条件が節点間に課せられる場合は、Lagrangeの未定乗数法を適用して式（８８）にその束縛条件を考慮する。時刻ｔ＝ｔ＋Δｔにおいて設定された変位、変位速度又は加速度に関する幾何的境界条件については、位置ベクトルの差すなわち変位ベクトルで与えられた場合は式（９２）により、速度ベクトルで与えられた場合は式（９３）により、加速度ベクトルで与えられた場合は式（９４）により、時刻ｔ＝ｔ＋θΔｔでの加加速度（躍度）ベクトルに換算して与える。【０１３１】
【数６２】
2014002977.tif 000064
【０１３２】
【数６３】
2014002977.tif 000065
【０１３３】
【数６４】
2014002977.tif 000066
【０１３４】
【数６５】
2014002977.tif 000067
【０１３５】
続いて、時刻ｔ＝ｔ＋θΔｔにおける加速度場、速度場、座標、各種状態量を計算すると共に各Gauss点の負荷状態を判定する（ステップＳ２４０）。時刻ｔ＝ｔ＋θΔｔにおける全体加速度関連増分ベクトル、全体速度関連増分ベクトルおよび座標ベクトルは、式（８７）で得られた時刻ｔ＝ｔ＋θΔｔにおける全体加加速度（躍度）関連増分ベクトルと時刻ｔ＝ｔにおける幾何的諸量から、式（８９）〜式（９１）より得られる次式（９５）〜式（９７）により算出することができる。時刻ｔ＝ｔ＋θΔｔにおける有効応力、構造の程度R*、過圧密の程度R、異方性などの各状態量Aについては、式（９８）により算出することができる。例えば、時刻ｔ＝ｔ＋θΔｔにおける要素内の各Gauss点での有効応力速度の場合は、式（９６）より得られた全体速度関連増分ベクトルのうち、当該要素が有する節点の速度ベクトルから式（９９）を用いて算出した後、式（９８）を適用して有効応力を算出する。なお、構造の程度、過圧密の程度、異方性などの速度は、Gauss点が弾塑性状態にある場合は式（６４）、式（６５）、式（６８）などを用いて算出し、弾性状態にある場合はそれらの速度をゼロとする。時刻ｔ＝ｔ＋θΔｔにおける「塑性体積ひずみ」は、下負荷面・上負荷面・回転硬化付き修正カムクレイモデルの場合、Gauss点が弾塑性状態にあるとき式（１０３）より算出するが、弾性状態にあった場合は、最後に弾塑性状態にあった時点での値のままとする。したがって、この弾性状態のとき、時刻ｔ＝ｔ＋θΔｔにおける過圧密の程度Rは、式（１０３）をRについて解いた式を用いて算出することができる。また、時刻ｔ＝ｔ＋θΔｔにおける間隙水圧および間隙空気圧は式（８７）の全体接線剛性方程式より直接決まるが、間隙水圧及び間隙空気圧の速度は式（９８）を逆に用いて求める。また、土要素の各Gauss点が弾塑性状態にあるか弾性状態にあるかについては、土要素の各節点から計算されるGauss点のストレッチングと有効応力などを用いて塑性乗数の分子から得られる次式（１０４）に示す負荷基準を満たすときに、そのGauss点の状態は弾塑性状態にあると判定することができる。なお、ストレッチングの列ベクトル表現したものは式（１０５）より計算されるが、時刻ｔ＝ｔ＋θΔｔにおけるその列ベクトルの成分表示は、２次元平面ひずみ条件の場合は式（１０６）により、軸対称条件の場合は式（１０７）により、３次元条件の場合は式（１０８）により表される。【０１３６】
【数６６】
2014002977.tif 000068
【０１３７】
【数６７】
2014002977.tif 000069
【０１３８】
【数６８】
2014002977.tif 000070
【０１３９】
【数６９】
2014002977.tif 000071
【０１４０】
【数７０】
2014002977.tif 000072
【０１４１】
そして、各要素の各Gauss点の有効応力を用いて収束判定を行なう（ステップＳ２５０）。第１実施例では、全要素内の全Gauss点で次式（１０９）を充足するか否かにより収束を判定することとした。ここで、式（１０９）中の「ε」は十分に小さい正の値である。収束が判定されないときには、ステップＳ１５０の全体HWマトリクスと全体HAマトリクスの計算に戻り、ステップＳ１５０〜ステップＳ２５０までの処理を繰り返す。【０１４２】
【数７１】
2014002977.tif 000073
【０１４３】
ステップＳ２５０で収束したと判定されると、時刻ｔ＝ｔ＋Δｔにおける加速度場、速度場、座標、各種状態量を計算すると共に各Gauss点の負荷状態を判定する（ステップＳ２６０）。時刻ｔ＝ｔ＋Δｔの座標などは、時刻ｔ＝ｔ＋θΔｔにおける全体加加速度（躍度）関連ベクトルと時刻ｔ＝ｔにおける全体加速度関連増分ベクトル、全体速度関連増分ベクトルおよび座標ベクトルから、次式（１１１）〜式（１１４）により算出することができる。また、時刻ｔ＝ｔ＋Δｔの各状態量とその速度については、時刻ｔ＝ｔと時刻ｔ＝ｔ＋θΔｔにおける状態量から式（１１５）および式（１１６）により計算することができる。そして、時刻ｔ＝ｔ＋Δｔの各状態量と式（１０４）と式（１０５）を用いて、時刻ｔ＝ｔ＋Δｔでの負荷状態の判定を行う。【０１４４】
【数７２】
2014002977.tif 000074
【０１４５】
【数７３】
2014002977.tif 000075
【０１４６】
そしてステップＳ１００で入力した計算ステップ数に至ったか否かを判定し（ステップＳ２７０）、計算ステップ数に至っていないときには、次の計算ステップとしてステップＳ１３０の時刻ｔ＝ｔ＋θΔｔにおける各種状態（有効応力と間隙水圧・間隙空気圧、構造・過圧密・異方性、等）の予測処理に戻る。一方、計算ステップ数に至ったときには、計算結果を指定した出力ファイルに出力して（ステップＳ２８０）、空気と水と土骨格の連成計算プログラム３０を終了する。【０１４７】
出力ファイルの一例を図１２に示す。この例では、第１行目は、左から順に、要素数、節点数、分布荷重速度作用要素数、節点数、節点数×次元数、節点数／要素×次元数、指定計算ステップ数、通算実施計算ステップ数、次元数、Gauss点の数／要素、弾性要素数である。２行目は、左から順に、計算条件、変形理論、有限変形理論時に使用する有効応力速度、収束の有無、構成式、状態量出力データの方法、出力データの形式、自重考慮の有無、に関する選択である。３行目から１１行目までは、左から順に、節点番号、ｘ方向の境界条件の種別、その値、ｙ方向の境界条件の種別、その値、である。１２行目から１５行目は、左端が要素番号、要素番号の右から順に各要素における局所節点番号に対応する全体における節点番号である。１６行目は、反力の出力を必要とする節点数であり、１７行目は、反力の出力を必要とする節点における方向込みの値として示される節点番号である。１８行目および１９行目は、左から順に、分布荷重速度作用要素、作用辺（作用節点間）、その左側節点の分布荷重速度のｘ方向の値、その左側節点の分布荷重速度のｙ方向の値、その右側節点の分布荷重速度のｘ方向の値、その右側節点の分布荷重速度のｙ方向の値である。２０行目から２３行目までは水理境界条件、２４行目から２７行目までは空気境界条件である。図１２では、「中略」として略しているが、各要素における各Gauss点の材料定数、各節点の成分の値、各要素の間隙水圧の値、各要素の間隙空気圧の値、各要素における各Gauss点の有効応力の成分の値、構造の程度の値なども出力される。「中略」の後の１行目は、２節点間の長さ不変の条件数であり、条件数が値１以上であればその条件番号や当該条件に関する節点の指定などがその下の行に出力される。「中略」の後の２行目は、３節点間の角度不変の条件数であり、次の３行目は、この３節点間の角度不変の条件における条件番号、当該条件に関する節点の指定（３箇所）である。「中略」の後の４行目は、前２節点の相対位置ベクトルと後２節点で指定された２節点の相対速度の方向が不変の条件数であり、その次の５行目は、その条件番号、該当条件に関する節点の指定である。「中略」の後の６行目は、等変位（速度）の条件数であり、等変位（速度）の条件があるときにはその節点番号などがその次の行に出力される。「中略」の後の７行目は、各要素における土の種類である。中略の後の８行目は、規則波の入力数と不規則波の入力数であり、規則波については９行目と１０行目に規則波振動の種類、作用する節点数、振幅、周期、規則波が作用する方向込みの節点番号が出力され、不規則波の入力数があるときには同様にその次の２行の不規則波のファイル名や作用する節点数、不規則波が作用する方向込みの節点番号などが出力される。そして、最終行は計算開始からの通算時間数である。【０１４８】
以上説明した第１実施例の空気と水と土骨格の連成計算プログラム３０では、土の種類を複数入力して要素に設定することにより、室内で人工的に作成するような練り返し土だけでなく、自然に堆積し、構造が発達した過圧密な地盤あるいは供試体を対象にして、動的・静的を問わない各種載荷条件に対する空気と水と土骨格の連成場での土骨格の弾塑性変形挙動を求めることができる。このため、第１実施例の空気と水と土骨格の連成計算プログラム３０では、供試体の変形挙動（分岐挙動、ポストピーク挙動を含む）や自重圧密計算、道路盛土など荷重載荷に伴う地盤の圧密沈下・側方変位挙動（「２次圧密」的な長期継続沈下挙動を含む）、掘削時の地盤変形挙動、揚水による地盤沈下問題、地盤の支持力問題、地震などによる繰返し負荷時の液状化と液状化後の沈下（揺すり込み沈下）挙動、浸透流による安定・不安定問題、砂杭等による円筒拡径を伴う締固め地盤改良問題、弾性構造物との相互作用問題（パイルドラフト敷設地盤の変形挙動）、補強土の変形挙動、マスパーミアビリティ改善による地盤改良問題、弾塑性盛土構造物の安定問題などを扱うことができる。さらに、対象範囲を飽和土から不飽和土まで拡張したことで、盛土や斜面の降雨浸透に伴う変形・破壊挙動や、豪雨と地震の複合時の盛土や斜面の変形・破壊挙動、河川堤防の耐震性評価、泥岩のスレーキング現象なども扱うことができる。即ち、土骨格の間隙に水が飽和している飽和土だけでなく、土骨格の間隙に水の他に空気が存在する不飽和土に対して、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。なお、これらの問題を扱う場合、必要に応じて図２に示した空気と水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて空気と水と土骨格の連成計算プログラムを複数回に亘って実行する必要がある。この場合、図１３のフローチャートに示すように、空気と水と土骨格の連成計算プログラムを実行した後に（ステップＳ３００）、計算を続行するか否かを判定し（ステップＳ３１０）、計算を続行する場合には計算条件や初期状態を変更して（ステップＳ３２０）、ステップＳ３００の空気と水と土骨格の連成計算プログラムの実行に戻る。一方、計算を続行しない場合には、それで計算を終了する。【０１４９】
図１４に、初期飽和度を最大飽和度に設定して、非排気・非排水三軸圧縮試験をシミュレートした際の軸差応力〜せん断ひずみ関係（図１４（ａ））と軸差応力〜平均有効応力関係（図１４（ｂ））とを、水と土骨格からなる飽和土の解析にて非排水三軸圧縮試験をシミュレートした場合と比較して示す。また、図１５に、初期飽和度を最大飽和度に設定して、排気・排水三軸圧縮試験をシミュレートした際の軸差応力〜せん断ひずみ関係（図１５（ａ））と軸差応力〜平均有効応力関係（図１５（ｂ））を、水と土骨格からなる飽和土の解析にて排水三軸圧縮試験をシミュレートした場合と比較して示す。両者において、空気と水と土骨格からなる不飽和土の解析による結果と水と土骨格からなる飽和土の解析による結果が一致していることから、本発明の空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラムは不飽和土から飽和土までを連続的に扱うことが可能であることがわかる。【０１５０】
図１６に、初期飽和度をそれぞれ変えた５つの場合に対して行なった非排気・非排水三軸圧縮試験のシミュレーション結果として、軸差応力〜せん断ひずみ関係（図１６（ａ））と軸差応力〜平均有効応力関係（図１６（ｂ））を示す。初期飽和度が小さくなるにつれて、強度と剛性が大きくなる、有効応力経路が右に進む、という不飽和土の典型的な力学挙動を表現できていることがわかる。【０１５１】
以上説明した第１実施例の空気と水と土骨格の連成計算装置２０によれば、全体Lマトリクス，全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体HWマトリクス，全体HAマトリクス，全体LBWマトリクス，全体LBAマトリクス、全体Mマトリクス，全体Kマトリクス、全体LWTマトリクス，全体LATマトリクス，全体SWAマトリクス、全体SACマトリクスを用いた全体接線剛性方程式（連立１次方程式）により変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件、間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件、間隙空気圧又は間隙空気の流量に関する空気境界条件を与えて未知の「加加速度（躍度）場」、間隙水圧場と間隙空気圧場の求解を行なうことにより、地盤等の土の幾何形状変化を考慮することができ、変形から破壊まで、あるいは破壊後の挙動を対象とした土の体積変化の正確な計量を必要とする大変形解析を行なうことができる。即ち、土骨格の間隙に水が飽和している飽和土だけでなく、土骨格の間隙に水の他に空気が存在する不飽和土に対して、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。しかも、土の種類を設定することができると共に要素毎に土の種類を設定することができるから、各種骨格構造（構造・過圧密・異方性）を有する砂から中間土、粘土までの連続的な土、あるいはピートやシラスなどの特殊土や問題土、もしくは砂や粘土，中間土が層をなす土など、広範な種類の土骨格の力学挙動を計算することができる。さらに、対象範囲を飽和土から不飽和土まで拡張したことで、盛土や斜面の降雨浸透に伴う変形・破壊挙動や、豪雨と地震の複合時の盛土や斜面の変形・破壊挙動、河川堤防の耐震性評価、泥岩のスレーキング現象なども扱うことができる。また、全体接線剛性方程式は、加速度の影響が小さいような応答が得られるときは、静的解析での全体剛性方程式と一致するから、静的・動的を問わないで、土骨格の力学挙動を計算することができる。すなわち、地震時の液状化後の圧密変形挙動だけでなく、圧密変形中の地震時変形挙動などのように、動的な解析に連続して静的な解析を行なったり、逆に静的な解析に連続して動的な解析を行なうことができる。【０１５２】
第１実施例では、平均化骨格応力（Bishopの有効応力式においてχ=swとした式）を用いたが、χ=swでなくても構わないし、他の有効応力式を用いても構わない。【０１５３】
第１実施例では、間隙水を非圧縮性流体、間隙空気を圧縮性流体とした場合を示しているが、間隙水を圧縮性流体として扱っても構わない。【０１５４】
本第１実施例では、温度変化は考慮していないが、温度変化を考慮しても構わない。【０１５５】
本第１実施例では、三相（固相、液相、気相）間の相変化は考慮していないが、各相間の相変化を考慮しても構わない。【０１５６】
第１実施例の空気と水と土骨格の連成計算装置２０では、土の種類を設定すると共に要素毎に土の種類を設定するものとしたが、土の種類を限定するものとしてもよいし、要素毎ではなく要素のブロックごとに土の種類を設定するものとしても構わない。【０１５７】
第１実施例の空気と水と土骨格の連成計算装置２０では、全体Lマトリクス，全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体HWマトリクス，全体HAマトリクス，全体LBWマトリクス，全体LBAマトリクス、全体Mマトリクス，全体Kマトリクス、全体LWTマトリクス，全体LATマトリクス，全体SWAマトリクス、全体SACマトリクスを各要素マトリクスとして厳密に計算するものとしたが、近似的に計算するものとしても構わない。【０１５８】
第１実施例の空気と水と土骨格の連成計算装置２０では、全体Lマトリクス，全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体ＨWマトリクス，全体HAマトリクス，全体LBWマトリクス，全体LBAマトリクス、全体Mマトリクス，全体Kマトリクス、全体LWTマトリクス，全体LATマトリクス，全体SWAマトリクス、全体SACマトリクスを計算し、これら複数種類のマトリクスのすべてを用いて全体接線剛性方程式（連立１次方程式）を連立するものとしたが、全体Lマトリクス，全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体HWマトリクス，全体HAマトリクス，全体LBWマトリクス，全体LBAマトリクス、全体Mマトリクス，全体Kマトリクス、全体LWTマトリクス，全体LATマトリクス，全体SWAマトリクス、全体SACマトリクスの一部を用いて全体接線剛性方程式（連立１次方程式）を連立するものとしてもよい。この場合、これら複数種類のマトリクス以外の他の要素を用いるものとしても構わない。【０１５９】
第１実施例の空気と水と土骨格の連成計算装置２０では、有限要素法を用いて土骨格における変形解析を行なうものとしたが、有限要素法を用いずに土骨格における変形解析を行なうものとしても差し支えない。【０１６０】
第１実施例の空気と水と土骨格の連成計算装置２０では、計算条件の設定の入力としては、図５および図６の具体例に示すように、計算条件として「平面ひずみ」であるか「軸対称」であるか「３次元」であるかの設定、変形理論として「微小変形理論」を選択するか「有限変形理論」を選択するかの設定、有限変形理論時に使用する客観性のある有効応力速度として「有効Cauchy応力のJaumann速度」であるか「Green-Nagdhi速度」であるかの設定、収束の有無の設定、構成式の選択として「下負荷面・上負荷面・回転硬化付きオリジナルカムクレイモデル」を選択するか「下負荷面・上負荷面・回転硬化付き修正カムクレイモデル」を選択するかの設定、状態量出力データの方法として「Gauss点の要素平均」であるか「ある特定のGauss点の値」であるかの設定、出力データの形式として「Text」であるか「Binary」であるかの設定、自重考慮の有無の設定、現在のプログラム使用番号、計算ステップ数、１計算ステップ当たりの計算時間DT(sec/step)、盛土・掘削時の荷重分割数、分布荷重作用時の分布荷重分割数、土の材料定数と初期値入力ファイル、メッシュデータ入力ファイル、次のプログラム時に用いる計算初期条件用出力ファイル指定を入力するものとしたが、計算条件の設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。【０１６１】
第１実施例の空気と水と土骨格の連成計算装置２０では、土の水分特性式としてvan Genuchtenの式、透水係数・透気係数のモデルとしてMualemモデルを用いる場合を考えているが、その他のヒステリシス、間隙比変化、温度依存性などを考慮したり、体積含水率や含水比などを用いる土の水分特性に関する式や、その他の透水係数・透気係数モデルを用いて行なうものとしても良い。【０１６２】
第１実施例の空気と水と土骨格の連成計算装置２０では、地盤を構成する各要素に対する粘土・中間土・砂などの土の設定として、図７の具体例として示したように、弾塑性体である土の種類の数や弾性体の種類の数、材料番号、過圧密の発展則に用いる塑性変形の進展尺度の指定、構造の発展則に用いる塑性変形の進展尺度の指定、練返し土の等方正規圧密線の切片N、限界状態定数M、圧縮指数λ、膨潤指数κ、ポアソン比ν、飽和透水係数ksw、土粒子の比重Gs、初期側圧係数K0、初期過圧密比1/R0、初期構造の程度1/R*0、正規圧密土化指数m、U*の形状、構造劣化指数a、構造劣化指数b、構造劣化指数c、構造劣化指数cs、回転硬化指数br、回転硬化限界面mb、初期異方性の程度、初期間隙比、乾燥透気係数kda、水分特性曲線のパラメータswmax，swmin，αv，mv，nv、透水係数・透気係数モデルのパラメータξ，γ、空気の温度、地盤の初期上載圧、セル圧(引張：正)、空気圧、サクション、空気の気体定数、単位系の設定（時間、力、長さ）、粘性境界の密度、圧縮波速度、せん断波速度などを入力するものとしたが、土の設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。【０１６３】
第１実施例の空気と水と土骨格の連成計算装置２０では、地盤立体構造モデルの設定としては、図８の具体例に示すように、第１行目の地盤立体構造モデルの要素数、節点数、分布荷重速度作用要素数、１要素中の節点数、計算ステップ数、弾性体の要素数、第２行目から第１０行目までの節点番号、ｘ方向の境界条件の種別、ｘ方向の境界条件の値、ｙ方向の境界条件の種別、ｙ方向の境界条件の値、第１１行目から第１４行目までの要素番号、要素番号の右から順に各要素における局所節点番号に対応する全体における節点番号、第１５行目から第２３行目までの節点番号、ｘ方向の座標、ｙ方向の座標、第２４行目および第２５行目の反力の出力を必要とする節点数、その節点番号、第２６行目および第２７行目の分布荷重速度作用要素、作用辺（作用節点間）、その左側節点の分布荷重速度のｘ方向の値、その左側節点の分布荷重速度のｙ方向の値、その右側節点の分布荷重速度のｘ方向の値、その右側節点の分布荷重速度のｙ方向の値、第２８行目から第３１行目までの要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号（水理境界の設定）、第３２行目から第３５行目までの要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号（空気境界の設定）、第３６行目の２節点間の長さ不変の条件を与えて計算を行う場合の条件数（とその条件の設定）、第３７行目および第３８行目の３節点間の角度不変の条件数とその条件の設定、第３９行目および第４０行目の前２節点の相対位置ベクトルと後２節点で指定された２節点の相対速度の方向が不変の条件数とその条件の設定、第４１行目の等変位（速度）の条件数（とその条件の設定）、第４２行目の土の種類、第４３行目のモデルに入力される規則波の入力数、モデルに入力される不規則波の入力数、第４４行目のモデルに入力される規則波振動の種類、作用する節点数、振幅、周期、第４５行目の規則波振動が作用する節点番号、第４５行目の下に新たな行を設けて指定する不規則波のファイル、などを入力するものとしたが、地盤立体構造モデルの設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。【０１６４】
第１実施例の空気と水と土骨格の連成計算装置２０では、空気と水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて空気と水と土骨格の連成計算プログラムを複数回に亘って実行することができるものとしたが、空気と水と土骨格の連成計算プログラムを複数回に亘って実行する回数を制限するものとしてもよいし、空気と水と土骨格の連成計算プログラムを一回のみ実行するものとしても差し支えない。【０１６５】
第１実施例の空気と水と土骨格の連成計算装置２０では、要素Ｋマトリクスを、式（５１）に示すように、特定された土材料（砂や粘土など）の力学挙動を与える接線剛性を示す項（右辺第１項）と、刻々の土の形状変化による接線剛性への寄与を表わす項（右辺第２項）と、衝撃荷重など土骨格に加速度が顕著に働く条件下での接線剛性への寄与を示す項（右辺第３項）と、物体力作用条件下で生じる接線剛性への寄与を示す項（右辺第４項）とにより計算するものとしたが、解析条件によってはこれらの一部の項を加味せずに要素Ｋマトリクスを計算するものとしてもよい。【０１６６】
第１実施例の空気と水と土骨格の連成計算装置２０では、２次元平面ひずみ条件や軸対称条件、３次元条件のいずれについても解析できるものとしたが、２次元平面ひずみ条件だけしか解析できないものや２次元平面ひずみ条件と軸対称条件だけしか解析できないものとしても構わない。【０１６７】
第１実施例の空気と水と土骨格の連成計算装置２０では、供試体の変形挙動（分岐挙動、ポストピーク挙動を含む）や自重圧密計算、道路盛土など荷重載荷に伴う地盤の圧密沈下・側方変位挙動（「２次圧密」的な長期継続沈下挙動を含む）、掘削時の地盤変形挙動、揚水による地盤沈下問題、地盤の支持力問題、地震などによる繰返し負荷時の液状化と液状化後の沈下（揺すり込み沈下）挙動、浸透流による安定・不安定問題、砂杭等による円筒拡径を伴う締固め地盤改良問題、弾性構造物との相互作用問題（パイルドラフト敷設地盤の変形挙動）、補強土の変形挙動、マスパーミアビリティ改善による地盤改良問題、弾塑性盛土構造物の安定問題、さらには盛土や斜面の降雨浸透に伴う変形・破壊挙動や、豪雨と地震の複合時の盛土や斜面の変形・破壊挙動、河川堤防の耐震性評価、泥岩のスレーキング現象などを扱うものとしたが、これらのすべてを扱うものとせず、これらの一部だけを扱うものとしたり、これらのうちの一つだけを扱うものとしても構わない。【０１６８】
第１実施例では、空気と水と土骨格の連成計算プログラム３０がインストールされた空気と水と土骨格の連成計算装置２０として説明したが、空気と水と土骨格の連成計算プログラム３０の形態や空気と水と土骨格の連成計算プログラム３０のうちの計算モジュール３４の形態としても差し支えない。 
【実施例２】
【０１６９】
次に、本発明の第２実施例としての空気と水と土骨格の連成計算装置２２０について説明する。第１実施例では、間隙水を非圧縮として定式化したが、第２実施例では、間隙水の圧縮性を考慮すると共に飽和度についてもより詳細に考慮するものとした。第２実施例の空気と水と土骨格の連成計算装置２２０は、図１に示す第１実施例の空気と水と土骨格の連成計算装置２０と同一のハード構成をしており、一般的な汎用コンピュータ２２２にアプリケーションソフトウエアとしての空気と水と土骨格の連成計算プログラム２３０がインストールされたものとして構成されている。第２実施例の空気と水と土骨格の連成計算プログラム２３０は、第１実施例の空気と水と土骨格の連成計算プログラム３０と同様に、計算条件の設定や土の性状の設定、地盤立体構造モデルの設定などを行なうためにデータを入力するデータ入力モジュール２３２と、入力されたデータを用いて不飽和土の力学的挙動を弾塑性構成式を用いて計算する計算モジュール２３４と、計算した結果を出力する結果出力モジュール２３６と、から構成されている。なお、第２実施例の空気と水と土骨格の連成計算装置２２０は、符号が２００だけプラスされているだけで図１と同一であるため、その図示は省略した。【０１７０】
図１７は、第２実施例の空気と水と土骨格の連成計算プログラム２３０における処理の一例を示すフローチャートである。第２実施例でも第１実施例と同様に、空気と水と土骨格の連成計算プログラム２３０の処理を説明するに当たり、説明を簡易にするために、解析対象としての地盤を２行２列のマトリクスの４要素（図３、図４参照）とするときを具体例として用いる。要素番号や節点番号、節点間番号については同様である。【０１７１】
第２実施例の空気と水と土骨格の連成計算プログラム２３０が実行されると、まず、計算条件の設定の入力処理が行なわれる（ステップＳ１１００）。ここで、計算条件の設定の入力としては、第１実施例と同様に、差分法に基づいた時間数値積分方法の選択、時刻の更新回数（計算ステップ数），２次元平面ひずみや軸対称，３次元などの計算条件を入力する。計算条件の設定を入力する具体例としてのファイルの一例として、第１実施例と同様に図５、図６を用いる。図５、図６に示す計算条件の設定の詳細については上述した。【０１７２】
続いて、地盤の各要素に対する粘土・中間土・砂などの土の設定を入力する（ステップＳ１１１０）。この土の設定では、第１実施例と同様に、材料定数や初期状態が入力される。土の設定を入力する具体例としてのファイルの一例として、第１実施例と同様に図７を用いる。図７に示す計算条件の設定の詳細については上述した。【０１７３】
そして、地盤立体構造モデルの設定を入力する（ステップＳ１１２０）。この地盤立体構造モデルの設定では、地盤モデルの要素に対する条件などを入力する。ここでも、第１実施例と同様に、図３および図４に示す具体例としての地盤立体構造モデルの設定を入力するファイルの一例として図８のものを用いる。図８に示す地盤立体構造モデルの設定の詳細については上述した。【０１７４】
こうしてデータの入力が終了すると、入力したデータを用いて計算処理を行なう（ステップＳ１１３０〜Ｓ１２７０）。計算処理は、第２実施例では、時間積分の方法としてWilsonのθ法を用いて行なったが、Newmarkのβ法などを用いて行なうものとしてもよい。また、第２実施例では、空間的離散化に対して、有限要素法を用いて間隙水圧、間隙空気圧および飽和度を要素中心に割り当てる田村流・Christian流に倣って行ったが、節点に間隙水圧を割り当てるSandhu流に基づいて行う方法やあるいはメッシュフリー法などの方法を用いて行ってもよい。以下、Wilsonのθ法と田村流に基づく有限要素法を用いた処理について説明する。まず、時刻ｔ＝ｔ＋θΔｔにおける「全体位置ベクトル」、「全体速度関連増分ベクトル」ならびに「全体加速度関連増分ベクトル」と、各種状態（有効応力と間隙水圧・間隙空気圧・飽和度、構造・過圧密・異方性等）を次式（２０１）〜（２０５）を用いて予測する（ステップＳ１１３０）。この予測計算処理は、短時間に収束させるために行なうものであり、必ずしも必要な処理ではない。予測計算処理を行わない場合は、時刻ｔ＝ｔで最終的に決められた値、ただし、初期時刻の場合は初期値、を時刻ｔ＝ｔ＋θΔｔにおける値として用いる。ここで、式（２０１）の左辺や右辺の各項のように全節点のあるベクトル量をまとめて列ベクトル表現するときには、式（２０６）に示すように成分表示される。また、式（２０４）の左辺や右辺の各項のように全要素のあるスカラー量をまとめて列ベクトル表現するときには、式（２０７）に示すように成分表示される。【０１７５】
【数７４】
2014002977.tif 000076
【０１７６】
【数７５】
2014002977.tif 000077
【０１７７】
【数７６】
2014002977.tif 000078
【０１７８】
次に、節点の変位速度を土骨格の時間的体積変化率に変換するマトリクスとして、次式（２０８）により全体Lマトリクス（全体体積変化率マトリクス）を計算する（ステップＳ１１４０）。全体Lマトリクスの各要素は式（２０９）により示される。全体LマトリクスはＮＥ行（ＮＰ×ｒ）列のマトリクス（ｒ：一節点あたりの成分の数）で、地盤立体構造モデルの設定時に指定した全体節点番号と局所節点番号の対応に従うとき、式（２０９）の要素ｉにおける要素Lマトリクスにおいて、ｍ番目の節点が全体の節点番号と一致しない列成分に対しては、数値ゼロを代入する。要素Lマトリクスを構成するBvマトリクスは、２次元平面ひずみ条件の場合は式（２１０）により、軸対称条件の場合は式（２１１）により、３次元条件の場合には式（２１２）により表わされる。【０１７９】
【数７７】
2014002977.tif 000079
【０１８０】
【数７８】
2014002977.tif 000080
【０１８１】
式（２０９）などの積分はいずれも、式（２１３）のように、全体座標系における解析対象領域を局所座標系に変換した後、Gaussの数値積分法を用いて行われる。ここに、式（２１３）は３次元条件の場合における体積積分を示す。なお、２次元平面ひずみ条件の場合は解析対象面に垂直な方向に単位厚さを考慮して、軸対称条件の場合は軸の周りの回転対称性を考慮して積分が行われる。【０１８２】
【数７９】
2014002977.tif 000081
【０１８３】
続いて、土の透水性を表す全体HWマトリクス（全体透水性マトリクス）を式（２１４）により、土の透気性を表す全体HAマトリクス（全体透気性マトリクス）を式（２１５）により計算する（ステップＳ１１５０）。全体HWマトリクスと全体HAマトリクスはＮＥ行ＮＥ列のマトリクスで、地盤立体構造モデルの設定時に指定した当該要素の節点間番号と全体要素番号の対応に従うとき、次式（２１４）と次式（２１５）で示す要素ｉの要素HWマトリクスと要素HAマトリクスにおいてｍ番目の節点間番号が全体の要素番号と一致しない全体HWマトリクスと全体HAマトリクスの成分に対しては、数値ゼロを代入する。なお、式（２１６）と式（２１７）の要素ｉにおける要素HWマトリクスと要素HAマトリクスにおける最後の成分は、要素ｉ番目（すなわち、ｉ行目）のｉ番目の列成分に代入する。式（２１６）の右辺の各成分は、式（２１８）により示される。式（２１７）の右辺の各成分は、式（２１９）により示される。ここで、排水・非排水等の水理境界条件と排気・非排気等の空気境界条件は、それぞれ式（２１８）と式（２１９）の透水係数と透気係数を変化させることで与える。例えば、非排水条件の場合は、この条件が設定されている節点間に対応する透水係数の値をゼロとする。非排気条件の場合は、この条件が設定されている節点間に対応する透気係数の値をゼロとする。各記号は、第１実施例と同様に、２次元平面ひずみ条件の場合は図９を、軸対称条件の場合は図１０を、３次元条件の場合は図１１を参照することになる。【０１８４】
【数８０】
2014002977.tif 000082
【０１８５】
【数８１】
2014002977.tif 000083
【０１８６】
次に、全体LBWマトリクス（全体水流量変換マトリクス）を式（２２０）により、全体LBAマトリクス（全体空気流量変換マトリクス）を式（２２５）により計算する（ステップS１１６０）。全体LBWマトリクスの各要素は式（２２１）により、全体LBAマトリクスの各要素は式（２２６）により示される。全体LBWマトリクスまたは全体LBAマトリクスはＮＥ行（ＮＰ×ｒ）列のマトリクスで、地盤立体構造モデルの設定時に指定した全体節点番号と局所節点番号の対応に従うとき、式（２２１）または式（２２６）の要素ｉにおける要素LBWマトリクスまたは要素LBAマトリクスにおいて、ｎ番目の節点が全体の節点番号と一致しない列成分に対しては、数値ゼロを代入する。要素LBWマトリクスまたは要素LBAマトリクスを構成するLBマトリクスは式（２２２）により表される。要素LBWマトリクスの各成分は式（２２３）により表され、βjwiの値を式（２２４）に示す。要素LBAマトリクスの各成分は式（２２７）により表され、βjaiの値を式（２２８）に示す。ここで、排水・非排水等の水理境界条件と排気・非排気等の空気境界条件は、それぞれ式（２２４）と式（２２８）の透水係数と透気係数を変化させることで与える。例えば、非排水条件の場合は、この条件が設定されている節点間に対応する透水係数の値をゼロとする。非排気条件の場合は、この条件が設定されている節点間に対応する透気係数の値をゼロとする。各記号は、第１実施例と同様に、２次元平面ひずみ条件の場合は図９を、軸対称条件の場合は図１０を、３次元条件の場合は図１１を参照することになる。【０１８７】
【数８２】
2014002977.tif 000084
【０１８８】
【数８３】
2014002977.tif 000085
【０１８９】
そして、次式（２２９）により時刻ｔ＝ｔ＋θΔｔにおける全体荷重増分ベクトルを、式（２４４）より間隙水の全体流量増分ベクトルを、さらに式（２４６）より間隙空気の全体流量増分ベクトルの計算を行なう（ステップＳ１１７０）。式（２２９）の右辺第１項は公称応力速度ベクトルに関する項で式（２３０）から構成され、右辺第２項はGreen-Naghdi速度を用いる場合に現れる項で式（２３１）から構成される。ここで、式（２３０）の右辺のNマトリクスと式（２３１）の右辺のBマトリクスは、２次元平面ひずみ条件については式（２３２）および式（２３３）で表わされ、軸対称条件については式（２３４）および式（２３５）で表わされ、３次元条件について式（２３６）および式（２３７）で表わされる。また、式（２３０）の公称応力速度ベクトルに関する項は、式（２３８）に示すように、右辺第１項に示す境界を通じて外部から与えられる表面力速度ベクトルと右辺第２項に示すこの境界の幾何形状変化に伴って現れる表面力ベクトルの増分の和からなる。このとき、右辺第１項の時刻ｔ＝ｔ＋θΔｔにおける表面力速度ベクトルは、Wilsonのθ法に従う式（２３９）により、時刻ｔ＝ｔ＋Δｔにおいて境界条件として設定された表面力速度ベクトルを換算して与える。式（２３８）右辺第２項に含まれる応力又は応力速度に関する力学的境界に立てた単位外向き法線ベクトルからなるnマトリクスは、２次元平面ひずみ条件または軸対称条件の場合は次式（２４０）で与えられ、３次元条件の場合は式（２４１）で与えられる。また、式（２３１）の右辺第１項の有効応力成分を含む列ベクトルは、２次元平面ひずみ条件または軸対称条件の場合は式（２４２）で与えられ、３次元条件の場合は式（２４３）で与えられる。【０１９０】
【数８４】
2014002977.tif 000086
【０１９１】
【数８５】
2014002977.tif 000087
【０１９２】
【数８６】
2014002977.tif 000088
【０１９３】
【数８７】
2014002977.tif 000089
【０１９４】
【数８８】
2014002977.tif 000090
【０１９５】
【数８９】
2014002977.tif 000091
【０１９６】
【数９０】
2014002977.tif 000092
【０１９７】
【数９１】
2014002977.tif 000093
【０１９８】
また、式（２４４）に示す時刻ｔ＝ｔ＋θΔｔにおける間隙水の全体流量増分ベクトルは、式（２４５）で与える全体流量速度ベクトルにθΔｔを乗じたものとして与える。全体流量速度ベクトルの成分となっているζeimの値を式（２４６）に示す。ここで、排水・非排水等の水理境界条件は、式（２４６）の透水係数を変化させることで与える。例えば、非排水条件の場合は、この条件が設定されている節点間に対応する透水係数の値をゼロとする。各記号は、２次元平面ひずみ条件の場合は図９を、軸対称条件の場合は図１０を、３次元条件の場合は図１１を参照することになる。さらに、式（２４７）に示す時刻ｔ＝ｔ＋θΔｔにおける間隙空気の全体流量増分ベクトルは、式（２４８）で与える間隙空気の全体流量速度ベクトルにθΔｔを乗じたものとして与える。【０１９９】
【数９２】
2014002977.tif 000094
【０２００】
【数９３】
2014002977.tif 000095
【０２０１】
次に、全体Mマトリクス（全体質量マトリクス）を式（２４９）により計算する（ステップＳ１１８０）。ここで、式（２４９）右辺の要素Mマトリクスは、式（２５０）で表される。なお、全体Mマトリクスは（ＮＰ×ｒ）行（ＮＰ×ｒ）列のマトリクスとなる。【０２０２】
【数９４】
2014002977.tif 000096
【０２０３】
そして、全体Kマトリクス（土骨格の全体接線剛性マトリクス）を式（２５１）により計算する（ステップＳ１１９０）。ここで、全体Kマトリクスは（ＮＰ×ｒ）行（ＮＰ×ｒ）列のマトリクスとなる。全体Kマトリクスを構成する要素Ｋマトリクスは式（２５２）で示される。この式（２５２）中、右辺第１項は特定された土材料（砂や粘土など）の力学挙動を与える接線剛性を示し、右辺第２項は刻々の土の形状変化による接線剛性への寄与を表わし、右辺第３項は衝撃荷重など土骨格に加速度が顕著に働く条件下での接線剛性への寄与を示し、右辺第４項は物体力作用条件下で生じる接線剛性への寄与を示す。材料定数の中で、特にこれらの「構造低位化指数」、「正規圧密土化指数」および「回転硬化限界定数」に対して相対的な差異を与えることによって、砂から中間土、粘土までを連続的に取り扱うことができる。式（２５２）の右辺第１項に含まれるDepマトリクス（弾塑性マトリクス）は、弾塑性状態（負荷状態）では式（２５３）で与えられ、さらに式（２５３）は２次元平面ひずみ条件で式（２５４）と式（２５５）より与えられ、軸対称条件で式（２５６）と式（２５７）より与えられ、３次元条件で式（２５８）と式（２５９）により与えられる。また、弾性状態では式（２６０）で与えられる。式（２５４）〜式（２５９）中に現れる記号は式（２６１）より与えられ、式（２６１）中に現れる記号は式（２６２）を一個ないしは複数個利用して与えられる。【０２０４】
【数９５】
2014002977.tif 000097
【０２０５】
【数９６】
2014002977.tif 000098
【０２０６】
【数９７】
2014002977.tif 000099
【０２０７】
【数９８】
2014002977.tif 000100
【０２０８】
【数９９】
2014002977.tif 000101
【０２０９】
【数１００】
2014002977.tif 000102
【０２１０】
【数１０１】
2014002977.tif 000103
【０２１１】
【数１０２】
2014002977.tif 000104
【０２１２】
p'〜q応力空間での塑性圧縮と塑性膨張の敷居線の傾きと硬化と軟化の敷居線の傾きに対する式（２６３）と式（２６４）の表記は、土の骨格構造のうちで、構造R*（０＜R*≦１、R*が０に近いほど構造は高位であることを示す）の低位化の仕方と過圧密の程度Ｒ（０＜R≦１、Rが０に近いほど過圧密であることを示す）の解消の仕方を与える発展則において、次式（２６５）および（２６６）のように、塑性測度として塑性ストレッチングのノルムを用いた場合を示している。土の種類に応じて塑性ストレッチングのせん断成分や有効応力テンソルと塑性ストレッチングの内積で表す塑性仕事率などで与えることもできる。なお、微小変形解析を用いる場合、塑性ストレッチングは塑性ひずみ速度で表される。また、式（２６５）に現れる「Ｕ*」および式（２６６）に現れる「Ｕ」はここでは非負の関数とし、例えば、次式（２６７）と（２６８）で与える。また、異方性を表す回転硬化変数テンソルの発展則については、次式（２６９）を用いている。【０２１３】
【数１０３】
2014002977.tif 000105
【０２１４】
【数１０４】
2014002977.tif 000106
【０２１５】
式（２５２）の右辺第２項を構成するＮ’マトリクスとT1マトリクスは、２次元平面ひずみ条件では式（２７０）と式（２７１）より与えられ、軸対称条件では式（２７２）と式（２７３）より与えられ、３次元条件では式（２７４）と式（２７５）により与えられる。【０２１６】
【数１０５】
2014002977.tif 000107
【０２１７】
【数１０６】
2014002977.tif 000108
【０２１８】
【数１０７】
2014002977.tif 000109
【０２１９】
次に、全体LWTマトリクス（全体間隙水圧分荷重変化率変換マトリクス）、全体LATマトリクス（全体間隙空気圧分荷重変化率変換マトリクス）および全体LSTマトリクス（全体飽和度分荷重変化率変換マトリクス）を式（２７６）、式（２７８）および式（２８０）により計算する（ステップＳ１２００）。ここで、式（２７６）右辺の要素LWTマトリクスは、式（２７７）で表され、式（２７８）右辺の要素LATマトリクスは、式（２７９）で表され、式（２８０）右辺の要素LSTマトリクスは、式（２８１）で表される。なお、全体LWTマトリクス、全体LATマトリクスおよび全体LSTマトリクスは（ＮＰ×ｒ）行ＮＥ列マトリクスとなる。【０２２０】
【数１０８】
2014002977.tif 000110
【０２２１】
【数１０９】
2014002977.tif 000111
【０２２２】
【数１１０】
2014002977.tif 000112
【０２２３】
全体SWCマトリクス（全体水圧縮性マトリクス）、全体SACマトリクス（全体空気圧縮性マトリクス）および全体SSWマトリクス（全体飽和度変化率マトリクス）をそれぞれ式（２８２）、式（２８３）および式（２８４）より求める。また、全体SCマトリクス（全体サクション寄与分水分特性マトリクス）と全体SWマトリクス（全体飽和度寄与分水分特性マトリクス）を式（２８５）、式（２８６）により求める。ここではvan Genuchtenの式を用いる場合を例として示している。上記のマトリクスは全て、ＮＥ行ＮＥ列となる。（ステップＳ１２１０）【０２２４】
【数１１１】
2014002977.tif 000113
【０２２５】
【数１１２】
2014002977.tif 000114
【０２２６】
全体透水修正原Lマトリクス、全体透気修正原Lマトリクス、全体飽和度修正原Lマトリクス、全体空気間隙比修正原Lマトリクス、全体間隙水修正原Lマトリクスおよび全体間隙空気修正原Lマトリクスをそれぞれ式（２８７）、式（２８８）、式（２８９）、式（２９０）、式（２９１）および式（２９２）より求める（ステップＳ１２２０）。全体Lマトリクスと同様にして作られるＮＥ行（ＮＰ×ｒ）列のマトリクスである。ただし、全体透水修正原Lマトリクスと全体透気修正原Lマトリクスは、全体接線剛性方程式（連立１次方程式）を構築する上では、それぞれ全体間隙水修正原Lマトリクスと全体間隙空気修正原Ｌマトリクスに包含されるため、全体透水修正原Lマトリクスと全体透気修正原Lマトリクスは実質計算する必要はない。【０２２７】
【数１１３】
2014002977.tif 000115
【０２２８】
【数１１４】
2014002977.tif 000116
【０２２９】
次式（２９３）に示す全体接線剛性方程式（連立１次方程式）より、変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件並びに間隙空気圧又は間隙空気の流量に関する空気境界条件を与えて未知の「加加速度（躍度）場」、間隙水圧場、間隙空気圧場と飽和度場の求解を行なう（ステップＳ１２３０）。この式（２９７）は、速度型運動方程式の弱形式の有限要素離散化と慣性項の影響を考慮した空気と水と土骨格連成式のモデル化によって最終的に得られる連立２階常微分方程式（式（２９４））において、時間微分項を陰解法（差分法）に基づいて処理したもので、変形場に対してはWilsonのθ法を用いて、間隙水圧、間隙空気圧および飽和度に対しては台形公式を利用した場合の（最終的に）解くべき連立１次方程式である。全体荷重増分ベクトル、間隙水の全体流量増分ベクトル及び間隙空気の全体流量増分ベクトル並びに全体Mマトリクス及び全体Kマトリクス並びに全体HWマトリクス及び全体HAマトリクス並びに全体LBWマトリクス及び全体LBAマトリクス並びに全体LWTマトリクス、全体LATマトリクス及び全体LSTマトリクス並びに全体SWCマトリクス、全体SACマトリクス及び全体SSWマトリクス並びに全体SCマトリクス及び全体SWマトリクス並びに全体LW2マトリクス及び全体LA2マトリクス並びに全体LθnWマトリクス及び全体LθnAマトリクスは時刻ｔ＝ｔと時刻ｔ＝ｔ＋θΔｔとの間の時間ステップ内で繰り返し計算を行ないながら更新する。なお、式（２９４）では、運動方程式に対しさらに土骨格から見た物質時間微分を施しているため、説明変数である速度ベクトルの時間に関する２階微分の項、すなわち、「加加速度（躍度とも称する）」項を有している点がもっとも特徴的である。これは、土骨格の構成式を速度型で与えるためだけでなく、地盤等の土の幾何形状変化を考慮するためで、これによって、変形から破壊まで、あるいは破壊後の挙動を対象にした土の体積変化の正確な計量を必要とする大変形解析を可能にするためである。またこの加加速度（躍度）導入のため、Wilsonのθ法の基礎となる線形加速度法の考えに基づいて、変形に対して式（２９５）〜式（２９７）を誘導している。なお、加速度の影響が小さいような応答が得られるときは、式（２９４）は、静的解析で得られる連立１階常微分方程式と一致する。このことは、動的および静的解析の使い分けを必要としないことを意味する。また、例えば、二節点の間に距離が不変あるいは三節点がなす角度が不変の条件などの束縛条件が節点間に課せられる場合は、Lagrangeの未定乗数法を適用して式（２９４）にその束縛条件を考慮する。時刻ｔ＝ｔ＋Δｔにおいて設定された変位、変位速度又は加速度に関する幾何的境界条件については、位置ベクトルの差すなわち変位ベクトルで与えられた場合は式（２９８）により、速度ベクトルで与えられた場合は式（２９９）により、加速度ベクトルで与えられた場合は式（３００）により、時刻ｔ＝ｔ＋θΔｔでの加加速度（躍度）ベクトルに換算して与える。【０２３０】
【数１１５】
2014002977.tif 000117
【０２３１】
【数１１６】
2014002977.tif 000118
【０２３２】
【数１１７】
2014002977.tif 000119
【０２３３】
【数１１８】
2014002977.tif 000120
【０２３４】
続いて、時刻ｔ＝ｔ＋θΔｔにおける加速度場、速度場、座標、各種状態量を計算すると共に各Gauss点の負荷状態を判定する（ステップＳ１２４０）。時刻ｔ＝ｔ＋θΔｔにおける全体加速度関連増分ベクトル、全体速度関連増分ベクトルおよび座標ベクトルは、式（２９３）で得られた時刻ｔ＝ｔ＋θΔｔにおける全体加加速度（躍度）関連増分ベクトルと時刻ｔ＝ｔにおける幾何的諸量から、式（２９５）〜式（２９７）より得られる次式（３０１）〜式（３０３）により算出することができる。時刻ｔ＝ｔ＋θΔｔにおける有効応力、構造の程度R*、過圧密の程度R、異方性などの各状態量Aについては、式（３０４）により算出することができる。例えば、時刻ｔ＝ｔ＋θΔｔにおける要素内の各Gauss点での有効応力速度の場合は、式（３０２）より得られた全体速度関連増分ベクトルのうち、当該要素が有する節点の速度ベクトルから式（３０５）を用いて算出した後、式（３０４）を適用して有効応力を算出する。なお、構造の程度、過圧密の程度、異方性などの速度は、Gauss点が弾塑性状態にある場合は式（２６５）、式（２６６）、式（２６９）などを用いて算出し、弾性状態にある場合はそれらの速度をゼロとする。時刻ｔ＝ｔ＋θΔｔにおける「塑性体積ひずみ」は、下負荷面・上負荷面・回転硬化付き修正カムクレイモデルの場合、Gauss点が弾塑性状態にあるとき式（３０９）より算出するが、弾性状態にあった場合は、最後に弾塑性状態にあった時点での値のままとする。したがって、この弾性状態のとき、時刻ｔ＝ｔ＋θΔｔにおける過圧密の程度Rは、式（３０９）をRについて解いた式を用いて算出することができる。また、時刻ｔ＝ｔ＋θΔｔにおける間隙水圧、間隙空気圧及び飽和度は式（２９３）の全体接線剛性方程式より直接決まるが、間隙水圧、間隙空気圧及び飽和度の速度は式（３０４）を逆に用いて求める。また、土要素の各Gauss点が弾塑性状態にあるか弾性状態にあるかについては、土要素の各節点から計算されるGauss点のストレッチングと有効応力などを用いて塑性乗数の分子から得られる次式（３１０）に示す負荷基準を満たすときに、そのGauss点の状態は弾塑性状態にあると判定することができる。なお、ストレッチングの列ベクトル表現したものは式（３１１）より計算されるが、時刻ｔ＝ｔ＋θΔｔにおけるその列ベクトルの成分表示は、２次元平面ひずみ条件の場合は式（３１２）により、軸対称条件の場合は式（３１３）により、３次元条件の場合は式（３１４）により表される。【０２３５】
【数１１９】
2014002977.tif 000121
【０２３６】
【数１２０】
2014002977.tif 000122
【０２３７】
【数１２１】
2014002977.tif 000123
【０２３８】
【数１２２】
2014002977.tif 000124
【０２３９】
【数１２３】
2014002977.tif 000125
【０２４０】
そして、各要素の各Gauss点の有効応力を用いて収束判定を行なう（ステップＳ１２５０）。実施例では、全要素内の全Gauss点で次式（３１５）を充足するか否かにより収束を判定することとした。ここで、式（３１５）中の「ε」は十分に小さい正の値である。収束が判定されないときには、ステップＳ１１５０の全体HWマトリクスと全体HAマトリクスの計算に戻り、ステップＳ１１５０〜ステップＳ１２５０までの処理を繰り返す。【０２４１】
【数１２４】
2014002977.tif 000126
【０２４２】
ステップＳ１２５０で収束したと判定されると、時刻ｔ＝ｔ＋Δｔにおける加速度場、速度場、座標、各種状態量を計算すると共に各Gauss点の負荷状態を判定する（ステップＳ１２６０）。時刻ｔ＝ｔ＋Δｔの座標などは、時刻ｔ＝ｔ＋θΔｔにおける全体加加速度（躍度）関連ベクトルと時刻ｔ＝ｔにおける全体加速度関連増分ベクトル、全体速度関連増分ベクトルおよび座標ベクトルから、次式（３１７）〜式（３２０）により算出することができる。また、時刻ｔ＝ｔ＋Δｔの各状態量とその速度については、時刻ｔ＝ｔと時刻ｔ＝ｔ＋θΔｔにおける状態量から式（３２１）および式（３２２）により計算することができる。そして、時刻ｔ＝ｔ＋Δｔの各状態量と式（３１０）と式（３１１）を用いて、時刻ｔ＝ｔ＋Δｔでの負荷状態の判定を行う。【０２４３】
【数１２５】
2014002977.tif 000127
【０２４４】
【数１２６】
2014002977.tif 000128
【０２４５】
そしてステップＳ１１００で入力した計算ステップ数に至ったか否かを判定し（ステップＳ１２７０）、計算ステップ数に至っていないときには、次の計算ステップとしてステップＳ１１３０の時刻ｔ＝ｔ＋θΔｔにおける各種状態（有効応力と間隙水圧・間隙空気圧・飽和度、構造・過圧密・異方性、等）の予測処理に戻る。一方、計算ステップ数に至ったときには、計算結果を指定した出力ファイルに出力して（ステップＳ１２８０）、空気と水と土骨格の連成計算プログラム２３０を終了する。【０２４６】
出力ファイルの一例としては、第１実施例と同様に図１２となる。この図１２の詳細についても上述した。【０２４７】
以上説明した第２実施例の空気と水と土骨格の連成計算プログラム２３０では、土の種類を複数入力して要素に設定することにより、室内で人工的に作成するような練り返し土だけでなく、自然に堆積し、構造が発達した過圧密な地盤あるいは供試体を対象にして、動的・静的を問わない各種載荷条件に対する空気と水と土骨格の連成場での土骨格の弾塑性変形挙動を求めることができる。このため、実施例の空気と水と土骨格の連成計算プログラム３０では、供試体の変形挙動（分岐挙動、ポストピーク挙動を含む）や自重圧密計算、道路盛土など荷重載荷に伴う地盤の圧密沈下・側方変位挙動（「２次圧密」的な長期継続沈下挙動を含む）、掘削時の地盤変形挙動、揚水による地盤沈下問題、地盤の支持力問題、地震などによる繰返し負荷時の液状化と液状化後の沈下（揺すり込み沈下）挙動、浸透流による安定・不安定問題、砂杭等による円筒拡径を伴う締固め地盤改良問題、弾性構造物との相互作用問題（パイルドラフト敷設地盤の変形挙動）、補強土の変形挙動、マスパーミアビリティ改善による地盤改良問題、弾塑性盛土構造物の安定問題などを扱うことができる。さらに、対象範囲を飽和土から不飽和土まで拡張したことで、盛土や斜面の降雨浸透に伴う変形・破壊挙動や、豪雨と地震の複合時の盛土や斜面の変形・破壊挙動、河川堤防の耐震性評価、泥岩のスレーキング現象なども扱うことができる。即ち、土骨格の間隙に水が飽和している飽和土だけでなく、土骨格の間隙に水の他に空気が存在する不飽和土に対して、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。なお、これらの問題を扱う場合、必要に応じて図２に示した空気と水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて空気と水と土骨格の連成計算プログラムを複数回に亘って実行する必要がある。この場合、第１実施例と同様に、図１３のフローチャートに示すように、空気と水と土骨格の連成計算プログラムを実行した後に（ステップＳ３００）、計算を続行するか否かを判定し（ステップＳ３１０）、計算を続行する場合には計算条件や初期状態を変更して（ステップＳ３２０）、ステップＳ３００の空気と水と土骨格の連成計算プログラムの実行に戻る。一方、計算を続行しない場合には、それで計算を終了する。【０２４８】
図１８に、初期飽和度をそれぞれ変えた５つの場合に対して行なった非排気・非排水三軸圧縮試験のシミュレーション結果として、軸差応力〜軸ひずみ関係（図１８（ａ））と応力経路（図１８（ｂ））を示す。飽和度の違いに応じた不飽和土の力学挙動を表現できていることがわかる。【０２４９】
図１９に、ほぼ飽和状態にある粘土地盤上に、不飽和状態にある盛土を載荷したときの飽和度分布図を示す。盛土載荷による沈下により、盛土下部が飽和化している様子を表現できている。図２０に、水位上昇時の盛土内の飽和度分布の様子を示す。水位上昇時の盛土内の浸潤面生成過程を表現できている。【０２５０】
以上説明した第２実施例の空気と水と土骨格の連成計算装置２２０によれば、全体Lマトリクス、全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体HWマトリクス、全体HAマトリクス、全体LBWマトリクス、全体LBAマトリクス、全体Mマトリクス，全体Kマトリクス、全体LWTマトリクス、全体LATマトリクス、全体LSTマトリクス、全体SWCマトリクス、全体SACマトリクス、全体SSWマトリクス、全体SCマトリクス、全体SWマトリクスを用いた全体接線剛性方程式（連立１次方程式）により変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件、間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件、間隙空気圧又は間隙空気の流量に関する空気境界条件を与えて未知の「加加速度（躍度）場」、間隙水圧場、間隙空気圧場及び飽和度の求解を行なうことにより、地盤等の土の幾何形状変化を考慮することができ、変形から破壊まで、あるいは破壊後の挙動を対象とした土の体積変化の正確な計量を必要とする大変形解析を行なうことができる。即ち、土骨格の間隙に水が飽和している飽和土だけでなく、土骨格の間隙に水の他に空気が存在する不飽和土に対して、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。しかも、土の種類を設定することができると共に要素毎に土の種類を設定することができるから、各種骨格構造（構造・過圧密・異方性）を有する砂から中間土、粘土までの連続的な土、あるいはピートやシラスなどの特殊土や問題土、もしくは砂や粘土，中間土が層をなす土など、広範な種類の土骨格の力学挙動を計算することができる。さらに、対象範囲を飽和土から不飽和土まで拡張したことで、盛土や斜面の降雨浸透に伴う変形・破壊挙動や、豪雨と地震の複合時の盛土や斜面の変形・破壊挙動、河川堤防の耐震性評価、泥岩のスレーキング現象なども扱うことができる。また、全体接線剛性方程式は、加速度の影響が小さいような応答が得られるときは、静的解析での全体剛性方程式と一致するから、静的・動的を問わないで、土骨格の力学挙動を計算することができる。すなわち、地震時の液状化後の圧密変形挙動だけでなく、圧密変形中の地震時変形挙動などのように、動的な解析に連続して静的な解析を行なったり、逆に静的な解析に連続して動的な解析を行なうことができる。【０２５１】
第２実施例の空気と水と土骨格の連成計算装置２２０では、平均化骨格応力（Bishopの有効応力式においてχ=swとした式）を用いたが、χ=swでなくても構わないし、他の有効応力式を用いても構わない。【０２５２】
第２実施例の空気と水と土骨格の連成計算装置２２０では、加加速度に加えて、間隙水圧、間隙空気圧および飽和度を未知数とした連立一次方程式を用いているが、これに限定されるものではなく、例えば飽和度を未知数にしないで方程式を構築するなど、その一部を未知数とした連立一次方程式を用いても構わない。【０２５３】
第２実施例の空気と水と土骨格の連成計算装置２２０では、温度変化は考慮していないが、温度変化を考慮しても構わない。【０２５４】
第２実施例の空気と水と土骨格の連成計算装置２２０では、三相（固相、液相、気相）間の相変化は考慮していないが、各相間の相変化を考慮しても構わない。【０２５５】
第２実施例の空気と水と土骨格の連成計算装置２２０では、土の種類を設定すると共に要素毎に土の種類を設定するものとしたが、土の種類を限定するものとしてもよいし、要素毎ではなく要素のブロックごとに土の種類を設定するものとしても構わない。【０２５６】
第２実施例の空気と水と土骨格の連成計算装置２２０では、全体Lマトリクス、全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体HWマトリクス、全体HAマトリクス、全体LBWマトリクス、全体LBAマトリクス、全体Mマトリクス，全体Kマトリクス、全体LWTマトリクス、全体LATマトリクス、全体LSTマトリクス、全体SWCマトリクス、全体SACマトリクス、全体SSWマトリクス、全体SCマトリクス、全体SWマトリクスを各要素マトリクスとして厳密に計算するものとしたが、近似的に計算するものとしても構わない。【０２５７】
第２実施例の空気と水と土骨格の連成計算装置２２０では、全体Lマトリクス、全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体HWマトリクス、全体HAマトリクス、全体LBWマトリクス、全体LBAマトリクス、全体Mマトリクス，全体Kマトリクス、全体LWTマトリクス、全体LATマトリクス、全体LSTマトリクス、全体SWCマトリクス、全体SACマトリクス、全体SSWマトリクス、全体SCマトリクス、全体SWマトリクスを計算し、これら複数種類のマトリクスのすべてを用いて全体接線剛性方程式（連立１次方程式）を連立するものとしたが、全体Lマトリクス、全体LW2マトリクス、全体LA2マトリクス、全体LΘnWマトリクス、全体LΘnAマトリクス，全体HWマトリクス、全体HAマトリクス、全体LBWマトリクス、全体LBAマトリクス、全体Ｍマトリクス，全体Kマトリクス、全体LWTマトリクス、全体LATマトリクス、全体LSTマトリクス、全体SWCマトリクス、全体SACマトリクス、全体SSWマトリクス、全体SCマトリクス、全体SWマトリクスの一部を用いて全体接線剛性方程式（連立１次方程式）を連立するものとしてもよい。この場合、これら複数種類のマトリクス以外の他の要素を用いるものとしても構わない。【０２５８】
第２実施例の空気と水と土骨格の連成計算装置２２０では、有限要素法を用いて土骨格における変形解析を行なうものとしたが、有限要素法を用いずに土骨格における変形解析を行なうものとしても差し支えない。【０２５９】
第２実施例の空気と水と土骨格の連成計算装置２２０では、計算条件の設定の入力としては、図５および図６の具体例に示すように、計算条件として「平面ひずみ」であるか「軸対称」であるか「３次元」であるかの設定、変形理論として「微小変形理論」を選択するか「有限変形理論」を選択するかの設定、有限変形理論時に使用する客観性のある有効応力速度として「有効Cauchy応力のJaumann速度」であるか「Green-Nagdhi速度」であるかの設定、収束の有無の設定、構成式の選択として「下負荷面・上負荷面・回転硬化付きオリジナルカムクレイモデル」を選択するか「下負荷面・上負荷面・回転硬化付き修正カムクレイモデル」を選択するかの設定、状態量出力データの方法として「Gauss点の要素平均」であるか「ある特定のGauss点の値」であるかの設定、出力データの形式として「Text」であるか「Binary」であるかの設定、自重考慮の有無の設定、現在のプログラム使用番号、計算ステップ数、１計算ステップ当たりの計算時間DT(sec/step)、盛土・掘削時の荷重分割数、分布荷重作用時の分布荷重分割数、土の材料定数と初期値入力ファイル、メッシュデータ入力ファイル、次のプログラム時に用いる計算初期条件用出力ファイル指定を入力するものとしたが、計算条件の設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。【０２６０】
第２実施例の空気と水と土骨格の連成計算装置２２０では、土の水分特性に関する式としてvan Genuchtenの式、透水係数・透気係数のモデルとしてMualemモデルを用いる場合を考えているが、その他のヒステリシス、間隙比変化、温度依存性などを考慮したり、体積含水率や含水比などを用いる土の水分特性に関する式や、その他の透水係数・透気係数モデルを用いて行なうものとしても良い。【０２６１】
第２実施例の空気と水と土骨格の連成計算装置２２０では、地盤を構成する各要素に対する粘土・中間土・砂などの土の設定として、図１８の具体例として示したように、弾塑性体である土の種類の数や弾性体の種類の数、材料番号、過圧密の発展則に用いる塑性変形の進展尺度の指定、構造の発展則に用いる塑性変形の進展尺度の指定、練返し土の等方正規圧密線の切片N、限界状態定数M、圧縮指数λ、膨潤指数κ、ポアソン比ν、飽和透水係数ksw、土粒子の比重Gs、初期側圧係数K0、初期過圧密比1/R0、初期構造の程度1/R*0、正規圧密土化指数m、U*の形状、構造劣化指数a、構造劣化指数b、構造劣化指数c、構造劣化指数cs、回転硬化指数br、回転硬化限界面mb、初期異方性の程度、初期間隙比、乾燥透気係数kda、水分特性のパラメータswmax，swmin，αv，mv，nv、透水係数・透気係数モデルのパラメータξ，γ、空気の温度、地盤の初期上載圧、セル圧 (引張：正)、空気圧、サクション、空気の気体定数、単位系の設定（時間、力、長さ）、粘性境界の密度、圧縮波速度、せん断波速度などを入力するものとしたが、土の設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。【０２６２】
第２実施例の空気と水と土骨格の連成計算装置２２０では、地盤立体構造モデルの設定としては、図８の具体例に示すように、第１行目の地盤立体構造モデルの要素数、節点数、分布荷重速度作用要素数、１要素中の節点数、計算ステップ数、弾性体の要素数、第２行目から第１０行目までの節点番号、ｘ方向の境界条件の種別、ｘ方向の境界条件の値、ｙ方向の境界条件の種別、ｙ方向の境界条件の値、第１１行目から第１４行目までの要素番号、要素番号の右から順に各要素における局所節点番号に対応する全体における節点番号、第１５行目から第２３行目までの節点番号、ｘ方向の座標、ｙ方向の座標、第２４行目および第２５行目の反力の出力を必要とする節点数、その節点番号、第２６行目および第２７行目の分布荷重速度作用要素、作用辺（作用節点間）、その左側節点の分布荷重速度のｘ方向の値、その左側節点の分布荷重速度のｙ方向の値、その右側節点の分布荷重速度のｘ方向の値、その右側節点の分布荷重速度のｙ方向の値、第２８行目から第３１行目までの要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号（水理境界の設定）、第３２行目から第３５行目までの要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号（空気境界の設定）、第３６行目の２節点間の長さ不変の条件を与えて計算を行う場合の条件数（とその条件の設定）、第３７行目および第３８行目の３節点間の角度不変の条件数とその条件の設定、第３９行目および第４０行目の前２節点の相対位置ベクトルと後２節点で指定された２節点の相対速度の方向が不変の条件数とその条件の設定、第４１行目の等変位（速度）の条件数（とその条件の設定）、第４２行目の土の種類、第４３行目のモデルに入力される規則波の入力数、モデルに入力される不規則波の入力数、第４４行目のモデルに入力される規則波振動の種類、作用する節点数、振幅、周期、第４５行目の規則波振動が作用する節点番号、第４５行目の下に新たな行を設けて指定する不規則波のファイル、などを入力するものとしたが、地盤立体構造モデルの設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。【０２６３】
第２実施例の空気と水と土骨格の連成計算装置２２０では、空気と水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて空気と水と土骨格の連成計算プログラムを複数回に亘って実行することができるものとしたが、空気と水と土骨格の連成計算プログラムを複数回に亘って実行する回数を制限するものとしてもよいし、空気と水と土骨格の連成計算プログラムを一回のみ実行するものとしても差し支えない。【０２６４】
第２実施例の空気と水と土骨格の連成計算装置２２０では、要素Kマトリクスを、式（２５１）に示すように、特定された土材料（砂や粘土など）の力学挙動を与える接線剛性を示す項（右辺第１項）と、刻々の土の形状変化による接線剛性への寄与を表わす項（右辺第２項）と、衝撃荷重など土骨格に加速度が顕著に働く条件下での接線剛性への寄与を示す項（右辺第３項）と、物体力作用条件下で生じる接線剛性への寄与を示す項（右辺第４項）とにより計算するものとしたが、解析条件によってはこれらの一部の項を加味せずに要素Kマトリクスを計算するものとしてもよい。【０２６５】
第２実施例の空気と水と土骨格の連成計算装置２２０では、２次元平面ひずみ条件や軸対称条件、３次元条件のいずれについても解析できるものとしたが、２次元平面ひずみ条件だけしか解析できないものや２次元平面ひずみ条件と軸対称条件だけしか解析できないものとしても構わない。【０２６６】
第２実施例の空気と水と土骨格の連成計算装置２２０では、供試体の変形挙動（分岐挙動、ポストピーク挙動を含む）や自重圧密計算、道路盛土など荷重載荷に伴う地盤の圧密沈下・側方変位挙動（「２次圧密」的な長期継続沈下挙動を含む）、掘削時の地盤変形挙動、揚水による地盤沈下問題、地盤の支持力問題、地震などによる繰返し負荷時の液状化と液状化後の沈下（揺すり込み沈下）挙動、浸透流による安定・不安定問題、砂杭等による円筒拡径を伴う締固め地盤改良問題、弾性構造物との相互作用問題（パイルドラフト敷設地盤の変形挙動）、補強土の変形挙動、マスパーミアビリティ改善による地盤改良問題、弾塑性盛土構造物の安定問題、さらには盛土や斜面の降雨浸透に伴う変形・破壊挙動や、豪雨と地震の複合時の盛土や斜面の変形・破壊挙動、河川堤防の耐震性評価、泥岩のスレーキング現象などを扱うものとしたが、これらのすべてを扱うものとせず、これらの一部だけを扱うものとしたり、これらのうちの一つだけを扱うものとしても構わない。【０２６７】
第２実施例では、空気と水と土骨格の連成計算プログラム２３０がインストールされた空気と水と土骨格の連成計算装置２０として説明したが、空気と水と土骨格の連成計算プログラム２３０の形態や空気と水と土骨格の連成計算プログラム２３０のうちの計算モジュール２３４の形態としても差し支えない。【０２６８】
以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において、種々なる形態で実施し得ることは勿論である。【０２６９】
この出願は、２０１２年６月２５日に出願された日本国特許出願第２０１２−１４１８８８号を優先権主張の基礎としており、引用によりその内容の全てが本明細書に含まれる。
【産業上の利用可能性】
【０２７０】
本発明は、空気と水と土骨格の連成計算装置を利用する製造産業などに利用可能である。
【図１】
2014002977.tif 000131
【図２】
2014002977.tif 000132
【図３】
2014002977.tif 000133
【図４】
2014002977.tif 000134
【図５】
2014002977.tif 000135
【図１３】
2014002977.tif 000136
【図６】
2014002977.tif 000137
【図７】
2014002977.tif 000138
【図８】
2014002977.tif 000139
【図９】
2014002977.tif 000140
【図１０】
2014002977.tif 000141
【図１１】
2014002977.tif 000142
【図１２】
2014002977.tif 000143
【図１４】
2014002977.tif 000144
【図１５】
2014002977.tif 000145
【図１６】
2014002977.tif 000146
【図１７】
2014002977.tif 000147
【図１８】
2014002977.tif 000148
【図１９】
2014002977.tif 000149
【図２０】
2014002977.tif 000150
【国際調査報告】
2014002977.tif 000151
2014002977.tif 000152