iRICの格子生成アルゴリズム"Simple Grid Generator"について

:::note
本テキストは公開されている資料を基に@starfieldKtが作成しています。
そのため開発者様が意図していることと異なる部分がある可能性があることをご注意ください。
また、本テキストに関することを開発者様へお問い合わせ頂かないようお願い申し上げます。
:::

その他格子生成ツールについては[こちら](https://qiita.com/starfieldKt/private/fcc0d0261bdfc169696c)

# 簡易直線・蛇行水路生成ツール *"Simple Straight and Meandering Channel Creator"*

|項目|内容                                                                                                  |
|:- |:-                                                                                                    |
|名称|*Simple Grid Generator*                        |
|開発者|Dr.Jonathan Mark Nelson (USGS)                                                      |
|説明|無し|



## 概要
iRICの発足人のひとりである、USGS(United States Geological Survey)のDr.Jonathan Mark Nelsonが作成した格子生成プログラム。
直線水路、蛇行水路の生成が可能。
~~日本語の辞書ファイルが同梱されていないのでパラメータ名等が英語になっている。~~
~~日本語表示にしたい場合は[iRICのソルバー開発者マニュアル](https://iric-solver-dev-manual-jp.readthedocs.io/ja/latest/02/04_dictionary.html)に辞書ファイルの作り方が載っているので自分で翻訳して辞書ファイルを作成すればよい。~~

~~翻訳精度等諸々を保証しませんが、私が自分用に作成した辞書ファイルは[こちら](https://drive.google.com/file/d/18ZqyWSMP7oTjNQukMdeYY7EUlSCSKGtW/view?usp=sharing)。~~

:::note
2022.9.6追記
私が作成した辞書ファイルをiRIC ver3のSimple Grid Generator(ver.1.0.2)に同梱していただけることになりました。
2022.9.6以降に本ソルバーをインストールまたは更新した場合日本語環境ではGUIが翻訳された状態になっていると思います。

※本記事は辞書ファイル同梱前に作成したものなので画像内のGUIが英語のままです。
※ソルバーの開発者と辞書ファイル作成者(@starfieldKt)は別人です、翻訳に関するお問い合わせは開発者へしないようにお願い致します。
:::


## 格子生成条件
格子生成条件の入力画面は以下のようになっている。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/f010e928-77de-28dd-3ced-f6618c07a74d.png)


### 水路の平面形

`水路延長(m)：Channel Length (m)`、`水路幅(m)：Channel Width (m)`で水路の延長、幅を設定できるほか、`蛇行角(度)：Channel Crossing Angle (degrees)`で水路の蛇行角の最大値を設定できる。
プログラムの計算に内容を見ていないため断言はできないが、$s$を水路中央線に沿った軸として、`水路延長(m)：Channel Length (m)`を$L$、`蛇行角(度)：Channel Crossing Angle (degrees)`を$\theta_0$とすると、s軸と水平方向のなす角度$\theta$は恐らく以下のコサイン関数で表されていると思われる。

```math

\theta=\theta_0cos\biggl(2\pi\frac{s}{L}\biggr)

```
なお、式が表しやすかったのでcos関数で示しているが、これは[格子生成アルゴリズム 簡易直線・蛇行水路生成ツールの記事](https://qiita.com/starfieldKt/items/62f440074cdbce59a023)で説明したサイン・ジェネレーテッド・カーブの位相をずらしたものと考えられる。

`蛇行波数：Number of Meander Wavelengths`は$0\leqq s \leqq L$の波形を何波分繰り返すかを設定できる。

横断方向の格子は設定した水路幅で水路中央線のs軸方向に直行するように生成される。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/f4b9bad0-578d-d324-6a59-d7a6a530b7f0.png)



### 格子の解像度
`計算格子：Computational Grid`の項目`流下方向格子点の数：Number of stream-wise nodes`と`横断方向格子点の数：Number of cross-stream nodes`では流下方向と横断方向の格子のサイズを設定できる。
ここで入力するのはそれぞれ上流端から下流端までの格子点の数と、右岸端から左岸端までの格子点の数であるので、水路を流下方向または横断方向にある分割数で分割したい場合は分割数＋1の数を入力してあげればよい。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/5893038f-6081-0357-bf8e-6a985bc130c8.png)

この時設定した格子点は上流端から流下方向に何個目かをi、右岸端から左岸方向に何個目かをjとしてそれぞれ格子点(i,j)として表すことにする。

### 水路中心線の河床高

各格子点の河床高は水路中心線の下流端の河床高と、`水路の勾配(m/m)：Channel Slope (m/m)`から水路中心線の初期河床高を求め、それを基準に各横断面で断面形状や蛇行洲の設定によって計算されると考えらえるが、パラメータ画面では水路中央部の下流端河床高を設定する項目がない。
試しに何パターンか作成してみたところ以下のような結果になった.
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/1dfed7f5-be70-a1ad-bbc4-40142ab1e38b.png)



結果から見ると、`平坦：Flat`を選択した時は下流端の水路中央の河床高は10mから`平坦河床の深さ(m)：Flat Bed Depth(m)`を引いたものであると推測できる。
`放物線状:Parabolic`についても`平坦：Flat`の時より約0.01m低い値になっているが、同様に10m(9.99m)から`放物線状河床の最大深さ(m)：Parabolic Maximum Depth(m)`を引いた値が下流端の水路中央部の河床高になると考えられる。



### 河床形状
`横断面形状：Cross-section Geometry`と`蛇行洲の形状：Point Bar Geometry`の項目では水路横断面の河床形状を設定できる。

まず初めに`蛇行洲高の振幅(0-1)：Point bar amplitude (0-1)`が0の場合を整理する。

`横断面形状の種類：Cross-section Type`で`平坦：Flat`を選択すると横断面の河床形状は平坦となり、水路中央部の標高(z軸)が一律で横断方向の格子点にも与えられる。

`放物線状：parabolic`を選択した場合は横断面の河床形状は水路中央部を頂点とした放物線の形状となり、`放物線状河床の最大深さ(m)：Parabolic Maimum Depth(m)`(おそらく`Parabolic Maximum Depth`なので以後そう示す。)を$h_{pmd}$、`水路幅(m)：Channel Width (m)`を$B$、格子点($i,j$)の標高を$z_{ij}$、$x,y$座標を$x_{ij},y_{ij}$、$i$断面の水路中央の標高を$z_{ic}$、$x,y$座標を$x_{ic},y_{ic}$とすると、作成される結果より以下の式で$z_{ij}$を表すことができると思われる。

```math
z_{ij}=z_{ic}+\frac{4h_{pmd}}{B^2}\sqrt{(x_{ij}-x_{ic})^2+(y_{ij}-y_{ic})^2}
```
これもソースコードを見ていないので正しいかはわかりませんが、iRICで作成した格子の断面$i$の格子点の座標を出力したものと、上記の式で算出した座標が以下の図のように一致していたので恐らく問題は無いと思います。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/10234f47-aff2-ff7f-0167-16337d3b7a2d.png)

次に、`放物線状の土手の高さ(m)：Parabolic Bank Height(m)`を設定した場合横断面の河床形状は、水路中央部を頂点とするのは同じですが、`放物線状河床の最大深さ(m)：Parabolic Maximum Depth(m)`を設定した時とz軸の方向に反転した放物線の形状となる。よって、`放物線状河床の最大深さ(m)：Parabolic Maximum Depth(m)`を$h_{pmd}$、`放物線状の土手の高さ(m)：Parabolic Bank Height(m)`を$h_{pbh}$とすると、$z_{ij}$は以下のようになる。

```math
z_{ij}=z_{ic}+\frac{4(h_{pd}-h_{ph})}{B^2}\sqrt{(x_{ij}-x_{ic})^2+(y_{ij}-y_{ic})^2}
```

イメージとしては以下のようなかたちになる。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/23a574e9-ee60-274d-b233-cbec4dbaefe6.png)


### 交互砂州

続いて`蛇行洲高の振幅(0-1)：Point bar amplitude (0-1)`を0以外に設定する場合についてですが、このパラメータを設定することで交互砂州のある水路形状を作成することができます。
このパラメータは0～1の倍率となっており、`平坦：Flat`を選択した場合は`平坦河床の深さ(m)：Flat Bed Depth(m)`の1.6倍に、`放物線状：parabolic`を選択した場合は`放物線状河床の最大深さ(m)：Parabolic Maximum Depth(m)`-`放物線状の土手の高さ(m)：Parabolic Bank Height(m)`にこの倍率をかけたものが、最大洗堀深と最大堆積厚の合計となる。

`蛇行洲高の振幅(0-1)：Point bar amplitude (0-1)`が1の時について交互砂州がどのように作成されるかを以下に示す。

#### 河床形状がflatの場合

`平坦：Flat`を選択し、交互砂州を設定した際は以下のように格子が作成される。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/c5f761f0-ed98-3c43-00f0-c04335c50635.png)


上の画像の平面図で最深洗堀深箇所と最大堆積圧を含む横断面と縦断面を示しているが、これらより交互砂州を設定した場合、横断方向は横断方向インデックスjを変数とした$0\leqq \theta \leqq \pi$の範囲でのcosカーブ、縦断方向については蛇行の1周期の長さ(`水路延長(m)：Channel Length (m)`)の半分が波長となり、上流端からちょうど蛇行の1周期の半分の地点で正と負が入れ替わるインデックスiを変数としたcosカーブで表せることが推測できる。

`蛇行波数：Number of Meander Wavelengths`が1の時の場合を考えると、格子点($i,j$)の標高$z_{ij}$は交互砂州を設定しない場合の各点の標高を$z_{ij0}$、`平坦河床の深さ(m)：Flat Bed Depth(m)`を$h_{fbd}$、`流下方向格子点の数：Number of stream-wise nodes`を$n_{i \\_ max}$、 `横断方向格子点の数：Number of cross-stream nodes`を$n_{j \\_ max}$として$ 1 \leqq i \leqq n_{i\\_ max} $、$ 1 \leqq j \leqq n_{j\\_max} $とすると以下のように推測できた。

$1 \leqq i \leqq \frac{n_{i \\_ max}-1}{2}+1$のとき
```math
z_{ij}=z_{ij0}+\frac{2h_{fbd}}{5}\Biggl(\cos\biggl(\frac{\pi(j-1) }{n_{j \_ max}-1}\biggr)\Biggr)\Biggl(\cos\biggl( \frac{2\pi(i-1)}{\frac{n_{i \_ max}-1}{2}}+\pi\biggr)+1\Biggr)
```
$\frac{n_{i \\_ max}-1}{2}+1 < i \leqq n_{i\\_ max}$のとき
```math
z_{ij}=z_{ij0}+\frac{2h_{fbd}}{5}\Biggl(\cos\biggl(\frac{\pi(j-1) }{n_{j \_ max}-1}\biggr)\Biggr)\Biggl(-\cos\biggl( \frac{2\pi(i-1)}{\frac{n_{i \_ max}-1}{2}}+\pi\biggr)+1\Biggr)
```

この式を基にExcelを用いて各格子点の標高を計算してみたところ、小数点以下の扱いの関係で100万分の1程度の誤差はあるが、格子生成ツールで作成し出力した結果の数値と一致していた。


#### 河床形状がparabolicの場合

`放物線状：parabolic`を選択し、交互砂州を設定した際は以下のように格子が作成される。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/3dbdd6a6-bb8a-9484-b3d4-e1ce5ec92b0f.png)

このままでは交互砂州を設定したことでどのようになったかがわかりづらいので、交互砂州有り無しでの各格子点の標高の差分を以下に示す。
以下の格子は上記の格子の交互砂州無しバージョンを作成し、上記格子と共にcsvに出力後、各格子点の差分を求めたcsvファイルを作成しiRICにインポートすることで作成しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/15dcb10e-bc42-f0fd-8601-6c955b4da613.png)

この格子についても平面図で最深洗堀深箇所と最大堆積圧を含む横断面と縦断面を示しているが、これらより交互砂州を設定した場合、横断方向は横断方向インデックスjを変数とした$0\leqq \theta \leqq 2\pi$の範囲でのsinカーブ、縦断方向については'Flat'の時と同様に、蛇行の1周期の長さ(`水路延長(m)：Channel Length (m)`)の半分が波長となり、上流端からちょうど蛇行の1周期の半分の地点で正と負が入れ替わるインデックスiを変数としたcosカーブで表せることが推測できる。

`蛇行波数：Number of Meander Wavelengths`が1の時の場合を考えると、格子点($i,j$)の標高$z_{ij}$は交互砂州を設定しない場合の各点の標高を$z_{ij0}$、`放物線状の土手の高さ(m)：Parabolic Bank Height(m)`を$h_{pbh}$、`放物線状河床の最大深さ(m)：Parabolic Maximum Depth(m)`を$h_{pmd}$、`流下方向格子点の数：Number of stream-wise nodes`を$n_{i \\_ max}$、 `横断方向格子点の数：Number of cross-stream nodes`を$n_{j \\_ max}$として$ 1 \leqq i \leqq n_{i\\_ max} $、$ 1 \leqq j \leqq n_{j\\_max} $とすると以下のように推測できた。

$1 \leqq i \leqq \frac{n_{i \\_ max}-1}{2}+1$のとき
```math
z_{ij}=z_{ij0}+\frac{h_{pmd}-h_{pbh}}{4}\Biggl(\sin\biggl(\frac{2\pi(j-1) }{n_{j \_ max}-1}\biggr)\Biggr)\Biggl(\cos\biggl( \frac{2\pi(i-1)}{\frac{n_{i \_ max}-1}{2}}+\pi\biggr)+1\Biggr)
```
$\frac{n_{i \\_ max}-1}{2}+1 < i \leqq n_{i\\_ max}$のとき
```math
z_{ij}=z_{ij0}+\frac{h_{pmd}-h_{pbh}}{4}\Biggl(\sin\biggl(\frac{2\pi(j-1) }{n_{j \_ max}-1}\biggr)\Biggr)\Biggl(-\cos\biggl( \frac{2\pi(i-1)}{\frac{n_{i \_ max}-1}{2}}+\pi\biggr)+1\Biggr)
```

この式を基にExcelを用いて各格子点の標高を計算してみたところ、小数点以下の扱いの関係で100万分の1程度の誤差はあるが、格子生成ツールで作成し出力した結果の数値と一致していた。

#### 交互砂州の位相を変化させる

`蛇行洲と蛇行頂点の位相差(度)(流下方向が正)：Point bar - bend apex phase difference (degrees) (+downstream)`に角度(degree)を与えると交互砂州の縦断方向の位相を移動することができる。
位相の変化距離は流下方向を正として`蛇行洲と蛇行頂点の位相差(度)(流下方向が正)：Point bar - bend apex phase difference (degrees) (+downstream)`を$\theta_{apd}$、蛇行の1周期の長さ(`水路延長(m)：Channel Length (m)`)を$L$とすると$\frac{\theta_{apd}}{720}L$となる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/32b5a687-5b44-d91e-2e15-ea0819f86a7b.png)

## 簡易直線・蛇行水路生成ツール "Simple Straight and Meandering Channel Creator"との違い
格子生成のパラメータ等を整理してみると、この格子生成ツールは[簡易直線・蛇行水路生成ツール](https://qiita.com/starfieldKt/items/62f440074cdbce59a023)と同じような格子を生成することが可能であるとわかると思うが、若干の違いがあるため整理する。

### 蛇行の位相
簡易直線・蛇行水路生成ツールとSimple Grid Generatorでは、蛇行水路の平面形状について、蛇行が開始する位相が異なる。
x軸方向に対するs軸の進行方向との角度は以下の式でそれぞれ表すことができる。

簡易直線・蛇行水路生成ツール
```math

\theta=\theta_0 \sin \biggl(2\pi\frac{s}{L}\biggr)

```

Simple Grid Generatorの場合

```math

\theta=\theta_0 \cos \biggl(2\pi\frac{s}{L}\biggr)

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/02d9e452-b979-dfb9-d514-da11b10363ce.png)

### 選択できる蛇行の種類
Simple Grid Generatorでは直線水路およびサイン・ジェネレーテッド・カーブ蛇行水路を作成できるが、簡易直線・蛇行水路生成ツールではこれに加えて歪度を考慮した木下蛇行水路も作成できる。

### 水路幅の縦断的変化
縦断的に水路幅を変化させる機能はSimple Grid Generatorに搭載されていない。

### 選択できる横断河床形状
簡易直線・蛇行水路生成ツールでも横断河床形状は平坦、砂州のある場合、放物線状と3選択することができた。
Simple Grid Generatorでは河床形状を平坦と放物線状から選択し、別途砂州について設定するので4パターンから選択することができる。

放物線状かつ砂州を考慮した水路を作成したい場合はSimple Grid Generatorでしか作成できない。

## 格子作成例

交互砂州のある直線水路
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/974bdd7c-ed82-8cd7-2a71-38fc859fe077.png)

蛇行水路（平坦な河床）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/4b6bc495-5cab-157c-2fd5-07329c2fbb05.png)

蛇行水路（放物線状の河床）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/4b60468f-7b6f-3b8e-a90b-4392e98e9b6c.png)

蛇行水路（砂州のある放物線状の河床）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2626673/3fe5ca3f-eb70-bbe6-78b7-8f83faaa7a02.png)