やり残し事項 2021

今年は人生で1，2番目の悪い年でした。

年末のダメ押しが決定的で、再起動まで少し時間がかかりました。

仕事のストレスも相まって、余暇にPCを触る気力も失せてしまい、どうしようかと思案。
胃痛のため食欲に意識を振ることができず、長く寝られるわけでもない。
結果、PCを閉じ、何もしないことにしました。

といっても本当に何もしないわけではなく、大掃除をしながら体を動かし、静かな曲を聴いたり小説や漫画を読んだりしました（技術者でなければ、情報収集や勉強をしないというのも普通の生活なのかもしれません）。

救いは思いがけないところにあるもので、 漫画の主人公が人の死に思い悩んで立ち止まりかけた時に他のキャラが吐いたセリフ。これにハッとさせられました（我ながら安いものです）。

それを契機とし、PCを開き計算結果を確認したり、新たな計算を始められるまでには回復。胃痛は治らないものの、心は明確に切り替わりました。まだ何かやりたいことがあるわけではありませんが、手を止めないということは決定事項。後悔や反省は忘れることなく、考えながら2022年を進もうと思います。

やり残し事項です。

優先度中：流体＋個体（不連続体＋連続体）＋振動
優先度低：PSInSAR
優先度低：Dtransu の GPU 対応
優先度低：地表流＋地下水＋移流拡散

中期目標5年目は失敗。ま、成功しても結果的には意味がない状況に陥っていました。こちらは当分目標ナシです。しばらく考えましょう。

個別要素法の接触力モデル

Catherine O'Sullivan 原著「粒子個別要素法」3章より、接触力を表現するモデルです。
これまで見たり接してきたりしたモデルが、どのような分類に属するのかを理解できました。この図書、良いかも。


法線方向の接触モデル
・線形弾性ばね
・簡易ヘルツ接触モデル
　材料特性（G, ν）の導入
・降伏を考慮した法線接触モデル
　・ウォルトン―ブラウン線形モデル（半ラチェットばねモデル、履歴減衰モデル）
　・ばね、ダッシュポット（線形ケルビンモデル）：塑性変形によるエネルギー散逸を表現

接線力計算
・ミンドリン―デレシビッツ接触力モデル
　・簡略版：ブクオックモデル
　・ソーントン―イェンモデル（荷重履歴依存性の導入）：非線形、再載荷の考慮

引張り力伝達モデル
・並列結合モデル

転がり抵抗
・岩下―小田の回転抵抗モデル
・ジアンらの回転抵抗モデル

時間依存性
レオロジー：ばねとダッシュポッの組み合わせで表現
・バーガーモデル（ケルビンとマックスウェルの直列）：クリープの表現
・マックスウェルモデル（直列）：経時変化による劣化

 

個別要素法での減衰

個別要素法の減衰の扱いを整理。テキストは以下の２つ。
Catherine O'Sullivan 原著「粒子個別要素法」
古川ほか「個別要素法を用いた構造物の動的解析における減衰のモデル化に関する基礎的検討」
https://www.jstage.jst.go.jp/article/jscejseee/70/4/70_I_89/_pdf

F=Mx''+Cx'+Kx

1.要素間に作用する減衰力
F=Cδ'+Kδ
1-1.ダッシュポット(衝突エネルギーを消散)
・反発係数e←→減衰定数h←→ダッシュポットの粘性減衰係数C（衝突する粒子の質量・ばねに依存）

1-2.剛性比例型減衰：C=βK=(2h/ω)K
・速度項に乗じて力になる→単位：力/速度
・β固定（周波数ωのh固定）→周波数が高いほど減衰（定数）大
・単独で主要なモードの減衰を再現でき，剛性低下に応じた減衰力の低下も表現し，かつ高振動数領域の減衰性能も確保できる
・計算ステップを細かくする必要あり

2.要素重心に作用する減衰力
F=Ma+Cv
2-1.質量比例減衰(主要なモードの減衰)：C=αM=2hω
・α固定（周波数ωのh固定）→周波数が高いほど減衰（定数）小→１-1の併用で高周波数帯をカバー
・相対的に速度の大きな領域で誤った物体力→破壊モードに影響

2-2.局所比例減衰力(主要なモードの減衰)：-α|Fp|sign(vp)
・加速させようとする運動のみを減衰→誤った減衰力が定常運動から生じない
・速度に依存しないので、履歴減衰と類似
・無次元のため、周波数に依存しない
・係数αを調整することで主要な固有モードの減衰定数を表現→減衰定数は高々数％→１-1の併用で高周波数帯をカバー

コロナ対策のバランス

忘年会があったようです。

私の周りの方々は参加されませんでしたが、他のフロアの方々を中心にバス2台でホテルへ向かわれていました。

勤め先の支店では、緊急事態宣言中も、その前後のテレワークが推奨されている期間も、基本、テレワークなしの出社体制を貫かれてきました。この支店のみのようなので理由は不明ですが、おそらく地域的な業務特性？を優先された結果だと想像します。

と思いきや、今月、昼食のための会議室利用が禁止さました。これまで使われていた方々は文句を言われていましたが、これは理解できます。が、今頃？という意見にも同感。
リモートはダメで全員出社を強いるのに、昼の少人数の集合はダメ、忘年会のような大人数会食はOKで、少人数会食はダメというように見えてしまう方針は、批判されても仕方ないでしょう。

なんだか、バランスが悪いですね。

メンテナンス その２

 引き続き、週末は計算を回しながら細々したメンテナンスを継続。

・自転車のブレーキ交換
　シュー、ワイヤー交換。快調です。

・自動車のタイヤ交換
　冬タイヤへ。ハブの錆が少し進んでいたので清掃。

・自転車チューブ交換
　パンクしたので開けてみると、異なるサイズのチューブが入っていました。即交換。

・大型ごみ搬出
　ヒーター、ストーブ、スキー、スノボ等々。

・図書整理
　配置換えとともに、不要なものを資源ごみへ。

・室外機のホース交換
　配管テープが破れていたので新品に巻きなおし。ついでに割れていたホース交換。

・Switch Joy-Con ステック交換
　代用品の DUALSHOCK3 がダメになり、壊れていた純正品を修理。パーツ交換で復活。

あとはバッテリーが気になるくらい。まだもつかな？

今年は仕事が暇で、土日を完全に休めています。このような年末は初めて。おかげで今までサボっていた身の回りの整理がはかどっています。晴耕雨読ではありませんが、日中動いて、夜文献整理。既に大掃除にも半分終わり。年末年始はゆっくりできそうです。

 

個別要素法でのパッキング

初めて3次元で DEM を扱っています。

理論は連続体に比べて単純。小さなモデルは特に問題なく流せました。
が、数万～数十万の粒子を扱いだしてから、パッキングで躓きました。結果が初期状態に大きく依存するのです。パラメータよりも、影響が大きいのかもしれません。

パッキングについて具体的なノウハウが書かれている図書、文献は見当たりません。手元の資料で関連する記述があったのは以下の通り。

伯野元彦「破壊のシミュレーション」よりDEM のパッキング方法です。
a） 規則的に配置する方法
b) ランダムに配置する方法
1. 棄却法
2. 局所移動法
3. 割込法
4. 落下法
5. 成長法

Catherine O'Sullivan 原著「粒子個別要素法」より
1. ランダム生成（半径拡大法、段階圧縮法）
2. 構築法
3. 三角分割
4. 重力堆積法

高橋ほか「個別要素法によるアスファルト混合物解析の要素パッキングに関する一検討」
https://www.jstage.jst.go.jp/article/journalpe1996/2/0/2_0_229/_article/-char/ja/
規則的にパッキングする場合、配置条件によって変形に方向性が現われやすい。
規則的な配置では、単一粒子径が望ましい。
粒子分割の粗密、すなわち粒子半径も解析結果に影響する。その程度は要素モデルによって異なる。
必要とする精度を満足する範囲で最も大きい粒子径を選定するのがよい。
要素モデルの組成(バネ要素とダッシュポット要素の組み合わせ)と、その物理定数のみではなく、パッキング方法も含めて個別要素の集合体モデルを構築しなければならない。

 

パッキング状態によって結果が大きく異なる事象は、DEM 界では良く知られているようです。この、「何とでもなる」という自由度が辛い。地盤を DEM で適切に模擬するのは、連続体で初期応力状態をつくるよりもかなり難しいと感じます。選択するパッキング方法、その組み合わせ、ノウハウは技術者の経験で培われているといえるでしょう。
が、逆に以下のような考え方もあるようです。

地盤工学のための個別要素法_7_パラメータの設定と土と地盤の作り方(その2)
「その作業ができるところがDEM解析の強みでもある。解析者はその地盤作成の深淵を感じながらも、速く、効率的に解析対象となる供試体や地盤をつくる必要がある。」

経験を積んで深淵を俯瞰できるようになりたいものです。

**************************************
20211217
粒子個別要素法の記述を追記。

PSInSAR ＆ Landslide その2

こちらは地すべり分布図と PSInSAR の組み合わせ。
前者はDEM＋RF。その精度を高めるのに後者を使用。

Andrea Ciampalini et al. (2016)
Landslide susceptibility map refinement using PSInSAR data
https://www.sciencedirect.com/science/article/pii/S0034425716302759

Landslide Susceptibility Map
The landslide susceptibility map (LSM) of Messina Province was produced by implementing a Random Forests (RF) algorithm in Matlab (Catani et al., 2013).
To avoid subjectivity in the choice of explanatory variables, several parameters were considered among the DEM-derived products: Aspect, Planar Curvature, Profile Curvature, Curvature s.s., Elevation, Flow Accumulation, Topographic Wetness Index (TWI), Log Flow Accumulation and Slope.
The choice of parameters depends on the map unit resolution (MUR) (Catani et al., 2013). Among the possible parameters, we used those suggested for a map unit resolution of 100 m.
Using this MUR, the expected Area Under the Curve (AUC) value is 0.81. A training set was created by randomly selecting 10% of the landslide database.
The DEM had a 20 by 20 m spatial resolution. In the LSM, the average value inside a 100 by 100 m cell was calculated.
Each pixel was classified using four susceptibility classes:
(i) low to null (0–0.3)
(ii) moderate (0.3–0.55)
(iii) high (0.55–0.75)
(iv) very high (0.75–1)

PSInSAR
SqueeSAR: TRE ALTAMIRA, detection of millimetre surface displacements, improving previous PSInSAR algorithm.
Xband COSMO-SkyMed (CSK) in Stripmap mode (40 × 40 km in the range and azimuth directions and 3 × 3 m ground resolution) in both ascending and descending geometry

Integration
The VSlope intervals were determined based on the standard deviation (σ = 7 mm/year) of VSlope for the whole PSI data set after merging the ascending and descending data sets. The matrix in Table 3 shows the correction factors for each considered case. The velocity intervals were determined to increase the susceptibility degree from level 1 to 4. For example, an area characterized by a susceptibility degree of 1 will increase by 1 degree if its VSlope falls within 1 σ to 2 σ of the stability threshold (7–14 in this case), 2 degrees if the VSlope is between 2 σ and 3 σ and so on. The higher the ground deformation velocity, the higher the susceptibility degree. A cell characterized by the maximum susceptibility degree cannot increase further.

Discussion
To improve the PS/DS detection, L-band SAR sensors can be used because their wavelength (30–15 cm) is more suitable to investigate densely vegetated areas.
Both the landslide susceptibility map and the ground deformation velocity map can be used in different ways to forecast landslide occurrence in a selected area. Combining them may improve the reliability associated with predicting these types of phenomena, particularly for slow moving landslides

PSInSAR ＆ Landslide

PSInSAR の活用のお話。
技術者が時間を割くのは PSInSAR 自体ではなく、その結果をどう問題に応用するかというところ。活用例として見本になるような文献でしたので要点を( ..)φメモメモ。

Kamila Pawluszek-Filipiak et al. (2021)
Multi-temporal landslide activity investigation by spaceborne SAR interferometry: The case study of the Polish Carpathians
https://www.sciencedirect.com/science/article/pii/S2352938521001658

However, due to the high temporal decorrelation resulting from vegetative cover and short wavelengths (X-band), it was only partially successful (Perski et al., 2009; Perski et al. (2011). This is mostly related to the low PS density due to the application of TerraSAR-X data. Therefore, the exploitation of C-band (Sentinel-1) and L-band (ALOS PALSAR) data can present more advantages, especially in rural areas.

The PSInSAR approach is based on the following steps:
(1)interferogram generation with respect to a common master image
(2) PS candidate selection based on amplitude dispersion index
(3) first estimation of atmospheric phase screen (atmospheric influence) and topographical and displacement components
(4) second estimation of interferometric components and final PS points selection
For a precise description of the algorithms applied within each step, see Ferretti et al. (2000; 2001).

The interferograms are then generated for each slave image always using the same master image.
Amplitude dispersion index of 0.25 was utilized as a threshold for PS candidate selection, and the final PSs were extracted with a coherence value higher than 0.75.

SARscape® software was used for SAR data processing (Sahraoui et al., 2006).

Persistent scattering (PS) postprocessing
R index: > 0.33
cosβ: >0.3
Moreover, we discarded PS points which showed Vslope>0 because positives values represent uphill movements and it is not representative for small landslide movements, even though positive values exist within landslides, especially in the toe area.

three categories based on PSI velocities :negligible, extremely slow, and very slow

It is worth emphasizing that the reality of the PSInSAR velocity estimation itself is not assessed due to the lack of external field measurements, but this measurement may be related to landslide activity. For this reason, field verification was performed.

In general, applying the PSI approach in rural areas is challenging. However, the results of this study prove that increasing the temporal sampling rate in view of launching the second Sentinel 1 satellite, provides a higher PS density and, therefore, this technique can deliver useful information for landslide activity assessment.
From another point of view, the application of two Sentinel datasets from ascending and descending orbits, using only one Sentinel satellite (C band), allows us to achieve a similar PS density in comparison to the one ascending dataset of ALOS PALSAR (L-band).

SAR の利用

勤め先でも、ようやくSAR の利用を考え 始めたようです。

といっても、担当はまだ本気ではないようで、
「使うための問題点や課題をこれから抽出します」
「社員の現状はこうでした」
といった報告をされていました。

SAR に関する私の要望はPSI一択です。実務なら数百万のソフトを買っていただく必要があります。手法の区別がつかない executives に承認権があるため、技術者も含めたそれらの方々の変化を待っている状態です。が、まだ時間がかかりそう。残念。

土木分野で見かける手法は以下の3種でしょうか。

強度画像（振幅のみ利用）
pros：RGB画像の合成のみなので、技術力や知識を必要としない。簡単。
cons：変位量が得られない。

DInSAR(differential interferometry SAR)
pros：広域の変位情報が得られる。
cons：土木構造物に対しては空間解像度が粗い。

MTInSAR (multi-temporal interferometric SAR)
SBAS (small baseline subset)
pros：広域の変位情報が得られる。変位の時間変化を把握できる。
cons：土木構造物に対しては空間解像度が粗い。
PSInSAR (Persistent Scatterer Interferometry)
pros：ピンポイントの変位情報を高精度で得られる。変位の時間変化を把握できる。
cons：欲しい位置の情報が必ずしもえられるとは限らない。

その他、TRE ALTAMIRA 社の SqueeSAR の利用を見かけます。精度が良いと謳われていますが、内容を知りません。パスコさんはじめ、国内でも提携する会社があるようですね。今後、報告が増えてくるでしょう。 楽しみです。

GMT

この年になると、新たなソフトの使い方を覚えるのが億劫になってきます。
できるだけシンプルに良いツールのみを使い続けられると楽なのですが、時間と人は止まってくれません。

今回は GMT。
OpenSWPCの結果が NetCDF で提供されており、図化に何を使用するか探っていました。

ArcGIS pro 2.8 から NetCDF を Basic でも扱えるようになりましたので、まずはこれで。
はい、読めませんでした。XYだと読めましたが、緯度経度ではダメ。後で中身を確認すると、後者の並びがダメでした。

で、Python で geotiff に変換。
これは Arc で読めましたが、全時間ステップを読み込んで動画を作成するのは手間でした。

Arc に持ち込まず、Folium 上にてプロットしようかと考えましたが、数が多いので却下。
https://phreeqc.blogspot.com/2020/02/folium.html

最終的には、Forium のかわりに PyGMT を使用。
地図、観測局、波形（NetCDF）、テキスト等を重ね合わせて図化。それを動画に変換して保存。計算後にスクリプトを動かすと動画が作成されるようにしました。
PyGMT は初めてだったのですが、Python っぽい書き方のためか違和感なく作業ができました。

今回は Ubuntu で環境を構築しましたが、Win 版も配布されていました。一昔前だと、Cygwin が必要だったように思います。ユーザーは増えているのかな？

得られた結果はイマイチ。
カラーマップを調整して何度かトライしましたがダメ。
OpenSWPC のソースを見ると、RGB の前2つに P,S が割り当てられていました。うーん、最初に見ておけば良かった。

差分近似の高精度化

微分の差分近似は、テイラー展開後に目的の項について整理するのみです。

簡単に高次精度化ができるようですが、今まで気にしたことがありませんでした。
反省も込めて、整理です。

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

1次精度‐1階微分
前進差分（テイラー展開）
\begin{align*}f\left(x+\mathrm{\Delta x}\right)=f\left(x\right)+\frac{\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x\right)}{\mathrm{\Delta x}}-\frac{{\mathrm{\Delta xf}}^{\prime\prime}\left(x\right)}{2!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x\right)}{\mathrm{\Delta x}}+O\left(\mathrm{\Delta x}\right)\end{align*}

 後退差分（テイラー展開）
\begin{align*}f\left(x-\mathrm{\Delta x}\right)=f\left(x\right)+\frac{-\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{-{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{{-\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x\right)-f\left(x-\mathrm{\Delta x}\right)}{\mathrm{\Delta x}}+O\left(\mathrm{\Delta x}\right)\end{align*} 

２次精度‐1階微分
中心差分（前進差分‐後退差分）
\begin{align*}f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)=\frac{2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{2\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)}{2\mathrm{\Delta x}}-\frac{{\mathrm{\Delta x}}^2f^{\prime\prime\prime}\left(x\right)}{3!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)}{2\mathrm{\Delta x}}+O\left({\mathrm{\Delta x}}^2\right)\end{align*}

２次精度‐2階微分
中心差分（前進差分＋後退差分）
\begin{align*}f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)=2f\left(x\right)+\frac{{2\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2{\mathrm{\Delta x}}^6f^{(6)}\left(x\right)}{6!}+\ldots\end{align*} \begin{align*}f^{\prime\prime}\left(x\right)=\frac{f\left(x+\mathrm{\Delta x}\right)-2f\left(x\right)+f\left(x-\mathrm{\Delta x}\right)}{{\mathrm{\Delta x}}^2}-\frac{{2\mathrm{\Delta x}}^2f^{(4)}\left(x\right)}{4!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+\mathrm{\Delta x}\right)-2f\left(x\right)+f\left(x-\mathrm{\Delta x}\right)}{{\mathrm{\Delta x}}^2}+O\left({\mathrm{\Delta x}}^2\right)\end{align*} 

2次精度‐1階微分
中心差分（2Δx）
\begin{align*}f\left(x+2\mathrm{\Delta x}\right)=f\left(x\right)+\frac{2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^2{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2^3{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2^4{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2^5{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f\left(x-2\mathrm{\Delta x}\right)=f\left(x\right)+\frac{-2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^2{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{-2^3{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{2^4\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{{-2}^5{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)=\frac{2^2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^4{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2^6{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)}{4\mathrm{\Delta x}}-\frac{{2^2\mathrm{\Delta x}}^2f^{\prime\prime\prime}\left(x\right)}{3!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)}{4\mathrm{\Delta x}}+O\left({\mathrm{\Delta x}}^2\right)\end{align*} 

4次精度‐1階微分
中心差分（Δx, 2Δx）
\begin{align*}f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)=\frac{2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)=\frac{2^2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^4{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2^6{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}8\left\{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)\right\}-\left\{f\left(x+2\mathrm{\Delta x}\right)-\ f\left(x-2\mathrm{\Delta x}\right)\right\}=12\mathrm{\Delta x}f^\prime\left(x\right)+O\left({\mathrm{\Delta x}}^5\right)\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{8\left\{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)\right\}-\left\{f\left(x+2\mathrm{\Delta x}\right)-\ f\left(x-2\mathrm{\Delta x}\right)\right\}}{12\mathrm{\Delta x}}+O\left({\mathrm{\Delta x}}^4\right)\end{align*} 

4次精度‐2階微分
中心差分（Δx, 2Δx）
\begin{align*}f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)=2f\left(x\right)+\frac{{2\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2{\mathrm{\Delta x}}^6f^{(6)}\left(x\right)}{6!}+\ldots\end{align*} \begin{align*}f\left(x+2\mathrm{\Delta x}\right)+f\left(x-2\mathrm{\Delta x}\right)=2f\left(x\right)+\frac{2^3{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2^5{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2^7{\mathrm{\Delta x}}^6f^{(6)}\left(x\right)}{6!}+\ldots\end{align*} \begin{align*}16\left\{f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)\right\}-\left\{f\left(x+2\mathrm{\Delta x}\right)+f\left(x-2\mathrm{\Delta x}\right)\right\}=30f\left(x\right)+\frac{24\mathrm{\Delta x}f^{\prime\prime}\left(x\right)}{2!}+O\left({\mathrm{\Delta x}}^6\right)\end{align*} \begin{align*}f^{\prime\prime}\left(x\right)=\frac{16\left\{f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)\right\}-30f\left(x\right)-\left\{f\left(x+2\mathrm{\Delta x}\right)+f\left(x-2\mathrm{\Delta x}\right)\right\}}{12{\mathrm{\Delta x}}^2}+O\left({\mathrm{\Delta x}}^4\right)\end{align*} 

 

LSTM の x_train データ

Keras ＋ TensorFlow で LSTM と Seq2Seq。

どちらも与えるデータの形は3次元。
[num_samples, num_times, num_features]

画像を与えるのと同じ感覚ですね。

Sony NNC では num_samples 個の2Dデータとなり、保存するのも削除するのも一苦労でした。が、２次元データを加工して３次元にするのであればハンドリングが容易。

ありがたい。

地下水観測手法と空間スケール

もう一つ、惹かれた文献。

Monitoring ground water storage at mesoscale using seismic noise 30 years of continuous observation and thermo-elastic and hydrological modeling
https://www.nature.com/articles/s41598-017-14468-9

地下水観測手法と対応する空間スケール

small
・traditional monitoring (piezometer networks)
medium (kilometric to decakilometric scale)
・continuous seismic monitoring (ambient seismic noise)
large (>200 000km²) or subsiding areas
・the GRACE gravimetric mission
・InSAR methods

内容

• Gräfenberg における30年間の地震観測結果（ambient seismic noise）を利用し、地下の力学的変化を検討。
• 1.2～10秒の周波数帯における速度変化を処理することで、年スケールでのレーリー波速度変化を0.01％オーダーでとらえた。
• 観測された速度変化は、局所的な帯水層（静水圧）、または地下全体での変化（熱弾性）によって生成される。
• dv/vの長期的な変化について熱弾性と静水圧のみを考慮するという限定的な近似条件の下で、両者の相対的な貢献度を推定。熱弾性効果を50％、静水圧効果を50％考慮した場合、観測された dv/v との間に最も高い相関係数（0.83）が得られた（Fig.5A）。
• 観測局グループ間の速度変化（Fig.3）は、水文学的な挙動の空間的なばらつきによるものであり、熱弾性的な影響はより均質であると予想される。
• 地質学的に安定しているこの地域では、熱弾性と水文学の両方の影響を、年単位の時間スケールと㎞単位の空間スケールでフォワードモデル化可能。
• 将来的には、深さ方向の分解能向上のため異なる周波数を用いて地下の4次元画像を解析することや、高密度地震アレイを使用して横方向の変動をとらえることが期待される。

GRACEは知りませんでした。
SAR の利用は最近聞いたり読んだりしました。海外では既に浸透しているのかもしれません。国内では山林が多いため、都市部を除き適用が難しい手法でしょうか。
微動を用いた観測手法は国内でも使えそう。というか、世界に誇るべき高密度の観測網が整備されていますので、積極的に研究に利用すべきでしょうね。F-net や Hi‐net のデータを取り入れて速度変化を表示・保存するシステム。技術的にハードルは高くないでしょう。サブダクションの影響がどこまで出てくるか？ん？影響が大きければ、そちら方面に使えるかな？

無償で公開されている地震波形や崩壊データは、近年、海外の研究者にも利用されています。日本の研究者が牽引する形で頑張っていただきたいところです。

 

東北地方太平洋沖地震後の速度低下

ambient noise を利用した cross correlation にて観測局間の表面波速度を求めている事例は海外の文献で多く見かけます。国内でもたまに見ますが、それらはすべて広帯域の低周波を利用しています。

昨日、はじめて Hi-net を利用した文献を見ました。
しかも、内容が興味深い。

Mapping pressurized volcanic fluids from induced crustal seismic velocity drops (science.org)

• 2011東北地方太平洋沖地震前後で地震波速度変化が認められた（Hi‐net、cross correlation で半年前と比較）。
• 火山地帯での低下率が大きい。
• 動的ひずみと動的応力を、PGV（KiK‐net）、レイリー波の平均位相速度（〜3km/s）、平均せん断弾性率（～30×109Pa）より求めた。
• 地震波による動的ひずみは、static coseismic strain （静的地震時ひずみ？）よりも1～2桁大きい。
• 速度低下率と動的応力の比を seismic velocity susceptibility （地震速度の感度）と定義。この分布には主要な火山地帯の分布との相関性が認められた。
• 地震速度の感度は、火山地帯における熱水やマグマ流体の加圧レベルの指標として使用できる 。
  
• 火山地帯、特に（静的応力の変化を無視できる）富士山麓での大きな速度低下の主な原因は、東北地方太平洋沖地震による動的応力であると結論づけられる。

使用しているのは Hi-net のデータだけなので、オンタイムでデータを受信し、日々の snapshot を残しておけば、地震の発生毎に全国レベルで速度変化を見比べることが可能かつ容易になります。水位変化や応力変化の把握にも使えるかな？

海外産ツールは多く出ていますが、トモグラフィーまで完走できていません。
https://phreeqc.blogspot.com/2021/05/passive.html
動かせるようにしておかないと。

********** 関連 *********
Static strain and stress changes in eastern Japan due to　the 2011 off the Pacific coast of Tohoku Earthquake, as　derived from GPS data
https://earth-planets-space.springeropen.com/articles/10.5047/eps.2011.06.049

Rockfall radar

正に ’cool’ という単語が当てはまるシステムです。
https://www.youtube.com/watch?v=o-oVXYkBwgw

使用されているのはスイスの GEOPREVENT 社の商品のようです。
Rockfall radar 
https://www.geopraevent.ch/technologies/rockfall-radar/?lang=en
・all-weather suitability (works in rain, fog, snowfall)
・works day and night

他にも魅力的な商品が多くあります。
Avalanche radar
https://www.geopraevent.ch/technologies/avalanche-radar/?lang=en
これは落石とアルゴリズムを変えているだけでしょうか。

GEORADAR
地上SARは国内でも実績ありですね。親会社？の商品でしょうか。

Deformation camera
https://www.geopraevent.ch/technologies/deformation-camera/?lang=en
カメラ画像から1日の動きをcmオーダーで判定。これは昼のみでしょう。
画像を用いた異常検知では機械学習を利用した文献を見かけます。学習してデプロイしてしまえばオンタイムで判定可能です。精度向上とともに追っかけの商品が出てくるでしょうね。

 
国内のシステムは子供のように見えます。Rockfall radarなど、一度は使ってみたい商品です。

メンテナンス

今週末は計算を回りっぱなしにしておいて、細々したメンテナンスを片づける。

・Microsoft Natural Ergonomic Keyboard 4000 パームレスト張替え
・自転車1のライト、ブレーキシュー、ブレーキワイヤー、Vブレーキパーツ手配
・自転車2のキズ消し
・自動車のフロントガラス　鱗取り

本皮のグレーの端切れを購入し、型を取って切ったところで表裏逆に気付く。ま、良いかとそのまま続行。なかなかの出来栄え。痛んだら、次は茶系にしましょう。

何を消そうと考えていたのか思い出せないほど前に購入していたソフト99「ツヤ仕上げキズ消しセット」。先週、道具を散らかしていたら出てきました。保管していた折り畳み自転車（白）に、拭くだけでは取れない黒いスリ汚れが多数あったことを思い出し、使用。初めて使いましたが簡単でよく落ちました。

来週末は別の自転車のブレーキ回り交換。シフト周りは調整のみでまだ持つかな。そうこうしていると冬タイヤが交換の季節。そのうち、大掃除、年末へ。
頑張ろう。

地表の地震動予測法

技法堂出版「入門数理地震工学」より

地表の地震動予測法は4つに大別される。

（１）経験的方法

• 観測地震動の最大加速度値等と地震規模、震源距離、地盤条件等の統計処理により地震動の特性値を予測
• 最大値等の震動特性値を簡単に評価できる
• 波形は得られない
• 水平多層弾性体にのみ適用可
• 逆解析による断層・地層構造の推定可

（２）半経験的方法（経験的・統計的グリーン関数法）

• 経験的：余震記録が伝播経路・サイト特性を表すと考え、観測記録を用いて大地震の波形を予測
• 統計的：観測波形の代わりに合成波形を利用
• 高精度
• 短周期(2Hz～)

（３）理論的方法（FDM,　FEM, 剛性行列法）

• 震源・伝播経路・サイト特性の全てをモデル化し、波動方程式から波形を予測
• 震源断層モデル：運動学的断層・動力学的断層
• 地盤モデル：水平多層弾性体・不正形多層弾性体・不均質不正形多層弾性体
• 長周期(～2Hz)

（４）ハイブリッド法（半経験的方法＋理論的方法）

• 短周期帯域を統計的グリーン関数法+長周期帯域を理論的方法（例えば剛性行列法）を用いて波形を予測
• 短周期～長周期をカバー

 

Ubuntu dist-upgrade

Ubuntu は手元に3台あるのですが、2台は 20.04、1台は18.04。

1台だけ環境が異なると支障の出る場合があり、最後の1台も上げましょうと重い腰を上げました。
「Software Updater」で 「upgrade」 をポチ、で、できるはずだったのですがボタンを押しても無反応。

$ sudo apt dist-upgrade
0 upgraded, 0 newly installed, 0 to remove and 1 not upgraded.
cuda-driver が引っ掛かっていました。

$ sudo apt install cuda-drivers
はい、確かにインストール不可→対処↓
https://forums.developer.nvidia.com/t/cuda-install-unmet-dependencies-cuda-depends-cuda-10-0-10-0-130-but-it-is-not-going-to-be-installed/66488/16
$ sudo apt clean
$ sudo apt update
$ sudo apt purge cuda-drivers
$ sudo apt purge nvidia-*
$ sudo apt autoremove
$ sudo apt install cuda-drivers

nvidia-* は迷いましたが、また入れたら良いでしょうと考え enter⏎。
これで完了。

最後に「Software Updater」 で upgrade ボタンを押すと反応しました。

Amplitude of the FTF

lsforce で距離を設定すると土量が少なく推定される原因２つ
https://pubs.geoscienceworld.org/ssa/srl/article-abstract/92/4/2610/596477/lsforce-A-Python-Based-Single-Force-Seismic?redirectedFrom=fulltext

The amplitudes of the FTF produced by lsforce will be lower than their true magnitudes for two principal reasons: the bandlimited nature of the input waveforms and the effects of the Tikhonov regularization scheme.

前者のみでは大きいなあと感じていたところです。
正則化でもそれぞれ影響があるといわれると、確かにそうかもしれないという程度の理解ですが、言われると納得できます。

PCを使う？

先日、以下のようなやり取りがありました。

私の質問

• 新しい共有PCを2か月くらい使用しようと考えている。
• 1日1回程度、計算状況を確認するだけで、席は空いている。
• そちらの浸透流解析で何コア使うのか？メモリーはどの程度使う予定か？影響しないように開けておく。

他部署の先輩の回答

• 浸透流用に導入予定のソフトが非力な PC では動かないケースを懸念している。
• 新しいPC はハイスペックなのでそれを利用する予定だった。それが使えなくなると非常に困る。

先輩にとっては難しい質問だったようです。知らないことには答えない、という典型的な回答でした。
どうも浸透流は初めてのようで、支配方程式や計算負荷を御存じない様子。そのためか計算はプロに任せ、結果をもらってソフトで表示して、報告書にまとめるだけのようでした。それなら、コアは1個でメモリも数㎇でしょう。ソフトの要件が「メモリ8GB以上」でしたのでどのような PC でも動くでしょうし、「1GB以上のビデオメモリを実装したOpen GL 4.5対応のグラフィックカード」ですので、古い Quadro で問題ないと思われます（新しいPC は GeForce です）。
まったく心配ないのですが、初めての知らないことで御心配だったのでしょう。2，3回やり取りしましたが「使ってくれるな」を繰り返されるだけでしたので、先輩ということで理詰めを控え、引きました（若い方からすると、こういうのを老害と思われるのでしょうか）。

ヒトは誰しも知らないことについて構えてしまうものですが、過大に備えて周りに迷惑をかけることは避けたいものです。今回、それを避けるには PC やソフトを知ることでした。普段から使用している道具の知識を update しておけば、周りに迷惑をかけることはなかったと思われます。無知は罪です。

道具に使われることのないよう、他山の石として肝に銘じておきましょう。

**************************************
ちなみに、このソフトの浸透流ソルバーは私が部分改変してコンパイルしています。Win11になればコンパイルしなおすべきなのですが、御存じないのでしょうね。ま、計算は御自分でされないみたいなので問題があること自体に気付かれなくても支障ないのですが。

20220123追記
ハイスペックPCを確認すると32コア。動かししたいと言われていたソフトは並列化されていません。計算されないものの、ソルバーも並列化非対応。
情報を仕入れず、知識をアップデートしないまま仕事をされると、迷惑なのは周りとお客様です。無知は罪です。

波の伝播

振動の伝播を解こうとして試行錯誤しています。

最初は簡単なモデルで試そうと考え、数十㎞の矩形領域の表面に点震源を与えました。

まずは CTC さんの SoilPlus でモデル化。しかし、計算しだすと完走までおよそ数か月かかることが判明。

これはダメだということで、Misdas さんの FEA NX に切り替え。ソルバーが GPU 対応です（なぜか GeForce 推奨）。
モデルを作って境界条件を張る際に気づいたのですが、側方の粘性境界が実装されていませんでした。伝達境界を張ろうとしましたが、コチラは３次元に対応していないようで晴れません。モデルを作る前に確認しておけばよかった。

で、FLAC 。
こちらは、node に 荷重を与えることが難しいようです。Fish で組み始めましたが、未だ完成せず。

最後にOpenSWPC。
これは以前、MPI が引っかかりコンパイルできなかったのですが、Ver5.20 になっていたためか今回はすんなりと通りました。サンプルを動かして動作確認OK。
実体力（物体力）を時系列で入力するのは左右対称のモーメント時間関数を利用し、点震源を重ね合わせることで対応できそうです。2D の例題で試してみました。

 

Boxcarで10秒間一定の力をかけた場合

 

Cosineで２秒毎（１秒間のラップ）の重ね合わせで10秒間

一定の実体力（Boxcar）をCosineで再現したのですが、内側の違い（９～10秒の関数形の相違）を除けば、一致しているるように見えます。振幅の比較は sac ファイルを確認ですね。ま、最終的にはコードの確認が必要でしょう。

こういった波の伝播に関しては、地震分野のプロがいくつかのコードを公開されていました。が、いずれも震源は断層をモデル化する方式で、時系列の力の入力はサポートされていません（当前なのですが）。OpenSWPC における重ね合わせの利用も、想定されていないでしょうね。

今のところ、FLAC が遅ければ OpenSWPC しか選択肢が残りません。これが正しく機能しないとわかったら、次に行きましょう。

Kratos

セミナーで 紹介されていた Kratos を試してみました。
https://github.com/KratosMultiphysics/Kratos

GiD を1か月フリーで使用できるとのことでしたので、それを利用。
Kratos に限らず、OpenFOAM なども扱えるようです。一通りのジオメトリファイルをインポート可能で、メッシュも nastran などに対応しているようでした。

MPM に関しては、Example に2D が含まれていました。
メッシュを切って run するだけでしたが、？でした。メッシュを切った時点で、解析領域と円の領域の接点が一致していません。円の方が粒子の発生させる領域なのかな？と思っていましたが、メッシュを切っていますし、よくわかりません。

うーん、どういうことなのでしょう。
うーん、寝かせておきましょう。時間が解決してくれそうな気がします。

Unless Making Progress

今日、セミナー動画を見ていると MPM が出てきました。
国内では流行っていないなあ、と思いながら検索。すぐに清水さん、大林さんの文献が引っ掛かりました。掲載が2012年ですから、意外と古い。私自身も2014年に書き残していますね。https://phreeqc.blogspot.com/2014/11/mpm-material-point-method.html

先日、他社さんの資料を見ていると、プレッシャーメーターからcφを求めている方がいらっしゃいました。これ、基準書に乗っているのですが、実務で報告さているのを初めて見ました。
https://phreeqc.blogspot.com/2011/01/blog-post_8311.html

 

前者の流行らなかった一因は PC スペックだろうと思われます。が、それも今は改善しつつあります。
後者は思考の停滞でしょう。組織慣性が働いているのかもしれません。

先達が示された通り、現状維持は後退と同義です。変化しないと進歩しません。
次の10年、新たなことに挑戦し続ける人が増えたら良いと期待します。

AMD + ifort

他支店で、AMD の32コアを積んだ PC を購入したとのこと。

まだ GPU に対応しておらず CPU 並列化のみのソフトは多いので、一つの選択肢かもしれません。が、OS が Windows10 のみ。以前コンパイルしたソフトを動かそうとしたが動かず困っているという連絡でした。

AMD Ryzen は手元に実機がありません。AMD + UBUNTU での並列化はよく聞くのですが、Winはナシ。手探りで始めることに。

エラー内容を見ると、どうやら Intel CPU用の命令が邪魔をしているようでしたので、まずそれを外すことに。そのままでは変更するオプションが多いので、コマンド打ちで０からオプションを追加する方針にしました。

何度か試行して、スタック領域を拡張し忘れていたことに気づき、追加。これで動きました。

今回、ifort ではコンパイルできたのですが、gfortran ではダメでした。何が違うのかわかりませんが、前者の速度と安定感は以前より感じています。

このあたり、いつかプロに御教授願いたいところです。

 

Animated Graph

 
SPH の動画を見ていて、ふと、グラフをどのように動かしているのか？と疑問に。

 

何で作られたのだろうと調べると、Matplotlib にその機能がありました（灯台下暗し）。

from matplotlib import animation
interval=100/3 #300 frames * interval 100/3ms = 10s

ani = animation.ArtistAnimation(fig, frames, interval,repeat=False)
ani.save('./Force.mp4',writer="ffmpeg",dpi=300)

for loop でグラフの作画を進めた後（frames リストに格納した後）、 animation.ArtistAnimation で集約。ani.save でmp4 として保存しました。

SPH の計算結果を動画にするのは blender。と思っていたのですが、これが落ちまくり。png で書き出してから動画にすることにしました。
apng だと PowerPointが対応していないため、これもm4へ。cv2.VideoWriter でOK。30 fps で10秒の動画とし、グラフと合わせます。

できた2つの mp4 を PowerPointに貼り付けたら、同時にスタート。
できました。

 

サンプリング周波数

空間周波数について調べていたのですが、サンプリング周波数として欲しい内容が書かれていました。空間周波数も同様ですのでメモしておきます。

https://www.keyence.co.jp/ss/products/recorder/lab/voltage/point.jsp

サンプリング周波数÷信号波形周波数	記録波形への影響
10倍以上	正確な波形表示・記録が可能
2倍以上10倍未満	波高値が小さくなります。波形形状が荒くなります。
2倍未満	エイリアシングが発生し、実際の波形とまったく異なる波形を表示。

関連→https://phreeqc.blogspot.com/2020/10/2_13.html

GEORAMA でボクセルメッシュ

GEORAMA v3.3
ソルバーの方の GEORAMA です。使用するのは初めてです。

複数の csv 標高データを SoilPlus へ取り込もうとしました。が、そのままでは薄層が多すぎてメッシングが困難になることが目に見えています。
で、事前に Civil3D＋GRORAMA でモデリングしようとしたのですが、なぜかLandXML を境界面として選択できませんでした。
結果、csv 読み込みが可能な GEORAMA 単体を選択。

メッシュの大きさに比べ薄い層が多くある場合、ソリッドやメッシュを切ることは難しいのが一般的です。この場合、割り切ってボクセル出力を選択せざるを得ません（解析の種類にもよりますが）。GEORAMA ならボクセルの重心付近の地質をあてはめたメッシュを自動で書き出してくれます。
また、異なる層が同標高となっている場合、他のソフトではある程度ダミーの層厚を持たせる配慮をしないと次の解析に持っていけませんでした。GEORAMA のボクセルなら、これも気にしなくて済みます。割り切ると便利です。

750,000メッシュの作成に1時間くらいでしょうか。その後は nastran 形式で作成し、Soilplus で読み込み（これも時間がかかります）。

無事にメッシュを作成できました。

***********************
20211019追記

ボクセルメッシュは重心位置の地質が適用される仕様だそうです（サポートさん談）。メッシュ内での各地質の体積は考慮されていません。残念です。

DesignSPHysics

セミナーで紹介されていた DesignSPHysics。
https://design.sphysics.org/

以前より公開されていることを知っていましたが、使っていませんでした。
今回、その機能を体験したのですが、これがなかなか素晴らしい。FreeCAD に 組み込むマクロ（Python）ですが、よくできていました。基本的なシミュレーションなら、これを通してジオメトリの取り込み、計算、可視化を容易に実施できます。

復習を兼ねて、手元にある3Dモデルに水を流してみました。ジオメトリは Civil3D で作成していましたので、まずはこれを取り込みます。

変換手順
1.Civil3D：ジオメトリ全体を iges 書き出し(部品毎の stl 書き出しは座標が飛ぶので不可)
2.FreeCAD：iges 読み込み → 部品毎に stl 書き出し（座標、スケールを保つ）
※正面図で書き出すこと。書き出したview の上方向がz軸となる
3.FreeCAD+DesignSPHysics：stl 読み込み（「Import GEO」 で1倍で取り込む）

 

読み込んだジオメトリに対し流体か boundary かを指定した後、諸々の計算条件を設定します。
「Execution Parameters」で人工粘性を組み込んだり、計算時間を決めたり。
「Object order」で粒子の発生順序（領域の重なる中で、どちらを後に発生させて上書きさせるか）を決めたり。計算を CPU で実施するか GPU で実施するかも、ここで指定できます。

 パラメータを決めたら、「Save Case」でモデルを保存。
「Run GenCase」の後、「Run」でDualSPHysics を実行します。

現行  Ver. では ｍDBC 未対応。書き出した xml を直接修正すればOKでした。DEM, Chrono も今後とのこと。これだけ簡単になると使い続けたくなります。楽しみですね。

計算後は「ComputeForce」で壁等にかかる力を抜き出せます。boundary の MK は＋10で指定すれば良いみたいですね。どこかに書かれているのかもしれませんが、わからなかったので VTK 書き出し後に MK を確認しました。
VTK書き出しは「PartVTK」、「IsoSurface」など。ParaView のパスを「Setup Plugin」で設定していたら、自動で立ち上げることも可能です。地味に便利な機能です。
VisualSPHysics を使う場合は、Blender でstl(boundary) + VTK(fluid)取り込み。Civil3Dではジオメトリを m単位 で作成していたため、stl を0.001倍指定で読み込むと VTK と座標が一致します。
残念ながら VisualSPHysics は開発が止まっているようで、 2.90LTS 未対応。2.83LTS に落としたら、問題なく動きました。

で、肝心の計算結果は、イマイチ。うーん、残念でした。

粒子法はメッシュを切らなくて済むので、手軽です。DesignSPHysics でより手軽になったかな。今後の開発が楽しみです。あらためて開発チームに感謝！

衝撃力のサンプリング

玉石がコンクリート壁に衝突する際の衝撃力を解析結果から抽出していました。

最初に衝突する6個の石の速度変化と一緒に図化すると、なぜか2個の衝撃力が取れていません。再計算しながら解析画像を書き出し、礫の衝突するタイミングや配置を追認してみましたが、きちんと壁に衝突しています。
1/1,000秒の書き出しではダメ？ということで1/10,000秒で書き出してみました。

はい、取り出せました。小さいのでΔtも小さいのかな。
1/10,000秒ではピーク値が取れていないように見える箇所があり、最終的には1/100,000秒で書き出して評価。この辺り、数値実験では容易です。

もし、実際に衝撃力を計るとすれば、このような高サンプリング周波数対応の応力測定器を探すだけで一苦労でしょう。ないかもしれませんね。

気をつけましょう。

************************************
20211019追記

https://www.jstage.jst.go.jp/article/sabo1973/32/1/32_1_40/_article/-char/ja/
衝撃力の時間に関し、「砂防ダムに対する土石流衝撃力算定とその問題点」で紹介がありました。内部の歪みは50μsec～1msec以内で消えるとのこと。今回は玉石なので、さらに短いのでしょう。1/1000秒間隔だと、拾えない衝撃力があるのも納得です。

Strontium isotopic composition

Strontium isotopic composition of hot spring and mineral spring waters, Japan
https://doi.org/10.1016/0883-2927(91)90053-R

温泉・鉱泉のストロンチウム同位体比が地質により異なるといった内容。
第四紀および新第三紀の火山地域の温泉水の87Sr/86Sr比は0.703-0.708の範囲を示し、花崗岩、堆積岩、変成岩地域のもの（0.706-0.712）よりも低いそうです。

SPH セミナー

昨年、パルマ大の SPH のセミナーを聞いていたのも、ちょうど今頃でした。

今年も開催されています。2週間の予定で、今日で前半戦が終わりました。
相変わらず英語が分かりません。説明は数式ベースなので何とかなるのですが、質疑応答がさっぱり。動画を見直せば何か引っかかるかもしれませんが、英語をクリアしても理解できない内容のような気がします。
参加されている皆さんの理解の速さに圧倒されています。否応なしに、立ち位置を確認させられています。

今年は dx と h 、それに伴う誤差の理解に焦点があてられているようです。これはありがたいですね。今まで、dx は計算速度から決定し、h＝1.3dx を当然のように使っていたのですが、エラーに関して気にしたことがありませんでした。求め方すら知りませんでした。以後、気を付けましょう。

理解促進のため、毎日コーディングの宿題が出ています。提出義務はないのですが、比較的容易な内容なので毎日実施しています。手は動かすもので、理解したつもりだった点、あいまいだった点が見つかります。
support domain = khと、カーネルに使用する h＝1.2dx の係数をあわせがちになっていたのですが、今回クリアになりました。初歩のガウスカーネルを選択して書いていたのですが、h が分散のところに入るので、カーネルの大きさに直に効くことにあらためて気づきました。セミナーでは、この大きさに誤差が関連することを説明されていました。
近傍探索アルゴリズムは、 以前、Fortran で組まれたコードを読んでいただけでした。今回、何も見ずに組んでみると迷いますね。あの部分、速く処理するためにどうしてたかな？などと。説明では、side effect としてメモリアクセス、キャッシュでの優位性を話されていました。が、さすがにこれを理解されていた受講生は少なかったような気がします。

翌朝の宿題に関する議論では、言語を問わないとあって FORTRAN, Julia, Pythonで書かれた方々が説明されていました。Julia が 出てきたのは面白い。グラフ化が必要ですからね。

あと4日。有意義な時間が続きますように。

df.plot()

Pandas 偉い。

index を datetime しておけば、

date_time	rain(mm/h)	Q(mm/h)
2009-10-07 01:00:00	0.019231	0.018576
2009-10-07 02:00:00	0.000000	0.018576
2009-10-07 03:00:00	0.019231	0.018576
2009-10-07 04:00:00	0.000000	0.018576
2009-10-07 05:00:00	0.038462	0.018576
2009-10-07 06:00:00	0.019231	0.018576
2009-10-07 07:00:00	0.019231	0.018576
2009-10-07 08:00:00	0.076923	0.018576
2009-10-07 09:00:00	0.057692	0.018576
2009-10-07 10:00:00	0.096154	0.021672
2009-10-07 11:00:00	0.076923	0.027864
2009-10-07 12:00:00	0.019231	0.034056
2009-10-07 13:00:00	0.038462	0.040248
2009-10-07 14:00:00	0.038462	0.043344
2009-10-07 15:00:00	0.038462	0.043344
2009-10-07 16:00:00	0.134615	0.043344

この呪文で

years=[2010,2011,2012,2013,2014,2015,2016,2017]
for year in years:
	df[df.index.year==year].plot()

こうなる。

 

指定した年毎に一気に書いてくれます。 

plot() の中に matplotlib の制御キーワードを渡せます（matplotlib の wrapper）。

Pandas 偉い。 

タンクモデルのパラメータ最適化

先輩＆後輩が、SCE-UA でタンクモデルのフィッティングを行っていました。

数年前にSCE-UAを利用した際、その恣意的な合わせこみ方法に引っかかっており、もっと一般的な最適化手法を適用すればどうなるの？と考えていました（手は動かしませんでしたが）。

良い機会ですので、何かコードが転がっていないか？と探し始めたところ、すぐに引っかかりました。

Python 流出解析用タンクモデルのパラメータ検索に挑戦（日流量解析）
https://qiita.com/damyarou/items/63480fb3e49ff2496a76

「パ－ソナル・コンピュ－タのためのタンク・モデル・プログラムとその使い方（第2報）」に沿って作られたようですね。オリジナルに乱数による範囲設定は入っていたかな？

一般的といえるかどうかはわかりませんが、タンクモデル本家の最適化アルゴリズムなので、ひとまず動くかどうか試してみることにしました。

動かないところを直しながら、少し手を加えて完成。動作確認もOK。結果が表示されました。Python だと計算はやや遅いのですが、図化まで一気にできるのが利点です。図化までのtotal 時間を考えると、小さいプログラムならこちらの方が速いので、手放せません。

で、結果はイマイチ。使われていたデータを眺めた段階で「ダメだろう」と思われたので、致し方ない結果。

ま、動きましたので、第2報の内容を読み返してみましょう。

大量処理

マージすると5億行になるデータを  Python で加工しいました。

for loop の速度が遅いことは有名ですが、まずは素直に書いて試行。
が、処理が終わるのに140日ほどかかる見込み。これはダメ。

joblib で parallel にするとともに、アルゴリズムを改良。これで、4時間程度に納まる見込み。
でしたが、メモリ不足を吐きます。SSD 500GB を仮想メモリに使うも足りず、HDD ２TB をたせば、、、フリーズ。頑張れ、Windows！

結局、小分けして処理。
が、これも夜中にフリーズしていました。Windows 嫌い。

最終的には９時間かかりました。計算は速いのですが、書き込みがネックになりますね。

他にも2週間計算していたPCが途中で再起動していたり、1週間かけて作っていた数億個のデータも途中でストップしていたり。

大きな処理は、PC ではなくクラウドが良いのでしょうか？データ転送がネックになるので、悩むところです。プロはどうされているのでしょう？ 

ArcGIS + ゼンリン

 ESRI さんが ArcGIS online でゼンリンさんの住宅地図を扱えるようにしたと案内が来ました。

これは嬉しい、と価格を問い合わせました。
結果、現ゼンリンさんとの契約額と1桁以上の差がありました。これは手が出ません。
ゼンリンさん提供の web 経由サービスでなく、数十倍の金額を払っても ArcGIS 経由で利用したい、という業種やニーズが思い浮かびません。123 で背景に使えると便利、という以上のメリットとなると、どのような利用法なのでしょう？

ArcMap が廃止予定となって実質サブスク化したり、サポートのレスポンスが1~2日程度かかるようになったり。ESRIサンには良い印象を持てなくなりつつあった時期に、このサービス提供。まだ優勢的地位は崩れず勢いがあるのでしょう。

numpy.reshape()

numpy で3次元配列を２次元配列にした場合の動作確認。

arr=np.array([[[1,2],[3,4],[5,6]],
              [[10,20],[30,40],[50,60]],
              [[100,200],[300,400],[500,600]],
              [[1000,2000],[3000,4000],[5000,6000]]])
              
print(arr.shape)
(4, 3, 2)

print(arr.reshape([4,-1]))
[[   1    2    3    4    5    6]
 [  10   20   30   40   50   60]
 [ 100  200  300  400  500  600]
 [1000 2000 3000 4000 5000 6000]]

転置すると、i,j,k が k,j,i になります。簡単。

 

GemPy

今日は Geoenvironmental Sciences に関する webinar を見ました。

1,2回目はライブで見ていたのですが、十分に理解できるわけもなく、3回目は Youtube にしました。
個人的には久しぶりとなる、地質に関する話でした。休日の穏やかな時間で見ていたこともあり、リラックスして聞けました。仕事では地質にかかわるヒトたちが CIM によるモデリングに迫られつつも見ぬふりをしていますが、ぜひ動画（で紹介されているようなツール）を見てほしいですね。

紹介されていたツールは GemPy。
https://www.gempy.org/
補間の計算には Theano が利用されています。PyMC3と組み合わせることで、確率論的手法（モンテカルロシミュレーション、ベイズ推論など）をサポートしていることが特徴とのことでした（そういえばPyMC3は動かなかったですね）。https://phreeqc.blogspot.com/2020/11/mcmc.html
textベースの入力ですので EVS と同様です。地質の３次元可視化という点では GEORAMA の方が操作性は良いと思われます。が、PyMC3が動けば強力な統合ツールになりそうです。地質って、不確実性もセットで表示すべきですからね。

webinar で、変化の速さを感じたのがココ↓。
https://youtu.be/dfKPxjSge90?t=2430
出典はわかりませんが、概ね理解できますし、利用してきたライブラリ等です。偶然、個人で選択していたツールが他者からも評価・整理されていると、安心します。なんとか、変化にはついていけていたと。一方で、ギリギリだったと複雑な思いも。

今日は webinar を見ただけですが、文献も出ています。
https://gmd.copernicus.org/articles/12/1/2019/
一度、読んでみましょう。

 

σレベル

σレベルが出てきました。

(x-平均)/標準偏差

あれ？標準化ですよね。
調べてみると、標準得点の z score と呼ばれるそうで、偏差値の仲間でした。統計と機械学習で呼び方というか、使われ方が違うのかな？

確率年よりも平易で処理が楽そう、というのが第一印象。
Today's Earth のリスク値としても利用されているようです。

有名なのでしょうね。避けずに学んでおけばよかった。

スペクトルの平滑化

連休中に、スペクトルの平滑化プログラムを書いておこうと思いつきました。

大崎本にサブルーチンが載っていますので、それを利用すれば完成するでしょうと高を括っていました。

が、本を持ち帰っていませんでした。

どうしましょうか、と検索してみると、大崎総合研究所から6月にソースコードが公開されていました。（正誤表もあります。）
http://www.ohsaki.co.jp/activity/download/index.html

昔の有名な図書ですので、プロは既にアレンジしたコードをお持ちです。が、私のようなアマにはありがたい。感謝です。

サブルーチンをダウンロードして、プログラムに組み込み、動かしてみました。が、正しく動きません。どこだ？と探しているうちに、連休が終わりました。

そういえば以前、倍精度に変更したサブルーチンをプロからいただいていました（今回も、連休中に頂きました）。それに組み替えると、改善。さらに、引数を変数に変えると正しく動きました。後者はともかく、前者はあるあるでしたね。感謝。

昨日、お礼がてら雑談していると、周波数によってバンド幅を変える方法もあると教えていただきました。バンド幅を決めて平滑化すると、高周波側と低周波側で見え方が変わってきます。それが改善されるのでしょう。良いですね。

といっても、それを使うのに妥当なスペクトルか判断できない知識レベルですので、今回は通常の平滑化で完成。バグ取りに数日かかりましたが、喰わせるデータや単位のチェックもあわせてできたので良しとしましょう。

天然ダムの側岸侵食

https://doi.org/10.1016/j.enggeo.2021.106322
Large-scale field model tests of landslide dam breaching

天然ダムの越流決壊を実験で確認した文献でした。

そういえば、決壊した後の形状はよく見ますが、越流初期からのプログレッシブな形状を追って見たことがありません。再現計算では流量と最終形状で評価しますし、そもそも映像がないので確認できないというのが実情でしょう。

この文献を読んでいて気づいたのが、側壁の侵食形状。水位から上で側岸崩落が発生しています。途中はそうなのかもしれませんね。

側方侵食速度、崩落角度も合わせこみパラメータの一つです。このような実験が増えることで、また一つ、この分野の経験式の積み上げが増えていくのでしょう。

トリチウムの利用

地下水の年代測定において、トリチウムの有用性が低下してから、ずいぶん経過しました。
https://phreeqc.blogspot.com/2017/04/blog-post_10.html

その間、他の同位体を利用していたため、動向のチェックを怠っていました。
久しぶりに実務で利用された方のグラフを見ると、壊変前の値でプロットされておらず、壊変後の値でプロットされえていました。ま、そういう方法もあるよね、と思いつつスルーしていました。が、先日、文献を読んでいて同様のグラフを見かけました。

https://doi.org/10.1016/j.apgeochem.2021.105070
Insight into watershed hydrodynamics using silica, sulfate, and tritium: Source aquifers and water age in a mountain river

出典はUSGS
https://doi.org/10.3133/sir20195090
Tritium as an indicator of modern, mixed, and premodern groundwater age
Scientific Investigations Report 2019-5090

年代を決定できなくても、mixed か premodern かがわかるだけでも価値があります。今後はこのような使い方が主流になるのでしょう。

Rayleigh damping

レイリー減衰でα、βを決定する考え方は、道路橋示方書に書かれています。

考慮すべき固有モードに対して安全側になるよう周波数、減衰比（減衰定数）を定めるような流れです。単純な構造、均一の部材であれば動解まで必要ないのですが、多様な形状・材料で橋が構成されているため、計算が必要になると理解しています。

これ、指定する周波数によっては、α、βの値が大きく変化します。結果にどの程度影響するのでしょう？SoilPlus で単純なモデルを作って計算してみました。

 

同じ減衰比0.01でも、狙う周波数帯によって、時刻歴が大きく異なります。今回のモデルでは質量比例減衰より、剛性比例減衰が効いています。 

 

 

 

周波数領域では、剛性比例のβが大きくなると、滑らかな形状になっています。値が大きくなると高周波側、低周波側の両方で効いています。
質量比例のαが加わると、低周波領域で振幅がやや低下しています。あまりに小さいと（注目する周波数帯が低いと）ほとんど効かないようですが。

 

思った以上に時刻歴波形に影響が出ていますし、思った以上に周波数領域で影響が出ていません。周波数領域で減衰力を変えて減算しているわけではないですからね。

 

もう少しケースを増やしてみないと感度解析にはなりませんが、参考にはなります。覚えておきましょう。

 

ArcGIS のライセンス統合

ESRI Japan さんの旧サポートサイトが12月でクローズする模様。

新しいサポートサイト MyESRI へ乗り換えないといけないのですが、紆余曲折を経て、ESRIの営業さんにお願いすることになりました。

その過程で、ライセンスの統合を進めました。驚くことに、2本目以降の年間保守料が半額近くまで下がりました。本体は微々たるものの、エクステンションの値下げ幅が大きい。もっと早く気付けばよかった。

世間のレベル

Python、なにそれ？とは異なる方との協働の話。

自己紹介もそこそこに協働することになった社外の方。メインの仕事とは外れる部分で、Python スクリプトをやり取りしたり、Fortran ソース、バッチファイル をいただいたり。
事前に確認したわけではありませんが、ある程度の言語スキル保有が前提でのお仕事でした。

先日のように社内の方をお世話するよりも、効率的かつストレスフリーでした。考えさせられます。

 

Python exe 変換 その2

数百の xml から 数項目の情報を抜き出して表にしたい、そんなプログラムを作ってくれるプログラマーがいないか？と他部署の先輩が探されていました。

それなら、先日の CIM スクリプトをダウングレードするだけだなあと思い、手を上げました。

30分ほどで出来上がったのですが、問題はPython（を使えない先輩とその部下）。

なので、重いし遅いのですが、実行ファイルに変換することにしました。
https://phreeqc.blogspot.com/2017/12/python-exe.html

そのまま変換すると 300MB の exe が出来上がりました。
これはダメ、軽くならないかと調べてみると、ありました。

【悲報】PyInstallerさん、300MBのexeファイルを吐き出すようになる
https://lets-hack.tech/programming/languages/python/pyinstaller-issue/

• コマンドを出した環境のPathにあるモジュールは根こそぎバンドルしに行くという仕様が原因。
• Pandas を入れる際はpip。 conda だとMKL も入ってしまう。

なるほど、全部入りでしたか。それなら仮想環境を毎回作るしかなさそうです（面倒ですが）。

手順
Jupyter の file - save as で .py 保存。

仮想環境確認
conda info -e

仮想環境作成
conda create -n pyinstaller
conda activate pyinstaller

インストール
conda install pyinstaller

今回利用するライブラリのインストール
pip install ●●

コンパイル
pyinstaller 〇〇.py --onefile --clean

変換された exe は 30MB。確かに軽くなりましたが、それでも重い。しかも遅い。

やはり、皆さんに覚えてもらうしかないでしょう。いえ、Python でなくてもよいのですが、何かしら一つくらい言語を習得しておかないと困ると思うのですが。手を貸すからダメなのかな。うーん。

個別要素法のパラメータ

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

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

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

先日、シュミットハンマーの反発度と反発係数を関連付けようとしていた文献を見かけました。この続きを読みたかったのですが、探せず。この先、どうなっているのか期待されます。
藤村ほか（１９９０） 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 を利用した現地踏査が始まっています。
「携帯端末を現地調査に使わない（使えない）」は、とっくに卒業すべき時期を迎えていたようです。