論 文
確率を考慮した3次元粒状体モデルのせん断変形特性の解析
土倉泰 岸野佑次 佐武正雄
（平成3年8月31日受理）
3D Shear Deformation Analysis of a Model of Granular
Material Considering Probability
TohruTSUCHIKURA YujiKISHINO MasaoSATAKE Abstract In an attempt to analyse 3D shear deformation of granular materials， calculation was performed in which the concept of probability was introduced on the condition of interparticle sliding obtained for a regularly packed spheres． The results of the calculation considering probability showed a reasonable agreement with the behavior of actual soils at small deforma− tion level． 1．ま え が き 土は不連続な多数の粒子で構成されている材料であ る。個々の粒子は固体の結晶のように互いに強く結合 されているものではなく，比較的自由に動くことがで きる。また，土を構成する粒子は回りの粒子により拘 束されているので，流体において考えられる構成要素 のように容易に動き回れるものではない。したがって， 土の力学は固体力学や流体力学とは区別して考えられ るべきものであり，実際，粒子集合体としてとらえる ことによって土の挙動をよく説明できることがわかっ ている。 粒子の集合体のことを粒状体といい，これを扱う力 学を粒状体力学という1）。粒状体力学の中でも，粒子性 に基づく微視構造に着目するマイクロメカニックス は，近年大きく発展している分野である2）。一方，土を 扱う力学としてまとめられているのが土質力学である が，現在の土質力学は土質材料に固有な力学体系とし て見通しよく構築されているとはいいがたいという指 摘がある。これは現在の土質力学において，土が粒子 集合体である事実を直接的に扱うことにより土に関す ＊土木環境工学科，Department of Civil and Environmental Engineering ＊＊ 喧k大学，Tohoku University ＊＊＊喧k学院大学，Tohoku Gakuin University る個々の現象を有機的に結び付けるということが行わ れていないためである。土質力学を，本来の土の特徴 である粒子集合体の力学として扱う，すなわち粒子レ ベルで土の変形を考えることによって見通しのよいも のとするために，粒状体のマイクロメカニックスに関 する研究のますますの発展が期待される。 ところで，粒状体の最も基本的な変形機構は隣接粒 子間の接触点における滑りと考えられる。著者らは三 主応力下でのせん断に対する粒状体の挙動を粒子レベ ルで把握するために，等方構造をもつ球の規則配列モ デルを用いた解析を行った3）。モデルの示す挙動には 現実の土の変形特性と共通する部分があり，これを粒 子間の滑りの発生などと関連させて議論できることを 明らかにしている。 規則配列では粒子間の滑りの発生条件を完全に把握 することができる。しかし，現実の粒状体の配列は不 規則であり，また接触点数は膨大であるので，各接触 点における滑りの発生条件をとらえることは到底でき ず，変形を確定論的に扱うことは困難である。そこで， 規則配列によって把握した滑りの発生条件に確率を考 慮することによって，不規則配列でみられるような変 形を表現することを考えた。 本報告は，粒状体のマイクロメカニックスの基礎と なる粒子間の滑りの発生条件に着目して規則配列の変 形をもとに不規則配列の変形を表現し，新しい視点か らの粒状体の三主応力下での挙動の解析を試みたもの
である。 2．単純な粒状体モデルの三主応力下のせん断 土などの粒状体内部の変形機構を実物実験において とらえることは困難であるので，個々の粒子レベルで の挙動の観察のために，現実の粒状体の代用として， いくつかのモデル粒子を積み重ねて作った粒子集合体 を用いることが多い。この場合，粒子数は制限される。 ここで注意しなけれぽならないのは，粒子数に限りの ある粒子集合体を用いて変形を解析しようとする場 合，粒子配列の不規則性が結果に多大な影響を及ぼし かねないという点である。 いま，粒状体のせん断変形特性に関する中間主応力 の影響を粒子レベルで調べるものとする。純粋に中間 σf
z
図一1 規則配列モデル 主応力の影響のみを調べようとするならば，解析する 粒状体自体は等方的な構造をもつ必要がある。ところ が，有限個の粒子モデルを不規則に配列して等方体を 作ることは難しい。そこで，図一1に示す規則配列モ デルを考えることとする。6粒子は粒径の等しい球で ある。応力主軸はx，y， z軸に合わせる。このよう に軸を固定することによって，モデルを三主応力方向 に滑り変形の可能な等方的モデルとみなすことができ る。なお，規則配列を扱う利点としては，後述のよう に変形のメカニズムを確実に把握できる点を挙げるこ とができる。 図一2は図一1に示したx，g， z軸のうちの2軸 を含んだ面での断面図である。解析領域は両図に破線 で示す立方体とし，これを単位領域と呼ぶことにする。 応力やひずみは単位領域において算定される。この単 位領域を一要素として無限に続く規則配列を考えると 面心立方構造となる。図一3には面心立方構造の中で モデルの占める位置を斜線によって示している。 応力主軸をx，y， z軸に一致させたせん断を考え た場合，モデルには3つの基本的な滑りのモードがあ る。先に断面図としてみた図一2にはモデル内の12の 滑り面と3つの基本的な滑りのモードを示している。 6、とσ、の差によるせん断応力によって，x軸と平行な 滑り面で滑りが生じる。この滑りのモードをAモード と呼ぶことにする。またO、と62の差によるせん断応力 によって，y軸と平行な滑り面で滑りが生じる。この 滑りのモードをBモードと呼ぶことにする。さらにσ2 と63の差によるせん断応力によって，z軸と平行な滑 （a） Aモード （b） Bモード （c） Cモード 図一2 滑りのモード図一3 面心立方構造 り面で滑りが生じる。このモードをCモードと呼ぶこ とにする。A， B， Cの3つのモードの滑りは，一般 には同時に発生せず，段階的に生じるものであること に注意する。 上述のモデルに対して，粒状要素法3）・4）を用いた計算 を行うこととする。粒状要素法とは，粒状体を形成す る個々の粒子を剛体要素と仮定し，それぞれの粒子自 体に移動・回転を許して粒子集合体の変形を解析する 数値シミュレーション法の1つである。粒子には一般 に円や球を用いる。隣接する粒子間にはごく小さな重 なりを許す。この重なりの生じた場合に2粒子は接し ていると判定し，粒子間に仮定した剛性に基づいて接 触力を与える。DEM（Distisct Element Method）5）は同 様のシミュレーション法であるが，両方法では各粒子 の移動・回転量の算出方法が異なっている。DEMでは これらを運動方程式に基づいて算出するのに対し，粒 状要素法では各粒子の剛性に基づいて算出する。また， 粒状要素法は載荷ステップごとに必ず力の平衡状態を 得る点をDEMと異なる特徴としており，この点で同 手法は準静的挙動の解析に適する手法と考えられる。 実際にモデルをせん断するに当たり，用いる定数や 条件は，現実の土を想定して定める。粒径は0．2mm， 粒子間にはクーロンの摩擦則を仮定し粒子間摩i擦角 φμは25°とする。初期状態は5kgf／cm2（4．9×105Pa） で等方圧縮した状態とする。また粒子間の接触点には 線形剛性を仮定し，法線方向及び接線方向の接触剛性 をそれぞれk。＝142kgf／cm（1．395×105N／m）， kt＝ 99．5kgf／cm（9．762×104N／m）とする。なお，ここでは 線形接触剛性を仮定しているが，厳密にはMindlinの 接触理論解に基づく非線形接触剛性を用いるべきであ る。しかし，同一のモデルに対して線形剛性と非線形 剛性を用いて単調載荷のせん断を行って比較したとこ ろ，それぞれから得られる応力比とせん断ひずみの関 係に大差のないことが確認された6）。そこで，変形機構 をより単純化してとらえるために，線形接触剛性を用 いることとした。 以下には記号の定義式を示す。 P−÷（・1＋・2＋・3）…・……・………・…・・………・（1） 玩一 ?o（σ、一ヵ）2＋（σ、一ヵ）2＋（σ，一ヵ）2｝…・…・…・（2） μ・−22≡篇一1・…・……・……・…………・・…・・…・（3） zノ＝ε1」一ε2十ε3 ・・・・・・・・・・・・・・・・・・・・・・・・・・… 一・一・・・・・・・・… （4） r・・t−2t｛（Ei−−9−）2＋（e2一丁）2＋（烏一Ty｝…（5） dε2 一 dε3 μdε＝2 −1 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・… （6） dEi 一 dε3 ここに，σ1，62，d3は主応力，ε、，ε2，ε3は主ひずみ， dε1，dE，， dE、は主ひずみ増分である。 、平均応力pを一定とし，式（3）に示すLod6のパラ メーターμσの値を9通りに固定して応力制御で単調 載荷する。応力増分は平均応力と比較して十分に小さ い値とする。 ここに挙げた問題では，個々の粒子の剛性を基にモ デル全体の剛性を構成することによって，変形量を簡 単な手計算から求めることもできる。しかし，そのた めには滑りの発生の仕方によって何通りもの剛性行列 を予め作成しておく必要があり，その場合分けは煩雑 を極める。そこで，計算機を用いて，予め滑りの発生 の仕方を考慮する必要のない，逐次計算を行うことと した。逐次計算の内容を簡単に説明すると次のように なる。モデル全体の剛性は作成せず，各粒子の剛性を 他の粒子を固定した状態で求め，これを用いて各粒子 に与えるべき移動・回転量を算定する。算定した移動・ 回転量に1より小さい正数（ここでは0．75とする）を 乗じた値によって粒子を動かす。粒子の移動・回転に よって接触力は変化することになるが，接線方向の接 触力が法線方向接触力から定まる滑り発生の限界値を 越えた場合，超過分は力として与えないものとする。 以上の操作を，所期の応力状態のもとで各粒子が力の 平衡状態に落ち着くまで繰り返す。ところで，簡単な 例についてはモデル全体の剛性を用いて手計算から直 接求めた変形量と逐次計算から得られた変形量は一致 することを確認している。 図一一一 4が得られた結果であり，各載荷径路における 応力比とひずみの関係を示している。変形過程におい ては，まず前述のAモードの滑りが発生し，続いてB， Cモードの滑りが段階的に発生する。これらの滑りの 発生により，図中の直線が折れ曲がる。また，変形が 大きくなるとOl軸と直交する面内の接触点が分離す
ミ § ．2 お 江 器
2
拐 1．O 0．8 O．6 O．4 0．2 μσ 0．O ．0 0．25 ・−O．75−一一一 〇．5 0．75 ．0 O O．G O．2 0．3 0ctahednal Shean Stnain γ・ct（％） 図一4 規則配列モデルの応力比∼ひずみ関係 る。接触点の分離によっても直線は折れ曲がる。同図 は三軸圧縮（μσ＝−1）よりも三軸伸張（μσ＝1）に近 いほど同一の応力比に対する変形量の大きいことを示 しているが，これは現実の土の変形とよく対応してい る。なお，変形が進み接触点が分離した状態でA，B 両モードの滑りが生じると図一5に示すようにσ、方 向の2粒子が急激に近づく不安定な状態となり，その 結果それまでのひずみと比べると非常に大きなひずみ が発生する。 このように本モデルでは3つのモードの滑りと接触 点の分離とによっていくつかの変形モードが現れる。 その中で応力とひずみ増分に関するLod6のパラメー ターμσとμ4。を比較検討した場合に，両パラメーター間 のずれの傾向が現実の土の変形と合致するのは，A， C両モードの滑りの発生している変形モードであると いう結論を得ている3）。 3．確率分布の導入 上述の解析では変形モードの変化は接触点での滑り の発生がその主要因となっているが，この滑りは規則 配列全体で一斉に生じる滑りとみなされる。すなわち， 本モデルで滑りが生じている時には無限に続く配列の 中での対応する接触点においても同時に滑りが生じて いる。しかし，現実の粒状体は本モデルのように同一 粒径ではなく，不規則な構造をもっており，全体で一 斉に生じるような明瞭な滑りの発生はみられない。そ こで，ここでは1つの試みとして，不規則配列中の各 接触点で滑りが発生し始める応力比を，規則配列にお いて求められた応力比をもとに人為的に定あることを 考える。 用いる規則配列の変形においては，2．で述べたよ うに変形が進むと6、軸と直交する面内の接触点が分 離するが，ここでは変形初期に着目し，接触点の分離 する直前（後に示す計算結果では図一4のA，B， C→
図一5 不安定状態 点）までを扱うこととする。 まず，いくつかの変形モードの中から，まったく滑 りの生じていない弾性変形状態と，塑性変形に対応す ると考えられる，滑りの発生している状態の2状態を 考える。後者の状態としては，A， C両モードの滑り が同時に発生している状態を選択する。これは，2． に述べたとおり，A， C両モードの滑りが同時に発生 している状態におけるモデルの変形が，現実の土の変 形とよい対応を示すからである。せん断応力τb，tと平 均応力pの比をrで表し，滑りの生じ始める応力比を riとすれば，図一6（a）のような応力比∼ひずみ関係 が得られる。ここに，εはひずみの成分に対応するε、， ε、，ε，を代表して示すものである。図一6（b）はriが変 わることによって応力比∼ひずみ関係が変化する様子 を示す。 いま，配列に不規則性を考慮すれぽ，滑りの生じ始 める応力比rlは各接触点で異なると考えられるので， γ凌人為的に変化させることとする。riは確率的に分 布しているととらえる。本報告では，確率分布関数f 伍）には図一7に示すものを与えることとする。同分 布関数は，ri＝0， rfにおいて零となり， r、＝ryにおい て最大値を与える最も簡単な関数である。ここに，rf は応力比の最終的な値を表し，本計算では接触点の分 離の起こる応力比（μσ＝1では接触点の分離は起こら ないのでσ3＝0となる応力比）を表すものとする。ま た，f（ri）の最大値を与えるryは降伏点に相当する応力 比と考えることができる7）。r。は弾性変形状態とA， C 両モードの滑りの発生している状態との境界という意 味から，Aモードの滑りの発生し始める応力比とC モードの滑り（本解析条件下では，Cモードの滑りの 発生している場合にはAモードの滑りも同時に発生し ている）の発生し始める応力比の中間点（Ptσ ・−1では 後者の滑りは発生しないので前者の滑りの発生し始め る応力比を採用）にとるものとする。 与えられたrに対するε、，ε、，ε、の各ひずみ成分は， 上述の2状態における応力比∼ひずみ関係を用いて次 式より算定できる。 ・（・）−c・（1一元ン（ri）dri）r
（a）滑りの生じ始める応力比riE
r
（b）riによる応力比∼ひずみ関係の変化E
図一6 2つの状態で代表させた応力比∼ひずみ関係 rr 図一7 riの確率分布 2／rf f（ri） 1．O ミ§o．8
．2 富 匡 器2
拐 0．6 0．4 0．2 O．O O．0 ・・…@ regular packing−
obtained by eq （7） O．1 0．2 0．3 Octahedral Sheap Stpain γ。ct（％） 図一8 確率の導入によって生じる残留ひずみ 頂｛Cri＋G［r−・1）｝f（ri）］dri…・…・…・（7） ここに，C， Gは滑り発生前後におけるεのrに対す る変化の割合である。C， Gは規則配列の結果から与 えられるものであり，μσの値やひずみ成分により異な る値となる。4．計算結果
図一8は式（7）を用いて得られる三軸状態の応力比と せん断ひずみの関係で，点線は2．に述べた規則配列 で得られたものを示している。確率を導入した計算結 果は，規則配列で得られる折れ線とは異なって滑らか な曲線となっており，不規則配列の示す変形を表現で きていることがわかる。ところで，除荷曲線は規則配 列で得られる折れ線のうち弾性状態に相当する原点を 通る直線と平行となり，矢印で示すひずみ径路をとる ものと考えられる。規則配列の場合にはτ。ct／p＝0．4は 滑り発生条件に達しておらず弾性変形の範囲内であ り，この状態から除荷しても残留ひずみは生じない。 ところが，確率を導入した計算では残留ひずみが生じ ている。しかも，同図より三軸伸張（μσ＝1）の方が三 軸圧縮（μσ＝−1）より残留ひずみが大きく出ており， 不規則配列において一般にみられる変形と傾向がよく 一致する。 図一9はμσ＝−0．5の場合のひずみ増分に関する Lod6のパラメーターμd．の変化である。規則配列では 滑りの発生とともに急激に変化したμd，が，確率を導入 することにより徐々に変化するようになることが示さ れている。なお，7。、t≒0．17％において点線が急落する のは，接触点の分離によるものである。確率の導入に 際し，接触点の分離後は考慮しないこととしたので， 実線ではこの部分に対応する変化は表現されていな い。jne de 1．O 0．5 O．O 一〇．5 一1．0 0．0 0．1 0．2 0．3 Octahedral Sheap Stnain γoct（％） 図一9 μσニー0．5の場合のde、の変化 しかしながら，解析は微小変形の範囲にとどまった。 それ以降については接触点の分離や新たな接触点の生 成といった内部構造の変化を考慮する必要があると考 えている。なお，式（7）に示した計算式を発展させた形 で接触点の分離後までを考慮した計算を行うには次の 2点が障害となる。まず，2．で示した等方体モデル では三軸伸張状態でのみ接触点の分離が生じず，三軸 伸張状態を他の状態と同等に評価できない点である。 いま1つは，同モデルにおいて接触点の分離と同時に 不安定状態となる場合があり，この際，分離後に式（7） のGが無限大となるためにεが無限大となってしま い，変形量を有限値として特定できない点である。ま た，ここではせん断変形量のみを取り扱ったが，体積 ひずみに関しては別途検討する必要がある。
