３．材料非線形問題の解析法


３－１　はじめに

 土や岩盤の応カーひずみ関係を調べてみると、ほとんどの場合、非線形挙動を示しており、このことが地盤の応力解析を複雑なものとしている原因の一つとなっている。応カーひずみ関係が非線形であれば、それを利用して作成される剛性方程式も、当然、非線形となるため、何らかの近似解法が必要となる。このような材料非線形問題の有限要素法による計算法についてはこれまで数多くの研究が発表されてきた。基本的には代数方程式の解法であるため、どのような方法を用いても解析可能であろうが、それぞれに長所、短所があり、常に好ましい結果が得られるとは限らない。解法によってはある特定の構成式に適したものもあるであろうし、使用するモデルによってはあまり得策ではない場合もある。また、同じ非線形問題であったとしても非線形性の強い場合と弱い場合とでは扱いが異なる。当然、非線形代数方程式の近似解法であるため各種の誤差も発生するであろう。本章ではこのような材料非線形問題の解法について説明を行うが、その前に、前提となる応カーひずみの関係を整理しておこう。
 図3.1は代表的な3タイプの応カーひずみ関係を示した図である。
 
$\hspace{1em}$(a)剛塑性タイプ$\hspace{2em}$(b)完全弾塑性タイプ $\hspace{1em}$(c)一般タイプ
$\hspace{4em}$図3.1　モデル化された応力ーひずみ関係
 図(a)に示した剛塑性タイプの関係は塑性解析に利用されており、変形係数を無限大と仮定している。このため、有限要素法やRBSMのように変形から応力を求める変位型の解法ではこのような関係を単純に利用することはできない。図(b)は完全弾塑性タイプの場合で、とりあえず図中の点線で示されているような直線で表されるハードニングやソフトニングも含まれるものとする。数値計算上は、図3.2に示されるような折れ線によって表される応カーひずみ関係も後述する理由により、この図(b)に含まれるものとする。
 
$\hspace{0em}$図3.2　折れ線で近似された応力ーひずみ関係
 図(b)の関係は、ある意味で極端に理想化された応カーひずみ関係とみなすことができるが、土のような粒状体から構成される材料では、マクロ的に眺めた場合、図(C)のように明確な降伏点が得られないのが一般的である。図(b)と図(C)とでは、非線形解法の取り扱いに相違が生ずる場合もあるため、ここでは、初めに図(C)の一般的タイプに対する解法について説明し、次いで図(b)のようなタイプに対する追加説明を行うことにする。

３－２　非線形解析法の分類

 非線形解析法を大まかに整理すると、以下に示すの3つの基本的な方法に分類できる。
ただし、ここで 使用した表現は代表的なものであり、分類法そのものについても若千の相違はあるかもしれない。例えば、可変剛性法と等価荷重法といった2つに分類しているものもある。
 荷重増分法は目的とする荷重を適切な増分荷重に分割し、各荷重増分段階毎に剛性を変化させながら計算する方法で、増分荷重が作用している間は剛性を一定として線形計算を行う。この方法には
（荷重増分法）
などがある。
 一方、反復法は目的とする全荷重を作用させ、解を修正してゆく方法であり、剛性を変化させる方法とそうでなものがある。代表的な方法には
（反復法）
などかある。初めの2つが剛性を修正する方法で、3番目が初期剛性を利用する方法である。2章で説明した応力遷移法は、この反復法に含まれる。
 反復法では、はじめに仮定した応カーひずみ関係と計算による値を比較し不適切な応力もしくはひずみを修正しなければならない。この方法として以下の 2つの考え方がある。
（反復法における修正法）
ただし、分類の方法が異なれば、その位置付けも異なってくる。
 また、混合法は荷重増分法と反復法の特徴を組み合わせた方法であり、step-iteration法とも呼ばれている。代表的な方法として、修正Newton-Raphson法と組み合わせた方法が提案されている。

３－３　荷重増分法

 本節では荷重増分法の一般的な手順について説明し、次いでその長所と短所を整理し、RBSMによる離散化極限解析に適用する上での注意事項を述べる。
 一般的な荷重増分法では、与えられた応カーひずみ関係の非線形性の強さに応じて全体の荷重を幾つかの小さな増分荷重に分割する。次に、この分割された増分荷重を作用させ、線形解析を行う。その結果得られた変位や応力を前回までの変位や応力に次のように重ね合わせる。
\[{\rm (3.1)}   \left. \begin{array}{lcr} \boldsymbol{u}^{n+1} &=& \Delta \boldsymbol{u}^{n+1} \displaystyle + \sum_{i=1}^{n} \Delta \boldsymbol{u}^{i} \\[0.5em] \boldsymbol{\sigma}^{n+1} &=& \Delta \boldsymbol{\sigma}^{n+1} \displaystyle + \sum_{i=1}^{n} \Delta \boldsymbol{\sigma}^{i} \end{array} \right\} \]
ここで、$\boldsymbol{u}^{n+1}$ あるいは $\boldsymbol{\sigma}^{n+1}$ は(n + 1)番目の増分荷重が作用し終わったときの全変位および全応力を示しており、$\Delta \boldsymbol{u}^i$ あるいは $\Delta \boldsymbol{\sigma}^i$ は $i$ 番目の増分荷重が単独で作用したときの増分変位あるいは増分応力である。
 荷重増分法では、各増分荷重の作用中に剛性を変化させないという前提により線形解析を可能にしており、非線形性の導入は、次の増分荷重による解析に移る直前に行われる。ところが、 図3.1(c)のように応カーひずみ関係が曲線で表される材料の場合、剛性が連続的に変化しているため、剛性決定にあたり、何らかの近似が必要となる。
 この剛性の決定法として、一般的によく利用されているのが、接線剛性を用いる方法である。図3.3は各応カレベルでの接線剛性を示した図である。
 
$\hspace{4em}$図3.3　接線剛性
この図を眺めてみると、荷重増分法により一見正しく解が求められるように見受けられるが、実際には、図3.4に示すよう、荷重増分回数の増加とともに誤差が累積してゆく。
 
$\hspace{2em}$図3.4　荷重増分法の応力ーひずみ曲線
 図からも理解できるように、荷重増分法ではA点の接線剛性を用いて荷重増分1を作用させることで線形解析を行う。したがって、B点の応力とひずみではなくA'点の応力とひずみが得られる。同様に、荷重増分2の段階ではB点の接線剛性を用いるため、B点を基準に考えるとC点ではなくB'点の値が求まる。しかし、実際には式(3.1)のような重ね合わせを行うため、荷重増分2が作用し終わった段階ではB”点での応力とひずみが得られる。したがって、荷重増分量を無限に小さくしてゆけば正解に近づくが、数値計算上それは現実的ではないため、誤差が累積する。この方法は接線剛性を利用することから、接線剛性法あるいは接線係数法とも呼ばれている。
 オイラー法は、このようにして決定された性を利用し、
\[{\rm (3.2)}   \boldsymbol{K}^{n} \cdot \Delta \boldsymbol{u}^{n+1} = \Delta \boldsymbol{F}^{n+1} \] を解く最も単純な方法である。ここで、$\boldsymbol{K}^{n}$ は前回の応カレベルにおける接線剛性を用いていることを示している。
 一方、誤差の累積を減らすため、図3.5に示すようなRunge-Kutta法を利用した荷重増分法も提案されている。
 
$\hspace{4em}$図3.5　Runge-Kutta法による増分法
この方法は、ある荷重増分段階で、一旦、$\theta \cdot \Delta \boldsymbol{F}^{n+1} (0 \lt \theta \lt 1)$ なる仮の増分荷重を作用させ、A点の接線剛性より B'点の応力とひずみを求める。次にこの点B'に対応するB点の接線剛性により 再度A'点から計算をやり直してC'点の応力とひずみを求める。すなわち、
\[{\rm (3.3)}   \left. \begin{array}{lcr} \boldsymbol{K}^{n} \cdot \Delta \boldsymbol{u}^{n+\theta} & = & \theta \cdot \Delta \boldsymbol{F}^{n+1} \\[0.5em] \boldsymbol{K}^{n+\theta} \cdot \Delta \boldsymbol{u}^{n+1} & = & \Delta \boldsymbol{F}^{n+1} \end{array} \right\} \]
なる関係を解くことになる。この方法は、実際の点Cからは離れているが、オイラーの方法によるC”点の値よりはC点に近づく。ただし、剛性方程式を2回解くため、1回当たりの荷重増分量に対する計算時間は増加する。
 以上のように、応力とひずみの関係が連続的に変化するような場合、荷重増分法による解法では誤差が累積してゆくが、もし、 図3.1(b)のように完全弾塑性タイプであった場合は、正しい荷重増分を選択することによって誤差を累積させずに再現することができる。これは、降伏点までの区間において応力ーひずみ曲線が線形（直線）で表されるため、前回の接線勾配が、今回の応力レベルまで変化しないということに起因する。このことを拡張して考えれば、折れ線で表されるような応カーひずみ関係を持つ材料に対しても同様なことを述べることができる。
 ただし、この場合、領域内の各点において応カレベルが異なるため、どの点かどの荷重レベルで降伏するかを予め想定し、荷重増分をデーターとして入力することは不可能に近い。そこで考えられたのが山田の方法である。山田の方法では、ある要素が降伏するのに必要となる荷重増分量を応力状態から計算により決定する方法であるため、荷重増分量を入力せずにすむ。
 次に、この山田の方法を用いた荷重増分法のアルゴリズムを説明する。山田の方法に関する計算アルゴリズムは以下の通りである。
この記述に基づいた川井モデルにおけるプログラム化の処理手順については「離散化極限解析プログラム」に記載されているので参考にしてほしい。
 さて、山田の方法では、荷重増分量を荷重増分率という形式で求めている。図3.6には、この荷重増分率に関する考え方が示されている。
 
$\hspace{0em}$図3.6　山田の方法における荷重増分率
図中のP点は前回の、すなわち、nステップ目の応力状態を表しており、R点は今回求まった増分応力 $\Delta \boldsymbol{\sigma}$ を前回の値に加え合わせた仮りの(n + 1)ステップ目における応力状態を表している。P点からR点までは当然、線形的に変化していることに注意してほしい。図からも理解できるように、今回の応力増分 $\Delta \boldsymbol{\sigma}$ は $\overline{PR}$ となっている。しかし、R点で明らかに初期降伏曲面上の点を越えていれば $\overline{QR}$ という実際には起こり得ない余分な力がその点に作用していることになる。そこで、とりあえず $\Delta \boldsymbol{\sigma}$ のうち $\overline{PQ}$ で表される弾性部分を取り出す。このとき、$\overline{PQ}$ と$\overline{PR}$ の比
\[{\rm (3.5)}   r = \frac{\overline{PQ}}{\overline{PR}} \]
を用いて、
\[{\rm (3.6)}   \boldsymbol{\sigma}^{n+1} = \boldsymbol{\sigma}^{n} + r \cdot \Delta \boldsymbol{\sigma} \]
とすれば(n + 1)ステップ目の応力は降伏曲面上に載る。これは、完全弾塑性の場合、降伏応力以下では弾性状態であることを利用しており、ここまでの計算は線形計算でよいことになる。式(3.5)で表されるような $r$ のことを山田の方法では荷重増分率と呼んでいる。この $r$ は1つの要素（川井モデルの場合、ばねに該当する）について求めたものであるから、これをすべての要素について求め、その最小値を今回の荷重増分率とすれば、すべての要素は弾性状態か少なくとも降伏曲面上にあることになる。
 このような荷重増分率より
\[{\rm (3.7)}   \Delta \boldsymbol{P}^{n+1} = r^{n+1}\boldsymbol{P} \]
とすることで荷重増分董を決定することができる。ここで、$\boldsymbol{P}$ は初めに与えた荷重値である。最終的に荷重増分率 $r_i$ の合計が1になれば、初めに設定した荷重がすべて作用したことになる。
 以上、各種の荷重増分法における解析の手順を中心に説明してきた。ここで、この荷重増分法の長所、短所を整理してみよう。
（長所）
（短所）
 さて、荷重増分法そのものは非線形解析法の一つであるため、有限要素法ばかりでなく川井モデルによる離散化極限解析に対しても同様な手順により適用することができる。川井モデルによる離散化極限解析では、応カーひずみ関係を完全弾塑性タイプ、もしくは折れ線で近似し、山田の方法による荷重増分法を適用することがよく行われている。この理由は、離散化極限解析の第一の目的が、極限耐力を求めることにあるためであり、この場合、応カーひずみ関係を直線で近似しても、極限耐力に与える影響は少ないと考えられるからである。
 このように考えると、荷重増分法、特に山田の方法は長所の箇所でも上げたように、より高い精度の崩壊荷重が得られるため、離散化極限解析に適した非線形解析法と考えることができる。また、荷重増分法は、すべり線の発達状況を把握することができるため、施工途中における各種の対策工法の検討にも利用できるものと思われる。
 山田の方法においてボイントとなるのが荷重増分率の計算方法であった。川井モデルではベクトル量としての表面力を扱うため、主応力を扱う有限要素法とは若干異なる。ここでは、破壊条件としてモール・クーロンの条件を用いた場合の荷重増分率の計算方法について説明する。図3.7はモール・クーロンの条件を用いた場合の応力経路を示したものである。
 
$\hspace{0em}$図3.7　モール・クーロンの条件における応力経路
 いま、前回までの表面力を $(\sigma_n,\tau)$ とし、新しい荷重増分による増分表面力を $(\Delta \sigma_n,\Delta \tau)$ とすれば、今回の表面力は
\[{\rm (3.8)}   ( \tau + r \cdot \Delta \tau)^{2} \leq \{ c - (\sigma_{n} + r \cdot \Delta \sigma_{n}) \tan \phi \}^{2} \]
なる関係を満足していなければならない。ここで、$r$ は該等するばねにおける未知の荷重増分率である。式(3.8)を $r$ について整理すると、
\[{\rm (3.9)}   A \cdot r^{2} + B \cdot r + C \leq 0 \] \[\hspace{6em} \begin{array}{l} A = \Delta \tau^{2} - (\Delta \sigma_{n} \tan \phi)^{2} \\[0.5em] B = 2 \{ \tau \Delta \tau - ( c - \sigma_{n} \tan \phi ) \Delta \sigma_{n} \cdot \tan \phi \} \\[0.5em] C = \tau^{2} - (c - \sigma_{n} \tan \phi)^{2} \end{array} \]
となる。ここで、等号が成立するのは、今回の表面力がちょうど降伏曲面上の点 $Q$ に載った場合である。山田の方法では、この等号が成立するための荷重増分率 $r$ を、すでに降伏しているばねを除き、すべてのばねについて求める。このようにして得られた値のうち、最小の値を系全体の荷重増分率とすれば、表面力は弾性状態か少なくとも降伏曲面上にあることになる。一旦、 降伏したばねについてはは、図3.7に示すよう降伏曲面上を移動するものとし、ばね剛性行列を1章で説明した塑性化後の剛性行列に変更する。
 他の破壊条件を利用する場合においても、ここで説明したように、合力によって表された式を利用すれば、川井モデルによる離散化極限解析においても、従来の有限要素法で利用されているような荷重増分法を適用することができる。

３－４　反復法と混合法

 地盤の応力解析における問題点の一つに現地盤の初期応力の推定がある。例えば、切土やトンネル安定問題などの掘削が伴う解析、あるいは盛土などのように既存の斜面に新たに荷重が加わるような解析、支持力解析など、どの問題を取り上げても初期応力の状態が解析結果に影響を与えていることがわかる。
 解析を行うに当たり、このような初期応力を計測などによって得られた結果から推定し、入力という手段で取り込むことも考えられるが、一般的には広範囲な領域で、且つ不均質な地盤を扱うことが多いため、自重を荷重として計算により求めることがよく行われている。
 計算により初期応力状態を求める場合、破壊がどのように進展してゆくかということよりも、むしろ、最終的な応力状態がどのようになるかが重要である場合が多い。これは地盤の応力解析が多くの場合、初期応力状態を基礎としてその後の施工による地盤の変化に重点を置いているためである。前節において説明した荷重増分法は、破壊の進展状況や崩壊荷重を求めるといった強非線形問題に適した解析法であり、初期応力の分布状態を知るといった弱非線形問題（自重により崩壊現象が発生しないという意味で）の場合には必ずしも適切な解法とはいいきれない。特に、応カーひずみ関係が前節の図3.1(c)で示したように曲線で表される場合は増分荷重の選び方により初期応力状態に誤差が含まれるため、注意が必要である。
 そこで、本節では荷重レベルに応じた破壊の進展状況は把握することができないが、応力状態をより簡単に求めることのできる反復法について説明する。反復法は計算によって得られた解を修正してゆく方法で、これに含まれる代表的な非線形解析法として以下の三つが上げられる。