JP4441693B2 - 水と土骨格の連成計算装置および水と土骨格の連成計算方法 - Google Patents

Description

本発明は、水と土骨格の連成計算装置および水と土骨格の連成計算方法に関し、詳しくは、土骨格の変形解析を行なう水と土骨格の連成計算を行なう装置および土骨格の変形解析を行なう水と土骨格の連成計算を行なう方法に関する。

従来、この種の土骨格解析プログラムとしては、３次元地盤変形解析や３次元圧密解析が可能なもの（例えば、非特許文献１参照）や地盤の静力学的挙動を表現した力学構成モデルを組み込んだもの（例えば、非特許文献２参照）などが提案されている。これらのプログラムは、透水、圧密、地盤内応力、せん断、土圧、斜面安定などの問題を精度よく解析することができる、とされている。
「地盤ＦＥＭ解析支援システム ＡＦＩＭＥＸ−ＧＴ」のカタログ，富士通エフ・アイ・ピー株式会社，No.0040708-2 「地盤解析汎用プログラム」のホームページ，財団法人 沿岸開発技術研究センター，http://www.ysk.nilim.go.jp/geofem/phgeo01j.html，２００５年８月１９日検索

しかしながら、上述の土骨格解析プログラムでは、土の性状としては粘土専用であったり、砂専用であったりするため、埋め立て人工地盤などのように、砂と粘土とが混在する中間土に対しては精度よく解析することができない。また、自然土の多くは、粒度が揃った砂や粘土ではなく、砂が多い粘土であったり粘土混じりの砂であったり様々であり、砂専用に開発されたプログラムによる解析結果や粘土専用に開発されたプログラムによる解析結果を適用すると、解析結果は実際の挙動と異なるものとなってしまう。さらに、上述のプログラムでは、粘土と砂とが層をなす地盤については解析することができない場合が多い。

また、砂は、ゆるい状態にあるときには、急速に非圧縮、非排水のまま繰り返し載荷されると液状化を示すが、ゆっくり排水を伴って繰り返し載荷されると締固めを示す。上述のプログラムでは、液状化を解析できるものの、締固めを解析できない結果、液状化専用のものとなってしまう。さらに、上述のプログラムでは、そのプログラムを使用するのに適した砂地盤あるいは粘土地盤であるか否かを別の指標を用いて判定した後に使用する必要がある。加えて、上述のプログラムでは、動的解析専用であったり、静的解析専用であったりするから、液状化後の揺すり込み沈下などのように動的な解析の後に連続して静的な解析を行なうことはできない。

本発明の水と土骨格の連成計算装置および水と土骨格の連成計算方法は、静的解析か動的解析かに拘わらずに土骨格の変形解析を行なうことを目的の一つとする。また、本発明の水と土骨格の連成計算装置および水と土骨格の連成計算方法は、砂から粘土に至るまでのあらゆる中間土の地盤についての解析を精度よく行なうことを目的の一つとする。さらに、本発明の水と土骨格の連成計算装置および水と土骨格の連成計算方法は、粘土や砂，中間土が層をなす地盤についても精度よく解析することを目的の一つとする。あるいは、本発明の水と土骨格の連成計算装置および水と土骨格の連成計算方法は、土の状態と載荷の状態に応じた解析を繰り返し可能にすることを目的の一つとする。加えて、本発明の水と土骨格の連成計算装置および水と土骨格の連成計算方法は、動的な解析に連続して静的な解析を可能にしたり、逆に静的な解析に連続して動的な解析を可能にすることを目的の一つとする。

本発明の水と土骨格の連成計算装置および水と土骨格の連成計算方法は、上述の目的の少なくとも一部を達成するために以下の手段を採った。

本発明の水と土骨格の連成計算装置は、土骨格の状態に関するデータと外力に関するデータとを入力するデータ入力手段と、前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧に関する連立速度型運動方程式を構築し、該連立方程式に変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて随時積分し、速度場および間隙水圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう解析手段と、該解析結果を出力する結果出力手段と、を備えることを要旨とする。

この本発明の水と土骨格の連成計算装置では、土骨格の加加速度項を含む土骨格の速度と間隙水圧に関する連立速度型運動方程式を変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて随時積分し、速度場および間隙水圧場の時刻歴応答を求めるから、加速度の影響が大きいときには動的解析による土骨格の変形解析を行なうことができ、加速度の影響が小さいときには静的解析による土骨格の変形解析を行なうことができる。したがって、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。

こうした本発明の水と土骨格の連成計算装置において、前記解析手段は、前記土骨格の状態に関するデータと前記外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と土の透水性に関する第２項目と質量に関する第３項目と土骨格の接線剛性に関する第４項目とのうちの少なくとも一つを含む複数の項目を計算し、該計算した複数の項目を用いて前記連立速度型運動方程式を構築する手段であるものとすることもできる。この場合、前記解析手段は、前記第１項目と前記第２項目と前記第３項目と前記第４項目とを用いて前記連立速度型運動方程式を構築する手段であるものとすることもできる。こうすれば、より適正に土骨格の変形解析を行なうことができる。

また、本発明の水と土骨格の連成計算装置において、前記解析手段は、有限要素法、差分法等の数値解析法を用いて前記連立速度型運動方程式を構築する手段であるものとすることもできる。こうすれば、計算を容易なものとすることができる。

さらに、本発明の水と土骨格の連成計算装置において、前記土骨格の状態に関するデータは、土骨格を構成する複数の要素における各要素の性状に関するデータおよび隣接する要素との関係に関するデータとが含まれるデータであるものとすることもできる。この場合、前記各要素の性状に関するデータは、土の性状に関するデータが含まれるデータであるものとすることもできる。こうすれば、土骨格を構成する複数の要素毎に要素の性状を設定することができると共に要素間の関係を要素間毎に設定することができる。したがって、砂から粘土に至るまでのあらゆる中間土の地盤についての土骨格における変形解析を精度よく行なうことができる。また、粘土や砂，中間土が層をなす地盤についての土骨格における変形解析を精度よく解析することができる。

あるいは、本発明の水と土骨格の連成計算装置において、前記外力に関するデータは、水と土骨格の連成系に与える荷重に関するデータが含まれるデータであるものとすることもできるし、前記外力に関するデータは、水と土骨格の連成系に与える振動に関するデータが含まれるデータであるものとすることもできる。こうすれば、荷重及び変位に対する土骨格における変形解析を行なうことができると共に振動に対する土骨格における変形解析を行なうこともできる。

また、本発明の水と土骨格の連成計算装置において、前記解析手段は、前記外力に関するデータを複数回変更して解析する手段であるものとすることもできる。この場合、前記解析手段は、前記外力に関するデータを変更するときには変更前の解析結果を初期状態の土骨格の状態に関するデータとして用いて解析する手段であるものとすることもできる。こうすれば、土の状態と載荷の状態に応じた土骨格における変形解析を繰り返し行うことができる。また、動的な解析に連続して静的な解析を行なったり、逆に静的な解析に連続して動的な解析を行なったりすることができる。

本発明の水と土骨格の連成計算方法は、土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧に関する連立速度型運動方程式を構築し、該連立方程式に変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて随時積分し、速度場および間隙水圧場の時刻歴応答を求めることにより土骨格の変形解析を行なうことを特徴とする。

この本発明の水と土骨格の連成計算方法では、土骨格の加加速度項を含む土骨格の速度と間隙水圧に関する連立速度型運動方程式を変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて随時積分し、速度場および間隙水圧場の時刻歴応答を求めるから、加速度の影響が大きいときには動的解析による土骨格の変形解析を行なうことができ、加速度の影響が小さいときには静的解析による土骨格の変形解析を行なうことができる。したがって、静的解析か動的解析かに拘わらず、土骨格における変形解析を行なうことができる。

こうした本発明の水と土骨格の連成計算方法において、前記土骨格の状態に関するデータと前記外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と土の透水性に関する第２項目と質量に関する第３項目と土骨格の接線剛性に関する第４項目とのうちの少なくとも一つを含む複数の項目を計算し、該計算した複数の項目を用いて前記連立方程式を求めるものとすることもできる。この場合、前記解析手段は、前記第１項目と前記第２項目と前記第３項目と前記第４項目とを用いて前記連立速度型運動方程式を構築する手段であるものとすることもできる。こうすれば、より適正に土骨格の変形解析を行なうことができる。

本発明の一実施例としての水と土骨格の連成計算装置２０の構成の概略を示す構成図である。 水と土骨格の連成計算プログラム３０における処理の一例を示すフローチャートである。 解析対象としての地盤を２行２列のマトリクスの４要素からなるときの具体例を示す説明図である。 各要素における局所節点を説明する説明図である。 計算条件の設定を入力する具体例としてのファイルの一例を示す説明図である。 計算条件の設定を入力する具体例としてのファイルの一例を示す説明図である。 土の設定を入力する具体例としてのファイルの一例を示す説明図である。 具体例としての地盤立体構造モデルの設定を入力するファイルの一例を示す説明図である。 ２次元平面ひずみ条件での間隙水の流れをモデル化した説明図である。 軸対称条件での間隙水の流れをモデル化した説明図である。 ３次元条件での間隙水の流れをモデル化した説明図である。 出力ファイルの一例を示す説明図である。 水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて水と土骨格の連成計算プログラムを複数回に亘って実行する際の処理の一例を示すフローチャートである。 境界での排水を許した条件で低周波の微小荷重振動を三軸砂供試体に載荷したときの計算結果の一例を示す説明図である。 境界での排水を許した条件で高周波の微小荷重振動を三軸砂供試体に載荷したときの計算結果の一例を示す説明図である。 盛土載荷に伴う自然堆積粘土地盤の遅れ沈下挙動（長期継続沈下挙動）とサンドドレーン工法などによるマスパーミアビリティ（地盤全体の透水性）改善との関係の一例を示す説明図である。 間隙水を圧縮性流体として計算する際の水と土骨格の連成計算プログラム３０における処理の一例を示すフローチャートである。 間隙水の圧縮性を考慮した場合の計算例としての盛土直下における地震動開始以降の沈下〜時間関係の一例を示す説明図である。

符号の説明

ｌ_ｍ：要素ｉの重心から要素ｍの重心への相対位置ベクトル
ｈ：要素ｉの重心における全水頭
ｈ_ｍ：要素ｍの重心における全水頭
ｈ_ｃ：ｌ_ｉと要素ｉと要素ｍとの境界面の交点における全水頭

ｒ_ａ，ｒ_ｂ：要素ｉと要素ｍとの境界面にある二つの節点の回転対称軸からの距離

次に、本発明を実施するための最良の形態を実施例を用いて説明する。

図１は、本発明の一実施例としての土骨格連成計算装置２０の構成の概略を示す構成図である。実施例の土骨格連成計算装置２０は、図示するように、一般的な汎用コンピュータ２２にアプリケーションソフトウエアとしての水と土骨格の連成計算プログラム３０がインストールされたものとして構成されている。水と土骨格の連成計算プログラム３０は、計算条件の設定や土の性状の設定，地盤立体構造モデルの設定などを行なうためにデータを入力するデータ入力モジュール３２と、入力されたデータを用いて飽和土の力学的挙動を弾塑性構成式を用いて計算する計算モジュール３４と、計算した結果を出力する結果出力モジュール３６と、から構成されている。

図２は、水と土骨格の連成計算プログラム３０における処理の一例を示すフローチャートである。以下、このフローチャートを用いて水と土骨格の連成計算プログラム３０の処理について説明する。ここで、水と土骨格の連成計算プログラム３０の処理を説明するに当たり、説明を簡易にするために、解析対象としての地盤を２行２列のマトリクスの４要素とするときを具体例として用いる。具体例の４要素を図３に示す。図示するように、４要素については下行左から要素１，要素２、上行左から要素３，要素４と要素番号を割り振り、節点については下段左から節点１，節点２，節点３、中段左から節点４，節点５，節点６、上段左から節点７，節点８，節点９と節点番号を割り振るものとした。また、各要素における局所節点については、図４に示すように、左下から左回りに局所節点１，局所節点２，局所節点３，局所節点４と局所節点番号を割り振り、各要素の節点間については図４に示すように下の節点間から左回りに節点間１，節点間２，節点間３，節点間４と節点間番号を割り振るものとした。

水と土骨格の連成計算プログラム３０が実行されると、まず、計算条件の設定の入力処理が行なわれる（ステップＳ１００）。ここで、計算条件の設定の入力としては、例えば、差分法に基づいた時間数値積分方法の選択、時刻の更新回数（計算ステップ数），２次元平面ひずみや軸対称，３次元などの計算条件を入力する。計算条件の設定を入力する具体例としてのファイルの一例を図５および図６に示す。図５のファイルでは、上から、計算条件として「平面ひずみ」であるか「軸対称」であるか「３次元」であるかの設定、変形理論として「微小変形理論」を選択するか「有限変形理論」を選択するかの設定、有限変形理論時に使用する客観性のある有効応力速度として「有効Cauchy応力のJaumann速度」であるか「Green-Nagdhi速度」であるかの設定、収束の有無の設定、構成式の選択として「下負荷面・上負荷面・回転硬化付きオリジナルカムクレイモデル」を選択するか「下負荷面・上負荷面・回転硬化付き修正カムクレイモデル」を選択するかの設定、状態量出力データの方法として「Gauss点の要素平均」であるか「ある特定のGauss点の値」であるかの設定、出力データの形式として「Text」であるか「Binary」であるかの設定、自重考慮の有無の設定、を入力する。図６のファイルでは、上から、現在のプログラム使用番号、計算ステップ数、１計算ステップ当たりの計算時間DT（以下Δｔと記載する場合がある。）(sec/step)、掘削時の掘削荷重分割数、分布荷重作用時の分布荷重分割数、土の材料定数と初期値入力ファイル、メッシュデータ入力ファイル、次のプログラム実行時に用いる計算初期条件用出力ファイル指定、座標の出力ファイル指定、状態分類(除荷、負荷(塑性圧縮/膨張、硬化/軟化))の出力ファイル指定、反力の出力ファイル指定、間隙水圧の出力ファイル指定、平均有効応力p'とせん断応力qの出力ファイル指定、過圧密の程度R(1/OCR、OCRは過圧密比)の出力ファイル指定、構造の程度R*の出力ファイル指定、硬化・軟化の敷居線の勾配の出力ファイル指定、塑性体積ひずみの出力ファイル指定、初期状態の出力ファイル指定、塑性圧縮・膨張の敷居線の勾配の出力ファイル指定、異方性の発達度合いの出力ファイル指定、経過時間の出力ファイル指定、file.dat(本データ)の出力、を入力する。

続いて、地盤の各要素に対する粘土・中間土・砂などの土の設定を入力する（ステップＳ１１０）。この土の設定では、材料定数や初期状態が入力される。土の設定を入力する具体例としてのファイルの一例を図７に示す。図７のファイルでは、弾塑性体である土の種類の数や弾性体の種類の数、構造の発展則に用いる塑性測度の種類の指定、練返し土の等方正規圧密線の切片N、限界状態定数M、圧縮指数λ、膨潤指数κ、ポアソン比ν、透水係数k、土粒子の比重Gs、初期側圧係数K0、初期過圧密比R0、初期構造の程度R*0、正規圧密土化指数m、U*の形状、構造劣化指数a、構造劣化指数b、構造劣化指数c、回転硬化指数br、回転硬化限界面mb、初期異方性の程度、初期間隙比、地盤／供試体の高さ、地盤の初期上載圧、供試体の初期圧密圧力、セル圧 (引張：正)、単位系の設定（時間、力、長さ）、などを入力する。

そして、地盤立体構造モデルの設定を入力する（ステップＳ１２０）。この地盤立体構造モデルの設定では、地盤モデルの要素に対する条件などを入力する。図３および図４に示す具体例としての地盤立体構造モデルの設定を入力するファイルの一例を図８に示す。図８のファイルの第１行目の１行６列のデータは、地盤立体構造モデルの要素数、節点数、分布荷重速度作用要素数、１要素中の節点数、計算ステップ数、弾性体の要素数である。第２行目から第１０行目までの９行５列のデータは、左から順に、節点番号、ｘ方向の境界条件の種別、ｘ方向の境界条件の値、ｙ方向の境界条件の種別、ｙ方向の境界条件の値である。ここで、境界条件の種別としては、座標制御（０）、荷重速度制御（１）、変位速度制御（−１）、変位加速度制御（−２）である。図８のファイルの第１１行目から第１４行目までの４行５列のデータは、左端が要素番号、要素番号の右から順に各要素における局所節点番号に対応する全体における節点番号である（局所節点番号と全体節点番号の対応の設定）。ファイルの第１５行目から第２３行目までの９行３列のデータは、節点番号、ｘ方向の座標、ｙ方向の座標である。ファイルの第２４行目および第２５行目については、第２４行目は反力の出力を必要とする節点数であり、第２５行目はその節点番号である。この場合、反力の出力は節点番号に対してｘ方向とｙ方向とがあるため、例えば節点番号７の場合のｘ方向については７×２−１＝１３として示し、ｙ方向については７×２＝１４として示している。ファイルの第２６行目および第２７行目の２行６列のデータは、左から順に、分布荷重速度作用要素、作用辺（作用節点間）、その左側節点の分布荷重速度のｘ方向の値、その左側節点の分布荷重速度のｙ方向の値、その右側節点の分布荷重速度のｘ方向の値、その右側節点の分布荷重速度のｙ方向の値である。ファイルの第２８行目から第３１行目までの４行６列のデータは、左から順に、要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号、土の種類である（節点間番号と全体要素番号の対応の設定と土の種類の設定）。ファイルの第３２行目のデータは、２節点間の長さ不変の条件を与えて計算を行う場合の条件数である。２節点間の長さ不変の条件数が存在する場合には、この第３２行目のデータの下に新たな行を設けて条件番号、当該条件に関する節点間の指定がなされる。第３３行目のデータは３節点間の角度不変の条件数であり、第３４行目のデータは、この３節点間の角度不変の条件に対して、左から順に、条件番号、当該条件に関する節点の指定である。第３５行目のデータは、前２節点の相対位置ベクトルと後２節点で指定された２節点の相対速度の方向が不変の条件数であり、第３６行目のデータはこの２節点の相対速度の方向が不変の条件に対して、左から順に、条件番号、当該条件に関する節点の指定である。第３７行目のデータは、等変位（速度）の条件数である。等変位（速度）の条件が存在する場合には、この第３７行目のデータの下に新たな行を設けて条件番号、当該条件に関する事項の指定がなされる。第３８行目のデータは、左から順に、モデルに入力される規則波の入力数、モデルに入力される不規則波の入力数であり、第３９行目のデータは、モデルに入力される規則波に対して、左から順に、規則波振動の種類、作用する節点数、振幅、周期である。第４０行目のデータは、規則波振動が作用する節点番号である。この場合、規則波振動が作用する方向はx方向とｙ方向とがあるため、例えば節点番号７の場合のｘ方向については７×２
−１＝１３として示し、ｙ方向については７×２＝１４として示している。ここで、モデルに入力する不規則波が存在する場合は、第４０行目の下に新たな行を設けて不規則波のファイルを指定する。

こうしてデータの入力が終了すると、入力したデータを用いて計算処理を行なう（ステップＳ１３０〜Ｓ２１０）。計算処理は、実施例では、時間積分の方法としてWilsonのθ法を用いて行なったが、Newmarkのβ法などを用いて行なうものとしてもよい。また、実施例では、空間的離散化に対して、有限要素法を用いて間隙水圧を要素中心に割り当てる田村流・Christian流に基づいて行ったが、節点に間隙水圧を割り当てるSandhu流に基づいて行う方法やあるいはメッシュフリー法などの方法を用いて行ってもよい。以下、Wilsonのθ法と田村流に基づく有限要素法を用いた処理について説明する。まず、時刻ｔ＝ｔ＋θΔｔにおける「全体位置ベクトル」、「全体速度関連増分ベクトル」ならびに「全体加速度関連増分ベクトル」と、各種状態（有効応力と間隙水圧、構造・過圧密・異方性、等）を次式（１）〜（５）を用いて予測する（ステップＳ１３０）。この予測計算処理は、短時間に収束させるために行なうものであり、必ずしも必要な処理ではない。予測計算処理を行わない場合は、時刻ｔ＝ｔで最終的に決められた値、ただし、初期時刻の場合は初期値、を時刻ｔ＝ｔ＋θΔｔにおける値として用いる。ここで、式（１）の左辺や右辺の各項のように全節点のあるベクトル量をまとめて列ベクトル表現するときには、式（６）に示すように成分表示される。また、式（４）の左辺や右辺の各項のように全要素のあるスカラー量をまとめて列ベクトル表現するときには、式（７）に示すように成分表示される。

次に、節点の変位速度を土骨格の時間的体積変化率に変換するマトリクスとして、次式（８）により全体Ｌマトリクス（全体体積変化率マトリクス）を、式（１０）により全体修正Ｌマトリクス（全体修正体積変化率マトリクス）を計算する（ステップＳ１４０）。全体Ｌマトリクスの各要素は式（９）により示される。全体ＬマトリクスはＮＥ行（ＮＰ×ｒ）列のマトリクスで、地盤立体構造モデルの設定時に指定した全体節点番号と局所節点番号の対応に従うとき、式（９）の要素ｉにおける要素Ｌマトリクスにおいて、ｍ番目の節点が全体の節点番号と一致しない列成分に対しては、数値ゼロを代入する。要素Ｌマトリクスを構成するＢｖマトリクスは、２次元平面ひずみ条件の場合は式（１２）により、軸対称条件の場合は式（１３）により、３次元条件の場合には式（１４）により表わされる。全体修正Ｌマトリクスも全体Ｌマトリクスと同様にして作られるＮＥ行（ＮＰ×ｒ）列のマトリクスである。

式（９）などの積分はいずれも、式（１５）のように、全体座標系における解析対象領域を局所座標系に変換した後、Gaussの数値積分法を用いて行われる。ここに、式（１５）は３次元条件の場合における体積積分を示す。なお、２次元平面ひずみ条件の場合は解析対象面に垂直な方向に単位厚さを考慮して、軸対称条件の場合は軸の周りの回転対称性を考慮して積分が行われる。

続いて、土の透水性を表す全体Ｈマトリクス（全体透水性マトリクス）を式（１６）により計算する（ステップＳ１４５）。ここで、排水・非排水等の水理境界条件は、式（１７）の要素Ｈマトリクスの成分を変化させることで与える。例えば、非排水条件の場合は、この条件が設定されている節点間に対応する要素Ｈマトリクスの成分の値をゼロとする。全体ＨマトリクスはＮＥ行ＮＥ列のマトリクス（ｒ：一節点あたりの成分の数）で、地盤立体構造モデルの設定時に指定した当該要素の節点間番号と全体要素番号の対応に従うとき、次式（１７）で示す要素ｉの要素Ｈマトリクスにおいてｍ番目の節点間番号が全体の要素番号と一致しない全体Ｈマトリクスの成分に対しては、数値ゼロを代入する。なお、式（１７）の要素ｉにおける要素Ｈマトリクスにおける最後の成分は、要素ｉ番目（即ち、ｉ行目）のｉ番目の列成分に代入する。式（１７）の右辺の各成分は、式（１８）により示される。各記号は、２次元平面ひずみ条件の場合は図９を、軸対称条件の場合は図１０を、３次元条件の場合は図１１を参照することになる。

そして、次式（１９）により時刻ｔ＝ｔ＋θΔｔにおける全体荷重増分ベクトルを、式（３４）より全体流量増分ベクトルの計算を行なう（ステップＳ１５０）。式（１９）の右辺第１項は公称応力速度ベクトルに関する項で式（２０）から構成され、右辺第２項はGreen-Naghdi速度を用いる場合に現れる項で式（２１）から構成される。ここで、式（２０）の右辺のＮマトリクスと式（２１）の右辺のＢマトリクスは、２次元平面ひずみ条件については式（２２）および式（２３）で表わされ、軸対称条件については式（２４）および式（２５）で表わされ、３次元条件について式（２６）および式（２７）で表わされる。また、式（２０）の公称応力速度ベクトルに関する項は、式（２８）に示すように、右辺第１項に示す境界を通じて外部から与えられる表面力速度ベクトルと右辺第２項に示すこの境界の幾何形状変化に伴って現れる表面力ベクトルの増分の和からなる。このとき、右辺第１項の時刻ｔ＝ｔ＋θΔｔにおける表面力速度ベクトルは、Wilsonのθ法に従う式（２９）により、時刻ｔ＝ｔ＋Δｔにおいて境界条件として設定された表面力速度ベクトルを換算して与える。式（２８）右辺第２項に含まれる応力又は応力速度に関する力学的境界に立てた単位外向き法線ベクトルからなるｎマトリクスは、２次元平面ひずみ条件または軸対称条件の場合は次式（３０）で与えられ、３次元条件の場合は式（３１）で与えられる。また、式（２１）の右辺第１項の有効応力成分を含む列ベクトルは、２次元平面ひずみ条件または軸対称条件の場合は式（３２）で与えられ、３次元条件の場合は式（３３）で与えられる。

また、式（３４）に示す時刻ｔ＝ｔ＋θΔｔにおける全体流量増分ベクトルは、式（３５）で与える全体流量速度ベクトルにθΔｔを乗じたものとして与える。

次に、全体Ｍマトリクス（全体質量マトリクス）を式（３６）により計算する（ステップＳ１６０）。ここで、式（３６）右辺の要素Ｍマトリクスは、式（３７）で表される。なお、全体Ｍマトリクスは（ＮＰ×ｒ）行（ＮＰ×ｒ）列のマトリクスとなる。

そして、全体Ｋマトリクス（土骨格の全体接線剛性マトリクス）を式（３８）により計算する（ステップＳ１６５）。ここで、全体Ｋマトリクスは（ＮＰ×ｒ）行（ＮＰ×ｒ）列のマトリクスとなる。全体Ｋマトリクスを構成する要素Ｋマトリクスは式（３９）で示される。この式（３９）中、右辺第１項は特定された土材料（砂や粘土など）の力学挙動を与える接線剛性を示し、右辺第２項は刻々の土の形状変化による接線剛性への寄与を表わし、右辺第３項は衝撃荷重など土骨格に加速度が顕著に働く条件下での接線剛性への寄与を示し、右辺第４項は物体力作用条件下で生じる接線剛性への寄与を示す。材料定数の中で、特にこれらの「構造低位化指数」、「正規圧密土化指数」および「回転硬化限界定数」に対して相対的な差異を与えることによって、砂から中間土、粘土までを連続的に取り扱うことができる。式（３９）の右辺第１項に含まれるＤepマトリクス（弾塑性マトリクス）は、弾塑性状態（負荷状態）では式（４０）で与えられ、さらに式（４０）は２次元平面ひずみ条件で式（４１）と式（４２）より与えられ、軸対称条件で式（４３）と式（４４）より与えられ、３次元条件で式（４５）と式（４６）により与えられる。また、弾性状態では式（４７）で与えられる。式（４１）ないし式（４６）中に現れる記号は式（４８）より与えられ、式（４８）中に現れる記号は式（４９）を一個ないしは複数個利用して与えられる。

p'〜q応力空間での塑性圧縮と塑性膨張の敷居線の傾きと硬化と軟化の敷居線の傾きに対する式（５１）の表記は、土の骨格構造のうちで、構造Ｒ＊（０＜Ｒ＊≦１、Ｒ＊が０に近いほど構造は高位であることを示す）の低位化の仕方と過圧密の程度Ｒ（０＜Ｒ≦１、Ｒが０に近いほど過圧密であることを示す）の解消の仕方を与える発展則において、次式（５２）および（５３）のように、塑性測度として塑性ストレッチングのノルムを用いた場合を示している。土の種類に応じて塑性ストレッチングのせん断成分や有効応力テンソルと塑性ストレッチングの内積で表す塑性仕事率などで与えることもできる。なお、微小変形解析を用いる場合、塑性ストレッチングは塑性ひずみ速度で表される。また、式（５２）に現れる「Ｕ＊」および式（５３）に現れる「Ｕ」はここでは非負の関数とし、例えば、次式（５４）と（５５）で与える。また、異方性を表す回転硬化変数テンソルの発展則については、次式（５６）を用いている。

式（３９）の右辺第２項を構成するＮ’マトリクスとＴ１マトリクスは、２次元平面ひずみ条件では式（５７）と式（５８）より与えられ、軸対称条件では式（５９）と式（６０）より与えられ、３次元条件では式（６１）と式（６２）により与えられる。

こうして全体Ｍマトリクスと全体Ｋマトリクスとを計算し、次式（６３）に示す全体接線剛性方程式（連立１次方程式）より、変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて未知の「加加速度場」と間隙水圧場の求解を行なう（ステップＳ１７０）。この式（６３）は、速度型運動方程式の弱形式の有限要素離散化と慣性項の影響を考慮した水から土の連成式のモデル化によって最終的に得られる連立２階常微分方程式（式（６４））において、時間微分項を陰解法（差分法）に基づいて処理したもので、変形場に対してはWilsonのθ法を用いて、間隙水圧に対しては台形公式を利用した場合の（最終的に）解くべき連立１次方程式である。全体荷重増分ベクトルおよび全体流量増分ベクトルと、全体Ｍマトリクスおよび全体Ｋマトリクスは、時刻ｔ＝ｔと時刻ｔ＝ｔ＋θΔｔとの間の時間ステップ内で繰り返し計算を行ないながら更新する。なお、式（６４）では、運動方程式に対しさらに土骨格から見た物質時間微分を施しているため、説明変数である速度ベクトルの時間に関する２階微分の項、即ち、「加加速度」項を有している点がもっとも特徴的である。これは、土骨格の構成式を速度型で与えるためだけでなく、地盤等の土の幾何形状変化を考慮するためで、これによって、変形から破壊まで、あるいは破壊後の挙動を対象にした土の体積変化の正確な計量を必要とする大変形解析を可能にするためである。またこの加加速度導入のため、Wilsonのθ法の基礎となる線形加速度法の考えに基づいて、変形に対して式（６６）〜式（６８）を誘導している。なお、加速度の影響が小さいような応答が得られるときは、式（６４）は、静的解析で得られる連立１階常微分方程式と一致する。このことは、動的および静的解析の使い分けを必要としないことを意味する。また、例えば、二節点の間に距離が不変あるいは三節点がなす角度が不変の条件などの束縛条件が節点間に課せられる場合は、Lagrangeの未定乗数法を適用して式（６４）にその束縛条件を考慮する。 時刻ｔ＝ｔ＋Δｔにおいて設定された変位、変位速度又は加速度に関する幾何的境界条件については、位置ベクトルの差即ち変位ベクトルで与えられた場合は式（６９）により、速度ベクトルで与えられた場合は式（７０）により、加速度ベクトルで与えられた場合は式（７１）により、時刻ｔ＝ｔ＋θΔｔでの加加速度ベクトルに換算して与える。

続いて、時刻ｔ＝ｔ＋θΔｔにおける加速度場、速度場、座標、各種状態量を計算すると共に各Gauss点の負荷状態を判定する（ステップＳ１８０）。時刻ｔ＝ｔ＋θΔｔにお
ける全体加速度関連増分ベクトル、全体速度関連増分ベクトルおよび座標ベクトルは、式（６３）で得られた時刻ｔ＝ｔ＋θΔｔにおける全体加加速度関連増分ベクトルと時刻ｔ＝ｔにおける幾何的諸量から、式（６６）ないし式（６８）から得られる次式（７２）ないし式（７４）により算出することができる。時刻ｔ＝ｔ＋θΔｔにおける有効応力、構造の程度R*、過圧密の程度R、異方性などの各状態量Aについては、式（７５）により算出することができる。例えば、時刻ｔ＝ｔ＋θΔｔにおける要素内の各Gauss点での有効応力速度の場合は、式（７３）より得られた全体速度関連増分ベクトルのうち、当該要素が有する節点の速度ベクトルから式（７６）を用いて算出した後、式（７５）を適用して有効応力を算出する。なお、構造の程度、過圧密の程度、異方性などの速度は、Gauss点が弾塑性状態にある場合は式（５２）、式（５３）、式（５６）などを用いて算出し、弾性状態にある場合はそれらの速度をゼロとする。時刻ｔ＝ｔ＋θΔｔにおける「塑性体積ひずみ」は、下負荷面・上負荷面・回転硬化付き修正カムクレイモデルの場合、Gauss点が弾塑性状態にあるとき式（８０）より算出するが、弾性状態にあった場合は、最後に弾塑性状態にあった時点での値のままとする。したがって、この弾性状態のとき、時刻ｔ＝ｔ＋θΔｔにおける過圧密の程度Rは、式（８０）をRについて解いた式を用いて算出することができる。また、時刻ｔ＝ｔ＋θΔｔにおける間隙水圧は式（６３）の全体接線剛性方程式より直接決まるが、間隙水圧の速度は式（７５）を逆に用いて求める。また、土要素の各Gauss点が弾塑性状態にあるか弾性状態にあるかについては、土要素の各節点から計算されるGauss点のストレッチングと有効応力などを用いて塑性乗数の分子から得られる次式（８１）に示す負荷基準を満たすときに、そのGauss点の状態は弾塑性状態にあると判定することができる。なお、ストレッチングの列ベクトル表現したものは式（８２）より計算されるが、時刻ｔ＝ｔ＋θΔｔにおけるその列ベクトルの成分表示は、２次元平面ひずみ条件の場合は式（８３）により、軸対称条件の場合は式（８４）により、３次元条件の場合は式（８５）により表される。

そして、各要素の各Gauss点の有効応力を用いて収束判定を行なう（ステップＳ１９０）。実施例では、全要素内の全Gauss点で次式（８６）を充足するか否かにより収束を判定することとした。ここで、式（８６）中の「ε」は十分に小さい正の値である。収束が判定されないときには、ステップＳ１５０の全体荷重増分ベクトルと全体流量増分ベクトルの計算に戻り、ステップＳ１５０〜ステップＳ１９０までの処理を繰り返す。

ステップＳ１９０で収束したと判定されると、時刻ｔ＝ｔ＋Δｔにおける加速度場、速度場、座標、各種状態量を計算すると共に各Gauss点の負荷状態を判定する（ステップＳ２００）。時刻ｔ＝ｔ＋Δｔの座標などは、時刻ｔ＝ｔ＋θΔｔにおける全体加加速度関連ベクトルと時刻ｔ＝ｔにおける全体加速度関連増分ベクトル、全体速度関連増分ベクトルおよび座標ベクトルから、次式（８８）ないし式（９０）により算出することができる。また、時刻ｔ＝ｔ＋Δｔの各状態量とその速度については、時刻ｔ＝ｔと時刻ｔ＝ｔ＋θΔｔにおける状態量から式（９１）および式（９２）により計算することができる。そして、時刻ｔ＝ｔ＋Δｔの各状態量と式（８１）と式（８２）を用いて、時刻ｔ＝ｔ＋Δｔでの負荷状態の判定を行う。

そしてステップＳ１００で入力した計算ステップ数に至ったか否かを判定し（ステップＳ２１０）、計算ステップ数に至っていないときには、次の計算ステップとしてステップＳ１３０の時刻ｔ＝ｔ＋θΔｔにおける各種状態（有効応力と間隙水圧、構造・過圧密・異方性、等）の予測処理に戻る。一方、計算ステップ数に至ったときには、計算結果を指定した出力ファイルに出力して（ステップＳ２２０）、水と土骨格の連成計算プログラム３０を終了する。

出力ファイルの一例を図１２に示す。この例では、第１行目は、左から順に、要素数、節点数、分布荷重速度作用要素数、節点数、節点数×次元数、節点数／要素×次元数、指定計算ステップ数、通算実施計算ステップ数、次元数、Gauss点の数／要素、弾性要素数である。２行目は、左から順に、計算条件、変形理論、有限変形理論時に使用する有効応力速度、収束の有無、構成式、状態量出力データの方法、出力データの形式、自重考慮の有無、に関する選択である。３行目から１１行目までは、左から順に、節点番号、ｘ方向の境界条件の種別、その値、ｙ方向の境界条件の種別、その値、である。１２行目は、反力の出力を必要とする節点数であり、１３行目は、反力の出力を必要とする節点における方向込みの値として示される節点番号である。図１２では、「中略」として略しているが、各要素における各Gauss点の材料定数、各節点の成分の値、各要素の間隙水圧の値、各要素における各Gauss点の有効応力の成分の値、構造の程度の値なども出力される。図１２中、「中略」の後の１行目は、要素番号１〜４に対して要素番号とその要素との関連要素の指定（４箇所）の繰り返しである。「中略」の後の２行目は、２節点間の長さ不変の条件数であり、条件数が値１以上であればその条件番号や当該条件に関する節点の指定などがその下の行に出力される。「中略」の後の３行目は、３節点間の角度不変の条件数であり、次の４行目は、この３節点間の角度不変の条件における条件番号、当該条件に関する節点の指定（３箇所）である。「中略」の後の５行目は、前２節点の相対位置ベクトルと後２節点で指定された２節点の相対速度の方向が不変の条件数であり、その次の６行目は、その条件番号、該当条件に関する節点の指定である。「中略」の後の７行目は、等変位（速度）の条件数であり、等変位（速度）の条件があるときにはその節点番号などがその次の行に出力される。「中略」の後の８行目は、各要素における土の種類である。「中略の後の９行目は、規則波の入力数と不規則波の入力数であり、規則波については１０行目と１１行目に規則波振動の種類、作用する節点数、振幅、周期、規則波が作用する方向込みの節点番号が出力され、不規則波の入力数があるときには同様にその次の２行の不規則波のファイル名や作用する節点数、不規則波が作用する方向込みの節点番号などが出力される。そして、最終行は計算開始からの通算時間数である。

以上説明した実施例の水と土骨格の連成計算プログラム３０では、土の種類を複数入力して要素に設定することにより、室内で人工的に作成するような練り返し土だけでなく、自然に堆積し、構造が発達した過圧密な地盤あるいは供試体を対象にして、動的・静的を問わない各種載荷条件に対する水と土骨格の連成場での土骨格の弾塑性変形挙動を求めることができる。このため、実施例の水と土骨格の連成計算プログラム３０では、供試体の変形挙動（分岐挙動、ポストピーク挙動を含む）や自重圧密計算、道路盛土など荷重載荷に伴う地盤の圧密沈下・側方変位挙動（「２次圧密」的な長期継続沈下挙動を含む）、掘削時の地盤変形挙動、揚水による地盤沈下問題、地盤の支持力問題、地震などによる繰返し負荷時の液状化と液状化後の沈下（揺すり込み沈下）挙動、浸透流による安定・不安定問題、砂杭等による円筒拡径を伴う締固め地盤改良問題、弾性構造物との相互作用問題（パイルドラフト敷設地盤の変形挙動）、補強土の変形挙動、マスパーミアビリティ改善による地盤改良問題、弾塑性盛土構造物の安定問題などを扱うことができる。なお、これらの問題を扱う場合、必要に応じて図２に示した水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて水と土骨格の連成計算プログラムを複数回に亘って実行する必要がある。この場合、図１３のフローチャートに示すように、水と土骨格の連成計算プログラムを実行した後に（ステップＳ３００）、計算を続行するか否かを判定し（ステップＳ３１０）、計算を続行する場合には計算条件や初期状態を変更して（ステップＳ３２０）、ステップＳ３００の水と土骨格の連成計算プログラムの実行に戻る。一方、計算を続行しない場合には、それで計算を終了する。

例えば、締固めと液状化後の揺すり込み沈下に伴う土骨格の剛性回復の計算例として、低周波の微小荷重振動を三軸砂供試体に載荷したときの計算結果を図１４に、高周波の微小荷重振動を三軸砂供試体に載荷したときの計算結果の一例を図１５に示す。図示するように、低周波の場合は振動を与えているときから締まるが、高周波の場合は振動を与えているときには締まらないで液状化し、振動を与えるのを止めると締まりだすことが解る。同じ繰り返し数のときには、圧縮量は高周波の微小荷重振動を与えた方が大きくなるのが解る。なお、液状化後の揺すり込み沈下については、初期状態に対して振動を与える条件で水と土骨格の連成計算プログラムを実行し、その後、振動を与える条件の結果を初期状態として振動を与えない条件で水と土骨格の連成計算プログラムを実行することにより得られる。

また、自然堆積粘土地盤の遅れ沈下の対策工の検討として、盛土載荷に伴う自然堆積粘土地盤の遅れ沈下挙動（長期継続沈下挙動）とサンドドレーン工法などによるマスパーミアビリティ（地盤全体の透水性）改善との関係の一例を図１６に示す。図中、マスパーミアビリティ「ｋ」が大きいほど透水性の改善を表わす。図示するように、透水性を改善しない場合は約５０年沈下が続き、盛土載荷後の残留沈下が大きいことが解る。また、透水性を改善した場合には、載荷終了後、早期に沈下量が収まり、残留沈下量も小さくなることが解る。

以上説明した実施例の水と土骨格の連成計算装置２０によれば、全体Ｌマトリクス，全体修正Ｌマトリクス，全体Ｈマトリクス，全体Ｍマトリクス，全体Ｋマトリクスを用いた全体接線剛性方程式（連立１次方程式）により変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて未知の「加加速度場」と間隙水圧場の求解を行なうことにより、地盤等の土の幾何形状変化を考慮することができ、変形から破壊まで、あるいは破壊後の挙動を対象とした土の体積変化の正確な計量を必要とする大変形解析を行なうことができる。しかも、土の種類を設定することができると共に要素毎に土の種類を設定することができるから、各種骨格構造（構造・過圧密・異方性）を有する砂から中間土、粘土までの連続的な土、あるいはシラスなどの特殊土、もしくは砂や粘土，中間土が層をなす土など、広範な種類の土骨格の力学挙動を計算することができる。また、全体接線剛性方程式は、加速度の影響が小さいような応答が得られるときは、静的解析での全体剛性方程式と一致するから、静的・動的を問わないで、土骨格の力学挙動を計算することができる。即ち、地震時の液状化後の圧密変形挙動だけでなく、圧密変形中の地震時変形挙動などのように、動的な解析に連続して静的な解析を行なったり、逆に静的な解析に連続して動的な解析を行なうことができる。

実施例の水と土骨格の連成計算装置２０では、土の種類を設定すると共に要素毎に土の種類を設定するものとしたが、土の種類を限定するものとしてもよいし、要素毎ではなく要素のブロックごとに土の種類を設定するものとしても構わない。

実施例の水と土骨格の連成計算装置２０では、全体Ｌマトリクス，全体修正Ｌマトリクス，全体Ｈマトリクス，全体Ｍマトリクス，全体Ｋマトリクスを各要素マトリクスとして厳密に計算するものとしたが、近似的に計算するものとしても構わない。

実施例の水と土骨格の連成計算装置２０では、全体Ｌマトリクス，全体修正Ｌマトリクス，全体Ｈマトリクス，全体Ｍマトリクス，全体Ｋマトリクスを計算し、これら複数種類のマトリクスのすべてを用いて全体接線剛性方程式（連立１次方程式）を連立するものとしたが、全体Ｌマトリクス，全体修正Ｌマトリクス，全体Ｈマトリクス，全体Ｍマトリクス，全体Ｋマトリクスの一部を用いて全体接線剛性方程式（連立１次方程式）を連立するものとしてもよい。この場合、これら複数種類のマトリクス以外の他の要素を用いるものとしても構わない。

実施例の水と土骨格の連成計算装置２０では、有限要素法を用いて土骨格における変形解析を行なうものとしたが、有限要素法を用いずに土骨格における変形解析を行なうものとしても差し支えない。

実施例の水と土骨格の連成計算装置２０では、計算条件の設定の入力としては、図５および図６の具体例に示すように、計算条件として「平面ひずみ」であるか「軸対称」であるか「３次元」であるかの設定、変形理論として「微小変形理論」を選択するか「有限変形理論」を選択するかの設定、有限変形理論時に使用する客観性のある有効応力速度として「有効Cauchy応力のJaumann速度」であるか「Green-Nagdhi速度」であるかの設定、収束の有無の設定、構成式の選択として「下負荷面・上負荷面・回転硬化付きオリジナルカムクレイモデル」を選択するか「下負荷面・上負荷面・回転硬化付き修正カムクレイモデル」を選択するかの設定、状態量出力データの方法として「Gauss点の要素平均」であるか「ある特定のGauss点の値」であるかの設定、出力データの形式として「Text」であるか「Binary」であるかの設定、自重考慮の有無の設定、現在のプログラム使用番号、計算ステップ数、１計算ステップ当たりの計算時間DT(sec/step)、掘削時の掘削荷重分割数、分布荷重作用時の分布荷重分割数、土の材料定数と初期値入力ファイル、メッシュデータ入力ファイル、次のプログラム時に用いる計算初期条件用出力ファイル指定、座標の出力ファイル指定、状態分類(除荷、負荷(塑性圧縮/膨張、硬化/軟化))の出力ファイル指定、反力の出力ファイル指定、間隙水圧の出力ファイル指定、平均有効応力p'とせん断応力qの出力ファイル指定、過圧密の程度R(1/OCR、OCRは過圧密比)の出力ファイル指定、構造の程度R*の出力ファイル指定、硬化・軟化の敷居線の勾配の出力ファイル指定、塑性体積ひずみの出力ファイル指定、初期状態の出力ファイル指定、塑性圧縮・膨張の敷居線の勾配の出力ファイル指定、異方性の発達度合いの出力ファイル指定、経過時間の出力ファイル指定、file.dat(本データ)の出力、を入力するものとしたが、計算条件の設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。

実施例の水と土骨格の連成計算装置２０では、地盤を構成する各要素に対する粘土・中間土・砂などの土の設定として、図７の具体例として示したように、土の種類の数や弾性体の種類の数、構造の発展則に用いる塑性測度の種類の指定、練返し土の等方正規圧密線の切片N、限界状態定数M、圧縮指数λ、膨潤指数κ、ポアソン比ν、透水係数k、土粒子の比重Gs、初期側圧係数K0、初期過圧密比R0、初期構造の程度R*0、正規圧密土化指数m、U*の形状、構造劣化指数a、構造劣化指数b、構造劣化指数c、回転硬化指数br、回転硬化限界面mb、初期異方性の程度、初期間隙比、地盤／供試体の高さ、地盤の初期上載圧、供試体の初期圧密圧力、セル圧 (引張：正)、単位系の設定（時間、力、長さ）、などを入力するものとしたが、土の設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。

実施例の水と土骨格の連成計算装置２０では、地盤立体構造モデルの設定としては、図８の具体例に示すように、図８のファイルの第１行目の地盤立体構造モデルの要素数、節点数、分布荷重速度作用要素数、１要素中の節点数、計算ステップ数、弾性体の要素数や、第２行目から第１０行目までの節点番号、ｘ方向の境界条件の種別、ｘ方向の境界条件の値、ｙ方向の境界条件の種別、ｙ方向の境界条件の値、第１１行目から第１４行目までの要素番号、各要素における局所節点番号に対応する全体における節点番号、第１５行目から第２３行目までの節点番号、ｘ方向の座標、ｙ方向の座標、第２４行目および第２５行目の反力の出力を必要とする節点数、節点番号、第２６行目および第２７行目の分布荷重速度作用要素、作用辺（作用節点間）、その左側節点の分布荷重速度のｘ方向の値、その左側節点の分布荷重速度のｙ方向の値、その右側節点の分布荷重速度のｘ方向の値、その右側節点の分布荷重速度のｙ方向の値、第２８行目から第３１行目までの要素番号、その要素の節点間１に隣接する要素番号、その要素の節点間２に隣接する要素番号、その要素の節点間３に隣接する要素番号、その要素の節点間４に隣接する要素番号、土の種類、第３２行目の２節点間の長さ不変の条件数、第３３行目の３節点間の角度不変の条件数、第３５行目および第３６行目の前２節点の相対位置ベクトルと後２節点で指定された２節点の相対速度の方向が不変の条件数、この２節点の相対速度の方向が不変の条件に対する条件番号、当該条件に関する節点の指定、第３７行目の等変位（速度）の条件数、第３８行目のモデルに入力される規則波の入力数、モデルに入力される不規則波の入力数、第３９行目のモデルに入力される規則波に対する規則波振動の種類、作用する節点数、振幅、周期、第４０行目の規則波振動が作用する節点番号、などを入力するものとしたが、地盤立体構造モデルの設定はこれに限定されるものではなく、上述の設定項目の一部を入力するものとしてもよいし、上述の設定項目以外の項目を入力するものとしても構わない。また、入力形式も如何なる形式としても構わない。

実施例の水と土骨格の連成計算装置２０では、間隙水を非圧縮の流体として計算しているが、間隙水を圧縮性流体として計算するものとしてもよい。この場合、上述した全体接線剛性方程式（式（６３））については、収束計算を行っている段階で要素iの時刻ｔ＝ｔ＋θΔｔの値として予測された間隙比eおよび式（９３）で与えられる間隙水の密度を用いつつ、修正を加えて解を求める。即ち、式（６３）の左辺の係数マトリクスについては、第一に、式（６３）中の全体Lマトリクス（全体体積変化率マトリクス）の転置の項（2箇所）をそのＬと式（９４）で与えられる全体Ｌcマトリクス（全体間隙水圧縮性マトリクス）との和の転置で置き換える。

第二に、式（６３）中の全体Hマトリクス（式（１６））を式（９６）の全体修正Hマトリクスで置き換える。また、式（６３）の右辺については、その下の項に式（９８）を加える。

したがって、最良の形態として用いた田村流・Christian流の場合のアルゴリズム、即ち、図２に例示した水と土骨格の連成計算プログラム３０における処理を示すフローチャートは、図１７に示すフローチャートのようになる。図１７のフローチャートでは、図２のフローチャートのステップＳ１４５の全体Ｈマトリクスの計算に代えて式（９６）の全体修正Ｈマトリクスを計算し（ステップＳ１４６）、その後に式（９４）の全体間隙水圧縮性マトリクス（全体Ｌｃマトリクス）を計算する（ステップ１４７）。また、ステップＳ１５０の全体荷重増分ベクトルと全体流量ベクトルの計算の中で式（９８）も計算する。なお、各ステップにおいて間隙水の密度を必要とする項（式（３５）、式（３６）など）のその値は、式（９３）を用いて、要素中心の水圧から計算した値を代入する。そして、ステップＳ１９０の収束判定までの繰り返しループをステップＳ１４０〜Ｓ１８０とする。ここで、土の設定を入力する具体例として示した図７などに、間隙水の圧縮性を表すために、間隙水の体積圧縮係数Kfに関する項目を加えるのが好ましい。また、透水係数については間隙比に応じて変化させるものとしてもよい。

間隙水の圧縮性を考慮した場合の計算例として、載荷盛土の下部領域のみについて間隙水の圧縮性を変化させたゆるい砂地盤（層厚20m、盛土高4ｍ）に、地盤底部から最大約100galの「地震動」を周期約１秒で与えた場合の、盛土直下における地震動開始以降の沈下〜時間関係の一例を図１８に示す。図示するように、間隙水が圧縮性を有する場合は、震動中に盛土下部で締め固めが顕著に発生し、非圧縮の場合と比べて盛土直下の沈下量が大きくなるが、その後の沈下量は小さくなることがわかる。なお、載荷盛土下部領域の間隙水の圧縮性変化時は、一般に地盤の各種状態が変化すると考えられるが、この計算例ではこの変化に伴って地盤状態に変化が生じないように、基準となる水圧u_w0における圧縮性間隙水の密度ρ_w0も変化させている。

実施例の水と土骨格の連成計算装置２０では、水と土骨格の連成計算プログラムにより計算された結果を新たな計算条件や初期状態として用いて水と土骨格の連成計算プログラムを複数回に亘って実行することができるものとしたが、水と土骨格の連成計算プログラムを複数回に亘って実行する回数を制限するものとしてもよいし、水と土骨格の連成計算プログラムを一回のみ実行するものとしても差し支えない。

実施例の水と土骨格の連成計算装置２０では、要素Ｋマトリクスを、式（３９）に示すように、特定された土材料（砂や粘土など）の力学挙動を与える接線剛性を示す項（右辺第１項）と、刻々の土の形状変化による接線剛性への寄与を表わす項（右辺第２項）と、衝撃荷重など土骨格に加速度が顕著に働く条件下での接線剛性への寄与を示す項（右辺第３項）と、物体力作用条件下で生じる接線剛性への寄与を示す項（右辺第４項）とにより計算するものとしたが、解析条件によってはこれらの一部の項を加味せずに要素Ｋマトリクスを計算するものとしてもよい。

実施例の水と土骨格の連成計算装置２０では、２次元平面ひずみ条件や軸対称条件、３次元条件のいずれについても解析できるものとしたが、２次元平面ひずみ条件だけしか解析できないものや２次元平面ひずみ条件と軸対称条件だけしか解析できないものとしても構わない。

実施例の水と土骨格の連成計算装置２０では、供試体の変形挙動（分岐挙動、ポストピーク挙動を含む）や自重圧密計算、道路盛土など荷重載荷に伴う地盤の圧密沈下・側方変位挙動（「２次圧密」的な長期継続沈下挙動を含む）、掘削時の地盤変形挙動、揚水による地盤沈下問題、地盤の支持力問題、地震などによる繰返し負荷時の液状化と液状化後の沈下（揺すり込み沈下）挙動、浸透流による安定・不安定問題、砂杭等による円筒拡径を伴う締固め地盤改良問題、弾性構造物との相互作用問題（パイルドラフト敷設地盤の変形挙動）、補強土の変形挙動、マスパーミアビリティ改善による地盤改良問題、弾塑性盛土構造物の安定問題などを扱うものとしたが、これらのすべてを扱うものとせず、これらの一部だけを扱うものとしたり、これらのうちの一つだけを扱うものとしても構わない。

実施例では、水と土骨格の連成計算プログラム３０がインストールされた水と土骨格の連成計算装置２０として説明したが、水と土骨格の連成計算プログラム３０の形態や水と土骨格の連成計算プログラム３０のうちの計算モジュール３４の形態としても差し支えない。

以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において、種々なる形態で実施し得ることは勿論である。

本出願は、２００５年１０月１７日に出願された日本国特許出願第２００５−３０１４１５号を優先権主張の基礎としており、本明細書にはその内容の全てが引用により含まれる。

本発明は、水と土骨格の連成計算装置を利用する製造産業などに利用可能である。

Claims (12)

1. 土骨格の状態に関するデータと外力に関するデータとを入力するデータ入力手段と、
   前記入力した土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧に関する連立速度型運動方程式を構築し、該連立方程式に変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて随時積分し、速度場および間隙水圧場の時刻歴応答を求めることにより土骨格の変形解析を行なう解析手段と、
   該解析結果を出力する結果出力手段と、
   を備える水と土骨格の連成計算装置。
2. 前記解析手段は、前記土骨格の状態に関するデータと前記外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と土の透水性に関する第２項目と質量に関する第３項目と土骨格の接線剛性に関する第４項目とのうちの少なくとも一つを含む複数の項目を計算し、該計算した複数の項目を用いて前記連立速度型運動方程式を構築する手段である請求項１記載の水と土骨格の連成計算装置。
3. 前記解析手段は、前記第１項目と前記第２項目と前記第３項目と前記第４項目とを用いて前記連立速度型運動方程式を構築する手段である請求項２記載の水と土骨格の連成計算装置。
4. 前記解析手段は、有限要素法、差分法等の数値解析法を用いて前記連立速度型運動方程式を構築する手段である請求項１ないし３いずれか記載の水と土骨格の連成計算装置。
5. 前記土骨格の状態に関するデータは、土骨格を構成する複数の要素における各要素の性状に関するデータおよび隣接する要素との関係に関するデータとが含まれるデータである請求項１ないし４いずれか記載の水と土骨格の連成計算装置。
6. 前記各要素の性状に関するデータは、土の性状に関するデータが含まれるデータである請求項５記載の水と土骨格の連成計算装置。
7. 前記外力に関するデータは、水と土骨格の連成系に与える荷重及び変位に関するデータが含まれるデータである請求項１ないし６いずれか記載の水と土骨格の連成計算装置。
8. 前記外力に関するデータは、水と土骨格の連成系に与える振動に関するデータが含まれるデータである請求項１ないし７いずれか記載の水と土骨格の連成計算装置。
9. 前記解析手段は、前記外力に関するデータを複数回変更して解析する手段である請求項１ないし８いずれか記載の水と土骨格の連成計算装置。
10. 前記解析手段は、前記外力に関するデータを変更するときには変更前の解析結果を初期状態の土骨格の状態に関するデータとして用いて解析する手段である請求項８記載の水と土骨格の連成計算装置。
11. 土骨格の状態に関するデータと外力に関するデータとに基づいて土骨格の加加速度項を含む土骨格の速度と間隙水圧に関する連立速度型運動方程式を構築し、該連立方程式に変位、変位速度又は加速度に関する幾何的境界条件及び応力又は応力速度に関する力学的境界条件並びに間隙水圧又は全水頭及び間隙水の流量に関する水理境界条件を与えて随時積分し、速度場および間隙水圧場の時刻歴応答を求めることにより土骨格の変形解析を行なうことを特徴とする水と土骨格の連成計算方法。
12. 前記土骨格の状態に関するデータと前記外力に関するデータとに基づいて、土骨格の時間的体積変化率に関する第１項目と土の透水性に関する第２項目と質量に関する第３項目と土骨格の接線剛性に関する第４項目とのうちの少なくとも一つを含む複数の項目を計算し、該計算した複数の項目を用いて前記連立方程式を求める請求項１１記載の水と土骨格の連成計算方法。

Images