個別要素法のパラメータ

昔、落石対策シミュレーションの結果を発表していた後輩君に、「結果の妥当性をどのように担保しているのか？」と聞いたことがあります。

通常、現地で落石実験を実施できないため、再現解析ができません。おそらく、文献値による予測のみの計算だったのでしょう。明確な返答はありませんでした。

それから十年以上経過しましたが、未だ落石シミュレーションを実務で積極的に利用しているという状況にはなっていません。というか、私も実施したことはありません。
試験ができない、再現計算ができない、では進みません。が、その上での設定方法やツールの特性、適用限界を知ろうとしない、知らないはダメでしょうね（反省）。

先日、シュミットハンマーの反発度と反発係数を関連付けようとしていた文献を見かけました。この続きを読みたかったのですが、探せず。この先、どうなっているのか期待されます。
藤村ほか（１９９０） DEM解析における要素定数の検討

反発係数がわかると、剛性に対する減衰を設定できます。「落石対策便覧に関する参考資料」で頻繁に出てきますので、実務利用可能でしょう。エコーチップも使えるかもしれません。こちらの方が原理は近いので。

ちょうど PFC を触っていますので、DEM のパラメータに関する文献の抜粋をメモしておきます。(今後、追加します）
**************************************
●ばね定数
日本道路協会（２００２）落石対策便覧に関する参考資料ー落石シミュレーション手法の調査研究資料ー

固定するバネ係数の例としては、落石の辺の大きさが0.5m程度の場合、前述の計算例にしたがい 3,000kN/m-100,000kN/m程度の値が参考となるが、低い値を用いると貫入量が落石径の10%を越え、Hertzの接触理論式の仮定とも矛盾するので、20,000kN/m以上の値が適当と考えられる。例えば本報告書の第4章におけるDEMの計算例では、バネ係数を50,000kN/mに固定して計算している。今後は、バネ係数の影響を詳細に分析するとともに、衝突時の圧縮力、せん断力、接触位置の変位を実験により分析して衝突時の接触モデルの開発を継続することが望まれる。

地盤工学会誌
初級講座 地盤工学のための個別要素法　3．一次元の個別要素法

多くの問題では多少ばね定数を下げても大きな問題にならないことが多い。例えば、筆者の経験に基づけば、落石や土砂流動問題に DEM を適用する場合，弾性
波速度から推定するようなばね定数を基準として、その値を 1 オーダー、場 合によっては 2オーダー下げたとしても結果に大きな影響を与えないことが多々ある。そのため 、この事実を知っていれば、実務レベルで DEM を利用する場合に必要以上に計算時間に悩まされずに済む。ただし、どの程度ばね定数を下げてよいかというのは問題に依存するため、例えば反発係数一定の条件でばね定数を変化させて、得られる結果の平均的な挙動が大きく変化しないことを確かめた上で低めのばね定数を用いるというのが理想的である。

地盤工学会誌
初級講座 地盤工学のための個別要素法　6．パラメータの設定と土と地盤の作り方 その1

線形ばね係数を試算する方法は、伯野5）によって、多質点一ばね連結系の一次元波動伝播速度の考察に基づいて、粒状体中の弾性波P波、S波の速度 Vp、Vs を用いることで、次式のように提案されている。
kn ＝1/4πρVp^2，ks＝l/4πρVs^2　・・・式（6．10）
なお、岩のように Vp=1000（m/s） を超えるような材料ではknを10^9N/m以上などの大きな桁数で設定することになる。しかし，knが大きくなると計算時間刻み（⊿t）を小さくせざるを得ないため、結果的に計算時間を増大させる原因となる。一方、後述するが、粒状体の変形・破壊挙動、流れ挙動，衝撃力の挙動を調べた結果、二次元ではばね定数が10^6N/m以下では粒子要素の重なりによる圧縮性が顕著になり、粒状体本来の挙動とはかけ離れることが分かっている。式（6．10）はあくまでも試算値で，一般に，この試算値は，かなり大きめのばね定数値を与え、計算時間刻みが小さくなり、計算時間を増大させる原因となる。一方で、粒状体の変形挙動は粒子構造が支配的であるので、ある程度以上のばね定数であれば、計算で得られる力学特性に与えるばね定数の影響が小さい。これに鑑みて，先の試算値よりも小さなばね定数を設定し、実効性を高めることができる。

●粘性係数、減衰定数
地盤工学会誌
初級講座 地盤工学のための個別要素法　6．パラメータの設定と土と地盤の作り方 その1

円形（球形）粒子を用いた解析（減衰定数以外にエネルギー損失メカニズムがない）において、解析対象を弾性衝突としてモデル化するのであればh＝0、接触部でのわずかな塑性変形が発生するような接触状態に相当するのはh＝0.2〜0.4程度であるといえる。また，接触部の塑性変形や隅角部の欠損が生じるような接触状態を考慮するにはh=0.4〜0.7程度，隅角部の破損や非球形に起因して並進運動が回転運動へと転換されることを考える場合にはh＝0.7〜1.0程度が妥当と考えることができる。従来から数値計算上の安定性確保から臨界減衰h=1.0がよく使用されているが，落体運動からみると，非円形（非球形）な粒子の接触部の塑性変形や隅角部の破損によるエネルギー損失を円形（球形）粒子で表現していたことに相当していたと言える。

第21回中部地盤工学シンポジウム
21.落石挙動のばらつきを考慮した堆積層の衝撃吸収効果

堆積層に落下させた場合には、hによる影響はほとんどないことが分かる。これは，落体衝突エネルギーが主に堆積層の変形に伴って逸散することによるところも大きいが、堆積層の場合では、粒子間の力のやり取りが非常に短期的な接触中でしか生じず、その効果が十分に発揮されるより前に接触が解放されてしまうためだと考えられる。そのため今回は、一般的に使われている臨界減衰h=1 を用いた。

●摩擦係数
地盤工学会誌
初級講座 地盤工学のための個別要素法　6．パラメータの設定と土と地盤の作り方 その1

DEM解析において，基本的な粒子要素としてよく用いられている円形（又は球形）粒子による解析を行うと，土塊は強い正のダイレイタンシー特性（せん断すると膨張する特性，圧縮性が小さい特性）を示す。これは，破壊の進行性に大きな影響をもたらす。また，実際の砂レベルの高い破壊時の内部摩擦角を得るために，粒子間摩擦係数を限りなく大きくしても内部摩擦角は一向に大きくならない（せいぜい30°である）。一方，粒子に少し凹凸を付けた粒子形状の要素を用いるだけで，十分な強度が発現するようになるが，圧縮量も増加し，堆積層の緩衝効果も高くなる（ただし，滑らかな楕円形状では余り強度は大きくならない）。非球形粒子を用いることで，計算時間が長くなるが，特に変形の局所化，破壊現象や流動後の堆積現象などを丁寧に精度良く扱いたい場合には粒子形状に関する検討が必要である。強度とダイレイタンシー挙動の関連を正しく理解することが重要である。

第21回中部地盤工学シンポジウム
21.落石挙動のばらつきを考慮した堆積層の衝撃吸収効果

μ（=tanφμ）をいくら高く設定しても，ある程度までのφfは表現できるがそれには限界があり，さらに高いφfの表現には，形状のモデル化が不可欠になることが分かる。このことから，Pm の不完全性を，Km の一部であるμ の変化で補完するには限界があることが分かる．ただし今回は，φf＝25(deg)程度で十分なため，粒子形状は円形でμ＝0.466 として解析を行っている。

●ボンド
第21回中部地盤工学シンポジウム
21.落石挙動のばらつきを考慮した堆積層の衝撃吸収効果

パラレルボンドの設定に必要なパラメータは，ボンドを付ける材料幅 rb、ボンドのばね係数 kb、引張強度 sbの 3 つである。ここで、ボンド強度sbについては，DEM による一軸圧縮試験を行い，そこで得られる強度と整合するよう設定している。
例えば、斜面を構成する岩盤強度は、試体サイズが 1(m)以下のような小さな供試体試験で得られる強度の 3 割程度まで低減することが報告されており，この寸法効果を考慮して最終的に sb を決定する必要がある。