個別要素法のパラメータ

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

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

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

先日、シュミットハンマーの反発度と反発係数を関連付けようとしていた文献を見かけました。この続きを読みたかったのですが、探せず。この先、どうなっているのか期待されます。
藤村ほか（１９９０） 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 を決定する必要がある。

 

 

landslide の低周波と高周波

以前、ある深層崩壊時の振動波形を見たときに、なぜ低周波が早く、高周波が遅れて発生しているのか疑問でした。

波動の知識不足もあり寝かせていたのですが、今日、似た事象を扱っている文献を見かけました。
Analysis of Broadband Seismic Recordings of Landslide Using Empirical Green's Function - Zhang - 2019 - Geophysical Research Letters - Wiley Online Library

すべり土塊が一体の時は低周波、それがバラバラになって高周波が発生するという解釈でした。根拠がわからないので利用できないのですが、そうなのかもと思える解釈です。

 

CIM への対応？

CIM 関連部署で話し合いが行われました。
CIM 対応に迫られているので、モノ・情報を継続して提供する立場でのお話でした。

一方、提供を受ける実務側のお話。
ある部署では３次元を避けてきたので、ノウハウも人もナシ。でも、対応に迫られている。モノも情報も欲しいけど、本当にすぐに欲しいのはできる人。ミスマッチです。

後者の executives が集まる会議では、以前より「人を育てないといけない」という意見で満場一致のようでした。でも、自分で手を動かすという意見は聞きません。自分は管理者であり、誰かに習得させる必要がある、でも対象の若手は限られている。困った。そうだ、できる人を雇おう！という発想なのでしょうね。総監の試験なら OUT です。（総監を持っていなくてもexecutivesになれるので、当然の帰結なのですが）。

決められたルールの中で効率よく物事を処理する能力は必要です。それよりも重要なのが、初めての事象に主体的に対応できる能力。俯瞰すると、前者に自信を持つ方々のみで議論を続けているように見えます。

わからなくても、手を動かしてみるべきでしょう。
多くの船頭のうち、何人かが気づくはずです。産む方が易かったと。

**************************************

20211003追記。

結局、雇える人がいなかったようです。 同時に手を動かしていたら、ある程度できるようになり方針が変わったかもしれません。

 

減衰の種類

減衰の種類について、復習です。
https://phreeqc.blogspot.com/2019/05/blog-post_29.html

 減衰定数（減衰比）h[無次元]
・臨界減衰係数㏄=2√(mk)で粘性減衰係数cを無次元化。h=c/cc=1では振動しない。
・理論的に求めることができない。観測に見合うよう導入。
・RC構造構造物：h=5%程度、鋼構造物：h=1～2%程度、地盤：h=5～30％程度

減衰(力)[N]の主な種類
・空気や液体の抵抗に起因する粘性減衰(viscous damping):F=cv、粘土の高い液体は2乗減衰:F∝v2
・構造・材料内部の摩擦に起因するクーロン減衰(coulomb damping)
・応力-ひずみループから求める履歴減衰(hysteric damping）
・構造物から地盤へ伝播する地下逸散減衰（radiation damping）

解析で利用される減衰
・ダッシュポット：粘性減衰係数マトリックスc[N・s/m=kgm/s2・s/m=kg/s]*速度ベクトル
・レーリー減衰：c[kg/s]=減衰係数α*質量マトリックス＋減衰係数β*剛性マトリックス
・複素剛性：構造減衰：iG*剛性マトリックス*変位

イマイチ、すっきりしていません。もう一歩なんですが。良い資料を見かけたら整理しましょう。

**************************************
20211218 追記
個別要素法での減衰を整理しました。少しすっきり。
https://phreeqc.blogspot.com/2021/12/dem_18.html

FLAC3D その3

FLAC3D の動的解析における境界条件です。

まずは「その2」の Local nonviscous damping を利用する方法↓。
A Local Damping Scheme For Nonreflecting Boundary Conditions In Surface Wave Modeling
https://www.earthdoc.org/content/papers/10.3997/2214-4609-pdb.190.sur03

次いで、Document の説明から2種。丁寧に説明されています。
Dynamic Modeling Considerations — FLAC3D 7.0 documentation (itascacg.com) 

動的解析の Quiet Boundaries（吸収境界）

In dynamic problems, however, such boundary conditions cause the reflection of outward propagating waves back into the model, and do not allow the necessary energy radiation. The use of a larger model can minimize the problem, since material damping will absorb most of the energy in the waves reflected from distant boundaries. However, this solution leads to a large computational burden. The alternative is to use quiet (or absorbing) boundaries. Several formulations have been proposed. The viscous boundary developed by Lysmer and Kuhlemeyer (1969) is used in FLAC3D. It is based on the use of independent dashpots in the normal and shear directions at the model boundaries. The method is almost completely effective at absorbing body waves approaching the boundary at angles of incidence greater than 30°. For lower angles of incidence, or for surface waves, there is still energy absorption, but it is not perfect. However, the scheme has the advantage that it operates in the time domain. Its effectiveness has been demonstrated in both finite-element and finite-difference models (Kunar et al. 1977). A variation of the technique proposed by White et al. (1977) is also widely used.

Dynamic analysis starts from some in-situ condition. If a fixed boundary is used while generating the static stress state, this boundary condition can be replaced by quiet boundaries; the boundary gridpoints will be freed, and the boundary reaction forces will be automatically calculated and maintained throughout the dynamic loading phase. However, changes in static loading during the dynamic phase should be avoided. For example, if a tunnel is excavated after quiet boundaries have been specified on the bottom boundary, the whole model will start to move upward. This is because the total gravity force no longer balances the total reaction force at the bottom that was calculated when the boundary was changed to a quiet one.

Free-Field Boundaries

The boundary conditions at the sides of the model must account for the free-field motion that would exist in the absence of the structure.
In order to apply the free-field boundary in FLAC3D, the model must be oriented such that the base is horizontal and its normal is in the direction of the z-axis, and the sides are vertical and their normals are in the direction of either the x- or y-axis. If the direction of propagation of the incident seismic waves is not vertical, then the coordinate axes can be rotated such that the z-axis coincides with the direction of propagation. In this case, gravity will act at an angle to the -axis, and a horizontal free surface will be inclined with respect to the model boundaries.

Vp,Vsの計算方法は SPH でも同様の計算をしていました。そうすると、初期定常の計算時から動的なポアソン比、剛性等を入力する必要があるのでしょう。静的な値で初期定常を行っていると、動的な値へ切り替えた際に定常状態が崩れます。

また、放射状の波を扱う際は、吸収境界が良さそうですね。張替えのサンプルはコチラ↓

http://docs.itascacg.com/flac3d700/flac3d/zone/test3d/Dynamic/ShearWaveFreeFieldBoundModel/ShearWaveFreeFieldBoundModel.html

; Set Initial Conditions
model gravity 10
; Set Boundary Conditions - Roller boundaries
zone face apply velocity-normal 0 range group 'West2' or 'East2'
zone face apply velocity-normal 0 range group 'Bottom'
; Static equilibrium
model dynamic active off
model solve convergence 1
; dynamic
zone face apply-remove range group 'Bottom'  ; Remove velocity boundary
zone face apply quiet range group 'Bottom' ; Add quiet boundary
zone dynamic free-field on ; Free field boundaries
; Solve to time 0.015
model dynamic active on
model solve time-total 0.015

FLAC3D その2

 FLAC3Dの続きです。

減衰は Local nonviscous damping が静的でデフォ。これ、初めて聞いた減衰だったのですが、古くから使われているようです。
↓の文献のように、動的では反射を吸収する境界として利用されているようです。手元のモデルで底に敷いたり外したり検討していたのですが、当たりでした。
A Local Damping Scheme For Nonreflecting Boundary Conditions In Surface Wave Modeling
https://www.earthdoc.org/content/papers/10.3997/2214-4609-pdb.190.sur03

The absorption coefficient α is a constant ranging from 0 to 1 and describes the degree of energy dissipation.
Therefore, local damping seems to be interpreted as energy dissipation by an internal friction mechanism. The magnitude of coefficient α in equation (2) is dimensionless. Therefore, it is independent of properties or boundary conditions and can be set to vary from one grid point to another
(Cundall, 1987).

As an empirical approach to find proper value of the coefficient, we compare results of the presented scheme for a given value of α with those obtained from viscous damping similar to Shin (1995) and from the spatial filtering scheme by Cerjan et al. (1985)

静的でもこの機械的減衰が必要なのは、動的な支配方程式を解いて収束結果を得ているため。大きめの減衰（静的では0.8がデフォ）を入れた方が早く収束するということです。

レーリー減衰については、FEM 等と扱いが同様です。
Dynamic Modeling Considerations — FLAC3D 7.0 documentation (itascacg.com) 

Thus, for many dynamic analyses that involve large-strain, only a minimal percentage of damping (e.g., 0.5%) may be required.

動的解析では、他にレーリー減衰、人口粘性、ヒステリティックな減衰が実装されていました。 前2つは Local nonviscous dampingと併用できるとありますが（実際、併用できていましたが）、最後は書かれていませんし試していません。

次は境界条件。

FLAC3D

Itasca 社 の FLAC3D の Document を読み直していました。

FLAC3D は 陽解法を用いた Lagrangian finite-volume プログラムです。
https://www.itascacg.com/software/flac3d 

陽的、Lagrangian。当然なのですが、計算の大まかな流れはSPH と同様でした。
http://docs.itascacg.com/flac3d700/flac3d/docproject/source/theory/theoreticalbackground/formulation.html?node2125

1. 速度（→速度勾配テンソル）→ひずみ勾配テンソル、スピンテンソルを算出
2. せん断剛性＋ひずみ速度テンソル→偏差応力テンソルの共回転微分（ヤウマン微分）＋スピンテンソル→コーシー応力のヤウマン応力速度を算出
3. 塑性化判定
4. 運動方程式を解いて加速度を算出。重力加速度、地震加速度、補正（レイリー減衰、人工粘性）→速度、変位を算出→ステップ1へ

ひずみ勾配テンソル等を求めるのに tetrahedral shape で 有限体積近似 を利用するのが FLAC3D、 粒子配置を利用するのが SPHと理解。文献では差分法と書かれていたのですが、FLAC（2D）のことでした。2D と3D で採用している空間の離散化手法が異なるようです。ただし、運動方程式の加速度の算出（時間方向）は3Dでも差分近似を利用しています。時間方向だけ差分というのは他でもよく見る手法です。

回転を考慮するヤウマン応力速度を利用するかどうかはスイッチで切り替え可能（large-strain mode）。SPH はタイムステップ毎に近傍粒子を更新しますが、FLAC3Dでは初期のコネクトを保存しており更新しません。ある程度 volume がつぶれると精度が悪化するのは変わらないようです。

陽解法なので全体剛性マトリックスを解く必要がありません。解自体は必ず得ることができます。1回の計算時間は少なくて済みますが、初期定常でも動的な挙動を解いて収束値を探しているため、案外、同じくらいになるのかも。並列化はしやすいと思います。

次は減衰。

スマート調査

昨年度から今年度にかけて、井戸調査を含む仕事が多く入っていたようです。

私は解析業務の管理に携わった程度ですが、現地担当は忙しそうでした。

台帳を見てみると、野帳に記入→EXCEL入力でした。20年以上前から変わっていません。GIS を利用していない方々でしたので、iPad 等での現地入力という選択肢をもっていなかったのかもしれません。

そういう私も、端末用のデジタルフォームを作っていないなあと思い、少し手を動かしてみました。
https://phreeqc.blogspot.com/2017/12/survey123-for-arcgis.html 

今回は、ArcGIS Survey123 Connect を利用。EXCELで項目を追加するだけで、入力フォーマットが作成されて行きます。非常に簡単です。

1～2時間ほどでプロトタイプが完成。

で、アップロード。 

ジオリファレンスは1個のみ、とエラーを吐かれたので、修正して再度アップロード。

実際にスマホで入力してみると、想像以上に使い勝手が良い（iPadよりも）！
1タップで選択できる項目を多く配置したので、入力が楽＆速い。これは後のフィルタリングで困らないためだったのですが、かなりの入力効率化につながりました。
リバースジオコーディングで位置情報から住所の入力も可能です（クレジットを消費します）。必須入力項目を指定しておけば、記入漏れを防ぐことができます。何より、帰社後の写真整理や位置図の貼り付けが不要になるのが嬉しい。しかも、入力したデータは既に GIS データになっています。これで帰社後の整理作業がほぼなくなりました。

入力の準備が終わったので、さあメンバーの追加！と思いきや、グループに現地調査員（仮）を追加できません。ん？

サポートに確認すると、デフォでは同種類のメンバーしか共有できないとのこと。組織アカウントでフォームを作成したなら、基本は同じ組織アカウントしかデータを入力・共有できない仕様。うーん、ここにきて使えない。ライセンス数分のアカウントを各班に割り当ててしまうと、他の GIS 作業がストップします。うーん。
不特定多数への完全公開は論外。残る手段は Field Worker ライセンスの購入だそうです。online では $350/yr で販売されていました。うーん、商売上手。私が担当なら買っちゃうでしょうね。

次は、入力したデータを含む台帳を印刷。
と思いましたが、ココでもクレジットを消費するそうです。これは VBA 等でスクリプトを組んだ方が自由に加工できるかも。今回はココまでかな。

他部署では、ソフトの作成・利用で現地調査の効率化・ペーパーレス化が進んでいます。設計部署でも AR を利用した現地踏査が始まっています。
「携帯端末を現地調査に使わない（使えない）」は、とっくに卒業すべき時期を迎えていたようです。

スポーツクライミング

オリンピックの新競技、スポーツクライミングをTVで観ています。

面白いですね。
難しい課題を登っていく選手。見ていてコチラまで腕が疲れそうになります。

三種目の順位の掛け算がポイントで少ないほうが優勝なんですね。知りませんでした。

競技人口、多くなるかな？

Today's Earth

https://www.eorc.jaxa.jp/water/index_j.html

"Today's Earth (TE)" は、JAXAが東京大学と共同で開発した陸上の水循環シミュレーションシステムです。河川流量や土壌水分量など、人間社会にとって極めて重要な水に関するデータの配布と可視化を通して、災害監視や水文学研究に役立てられています。 

全球と日本（詳細版）の2種類が web で公開されています。近年、陸水シミュレーションで全国を扱っている例を稀に見ますが、入力値に衛星‐データを利用しているのはココだけでは？さらに結果をオンタイムで web 上に公開されている点に驚き。先日の LiCSAR もそうですが、計算分野の処理技術には驚かされます。

データも公開されています。
wget で 14年分の daily データを３種類落としてみました。５時間程度かかりました（並列化しないといけないかな）。
データ形式は netCDF。データの並びについて説明がなかったのですが、これは気象業界の常識があるのかもしれません。
期待していたリスク値は含まれていませんでした。確率年等の数字は全国一律で表現できるので有用なのですが、残念。

1日1か所のデータを1行にまとめると、1年で5千万行を超えました。手元のPCはメモリーオーバーを吐きます。グリッド形式のまま扱うしかないかな？
大量のデータをオンタイムで提供していただけるようになったのはありがたいのですが、自身の update が追いついきません。新品に買い替えできたらよいのに。