趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

勉強会の議事メモ:2024/03/02

勉強会お疲れ様でした。


一般座標系平面二次元河床変動計算を書く

computational-sediment-hyd.hatenablog.jp

上記記事のパターン1をベースに検討を進める。

⇒ SIMPLE法を基本とした離散化を検討 ⇒ プリミティブ変数はv,hか?q,hでは対応できない。

小さなメッシュサイズおける浅水流計算の問題

computational-sediment-hyd.hatenablog.jp

⇒ close

スタッガードスキームは離散化方法の修正に対応可(別途書きます。)

AMR法(Adaptive Mesh Refinement、適合格子細分化法)

引き続き、文献収集を行う。

  • PARAMESHが参考になる?

備忘録:pythonでwebスクレイピング:プログラミングの筋トレ


導入

たまにやらないと忘れるので。簡単なスクレイピングの話です。

海上保安庁の潮汐推算のwebページスクレイピングです。

データ公開ページ(例えば、https://www1.kaiho.mlit.go.jp/TIDE/pred2/cgi-bin/TidePredCgi.cgi?area=1603&back=../TidePred/tide_pred/2.htm&btn=ForAreaWindow_Japanese )ですが、 なぜかボタンをクリックするとバグります。

pythonによるwebスクレイピング

pythonでのwebスクレイピングのコツは、

  • seleniumに頼らない。出来るだけbeautifulsoupを使う。
  • Chromeデベロッパーツールでwebサイトをしっかり分析する。

くらいじゃないでしょうか。

それだけ、上のhtml内にタグで数値が埋め込まれているような化石サイトでも、下記のような短いコードでスクレイピングできます。

import pandas as pd
import numpy as np
from bs4 import BeautifulSoup
import requests
def getHarr(Area,Area2,year,month,day):
    url = f'https://www1.kaiho.mlit.go.jp/TIDE/pred2/cgi-bin/TidePredCgi.cgi?area={Area}&back=../TidePred/tide_pred/{Area2}.htm&year={year}&month={month}&day={day}&btn=%E6%8E%A8%E3%80%80%E7%AE%97'
    res = requests.get(url)
    soup = BeautifulSoup(res.content, 'html.parser')
    text = [t.get_text() for t in soup.findAll("td")]
    data = text[14:26] + text[40:52]
    # マイナス値の処理
    data = [d.replace(' ','') for d in data]
    return np.array(data, dtype=float)
%%time
Area = 1603
Area2 = 5
sdate = '2023/01/01'
edate = '2023/12/31'

val = [getHarr(Area, Area2, d.year, d.month, d.day) for d in pd.date_range(sdate, edate, freq='D')]
val = np.array(val).flatten()

df = pd.DataFrame(val
             , columns=[Area]
             , index=pd.date_range(sdate, periods=len(val), freq='H')
                 )
CPU times: total: 2min 16s
Wall time: 3min 49s
  • export
df.to_csv('tmp.csv', index_label='time')
  • export
df.plot()

Gist

gistd8b965fd743e0a46716787187dccfc82

Google Colab

※Colabの解説記事はこちら

Open In Colab

一次元浅水流計算における合流による水位上昇量の計算方法

本記事はGitHub、nbviewer、Colabでも公開しています。

※Colabの解説記事はこちら

a nbviewer Open In Colab



導入

「河道計画検討の手引き」に示される合流による水位上昇の計算方法について解説します。 この手法は、室田 明, 多田 博登:一次元水面形解析における合流点モデルに関する研究が元になっています。

計算式

合流による水位上昇量の計算方法の概要を以下に示す。

手引と同様に合流後河道の流向軸の運動量保存則は以下のとおりとなる。 コントロールボリュームの考え方が若干複雑なため、詳細は元論文を参考にされたい(時間ができたときにまとめます)。

\begin{align}
& \beta \rho Q_2 \frac{Q_2}{B_2h_2}\cos \theta _2 + \beta \rho Q_1 \frac{Q_1}{B_1h_1}\cos \theta _1-\beta \rho Q_3 \frac{Q_3}{B_3h_3}\cos \theta _3 \\
& = B_3\frac{1}{2}\rho g {h_3}^2 - B'_2 \cos \theta _2 \frac{1}{2}\rho g {h_2}^2 - B'_1 \cos \theta _1 \frac{1}{2}\rho g {h_1}^2 
\end{align}

ここに、\rho:水の密度、g:重力加速度、Q:流量、B:水路幅、h:水深、\theta:合流角度(下図参照)、\beta:運動量補正係数、B':修正した川幅(下図参照)、添字1,2,3は河道1,2,3(下図参照、ここでは1:本川、2:支川とする)の諸量である。

元論文のとおりh_1=h_2を条件として式変形を行うと以下のとおりとなる。

\begin{align}
&\left(\frac{B'_2}{B_3}\cos \theta _2 + \frac{B'_1}{B_3}\cos \theta _1\right)\left(\frac{h_1}{h_3}\right)^3-(1+2\beta {Fr}^2)\frac{h_1}{h_3}\\
&\hspace{20mm}+2 \beta \left\{\left(\frac{Q_2}{Q_3}\right)^2 \frac{B_3}{B_2}\cos \theta _2+\left(\frac{Q_1}{Q_3}\right)^2 \frac{B_3}{B_1}\cos \theta _1\right\}{Fr}^2=0\\
\end{align}

ここに、フルード数{Fr}^2=\dfrac{{Q_3}^2}{g{B_3}^2{h_3}^3}とする。

\begin{align}
X&=\dfrac{h_1}{h_3}=\dfrac{h_2}{h_3}&\\
\alpha &= \left(\frac{B'_2}{B_3}\cos \theta _2 + \frac{B'_1}{B_3}\cos \theta _1\right)&\\
\gamma &= \left\{\left(\frac{Q_2}{Q_3}\right)^2 \frac{B_3}{B_2}\cos \theta _2+\left(\frac{Q_1}{Q_3}\right)^2 \frac{B_3}{B_1}\cos \theta _1\right\}
\end{align}

として整理すると、

\begin{align}
X^3-\frac{1+2\beta {Fr}^2}{\alpha}X+ \frac{2 \beta \gamma {Fr}^2}{\alpha}=0
\end{align}

のとおり、Xの三次方程式となる。これを解くことによって水位上昇量を計算する。

なお、元論文では\alpha=1(角度補正した合流前川幅の和が合流後川幅と等しい)を前提条件としているが、前出の手引には記載されていない。

本記事では簡易的に\alpha=1\beta=1として水位上昇量Xの計算を行った。

テスト計算

計算モデル

  • 三次方程式はnp.rootsで解く。
  • 上記の解のうち、実数かつ1以上(水位上昇量のため)のものを水位上昇量とする。
  • 正しい条件を与えると解は1または0個になる。

参考 三次方程式の計算方法:https://west-village-tech.com/cubic_equation/

import numpy as np
def caldhConfluence(Q1Q3 ,Q2Q3 ,B3B1 ,B3B2 ,Theta1 ,Theta2 ,Fr):
    
    g = (Q2Q3)**2*B3B2*np.cos(Theta2*np.pi/180) \
     +  (Q1Q3)**2*B3B1*np.cos(Theta1*np.pi/180)
    
    x = np.roots([1,0,-1-2*Fr**2,2*g*Fr**2])
    
    # 実数のみ抽出
    x2 = x[np.iscomplex(x)==False]
    # 水位上昇量なので1以上のみ
    x3 = x2[np.where(x2>=1.0)]
    
    if len(x3) == 0:
        ans = np.nan
    elif len(x3) == 1:
        ans = x3[0].real
    else:
        print('error')
        ans = np.nan
    return ans

合流角度および合流前流量比の感度分析

合流角度は以下の3ケースを設定する。

\begin{align}
\begin{array}{l:c:c}
\mathrm{Case} & \theta_1[\mathrm{deg}]    & \theta_2[\mathrm{deg}] \\ \hline
\mathrm{Case A}  &  10  &  10 \\ 
\mathrm{Case B}  &  45  &  45 \\ 
\mathrm{Case C}  &  0   &  90 \\ 
\end{array}
\end{align}

合流前流量比は以下の3ケースを設定する。また、川幅はRegime theory(B = 5 \sqrt{Q})により設定する。よって、式中の係数は以下のとおりとなる。

\begin{align}
\begin{array}{l:c:c:c:c:c}
\mathrm{Case} & Q_1:Q_2 & Q_1/Q_3 & Q_2/Q_3 & B_3/B_1 & B_3/B_2 \\ \hline
\mathrm{Case1} & 1:1 & 1/2 & 1/2 & \sqrt{2} & \sqrt{2}  \\
\mathrm{Case2} & 2:1 & 2/3 & 1/3 & \sqrt{3/2} & \sqrt{3}  \\
\mathrm{Case3} & 5:1 & 5/6 & 1/6 & \sqrt{6/5} & \sqrt{6}  \\
\end{array}
\end{align}

これらを組み合わせ9ケースについて合流後のフルード数を0.1~1.0まで変化させた場合の水位上昇量の計算結果を以下に示す。

これらより、以下の点が確認できる。

  • フルード数が大きい程水位上昇量が大きい
  • 合流角度が大きい程水位上昇量が大きい。(CaseA,B)
  • 合流流量比は1:1に近い程水位上昇量が大きい。
  • 直角合流CaseCの場合、合流流量が大きいCase1では水位上昇量は大きいが、Case3で相対的に水位上昇量が小さくなっている。
Fr = np.arange(0.1,1.01,0.1)

Q1Q3=1/2
Q2Q3=1/2
B3B1=np.sqrt(2)
B3B2=np.sqrt(2) 
dhoutA1 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=10 ,Theta2=10 ,Fr=f) for f in Fr]
dhoutA2 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=45 ,Theta2=45 ,Fr=f) for f in Fr]
dhoutA3 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=0 ,Theta2=90 ,Fr=f) for f in Fr]
Q1Q3=2/3
Q2Q3=1/3
B3B1=np.sqrt(3/2)
B3B2=np.sqrt(3) 
dhoutB1 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=10 ,Theta2=10 ,Fr=f) for f in Fr]
dhoutB2 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=45 ,Theta2=45 ,Fr=f) for f in Fr]
dhoutB3 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=0 ,Theta2=90 ,Fr=f) for f in Fr]
Q1Q3=5/6
Q2Q3=1/6
B3B1=np.sqrt(6/5)
B3B2=np.sqrt(6) 
dhoutC1 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=10 ,Theta2=10 ,Fr=f) for f in Fr]
dhoutC2 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=45 ,Theta2=45 ,Fr=f) for f in Fr]
dhoutC3 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=0 ,Theta2=90 ,Fr=f) for f in Fr]
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context='paper', palette="bright", style="whitegrid", font_scale=0.9)
%matplotlib inline
fig = plt.figure(figsize=(5, 4), dpi=150)
ax = fig.add_subplot() 
ax.plot(Fr,dhoutA1, label='$ \\theta_1 =10,\\theta_2 =10, Q_1:Q_2=1:1$'         , c='b', ls='-')
ax.plot(Fr,dhoutA2, label='$ \\theta_1 =45,\\theta_2 =45, Q_1:Q_2=1:1$'         , c='g', ls='-')
ax.plot(Fr,dhoutA3, label='$ \\theta_1 = 0,\,\,\,\, \\theta_2 =90, Q_1:Q_2=1:1$', c='r', ls='-')

ax.plot(Fr,dhoutB1, label='$ \\theta_1 =10,\\theta_2 =10, Q_1:Q_2=2:1$'         , c='b', ls='--')
ax.plot(Fr,dhoutB2, label='$ \\theta_1 =45,\\theta_2 =45, Q_1:Q_2=2:1$'         , c='g', ls='--')
ax.plot(Fr,dhoutB3, label='$ \\theta_1 = 0,\,\,\,\, \\theta_2 =90, Q_1:Q_2=2:1$', c='r', ls='--')

ax.plot(Fr,dhoutC1, label='$ \\theta_1 =10,\\theta_2 =10, Q_1:Q_2=5:1$'         , c='b', ls=':')
ax.plot(Fr,dhoutC2, label='$ \\theta_1 =45,\\theta_2 =45, Q_1:Q_2=5:1$'         , c='g', ls=':')
ax.plot(Fr,dhoutC3, label='$ \\theta_1 = 0,\,\,\,\, \\theta_2 =90, Q_1:Q_2=5:1$', c='r', ls=':')

ax.set_xlabel('Froude number')
ax.set_ylabel('$X=h_1/h_3=h_2/h_3$')
# ax.set_title('[tex:Q\_1:Q\_2=1:1]')
ax.legend()
plt.show()

合流前川幅の感度分析

次に川幅の影響について感度分析を行った。

ケースとして、合流前川幅がRegime theory(B = 5 \sqrt{Q})の1.2倍、0.8倍を設定した。

式中の係数は以下のとおりとなる。ここでは、Q_1:Q_2=1:1のみとした。

\begin{align}
\begin{array}{l:c:c:c:c:c}
\mathrm{Case} & B_1,B_2 &  Q_1:Q_2 & Q_1/Q_3 & Q_2/Q_3 & B_3/B_1 & B_3/B_2 \\ \hline
\mathrm{Case I}                                 & \mathrm{by\,Regime\,theory}\,\times\,1.0 & 1:1 & 1/2 & 1/2 & \sqrt{2} & \sqrt{2}  \\
\mathrm{Case I\hspace{-1.2pt}I}                 & \mathrm{by\,Regime\,theory}\,\times\,1.2 & 1:1 & 1/2 & 1/2 & \dfrac{5\sqrt{2}}{6} & \dfrac{5\sqrt{2}}{6}   \\
\mathrm{Case I\hspace{-1.2pt}I\hspace{-1.2pt}I} & \mathrm{by\,Regime\,theory}\,\times\,0.8 & 1:1 & 1/2 & 1/2 & \dfrac{5\sqrt{2}}{4} & \dfrac{5\sqrt{2}}{4}   \\
\end{array}
\end{align}

これらについて合流後のフルード数を0.1~1.0まで変化させた場合の水位上昇量の計算結果を以下に示す。

これらより、以下の点が確認できる。

  • 合流前川幅が大きいほど水位上昇量が大きくなる。これは流速が小さくなるためである。
  • 合流前川幅を小さくする場合、本計算条件では約0.7倍以下でXが1以下、つまり水位が上昇しない条件になる(本モデルではnanを返す)。これは、合流前の流れが射流になるためである。
Fr = np.arange(0.1,1.01,0.1)

Q1Q3=0.5
Q2Q3=0.5
B3B1=np.sqrt(2)
B3B2=np.sqrt(2)
dhout1 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=10 ,Theta2=10 ,Fr=f) for f in Fr]

B3B1=5/6*np.sqrt(2)
B3B2=5/6*np.sqrt(2)
dhout2 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=10 ,Theta2=10 ,Fr=f) for f in Fr]

B3B1=5/4*np.sqrt(2)
B3B2=5/4*np.sqrt(2)
dhout3 = [caldhConfluence(Q1Q3,Q2Q3,B3B1,B3B2,Theta1=10 ,Theta2=10 ,Fr=f) for f in Fr]
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context='paper', palette="bright", style="whitegrid", font_scale=0.9)
%matplotlib inline
fig = plt.figure(figsize=(5, 4), dpi=150)
ax = fig.add_subplot() 
ax.plot(Fr,dhout1, label='$\mathrm{by\,Regime\,theory}\, \\times \,1.0$')
ax.plot(Fr,dhout2, label='$\mathrm{by\,Regime\,theory}\, \\times \,1.2$')
ax.plot(Fr,dhout3, label='$\mathrm{by\,Regime\,theory}\, \\times \,0.8$')

ax.set_xlabel('Froude number')
ax.set_ylabel('$X=h_1/h_3=h_2/h_3$')
ax.set_title('[tex:Q\_1:Q\_2=1:1]')
ax.legend()
plt.show()

GitHub

github.com

参考文献

関連記事

computational-sediment-hyd.hatenablog.jp

備忘録:google Colabでcondaを使うためのライブラリ:condacolabの使い方

本記事はColabでも公開しています。

※Colabの解説記事はこちら

Open In Colab



https://pypi-camo.freetls.fastly.net/6cae3dec9a1e98ece1b184bd4863b26a0f73218d/68747470733a2f2f6769746875622e636f6d2f6a61696d657267702f636f6e6461636f6c61622f7261772f6d61696e2f636f6e6461636f6c61622e706e67

condacolabとは?

https://pypi.org/project/condacolab/(condacolab · PyPI]

Colabでcondaを使うためのライブラリです。

pythonではcondaによるインストールが推奨されるライブラリが多くあります。 Colabではpipが標準でcondaは準備されておりません。 Colabでconda環境を整備できなくはないですが、手間がかかる&安定しません。 そのため、現時点(2024/2/10)ではcondacolabの使用が最適と思われます。

condacolabの使い方

公式https://pypi.org/project/condacolab/のとおりですが簡単にまとめておきます。

インストール

pipでインストールします。

!pip install -q condacolab
import condacolab

メソッド:condacolab.install

condacolab.installでconda環境を構築します。

以下のメソッドが準備されています。

  • 現時点(2024/2/10)ではpython3.10の環境が準備されています。pythonのバージョンは指定できないようです。
  • たまにColabがクラッシュすることがありますが再接続で解決します。

condacolab.install_miniconda()

Minicondaディストリビューションをインストールします。

condacolab.install_miniconda()

condacolab.install_miniforge()

Miniforgeディストリビューションをインストールします。 Miniforgeディストリビューションはconda-forgeによって公式に提供されています。

condacolab.install_miniforge()

condacolab.install_mambaforge()

推奨。Miniforgeのようなものですが、mambaが含まれています。Mambaforgeディストリビューションはconda-forgeによって公式に提供されています。

condacolab.install_mambaforge()

condacolab.install_anaconda()

Anacondaのフルパッケージをインストール?(未確認)

condacolab.install_anaconda()

condacolab.install_from_url()

スクリプトのURLを指定してインストール?(未確認)

condacolab.install()

condacolab.install_mambaforge()と同じ。

メソッド:condacolab.check

condaが正しくインストールされていることを確認する。

condacolab.check()

ライブラリのインストール

ディストリビューションのインストールが完了したら、通常のcondaと同様にライブラリをインストールします。

!conda install [package]

Gist

gista7996070dcafd41c23e421be4639ea5a

備忘録:python shapelyの互換性の問題 2024/02/05時点

pythonのライブラリshapelyは、Version 2.0.1 (2023-01-30)、Version 2.0.2 (2023-10-12)にリリースされていますが、下記のとおり、cartopyでのlgeosのインポートエラーにより一部機能に不具合がでます。

参照:Version 2.x — Shapely 0 documentation

cartopyでのlgeosインポートエラー(2023/06/04)

2023年5月のアップデートに伴い、MacPortsでインストールされるshapelyのバージョンが1.8.5から2.0.1になりました。このバージョンでは、cartopyy-0.21.0が内部でインポートしているlgeosが含まれないようになったため、以下のようなImportErrorが発生します(Python3.10の場合のエラーメッセージ例)。このエラーは、cartopyインポート時に不可避な部分で発生しており、shapely-2.0.1のアップデートが適用されてしまった場合には、shapelyのバージョンを1.8.5に戻す必要があります。

引用元:気象データ解析のためのmatplotlibの使い方:更新・エラー情報

現時点では1.8.5が使うしかないです。2.0で色々な機能が追加されたので、早く修正されると助かります。

conda install shapely==1.8.5 --channel conda-forge

河川技術者向け基礎講座 準二次元不等流計算3/4:一般断面の不等流計算2-分割断面法(平均流速公式レベル2)

教科書の「不等流計算」は理解できるが、実務で汎用される「準二次元不等流計算」がなかなか理解し難い方に向けて、記事を書きました。 多分相当わかりやすいと思いますので、河川行政に関わる方やコンサルの方に読んで頂きたいです。

全4回の第3回目です。


本記事はGitHub、nbviewer、Colab、Marp(スライド形式)でも公開しています。

※Colabの解説記事はこちら

※Marpの解説記事はこちら

a nbviewer Open In Colab a



はじめに

  • 前回は一般断面で用いて断面内の流速が一様である条件での計算方法を示した。
  • 今回は下図の複断面形状(低水路+高水敷)にように、断面内の領域ごとに明らかな流速差が生じる場合の計算方法として分割断面法(平均流速公式レベル2)について説明を行なう。

A


  • 分割断面とは、粗度が異なる領域ではなく、流速が異なる領域で区分する。分割断面内で粗度が一定である必要はない。
  • 分割断面間で流れの干渉は生じないこと条件とする。
  • 分割断面法において、
    • 分割断面ごとに変わる:流速、河積、潤辺、粗度係数など
    • 分割断面ごとに変わらない:水位、エネルギー勾配

基礎式

平均流速公式:マニング則

前回示した流速が一様の領域の平均流速公式は以下である。

\begin{align}
 v &= \dfrac{1}{n'}i_e^{1/2}R^{2/3} \\
 Q &= K i_e^{1/2} \\
 K &\equiv \dfrac{A^{5/3}}{n'S^{2/3}} \\
 n' &= \left( \dfrac{  \displaystyle \sum_{i=1}^{imax-1} S_i n_i^{3/2}}{ \displaystyle \sum_{i=1}^{imax-1} S_i } \right)^{2/3} 
\end{align}

ここに、Q:流量、H:水位、A:河積、v:流速、R=A/S:径深、S:潤辺、K:通水能、n':マニングの粗度係数(合成粗度)とする。


分割断面法では、各分割領域で上式が成り立つことを条件に分割断面ごとの流量Q_nの総和が全体の流量Qとして次式を導出する。

\begin{align}
 Q &= \sum^{n} Q_n \\
  &=  \sum^{n} K_n i_e^{1/2} \\
 K &\equiv \dfrac{A_n^{5/3}}{n_n'S_n^{2/3}} \\
\end{align}

ここに、添字nは各分割断面の諸元を示す。

運動方程式の生成項が摩擦損失のみの場合、マニング則を用いると以下のとおりになる。

\begin{align}
 \dfrac{\tau}{\rho g A} = \dfrac{Q^2}{\left(\displaystyle \sum^{n} K_n \right)^2}
\end{align}

運動方程式

分割断面法における不等流計算の運動方程式は次式を用いる。

\begin{align}
& \dfrac{d}{dx}\left( \frac{\beta Q^2}{2gA^2} + H \right) = -\dfrac{\tau}{\rho g A} 
\end{align}

ここに、Q:流量、H:水位、A=\sum A_n:河積、\tau:コントロールボリュームに作用する力、\beta:運動量補正係数とする。

前回示した一般断面の運動方程式とほぼ同型である。 新たに追加された\beta:運動量補正係数について詳述しておく。


運動量補正係数

運動量補正係数は運動方程式の導出過程で、断面積分を行う際に必要となる係数である。

断面平均流速Vは局所流速u積分で表すことができる。 分割断面法では次式となる。

\begin{align}
& V = \dfrac{ \displaystyle \sum^n {u_n}A_n}{A}
\end{align}

一次元の運動方程式の導出過程で、\dfrac{ \displaystyle \sum^n {u_n^2}A_n}{A}が現れるが、分割断面下では次式は成立しない。

\begin{align}
\dfrac{ \displaystyle \sum^n {u_n^2}A_n}{A} \neq V^2 
\end{align}

そこで次のように示す。

\begin{align}
\dfrac{ \displaystyle \sum^n {u_n^2}A_n}{A} &= \beta V^2  \\
 \beta &\equiv \dfrac{ \displaystyle \sum^n {u_n^2}A_n}{V^2 A}
\end{align}

このとき、\betaを運動量補正係数と定義する。

分割断面法の運動方程式では、運動量補正係数を考慮する必要がある。


エネルギー補正係数

同様に、エネルギー保存則においても、エネルギー補正係数\alphaが必要となる。

\begin{align}
 \dfrac{d}{dx}\left( \dfrac{\alpha q^2}{2gh^2} + h + z_b \right) &= -i_e 
\end{align}

一次元のエネルギー保存則の導出過程では、\dfrac{ \displaystyle \sum^n {u_n^3}A_n}{A}が現れるが、分割断面下では次式は成立しない。

\begin{align}
\dfrac{ \displaystyle \sum^n {u_n^3}A_n}{A} \neq V^3 
\end{align}

そこで次のように示す。

\begin{align}
\dfrac{ \displaystyle \sum^n {u_n^3}A_n}{A} &= \alpha V^3  \\
 \alpha &\equiv \dfrac{ \displaystyle \sum^n {u_n^3}A_n}{V^3 A}
\end{align}

このとき、\alphaをエネルギー補正係数と定義する。

分割断面法のエネルギー保存則では、エネルギー補正係数を考慮する必要がある。


ここで、水面形方程式を整理すると次式となる。

\begin{align}
 \dfrac{d}{dx}\left( \dfrac{\alpha q^2}{2gh^2} + h + z_b \right) &= -i_e \\
- \dfrac{\alpha q^2}{gh^3}\dfrac{dh}{dx} + \dfrac{dh}{dx} + \dfrac{d z_b}{dx} &= -i_e \\
\dfrac{dh}{dx} &= \dfrac{i_b -i_e }{ 1 - \dfrac{\alpha q^2}{gh^3}} 
\end{align}

水面形方程式の分母に\alphaが含まれるため、前々回の説明したとおり、限界水深を計算する式が変更される(詳細は後述)。


補足:井田法による合成径深

参照:井田:広巾員開水路の定常流-断面形の影響について (わかりづらいので読まなくても良いです)

ここで、分割断面法において水深の代替として用いることが多い井田法による合成径深について説明しておく。

分割断面法による平均流速は次式で示される。

\begin{align}
 Q &= \left( \sum^{n} \dfrac{1}{n'_n}R_n^{2/3} A_n \right) i_e^{1/2}
\end{align}

断面全体においてもマニング則が成立するとして次式を定義する。

\begin{align}
 Q &= \dfrac{1}{N}R^{2/3} i_e^{1/2} A
\end{align}

ここに、N:断面全体のマニングの粗度係数、R:断面全体の径深を示す。


2式を連立すると次式が導出される。

\begin{align}
 R &= \left( \dfrac{N}{A} \sum^{n} \dfrac{1}{n'_n}R_n^{2/3} A_n  \right)^{3/2}
\end{align}

本式からわかるとおり、径深Rと粗度係数Nは独立していない。

井田は、水理計算上取り扱いやすいように、径深Rと粗度係数Nは独立させるために、 次式が成立するNN_cと定義した。

\begin{align}
  N_c \sum^{n} \dfrac{1}{n'_n}R_n^{2/3} A_n  \equiv  \sum^{n} R_n^{2/3} A_n
\end{align}

よって、

\begin{align}
  N_c \equiv   \dfrac{ \displaystyle \sum^{n} R_n^{2/3} A_n }{\displaystyle  \sum^{n} \dfrac{1}{n'_n}R_n^{2/3} A_n}
\end{align}

これを元のRの式に代入するとR_c:井田法による合成径深が得られる。

\begin{align}
 R_c &= \left( \dfrac { \displaystyle \sum^{n} R_n^{2/3} A_n }{A} \right)^{3/2}
\end{align}

R_cは水位のみで決まるため、分割断面法において水深の代替として用いられることが多い。

N_cは井田による等価粗度係数、合成粗度係数などと呼ばれる。前回定義した合成粗度係数と混同しないように注意すること。


限界水深、等流水深の定義

一般断面の基礎式を用いて、限界水深、等流水深を矩形断面と同様に定義で設定することは難しいため、便宜的に以下のとおりに設定する。


等流水深

全ての損失による水頭の勾配が河床勾配i_bと釣り合う状態を等流と定義してその水深を等流水深とする。 なお、一般断面では水深を用いないため、等流時の水位(以下、等流水位と定義)となる。

\begin{align}
 \dfrac{\tau}{\rho g A} = i_b
\end{align}

生成項が摩擦損失のみの場合は、

\begin{align}
 \dfrac{Q^2}{\left(\displaystyle \sum^{n} K_n \right)^2} = i_b
\end{align}

となり、これを満足する水位Hを反復法などにより求めれば良い。


限界水深

フルード数が1となる水位を限界流時の水位(以下、限界水位と定義)とする。

フルード数はエネルギー保存則の分母が0より、

\begin{align}
    1 &= \dfrac{\alpha q^2}{gh^3} \\
    Fr&=  \dfrac{V}{ \sqrt{\dfrac{gh}{\alpha}}} \\
    Fr&=  \dfrac{Q}{A\sqrt{\dfrac{gh}{\alpha}} }
\end{align}

となるが、平方根の中に水深hが含まれるため、一般断面ではそのままでは計算できない。

そのため、分割断面法では水深hの代替として井田の合成径深R_cを使用することが多い。その他には径深R=A/SA/Bを用いることもある。

参考:FORUM8ソフトウェア:等流・不等流の計算・3DCAD Ver.9 Q&A


離散化

離散化は次式となる。 なお、i:上流側、i-1:下流側とする。

\begin{align}
  \left(\frac{\beta_i Q^2}{2gA^2_i} + H_i  \right) 
-\left( \frac{\beta_{i-1} Q^2}{2gA_{i-1}^2} + H_{i-1}  \right) 
&= \dfrac{1}{2}\left(\dfrac{Q^2}{K_i^2} + \dfrac{Q^2}{K_{i-1}^{2}}\right)\Delta x \\
A&=\sum A_n \\
K&=\sum K_n
\end{align}

常流の場合、下流から逐次計算を行なうため未知数はH_iのみとなる。


数値計算方法

一般断面の場合、離散式の未知数H_iによる微分が難しいため、ニュートン法が使いづらい。そのため、二分法を使用することをおすすめする(もちろんニュートン法を使っても良い)。ただし、通常の二分法では安定的に計算することが難しいため、少し工夫が必要である。


分割断面法の応用

分割断面法では流速の異なる領域ごとに分割断面を設けるため、 高水敷と低水路という区分に限らない。 例えば、横断面内に樹林帯が存在し、著しく流速が低下する場合などはその領域に分割断面を設定しても良い。

import matplotlib.pyplot as plt
x=[0, 5,93,100,200,206,294,300]
y=[6, 3.5, 3.5,  0,  0,  3,  3,  6]
plt.plot(x,y, c='k')
plt.vlines([0,93,206,300], ymin=0, ymax=6, color='r')
plt.axvspan(xmin=30, xmax=60)
plt.show()                # 描画

png


分割断面法の必要性

  • なぜ分割断面法が必要性なのか?合成粗度による方法では駄目なのか?について考える。
  • 基礎式をふまえて正しく積分するという意味は当然ある。
  • 複断面河道の流れを適切に一次元で取り扱うこと自体に相当無理がある。
  • それでも、分割断面法を使う理由は河床抵抗を適切な評価できることである。
  • 次図に水位と通水能の関係を示すが、合成粗度では高水敷付近で不連続となるが、分割断面法では解消されている。これにより、安定的に水理計算が可能となる。

  • 合成粗度 A

  • 分割断面法 A


少し考えてみよう。

下に多自然川づくりの例でよく見る横浜市のいたち川の写真を示す。写真の下(施工後)の場合、どのように水理計算を行えば良いと思いますか?

合成粗度?分割断面法?

A 出典:(財)リバーフロント整備センター「多自然川づくり参考事例集」https://tenbou.nies.go.jp/science/description/detail.php?id=95


横断面内の流速分布の実態

  • 実験水路の計測結果を参照すると、横断面内の流速分布は水路幅と水深の比(アスペクト比)によって次図のようになる。
  • 流速分布は横断方法、鉛直方向の流れ(二次流:主流の数%程度の流れ)で大きさ決まる。
  • 現地スケールで、特に出水時の観測事例はほとんどない。
  • アスペクト比が小さい断面形状では、Velocity dipと呼ばれる最大流速位置が水面より少し下に沈むような現象が生じる。

冨永ら1985

A


冨永ら1990

A


禰津ら1993

A


まとめ

流速の異なる複数の領域を持つ河道横断形状を対象とした一次元(準二次元)不等流計算方法について説明を行った。

  • 一次元(準二次元)開水路流れの水理計算を行なう上ではここまで理解できれば十分。
  • 逆に、不等流計算に限らず、一次元不定流計算、一次元河床変動計算などで水位を正確に計算したい場合は、ここまでは考慮して欲しい。
  • 県管理河川の河道計画ではこれ以上の知識は不要。
  • 次回以降説明する平均流速レベル3は、直轄管理河川の河道計画に特化した方法で他で使うことほぼ無い。

断面特性を計算するプログラム例

ボリュームがあるため、別ページに作成

a nbviewer Open In Colab

不等流計算のプログラム例

Coming soon.

参考文献

関連記事

computational-sediment-hyd.hatenablog.jp

勉強会の議事メモ:2024/01/08

勉強会お疲れ様でした。


一般座標系平面二次元河床変動計算を書く

computational-sediment-hyd.hatenablog.jp

上記記事のパターン1をベースに検討を進める。

小さなメッシュサイズおける浅水流計算の問題

computational-sediment-hyd.hatenablog.jp

状況を共有。圧力項と重力項の分離等の対応方法を検討。

四分木法

AMR法(Adaptive Mesh Refinement、適合格子細分化法)も含めて、関連資料の収集を行う。 特に、2wayのやりとりについて離散化、コーディングも含めて収集する。