趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

河床粒子の抗力係数と揚力係数の比について:0.85って

新しいモデルを作る中で、河床粒子の抗力係数と揚力係数の比がふと気になったので調べてみました。

河床変動計算で揚力係数CL/抗力係数CDを使う場合、0.85を使うことが多いです。

これは、 池田先生の移動床河川の蛮曲部における二次流と動的横断平衡河床について から決まっており、W. S. ChepilさんのThe use of evenly spaced hemispheres to evaluate aerodynamic forces on a soil surfaceの実験結果を整理して導いたものです。

感覚的には結構大きなと感じます。というのは、しばしば、揚力の影響は無視する・考慮しないとありますが、こんなに大きいなら無視したら駄目じゃんと思うからです。

海外の文献もググってみましたがほとんど情報がないですね。

気が向いたらもう少し調べてみます。


イスコの髭がない

ArcGISの傾斜角、傾斜方向の計算方法

全然大した話じゃないですが、ArcGISの傾斜角、傾斜方向の計算方法について確認してみました。以下のページに書いてました。

傾斜角 (Slope) の仕組み—ヘルプ | ArcGIS for Desktop

傾斜方向の仕組み—ヘルプ | ArcGIS for Desktop

傾斜角、傾斜方向は普通に以下の式で計算していますが、問題はzの空間微分です。

 slope = \arctan{ \left( \sqrt{ \left( \frac{\partial z}{\partial x} \right)^2 + \left( \frac{\partial z}{\partial y} \right)^2}\right) }


aspect = \arctan{ \left( \dfrac{ \dfrac{\partial z}{\partial y} }{ -\dfrac{\partial z}{\partial x}} \right) }

結論としては、計算方法は周辺9点の座標を使って中央差分でこんな感じで計算してます。


\frac{\partial z}{\partial x} = \frac{1}{4}\left(\frac{z_{i+1,j+1}-z_{i-1,j+1}}{2\Delta x}
+ 2\frac{z_{i+1,j}-z_{i-1,j}}{2\Delta x} +\frac{z_{i+1,j-1}-z_{i-1,j-1}}{2\Delta x} \right)


で、ここからが本題というか面白い話ですが元のwebページにはこんな感じで書いてます。

aa

[dz/dx] = ((c + 2f + i) - (a + 2d + g) / (8 * x_cellsize)     
[dz/dy] = ((g + 2h + i) - (a + 2b + c)) / (8 * y_cellsize)

暗号みたいで気持ち悪い・・・。この書き方はないな。

砂防学会誌:新砂防の過去の有名な講座のリンク

流砂関連は書籍もあまり無いためなかなか勉強しづらいのですが、砂防学会誌内にまとめられた代表的な講座のリンクを備忘録的に貼っておきます。

伝説的な講座なので、これを理解できればかなりのレベルです。

講座「河床変動の数値計算法」(1997-1998)

これらは以下の書籍にまとめられておりますが絶版なので元論文を読むと良いかと思います。

山地河川における河床変動の数値計算法

山地河川における河床変動の数値計算法

  • 作者:
  • 出版社/メーカー: 山海堂
  • 発売日: 2000/09
  • メディア: 単行本

  1. 河床変動計算のススメ 水山高久(京都大学)
  2. 河床変動現象と基礎方程式 宮本邦明(鳥取大学)
  3. 急流河川の河床変動(その1) 岡部健士(徳島大学)
  4. 急流河川の河床変動(その2) 岡部健士
  5. 2次元河床変動の数値計算法 (その1) 檜谷 治(鳥取大学)
  6. 2次元河床変動の数値計算法 (その2) 檜谷 治(鳥取大学)
  7. 砂防ダムの堆砂計算法(その1) 藤田正治(京都大学)
  8. 砂防ダムの堆砂計算法(その2) 藤田正治(京都大学)
  9. 土石流の堆積 高浜淳一郎(砂防・地すべり技術センター)
  10. 数値計算の現状と今後の展開 宮本邦明

講座「土石流」(1992-1993)

土石流関連としては世界トップクラスでしょう。もしかすると書籍にまとまっているのかもしれませんがよくわかりません。

  1. 土石流概論 芦田和男
  2. 土石流の観測 諏訪 浩
  3. 土石流発生のメカニズム(1)太田猛彦
  4. 土石流発生のメカニズム(2)高橋 保
  5. 土石流流動のメカニズム(1) 宮本邦明
  6. 土石流流動のメカニズム(2) 宮本邦明
  7. 土石流流動のメカニズム(3) 宮本邦明
  8. 土石流の停止・堆積のメカニズム 江頭進治
  9. 土石流の停止・堆積のメカニズム 江頭進治
  10. 土石流の実験と相似則 水山高久

CLのドルトムント vs PSG凄かった。。。久しぶりにしびれました。あのインテンシティは、完全に照準を合わせてきた感じですね。

砂州の安定解析についてのメモ

黒木・岸の研究(https://www.jstage.jst.go.jp/article/jscej1969/1984/342/1984_342_87/_article/-char/ja/)による砂州の安定解析を分かりやすく説明して欲しいという無茶ぶりされたので無理ですと答えたのですがそうもいかずに。。。

導入部分だけ自分のメモ用に。

式はこんな感じですね。

 \tilde{\eta} = \hat{\eta} \cos(l\tilde{\eta}) e^{ik(\tilde{x} - c\tilde{t})}

 c = C_r + i C_i

 k = \dfrac{2 \pi H_0}{2L}

 l = \dfrac{m \pi H_0}{B}

mが砂州の列数です。t=0の想定する河床の形状(擾乱)を図化するとこんな感じです。

f:id:SedimentHydraulics:20200211132817p:plain

ここで、cは無次元複素移動速度で論文中で水深、河床高、流速xy方向、河床せん断力xy方向、流砂量xy方向の8変数と関数と定義されています。

あとはゴリゴリとcを解いてkcの正負で安定、不安定を判定する。

さらに論文内ではm=0,1,2で最も卓越するcの条件によって分類するようなことをやってます。 ちなみに河川用語では、m=0は砂州無し、m=1は交互砂州、m=>2は多列砂州と呼ばれます。

図化のソースは以下です。

Jupyter Notebook Viewer


3Dはholoviewsで上手くかけないので、久しぶりにmatplotlib使いました。

図化調整用のアフィン行列は以下を参照しました。

naga-tsuzuki.sblo.jp

論文紹介:河川流の三次元解析

https://www.jsnds.org/ssk/ssk_28_4_325.pdf

流体屋さんが河川の解析をやるとこうなるという感じが面白いです。勉強になります。

杉山先生の論文はよく読むんですが河川流を機械流体と同じ視点で見ているところが好きです。

面白いのはξ、η、ζ座標系でζ軸を水深方向に全然合わせてないところです。 河川屋さんはσ座標系とかソロバン格子とか使って意地でも水深方向に合わせるのにそれをしていない。 格子生成は面倒だと思いますが。

多列砂州河道の再現計算(失敗例)⇒ ちょっとだけ分析

先日の記事

computational-sediment-hyd.hatenablog.jp

をちょっとだけ分析として、局所勾配を計算してみました。

計算方法は、 傾斜角 (Slope) の仕組み—ヘルプ | ArcGIS for Desktop を使いました。

局所勾配はこんな感じです。

f:id:SedimentHydraulics:20200205194937p:plain

内部摩擦角約40度を超える部分が結構ありますね。

流砂量モデルにその影響を組み込まないと駄目ですね。

geoviewsの背景をgeotiffにする方法

これは結構質問がくるので書いておきます。

geoviewsでは背景画像にWMTSを使うことがよくありますが、 (例えば、

computational-sediment-hyd.hatenablog.jp

) 紙ベースにするとき(←止めればいいのに)にpng exportするとWMTSが表示されません。

そんなときの応急的対策です。

geotiffを作る。

背景画像のgeotiffを作ってください。tfw付きはtiffは駄目です。

tiff + tfw や jpeg + jgwからgeotiffの作り方は以下を参照ください。

qiita.com

geoviewsに読み込む

ただRGBで取り込むだけです。

import xarray as xr
import cartopy.crs as crs
import geoviews as gv

gv.extension('bokeh', logo=False)

#load geotiff
dstif = xr.open_rasterio('photo.tif')

r, g, b = dstif[0,:,:].values, dstif[1,:,:].values, dstif[2,:,:].values
dstif2 = xr.Dataset({'R': (['y','x'], r), 'G': (['y','x'], g), 'B': (['y','x'], b) }
                  , coords={'x': dstif['x'], 'y': dstif['y']}
               )
out = gv.RGB(dstif2, ['x', 'y'],['R','G','B'], crs=crs.epsg(xxxx)).opts(projection=crs.epsg(xxxx))

これでpng exportが可能です。 この上にコンター図とかベクトル図とか重ねれば良いと思います。

以下を参考にしました。

https://gitter.im/pyviz/pyviz?at=5d28deb53b41924785aa4d66

著作権に注意

WMTSはダウンロードして使うことを想定しておりませんので、使用の際は必ず配信元の仕様を確認下さい。

私がよく使う国土地理院のものはクレジットを記載して使用してください。