非線形スペクトル確率有限要素法の提案と断層問題への適用に関する基礎的研究

全文

(1)土木学会論文集A Vol.66 No.4,643-652,2010.11. 非線形スペクトル確率有限要素法の提案と断層問題への適用に関する基礎的研究堀宗朗1 ・中川英則2 1 正会員東京大学地震研究所教授災害科学系研究部門（〒. 113-0032 東京都文京区弥生 1-1-1） E-mail: [email protected] 2 正会員茨城工業高等専門学校准教授自然科学科（〒 312-8508 茨城県ひたちなか市中根 866） E-mail: [email protected]. 地表地震断層の挙動を予測することを念頭に，地盤の物性や破壊挙動など不確定性を含む非線形連続体問題に対し，スペクトル確率有限要素法に基づく新しい数値解析手法を提案した．提案する新しい数値解析手法では，大量の計算を行う必要のあるモンテカルロ法と比較して遥かに効率的に変位・歪・応力の確率関数を計算できる．この解析手法を使ったシミュレーションにより，横ずれ断層の進展に見られるリーデルせん断帯が再現できることが示された．解の収束性も吟味され，数値解析手法として基本的な妥当性が検討された．. Key Words : nonlinear spectral stochastic finite element method, surface earthquake fault, simulation of faulting. 1.. はじめに. 用することには無理がある．地表地震断層の脅威が認識された今日，何らかの対. 地表地震断層は震源断層の破壊過程が地表に達した. 策を検討することは不可避である．断層を避けるとい. 断層である．地表地震断層は地震の度に常に形成され. う対策は現実的ではなく，また，詳細な地質調査による. るとは限らず，また，被害も断層線に沿った地域のみ. 安全性の確認と同様，ずれ変位吸収の構造形式を設け. であるため，地震工学では地表地震断層の研究は限ら. ることはコストの点から慎重にならざるを得ない．構. れてきた．地震学でも地表地震断層そのものは研究対. 造物の重要度に応じた合理的な対応が必要であるが，合. 象にはなりづらい．地表地震断層は主に地質学で精力. 理性を確保するためには，地表地震断層の挙動をでき. 的に研究されてきた．その成果として，断層形状や履. るだけ正確に予測することが不可欠である．これは「震. 歴の同定方法である空中写真判読やトレンチ調査の手. 源断層が動いた時に，地表地震断層が該当地点に形成. 法が確立し，重要構造物の建設の際に断層の地質学的. されるか否か，また形成される場合にはどこにどのよ. 調査が欠かせぬものとなっている．1995 年の兵庫県南. うなずれ変位を生じさせるか」という予測である．. 部地震以来，都市直下型地震を引き起こす活断層に対する調査が加速されたこともあり，地表地震断層の研. 地質調査以外に断層挙動の予測手法が確立されてい. 究が活性化された．さらに 1999 年のトルコと台湾の地. ない現状をみると，対策の検討に有効な手法を提案・開. 震において，高架橋の落橋やダムの倒壊等，地表地震. 発することが必要である．地震動は震源断層の破壊過程. 断層のずれ変位によって構造物に甚大な被害を被った. に適切なシナリオを与えた地震波動伝播のシミュレー. ことが報告され，地表地震断層の対策を念頭に置いた. ションから予測されることもある．これを参考にする. 研究にも関心が高まっている．. と，シミュレーションが断層挙動の予測手法の候補と. 場所によっては 1 メートルを超えるずれ変位は構造. なる．すなわち，想定された震源断層のシナリオに対. 物には致命的であるが，前述のように，地表地震断層. して，適切な地下構造の力学モデルを構築し，破壊過. は必ずしも形成されるとは限らない1) ．形成されたとし. 程を計算するのである．. ても断層は帯状であり，ずれ変位の量は断層線に沿って. シミュレーションによる断層挙動の予測には 2 つの. も大きくばらつく．このため，カルフォルニア州の活断. 大きな課題がある．第一は，問題設定である．震源断層. 層法2) のように地表地震断層を避けることも考えられて. から地表までの破壊過程を追うことは，地殻構造のモ. いる．再現期間が短いプレート境界の地震. 3),4). には活. デル化の限界のため不可能に近い．完新統 (第四紀) の. 断層法は有効であろう．しかし，我が国にあるプレー. いわゆる未固結層内の破壊過程に限っても，一様でな. ト内の活断層は状況が異なり，活断層法をそのまま採. い地盤のモデル化は難しい．第二の課題は，ばらつき. 643.

(2) 土木学会論文集A Vol.66 No.4,643-652,2010.11. の評価である．地表地震断層の形成/非形成や位置や形. いて地盤材料を対象とした非線形解析を行い，実際の. 状，そして地表面のずれ変位量のばらつきである．地震. 地表地震断層で見られるような雁行状に並ぶ塑性ひず. 断層による地盤の破壊や変形に関する多くの研究5)∼12). み集中域を形成する分岐解の再現，分岐解と外力仕事. が行われてきたが，この 2 つの課題を完全に解決する. の関係，解の収束性等を検証し，その結果を議論する．. には至っていない．そのため著者らのグループでは，地表地震断層の挙. 2.. 非線形スペクトル確率有限要素法. 動予測を念頭に新しい数値解析手法を開発しそれを実際の断層に適用することで，地表地震断層の形成/非形. 本章では，準静的状態と微小変形の仮定に基づく NL-. 成や位置や形状，そして地表面のずれ変位量のばらつきの評価を行うことを目指して研究を進めてきている．. SSFEM の定式化を示す．未固結層内の破壊過程は準静的ではないが，脆性的な地殻と比較して軟弱な地盤の. 全体の構想を以下にまとめる．. 中の破壊の伝播速度は遅くなるため，慣性項を無視する．また，変形前と変形後の配置を明確に区別し，変. 第 1 段階の研究14) では，未固結層の確率モデルを対. 形を正確に記述する有限変形理論の考慮はとても大切. 象とした確率場のスペクトル展開を用いたスペクトル. であるが，複雑性は確実に増す．その上で，非線形性. 確率有限要素法13) (SSFEM, spectral stochastic finite. による変形局所化と確率過程問題を同時に扱うことに. element method) を基に，破壊過程という非線形現象を解析できるよう拡張した数値解析手法を提案し，その数値シミュレーションを行ってきた14) ．以下では，これ. なる．地盤材料の場合，有限変形の影響はある程度構成則のパラメータに含まれてくるため，有限変形の影響も当初としては無視する．数値解析法の提案段階に. を非線形スペクトル確率有限要素法 (NL-SSFEM) と呼. おいて，初めから現実的ではあるが複雑な状態を設定. ぶ．そこでは14) ，破壊過程という非線形現象において. するよりも，確率モデルの非線形解析という問題を解. 目的とする解が正しく拾えていることを，モンテカルロ. くための基本的な仮定を明らかにし，手法全体の見通. シミュレーションから得た確率密度関数と NL-SSFEM. しを良くすることを優先した．. で得た確率密度関数を比較することでその妥当性を示している．次の段階 (第 2 段階) として 3 次元化が図ら. (1). れ15) ，モンテカルロシミュレーションで得た確率密度. 連続体の確率モデルと確率変分問題. 関数との比較および破壊の進展に伴う確率密度関数の. 連続体の確率モデルとは，構造・形状や材料特性が正. 推移が検討されている．第 3 段階で「地表地震断層の. 確に分からない物体に対し，そのパラメータの不確か. 挙動予測を念頭にしたシミュレーションとその検討」を. らしさを確率的に記述したモデルである．空間的に非. 行い，最終的 (第 4 段階) にはおよび実断層データとの. 一様な未固結層の場合，材料特性のパラメータは，空. 比較・検証を計画している．未固結層のデータの質と. 間的に変化する確率変数として与えられる．以下，空. 量が限られていることを考えると，確率モデルは比較. 間座標をパラメータとする確率変数を，確率関数と呼. 的簡単なものとせざるを得ず，高度な数値解析手法を. ぶ．確率関数を決定するためには，空間における 1 点. もってしても予測精度には限界がある．地表地震断層. ごとの平均値や分散の他に，例えば 2 点間の相関関係. に関しては統計データに基づいた幾つかの経験則が見. が必要となる．確率モデルでは変位・ひずみ・応力も空. つかっており，この経験則を超えた予測を行うことが. 間的・確率的に変動する確率関数となる．このような. 数値解析手法の性能の目標となる．本論文は，この構. 物理場と材料特性のパラメータは確率的な相関がある. 想の中での第 3 段階に位置する「地表地震断層の挙動. ため，物理場の確率関数を求めることは容易ではない．. 予測を念頭にしたシミュレーションとその検討」にあ. 特に非線形の材料特性を持つ地盤材料では，正確に相. たる．. 関を考えて物理場を求めることは難問である．提案する NL-SSFEM では，この問題に適当な近似解を与える. 本論文の主な内容は，解析手法の定式化の概要，数. ことを目的としている．. 値計算のアルゴリズム，解析手法の適用性の検証，と. 簡単のため，線形弾性体を例として確率モデルと解. いう 3 点である．第 2 章において NL-SSFEM の定式化. 析理論を説明する．ヤング率 E が確率関数である連続. の概要を示す．確率モデルとその問題設定，スペクト. 体 B を考える．変位境界条件 ui が与えられた場合，変. ル展開を用いた離散化，非線形問題への適用について. 位 ui は次の確率境界値問題の解となる． { (E(x, ω)cijkl uk,l (x, ω)),i = 0 in B × Ω,. 説明する．第 3 章では NL-SSFEM の効率的な計算アルゴリズムを説明する．NL-SSFEM ではスペクトル展開を用いている分，解析領域を同じように離散化しても，. ui (x, ω) = ui (x). in ∂B × Ω. (1). 自由度は通常の FEM に比べて大きく増加し，大規模な. ここで ω は標本空間 Ω の標本点を表し，x は連続体 B 内. 数値計算が必要となるためである．最後に第 4 章にお. における点の位置座標を表している．cijkl は，ポアソン. 644.

(3) 土木学会論文集A Vol.66 No.4,643-652,2010.11. 比 ν によって cijkl =. として表し，通常の FEM を同様に，適当に形状関数を. ν 1 (1+ν)(1−2ν) δij δkl + 2(1+ν) (δik δjl +. (n). (n). δil δjk ) として与えられる．確率関数 E は，点 x を固定. 使って確率分布 ui. すると標本点 ω ごとに値が異なり，また標本点 ω を固. る変分から次のマトリクス方程式が導かれる． ∑ (hξ n+m i[K (0) ] + hξ n+m+1 i[K (1) ])[U (m) ] = [0] (5). 定すれば x の関数となる．正確には，E は，標本空間. Ω，そのボレル集合族 F ，確率 P が与えられた確率空. を離散化すると，各 δui. に関す. m. 間 (Ω, F, P ) で定義された確率関数が E である．. ここで [K (α) ] は E (α) に対応した剛性マトリクス，[U (n) ] (n). は各 ui. 確率変分問題は，確率境界値問題 (1) の弱形式から導. に対する離散化されたベクトルである．式 (4). かれる．弱形式は物理空間と確率空間を使って導くこ. の多項式を適当な次数で打ち切り，[U (n) ] の連立マトリ. とができる．すなわち，境界で 0 となる確率関数 δui を. クス方程式を解くことで確率境界値問題の近似解が計. 使って，δuj (Ecijkl uk,l )i を B × Ω で積分すればよい．. 算される．. この結果，次の汎関数が導かれる． ∫ 1 J(u) = E(x, ω)cijkl ui,j (x, ω)uk,l (x, ω), B×Ω 2. 式 (4) の形式の ui を汎関数に代入することが SSFEM. dxP (dω). 前の例で示されたように，式 (3) の形で E が与えられ，の本質である．SSFEM では，E の相関関数を用いて E. (2). を空間分布と確率変数の積で計算しているが，この計. 汎関数 J を停留させる ui は確率境界値問題 (1) の解と. 算には，確率変数が互いに独立な正規分布であること. 一致する．またヤング率 E が正値であり，また cijkl が. が仮定されている．これは E の KL(Karhunen-Loeve). 正定値テンソルである場合，J の停留値は最小値とな. 展開と呼ばれる13) ．さらに SSFEM では，確率変数の. る．したがって，離散化された ui の中で J を最小とす. 多項式を使うことで，正規分布をしない ui のばらつき. るものは確率境界値問題 (1) の最適な近似解となる．式. を計算する．これは ui の PC(Polynomial-Chaos) 展開. (2) を最小とする際には，ui と E の相関を直接考える必要がない．近似解は自動的により正しい相関を持つ. と呼ばれる13) ．この展開では，ui の多項式の項数の設. ようになるからである．なお，式 (2) の右辺から分かる. 度を上げるが，逆に，必要な数値解析の計算量も増加. ように，J の最小値は連続体 B に蓄えられる全ひずみ. させるというトレードオフがあるからである．. 定が重要である．なぜなら，項数の増加は近似解の精. E の平均は一様であることと，ばらつきは次の相関. エネルギーの期待値を与える．. 関数から決定されることを仮定する．. (2). 確率空間での関数展開. h(E(x, ω) − E (0) (x))(E(y, ω) − E (0) (y))i. 確率境界値問題 (1) の近似解を計算するために，式. = σ exp(−r/`). (2) の J を離散化する．スペクトル確率有限要素法 (SSFEM)13) では，物理空間と確率空間で E と ui を離. ここで r は点 x と点 y の距離であり，σ と ` は分散と相関距離である．この相関関数をスペクトル分解し，得. 散化する．この SSFEM は数理的にはやや抽象的であ. られた固有関数に対応する確率変数が平均 0 の正規分. るが，E を次のように展開した場合を考えると理解し. 布に従うと仮定すると，確率変数の分散を決定するこ. やすい．. とができる．スペクトル分解の展開を 2 で打ち切った. E(x, ω) = E (0) (x) + E (1) (x)ξ(ω) ここで E. (6). (0). は点 x での平均であり，E. (1). (3). 場合，対応する確率変数を ξ と η とすると，E を. と ξ がそれぞ. E = E (0) + E (1) ξ + E (2) η. れ物理空間と確率空間の関数，すなわちばらつきの空. と近似できる．ここで，E (0) は平均，E (1) と E (2) が固. 間分布と確率変数である．なお ξ の平均は 0 である．標. 有関数である．ui を 2 次の多項式. 本点 ω に応じて ξ の値が変わり，また位置 x に応じて. (00). ui = ui. E (1) が変わる．確率変数 ξ の多項式として，ui を次の. (10). + ui. (01). ξ + ui. (20) 2. η + ui. ∞ ∑. (n). ui (x)ξ n (ω). 多項式の係数. ξη. (02) 2. +ui (nm) とする．ui. (4). η (8). が未知の関数である．この関数を離散. 化して J に代入すると，次のマトリクス方程式が導か. n=0 (n) ui. (11). ξ + ui. ように離散化する．. ui (x, ω) =. (7). れる． ∑. は未知の関数であり，これを汎関数. J を使って決定する．確率変数の冪乗 ξ n の平均は必ず (0) しも 0 ではないため，式 (4) の第 1 項目 ui が平均を. 0. 0. 0. 0. (hξ n+n η m+m i[K 0 ] + hξ n+n +1 η m+m i[K ξ ]. n0 ,m0 0. 0. 0. 0. +hξ n+n η m+m +1 i[K η ])[U n m ] = [0]. 与えないことは注意すべきである．なお，ばらつきを摂. (nm) ui. (9). 動とする確率有限要素法は式 (4) を第 2 項で打ち切り，. ここで [K (α) ] と [U (nm) ] は，E (α) と. 第 1 項が常に平均となる．確率変数 ξ の平均を hξ i. 剛性マトリクスとベクトルである．なお ξ と η は独立. n. n. 645. に対応した.

(4) 土木学会論文集A Vol.66 No.4,643-652,2010.11. であることから hξ n η m i = hξ n ihη m i であり，さらに両者は正規分布に従うから奇数乗の平均は 0 になる．非線形問題への適用. (3). 未固結層の場合，簡単な弾塑性構成則として，Mohr-. Coulomb 型ないし Drucker-Prager 型の降伏関数を用いた関連ないし非関連な流動則が仮定される16)∼20) ．関連流動則を用いた場合，降伏関数を f とすると接線弾 ep 塑性テンソル Cijkl は次のように与えられる． ep Cijkl = Ecijkl −. ここで (∇f )ij =. ∂f ∂σij. Ecijpq (∇f )pq cklrs (∇f )rs cpqrs (∇f )pq (∇f )rs. 図–1 NL-SSFEM のアルゴリズムの概要. (10). である．応力も確率関数となるた. ep め，降伏関数 f も確率関数となる．このため Cijkl を確. と同様に，変位の増分に対する線形化方程式を使って. 率空間で離散化する際には，E と括弧の中の (∇f )ij の. 弾塑性材料からなる非一様かつ非線形の確率モデルの. 両方が離散化の対象となるが，(∇f )ij は応力の非線形. 応答を求めることが可能となる．特に，NL-SSFEM で. 関数であるため (∇f )ij の離散化は極めて難しい．特に，. は一つの増分解を求めるのに反復計算を要する静的陰. f と他の確率関数との相関を考慮することは極めて難. 解法を用いている．すなわち，式 (9) と同様，変位増分. しい． ep Cijkl の離散化を計算するため，摂動展開と bounding medium 理論21),22) を利用する．摂動展開は，応力 σij. (x) を物理空間で離散化するこ dui (x, ω) の係数 dui とでマトリクス方程式が導かれ，その結果として応力の平均値 hσij i が更新されるが，更新された hσij i が元. の期待値 hσij i からのずれを摂動として使って (∇f )ij を. の値と一致し (∇f )ij を変更する必要がなくなるまで繰. (nm). (∇f )ij (σ) = (∇f )ij (hσi) + O(σ − hσi). り返し計算を行い，変位増分を計算するのである．こ. (11). として展開するものである．第 2 項以下を無視すると，. (∇f )ij の確率的なばらつきはなくなる．摂動が小さい. れが NL-SSFEM のアルゴリズムの概要であり，そのフロー図を図–1 に示す．. 場合，すなわち，応力の分散が小さい場合，摂動展開は有効であり，分散が小さくなればなるほど hσij i を使って (∇f )ij を近似する誤差は小さくなる．通常，非線形弾塑性体の数値解析では，釣合い式と構成則を満た. 3.. 大規模数値計算のためのアルゴリズム. (1). マトリクスヤコビ法15). すよう変位増分を繰り返し計算で求めるが，式 (11) を. 前章で示されたように，変位増分を確率変数 ξ と η の. 利用する場合，hσij i の初期値が必要であり，bounding. 2 次の多項式として表すと，NL-SSFEM の 1 節点の自由度は 3 × 6 = 18 となる．元々自由度が大きい 3 次元問題に対し，自由度が多項式の次数に比例して大きく. medium 理論はこの初期値を与える．bounding medium 理論は，確率モデルに対し，次式に示すように，全ひずみエネルギーの期待値の上界と下界を与える 2 つの. なることはマトリクス方程式を解く際に大きな障害と. モデル (bounding medium) を与える． ∫ ∫ 1 − 1 − E cijkl u− Ecijkl ui,j uk,l dxi u dx < h i,j k,l 2 B B 2 ∫ 1 + + < E cijkl u+ (12) i,j uk,l dx B 2. なる．一般にマトリクス方程式を解く計算量は，直接法ではマトリクスの次元の 3 乗に比例し，また反復法では 1 反復あたりの計算量はマトリクスの次元の 2 乗に比例する．したがって連立したマトリクス方程式 (9) を，[U (00) ] から [U (22) ] を一つのベクトルにまとめたベ. ここで E + = hEi と E − = 1/h1/Ei であり，u± ij はヤ. クトル [U ] のマトリクス方程式として解くと，マトリ. ング率が E ± で与えられるモデルの解である．このモ. クス方程式の次元は元のマトリクスの次元の 6 倍とな. デルが bounding media である．未固結層の場合，E +. るため，例えば直接法だと計算負荷は 216 倍に比例し. の解を初期値とすると収束が早い．なお，分散が大き. た形で圧し掛かってくる．さらに，通常の FEM の剛性. い材料物性に対しても Hashin-Shtrikman の変分原理を. マトリクスは対角項やその付近の項のみが非 0 となる. 併用し上下界の幅を狭めることで，適切な応力の平均. 粗なマトリクスであるが，[K (α) ] を合成した [K] はこ. 値の初期値を与えることは可能である21),22) ．摂動展開. のような性質はもっていない．このため，例えば古典. と bounding medium 理論を用いることで，見かけ上，. 的ではあるが計算効率の良いスカイライン法等の FEM. ep Cijkl. ．この. に適した数値解法を NL-SSFEM に適用することができ. ため，NL-SSFEM では通常の非線形有限要素法23),24). ない．以上，1) 大きなマトリクスを扱うこと，2) マト. は E のみが確率的に変動することになる. 14). 646.

(5) 土木学会論文集A Vol.66 No.4,643-652,2010.11. リクスが粗とならないこと，の 2 点が NL-SSFEM のマ. 計算効率を上げるマトリクスヤコビ法を用いたアルゴ. トリクス方程式を解く際の問題点となる．. リズムは極めて重要である．. 1) の問題を解決するため，並列計算を使った大規模計算が必要となる．この際，2) の問題に対して，NL-. (2). SSFEM のマトリクス方程式を効率的に解くことが必要となる．並列計算機の各ノードに変位の係数のベクトル [U (nm) ] の計算を割り当てた上で，マトリクスヤコビ. 標準化された KL 展開15). 前章で示したように，NL-SSFEM の精度は変位増分の離散化に依存している．変位増分は確率関数である. 法という新しいアルゴリズムを開発した．なおヤコビ. ため物理空間と確率空間で離散化されるが，物理空間. 法はマトリクス方程式の解法の一つである．方程式を. の離散化による誤差は，通常の FEM と同様，形状関数. [A] [x] = [y] とすると，[A] の対角項が非対角項に比べて大きい場合に適用されるが，これは [A] のべき乗を計算すると，非対角項は対角項に比べて指数関数的に小. の数や次元によって決まる．一方，確率空間の離散化. さくなることを利用するからでる．左から逐次 [A] をマ. 散化の誤差は減少する．しかし確率変数はヤング率の. トリクス方程式にかけていくことで，対角項が大きい. 相関関数のスペクトル展開によって定義されているた. による誤差は確率変数 ξ と η を用いた多項式によって決定される．多項式の項数が増えれば確率空間での離. [A] を使った [A] [x] = [A] [y] に変換することが n+1 n できる．変換された [A] [x] = [A] [y] を解くことは元の [A] [x] = [y] を解くことに比べ容易である．これが. め，誤差の減少の仕方は複雑である．確率変数の個数. ヤコビ法の原理である．マトリクスヤコビ法は，ヤコ. すほうが高精度となり，逆に正規分布から離れている. ビ法を NL-SSFEM のマトリクス方程式に拡張したもの. 場合には M を大きくすると高精度となる．. n+1. n+1. n. と多項式の次数を N と M とすると，一般に正解の変位増分の確率分布が正規分布に近い場合には N を増や. 未固結層の問題では，ヤング率のばらつきに正規分. である．式 (9) は次のように書き直すことができる．  [ ] [ ]  K (0) · · · hξ 2 η 2 i K (2)   .. .. ..   (13) . . .   [ ] [ ] 2 (0) 4 (2) hη i K · · · hη i K  [  ]   U (00) [0]    .  .. = .  × .    .  [ (02) ] U [0]. 布が仮定されていても，変位増分の確率分布がガウス分布に近くなるとは限らない．実際，破壊近傍の点では変位増分が対数正規分布ないしワイブル分布に近くなる．このため多項式の次数 M を大きく取る必要があり，計算効率を確保するためには，N を小さくしなければならない．しかし，スペクトル展開を途中で打ち切って少ない固有関数でヤング率を表すことには無理 √ がある．これはヤング率の標準偏差 h(E − hEi)2 i を. 左辺のマトリクスは [K (α) ] から合成される．なお，合. 計算すると簡単に理解できる．フーリエ級数展開と同. 成の際には，ξ と η が独立の正規分布に従う確率変数で. 様，スペクトル展開は直交する固有関数を使うため，展. あるため，ともに偶数乗でなければ平均が 0 となるこ. 開の項数を増やすと標準偏差は単調に増加して真の値. とを利用する．左辺のマトリクスの対角には [K (0) ] の. に収束していく．逆に言えば，展開の項数が少ない場. 項が並び，かつ，確率変数 ξ や η のべき乗も他よりも. 合には標準偏差が過小評価されてしまうのである．実. 小さい．したがって，べき乗を計算すると，対角にあ. 際，N = 2 の場合では，ヤング率の標準偏差が 10% 以. るマトリクスに比べて非対角にあるマトリクスの項は. 上過小評価されることもある．確率モデルの不確から. 小さくなっていく．これがマトリクスヤコビ法の原理. しさ自体が過小評価されてしまうと，NL-SSFEM の精. である．次元の大きい [K] のべき乗は必要な計算量が. 度は当然のことながら下がってしまう．. 多いため，実際の数値解析では [K (α) ] の各々でべき乗. 以上を整理すると，計算量を増やさずに変位増分の. を計算する．並列計算のアルゴリズムの概略は次のように整理で. 離散化の精度を上げるためには高い次数の多項式を用. きる．計算ノードは 6 の倍数とし，各ノードに式 (13). いた展開が必要であるが，これはヤング率の離散化の. ] を左から同一の行のマトリクスにかけていく．この結果，対角. 誤差を増大させる．この点に対処するため，本研究の. のマトリクスは他の位置にあるマトリクスより次第に大. 展開を用いている．標準化 KL 展開とは，展開が打ち. きくなっていく．適当な回数，行単位のマトリクスの変. 切られた場合でも標準偏差が常に真の値と一致するよ. 換を行ってマトリクス方程式を書き直し，最後に共役勾. うに，固有値の値を増した展開である．式 (3) より展開 ∑N 項数が N の時の E の KL 展開は n=0 E (n) (x)ξ (n) (ω). の行を順に割り当てる．次に対角にある [K. (0). NL-SSFEM では通常の KL 展開の替わりに標準化 KL. 配法などの繰り返し計算によって未知の変位ベクトルを求めるのである．非線形解析を行うため，NL-SSFEM では繰り返しマトリクス方程式を解く必要があるので，. 647. (0) となる．ここで √∑ ξ (ω) = 1 である．この展開では標準 N (n) )2 に過小評価されるため，固有偏差 σ が n=1 (E.

(6) 土木学会論文集A Vol.66 No.4,643-652,2010.11. 表–1 シミュレーションに用いた材料パラメータ. es. fre rig id. urf. ace 0.05m. per. sid. e. ide ic s. iod. 5m. 0.2. 0.1 1m. e bas. z y. slip. x. 図–2 シミュレーションに用いた要素分割. 関数 E. (n). (x) を. E 0(n) (x) = √∑ N. σ. n=1. E (n) (x). (14). (E (n) )2. 略化された形状を使う方が断層挙動の予測には適している．. に取替え，次のようにヤング率 E を近似的に離散化するのである．. 上記の考察に基づき，簡単な幾何形状と構成則を用 N ∑. E(x, ω) ≈. E 0(n) (x)ξ (n) (ω). いた NL-SSFEM のシミュレーションが未固結層の破壊. (15). 過程に適用できるか否かを検証するため，地盤材料の. これが標準化 KL 展開である．簡単な処理ではあるが確. サンプルに対して破壊過程の簡単な 3 次元数値シミュ. 率モデルの不確からしさが正しく評価されるため，通. レーションを行った．検証のポイントは，地表地震断. 常の KL 展開を使う場合に比べて標準化 KL 展開を用. 層の特徴である「複雑な形状の面を作りながら断層が. いた NL-SSFEM がより正解に近い変位を計算すること. 進展すること」を定性的に再現することである．これ. n=0. が報告されている. 15). は，サンプルの底面に一様な横ずれが与えられた場合. ．. でも，破壊過程はサンプル内を一様に進展せず，表面. 4.. に雁行状のせん断帯が形成されることに対応する．サ. NL-SSFEM の地盤材料への適用と検証. ンプルの形状は板形として，弾塑性構成則として次の. 態であっても，ダイラタンシーや拘束圧依存性等の特異. Drucker-Prager 型の降伏関数を用いた． √ 1 f (σ) = ασ + sij sij − K (16) 2 ここで 1 sij = σij − σδij , σ = σii (17) 3 α と K はそれぞれ内部摩擦角 θ と粘着力 c に関係する. な挙動を示すことが知られている．このため，構成則の. 係数で，. 金属材料に比べ，未固結層を構成する砂や粘土のような地盤材料は複雑な弾塑性構成則を持つことが知られている．これは地盤材料が粒状体からなる固相と液相・気相から構成されることが原因である．固相と液相の連成は重要であるが，乾燥した固相と気相のみの状. モデルとして関連流動則を仮定する場合でも，降伏関数の材料パラメータは 10 を超える場合も少なくない．. 2 sin φ 6c cos φ α= √ , K=√ 3(3 + sin φ) 3(3 + sin φ). (18). 本研究の対象である地表付近の未固結層はさまざま. である．底面に横ずれに対応する変位境界条件を与え. な堆積層から構成される．各層の構成則は複雑である．. る．この境界条件では，塑性ひずみが集中するせん断. しかし，確率モデルを使って材料の不確からしさが引き. 帯には，1) 鉛直な一つの面，2) 周期的に並んだねじれ. 起こす挙動の不確からしさにを解析する際，構成則に. た面，という少なくとも二つの解があり，後者が実際の. 過度に精緻なモデルを使うことは得策ではない．地表. せん断帯に対応する4),25) ．数理的には，一様な解から. 地震断層の予測を目的とする場合，確率モデルには簡. 周期的な解に分岐することになる26)∼30) ．したがって，. 単な構成則を用いるほうが適している．第 2 章で説明したように，本論文では構成則のモデルとしては粘着力. NL-SSFEM によるシミュレーションが 2) のせん断帯を再現するのであれば，不安定な分岐解として 2) が自. と摩擦角の二つのパラメータで決まる Drucker-Prager. 動的に選ばれることを意味する．すなわち，破壊過程. 型を採用している．このパラメータは簡単な現場計測. とそれに伴う分岐現象が自然に再現されるのである．. や室内実験で決定することができる他，地盤の種類が. シミュレーションに用いた材料パラメータを表–1 に. 分かればある程度の精度で値を推定することが可能で. 整理する．なお，初期圧縮強度と初期引張強度は粘着. ある．確率モデルの形状に関しても，構成則と同様，簡. 力と内部摩擦角によって計算される．ヤング率の分散. 648.

(7) 土木学会論文集A Vol.66 No.4,643-652,2010.11. z=5 .0[cm]. z=5 .0[cm]. [%]. [%]. 10 8 6 4 2 0. 0.1. z=0.0[cm]. Y 0.0. −0.1 −0.1. 0.0 X. U =2 [mm]. 0.1. z=0.0[cm]. Y 0.0. −0.1 −0.1. U =4 [mm]. U =6 [mm]. U =8 [mm]. U =1 0[mm]. 0.0 X. 0.1. U =2 [mm]. U =1 2 [mm]. 3.0 2.4 1.8 1.2 0.6 0.0. 0.1. U =4 [mm]. U =6 [mm]. U =8 [mm]. U =1 0[mm]. U =1 2 [mm]. 図–3 hγi の分布. 図–4 σγ の分布. と相関距離も表–1 に載せている．横ずれ方向に対して. 10，12[mm] に達したときの 6 通りを取り上げた．図–3. 垂直な水平方向と平行な方向をそれぞれ x 軸と y 軸に，. には U が増えるにしたがって上の面に向かって hγi が. また鉛直方向に z 軸をとり，原点をサンプルの底面中. 増加していく様子が示されているが，底面 (z = 0[cm]). 央に置く．実際のモデル実験と同様，サンプルは y 方. では hγi は横ずれに沿って一様に分布していたが上面に. 向に十分長いものを想定するが，y 方向に周期性を仮. いくにしたがって斜めに非一様に分布するようになる．. ことで計算領域を直方体とすることができる．. これはせん断帯が周期的にねじた面を形成することを. 領域の x，y ，z 方向の各辺の寸法をそれぞれ Dx ，Dy ，. 示している．すなわち NL-SSFEM は，一様な垂直面の. Dz とし，Dx = 25[cm] と Dz = 5[cm] を固定し，Dy の値を変えてさまざまな周期性を設定する．領域の長. 代わりにこの周期的なねじれ面を解として選択した．こ. さ Dy が周期に対応し，せん断帯が一つの鉛直面とな. の基本的妥当性を裏付けるものである．同様に γ の標. る場合は勿論，長さ Dy の複数のねじれた面になる場. 準偏差 σγ を図–4 に示す．標準偏差の値は約 3.0[%] 近. 合の解を表すことができる．図–2 に Dy = 11[cm] とし. くにまで達しており，極めてずれ変位がばらつくこと. 15). 定する. の結果は断層のシミュレーションとしての NL-SSFEM. たの時の要素分割を示す．横ずれを与える位置を底面. が伺える．しかし σγ の分布は hγi と重なっており，形. の x = 0 の直線とすると，境界条件は次のようになる．. 状にさほど大きなばらつきはない結果となっている．せん断帯が地表に達するまでに必要な外力仕事 W を. (a) 表面はトラクションフリー (b) 底面に強制変位 U ，すなわち，x < 0 と x > 0 で. 計算する．底面に働く y-方向のトラクションの平均値. ux = U/2 と −U/2，そして uy = uz = 0 (c) x 方向の側面は底面と同じ. に定義される．. を hty i とすると単位長さあたりの外力仕事は次のよう. (d) y 方向の二つの側面は変位とトラクションが連続する周期境界. ∫ ∫ 1 U W (U ) = hty idU ds (20) L 0 S ここで S と L は解析領域の底面と長さ (Dy ) である．な. 底面に与えられる y 方向の強制変位 U を，0 から所. お dU は y-方向の分散 0 の変位増分であるから，この. 定の量まで一定の増分ステップで増加させる．与えられた U に対する底面のトラクションを計算し，強制変. W は外力仕事の期待値である．せん断帯が表面に達するときの U を求めるため，図–5 に U とサンプル表面. 位との積から外力仕事を計算する．周期 Dy を変える. 上の 1 点 ((x, y, z) = (0, 0, Dz )) における γ の平均 hγi. と外力仕事も変化するが，単位長さあたりに働く外力. の関係をプロットする．この点は表面に達したせん断. 仕事を最小とする周期 Dy の解が最も起こりやすい現. 帯の中央に位置する．図より U = 7[mm] 付近で hγi が. 象に対応することになる．. 大きく増加するようになり塑性変形が表面に達したこ. 最初に NL-SSFEM の基本的妥当性を検討するため. とが分かる．以降，この U を限界強制変位 Uc とする．. に，サンプル内の変形の様子を図–3 に示す．簡単のた. 比較のため図–5 に併せて U と γ の標準偏差 σγ の関係. め底面に平行な四つの面 (z = 0，1.7，3.3，5.0[cm]) で. をプロットする．平均と同様 Uc を境に σγ が大きく増. の最大せん断ひずみ √. 加することがわかる．せん断帯が鉛直面になる場合の. γ=. 1 (11 − 22 )2 + 212 4. 外力仕事を計算し，ねじれ面となる場合との比較を行っ. (19). た．各々の強制変位 U と W の関係を図–6 に示す．せ. の平均 hγi をプロットする．また U が，U = 2，4，6，8，. 649. ん断帯がねじれ面となる場合のほうが外力仕事は少な.

(8) me a n. Dy= 6.6 cm Dy= 8.8 cm Dy=11.0 cm Dy=13.2 cm Dy=17.6 cm. . S t a n d a r d D e v ia t io n. 2.10 2.05 2.00 1.95 1.90. Uc 7[mm]. base slip [mm]. 12.0. 1.85 11.5. 図–5 強制変位 U と表面中央の hγi および σγ との関係. 11.0. base slip [mm]. 10.5. work by external force / Dy [kN]. maximum shear strain [%]. 土木学会論文集A Vol.66 No.4,643-652,2010.11. 1.80 10.0. 2.10 work by external force / Dy [kN]. External work [kNm]. 図–7 周期 Dy と単位長さ当たりの外力仕事の関係. u n ifo r ms h e a r b a n d. R ie d e l s h e a r s. Uc 7[mm]. base slip [mm]. base slip 10 [mm]. 2.05 2.00. base slip 11 [mm]. 1.95 1.90. base slip 12 [mm]. 1.85. Dy= 6.6 cm Dy= 8.8 cm Dy=11.0 cm Dy=13.2 cm Dy=17.6 cm. . 1.80 6. 図–6 強制変位 U と外力仕事 W の関係. 8. 10 12 14 length of Dy [cm]. 16. 18. 図–8 周期 Dy と単位長さ当たりの外力仕事の比較. いことが分かる．Uc の時点で外力仕事 W は約 20%少なくなっていることが分かる．この結果は，鉛直面となるせん断帯よりも，周期的にねじれたせん断帯が形. てあった．実際の未固結層のシミュレーションでは解析. 成されやすいことを示している．図–6 に示された結果は，元々は長いサンプルに対し. 領域の幅を制限しなければならないが，x 方向の側面の. Dy = 11[cm] の周期性を仮定したものである．この Dy の周期性が実際に起こるか否かは，他の周期性を仮定. 境界条件の影響を除くためには幅をできるだけ大きく. した場合との比較が必要である．そこで Dy を変えてシ. 界条件の影響を受けないぐらいに十分大きいか否かを. ミュレーションを行い，せん断帯がサンプル表面に達す. 検討するため，幅を変えたシミュレーションを行い，分. るのに必要な単位長さあたりの外力仕事 Wc = W (Uc ). 岐解の選択への影響を検討した．図–9 にサンプルの幅. を計算し，この Wc を最小にする Dy が存在するかどう図–7 において U = 10, 11, 12[mm] となった 3 つの段階. Dx を 7.5[cm] から 28.0[cm] まで変えた時の U と単位幅あたりの外力 F の関係を示す．図から Dx = 24.0[cm] 以上とっていれば U と F の関係に有意な差がないこと. を取り上げ比較したものである．5 通りの Dy の結果を. がわかる．また γ の平均や標準偏差の分布についても. 内挿した曲線から Wc を最小にする Dy が 10∼12[cm]. ほとんど同様であった．したがって，解析領域の幅を. の間にあることがわかる．図–6 に示された鉛直面とね. 厚さの約 5 倍程度もとれば，側面に人工的な境界を設. じれた面の外力仕事の比較に加え，図–8 からせん断帯. 定したシミュレーションの結果は十分幅のある場合の. が 11[cm] 程度の周期性を持つねじれた面となることが. シミュレーションと顕著な差がないことが予想される．. とる必要がある．設定された Dx = 25[cm] が側面の境. かを調べてみた．結果を図–7 に示す．また，図–8 は，. 実際に起こりやすいことが分かる．. 変形局所化現象のシミュレーションでは，解が物理. 今までのシミュレーションでは，横ずれ方向と垂直の. 空間の要素分割に依存することが知られている．これ. 方向であるサンプルの幅は常に Dx = 25[cm] で固定し. はひずみが集中する面を形状関数で表すことが難しい. 650.

(9) 土木学会論文集A Vol.66 No.4,643-652,2010.11. 図–10 離散化を上げる場合の変位ノルムの収束性. 図–9 領域幅 Dx を変えた場合の強制変位 U と外力 F の関係. 5.. ことに起因する．通常は要素分割を細かくすればするほど，ひずみの集中の度合いが高くなり，ひずみが集. おわりに. 地表地震断層の挙動予測を目的とした NL-SSFEM に. 中する領域は薄くなっていく．比較的簡単な弾塑性構. 関して，確率モデルの問題設定と定式化・大規模数値. 成則を持つ金属材料では有限変形を仮定することで要. 計算のためのアルゴリズム・基本的な妥当性の検証を. 素分割依存性がある程度回避できる．複雑な構成則を. 整理した．第 4 章に示したように，破壊過程の分岐現. 持つコンクリート材料や地盤材料に関しては，せん断. 象が再現されたことは重要である．少なくとも著者ら. 帯の進展解析にはさまざまな理論と手法が提案されて. のグループが知る限りにおいて，特に数値計算上の工. いる．本論文で提案する NL-SSFEM では，要素分割依. 夫もなく，一様な変位境界条件から周期的な破壊過程. 存性を回避するために等価エネルギー法15),31)∼33) を用. が生じるという分岐現象が再現された FEM の解析例. いている．この方法は，元来，コンクリートの解析に. は見当たらない．NL-SSFEM は，横ずれとそれに伴う. 使われた方法であり，破壊の際に解放される単位体積. 雁行状断層の発生をも含む地表地震断層のシミュレー. 毎のひずみエネルギーの値を等価エネルギーとして設. ションにも適応可能であると判断される．. 定し，要素寸法に限界値を与えるという簡便なもので. 適用可能性が検証されたため，次の段階として地表. ある．等価エネルギーは破壊靭性値等から設定するこ. 地震断層のシミュレーションへの NL-SSFEM の妥当. とが可能であり，必要な材料パラメータが一つのみで. 性と有効性を検証することを考えている．最初に，詳. すむため実用的な手法である．. 細な観察や計測データが利用できるモデル実験の再現. 要素分割の影響を調べるため，同一の確率モデルに. を行い，実験結果とシミュレーション結果の比較から. 対し 9 通りの要素分割を行ってシミュレーションを行っ. NL-SSFEM の精度や限界を調べ，その妥当性を見極める．次に，実際の地表地震断層の再現を行い，観測データ等の比較からその有効性を検討する予定である．. た．限界強制変位 Uc が与えられた時に，表面の中央点 √ において，節点変位の平均値のノルム hui ihui i を要素分割数 (N ) の逆数に対してプロットしたものが図–10 である．1/N が 0.005 以下となるような要素分割数で. NL-SSFEM が行う確率的評価は，未固結層の不確か. あれば変位はほぼ同一と見なせ，収束しているようで. らしさに起因する挙動のばらつきであることを強調し. ある．せん断帯の幅に比べて要素の寸法が大きいため，. ておく．ばらつきの大きな原因である基盤の断層変位. この結果は変位の収束を厳密に示すものではない．ま. は，NL-SSFEM では入力される量である．したがって，. た，収束の度合いを検証するためには，変位のみではな. 例えば「基盤断層変位が 1[m] 以上となる確率は 5%」と. く，応力や歪の収束を確認する必要がある．膨大な計算. いうように与えられるのであれば，NL-SSFEM の確率. 資源が必要となるが，変形局所化現象を忠実に再現す. 評価を条件付確率として組み合わせることは容易であ. ることを目的としたシミュレーションでは，より細かい. る．すなわち，基盤断層変位が 1[m] 以上のときに地表. 要素分割が必要であることは自明である．NL-SSFEM. に断層が達する確率が 10% と評価されれば，地表に断. はばらつきを含む断層挙動の予測であるため，変形局. 層が出現する確率は 0.5% となる．. 所化現象を忠実に予測することは重要ではない．変位参考文献. の期待値が概ね収束したと判断される 1/N が 0.005 以. 1) 山崎晴雄: 地震断層の出現形態とその形成条件, 地質調査所月報, Vol. 32, No. 10, pp. 574-575, 1981. 2) 中田高: カルフォルニア州の活断層法「アルキスト-プリオ. 下の要素分割で十分であると判断している．. 651.

(10) 土木学会論文集A Vol.66 No.4,643-652,2010.11 18) 田中忠次, 鵜飼恵三, 河邑眞, 阪上最一, 大津宏康: 地盤の三次元弾塑性有限要素解析, 丸善, 2000. 19) 若井明彦，鵜飼恵三: 単杭の水平載荷挙動に関する模型実験と解析, 土木学会論文集，No. 517/III-31，pp. 159168, 1995. 20) 若井明彦，鵜飼恵三: 単杭の水平挙動解析に用いる砂の 3 次元構成則の検討, 土木学会論文集，No. 589/III-42，pp. 369-374, 1998. 21) Nasser, S. N. and Hori, M.: Micromechanics: Overall Properties of Heterogeneous Materials, second edition, North-Holland, London, 1998. 22) Munasinghe, H. M. S.: Analysis on Optimistic Pessimistic Behaviour of Heterogeneous Bodies with Statistically Varying Material Characteristics, Ph. D. thesis, University of Tokyo, 1997. 23) 久田俊明，野口裕久: 非線形有限要素法の基礎と応用, 丸善, 1995. 24) 日本塑性加工学会: 非線形有限要素法 −線形弾性解析から塑性加工解析まで−, コロナ社, 1992. 25) 小出仁，山崎春雄，加藤碵一：地震と活断層の本日本の地震断層・地震の発生・予知と活断層，(株) 国際地学協会，1979. 26) Ikeda, K., Murota, K. and Nakano, M.: Echelon modes in uniform materials, Int. J. Solids Structures, Vol. 31, No. 9, pp. 2709-2733, 1994. 27) Oguni, K., Hori, M. and Ikeda, K.: Analysis on evolution pattern of periodically distributed defects, Int. J. Solids Structures, Vol. 34, No. 25, pp. 3259-3272, 1997. 28) 池田清宏, 室田一雄: 構造系の座屈と分岐, コロナ社, 2001. 29) Ikeda, K. and Murota, K.: Imperfect Bifurcation in Structures and Materials, Engineering Use of GroupTheoretic Bifurcation Theory (Applied Mathematical Sciences), Springer, 2002. 30) 藤井文夫 , 大崎純, 池田清宏: 構造と材料の分岐力学 (計算工学シリーズ 3) , コロナ社, 2005. 31) 金刀督純, 船山哲, 吉川弘道: コンクリートのひび割れモデルと有限要素解析, PDF-distribution, 2003. 32) Feenstra, P. H. and de Borst, R.: A plasticity model and algorithm for mode-I cracking in concrete, Int. J. Numer. Meth. Eng., Vol. 38, pp. 2509-2529, 1995. 33) Feenstra, P. H. and de Borst, R.: A composite plasticity model for concrete, Int. J. Numer. Meth. Eng., Vol. 33, pp. 707-730, 1996.. ロ特別調査地帯法 (Alquist-Priolo Speacil Studies Zone Act)」と地震対策，地学雑誌 No. 437/I-17, pp. 1-18, 1991. 3) 金折裕司: 甦る断層, 近未来社, 1993. 4) 加藤碵一: 地震と活断層の科学，朝倉書店，第 5 刷，1996. 5) Bray, J. D., Seed, R. B. and Seed, H. B.: 1g smallscale modelling of saturated cohesive soils, Geotechnical Testing Jornal, ASTM, Vol. 16, No. 1, pp. 46-53, 1993. 6) Bray, J. D., Seed, R. B., Cluff, L. S. and Seed, H. B.: Earthquake Fault Rupture Propagation through Soil, Journal of Geotechnical Eng., Vol.120, No.3, pp. 543561, 1994. 7) Bray, J. D., Seed, R. B. and Seed, H. B.: Analysis of Earthquake Fault Rupture Propagation through Cohesive Soil, Journal of Geotechnical Eng., Vol.120, No.3, pp. 562-580, 1994. 8) Hori, M. and Oguni, K.: Bifurcation and Stability Analyses of Initiation and Evolution of Periodic Defects, Localization and Bifurcation Theory for Soils and Rocks, T. Adachi, F. Oka and A. Yashima (eds.), Balkema, Rotterdam, pp. 127-136, 1998. 9) Ikeda, K., Murota, K. and Nakano, M.: Echelon modes in uniform materials, Int. J. Solids Structures, Vol. 31, No. 9, pp. 2709-2733, 1994. 10) Lazarte, C. A. and Bray, J. D.: A study of strike-slip faulting using small-scall models, Geotechnical Testing Jornal, ASTM, Vol. 19, No. 2, pp. 118-129, 1996. 11) 谷和夫: ジョイント要素を用いた FEM による逆断層の模型実験のシミュレーション, 地盤の破壊と歪みの局所化に関するシンポジウム発表論文集, VIII-2, pp. 215-222, 1994. 12) 谷山尚，渡辺啓行: 逆断層運動に伴う砂質表層地盤の変形に関する研究, 土木学会論文集，No. 591/I-43，pp.313325, 1998. 13) Ghanem, R. G. and Spanos, P. D.: Stochastic finite elements: a spectral approach, Springer, Berlin, 1991. 14) Anders, M. S. and Hori, M.: Stochastic Finite Element Method for Elasto-Plastic Body, Int. J. Num. Meth. Eng., Vol. 46, pp. 1897-1916, 1999. 15) Anders, M. S. and Hori, M.: Three-dimensional stochastic finite element method for elasto-plastic bodies, Int. J. Num. Meth. Eng., Vol. 51, pp. 449- 478, 2001. 16) わかりやすい土質力学原論 [第一回改定版], (社) 地盤工学会, 1992. 17) 入門シリーズ 13 土の強さと地盤の破壊入門, 土質工学会, 1987.. (2009. 7. 29 受付). FUNDAMENTAL RESEARCH ON APPLICATION OF NON-LINEAR SPECTRAL FINITE ELEMENT METHOD FOR SURFACE EARTHQUAKE FAULT PROBLEMS Muneo HORI and Hidenori NAKAGAWA This paper proposes a new analysis method for a stochastic model of non-linear continuum, in order to predict surface fault behavior. Compared with Monte-Carlo simulation, the proposed method is capable of carrying out the numerical analysis of the stochastic model more efficiently. It is shown that the simulation that uses the proposed method reproduces the formation of Riedel shears for a lateral problem. The convergence of the numerical solution is examined, and the validity of the proposed method is discussed.. 652.

(11)