趣味で計算流砂水理

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

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

MENU

岸・黒木による実験式を修正する

スポンサーリンク

実験式修正シリーズの第2段として、前回の岩垣の式(岩垣による限界掃流力の式の不連続部を解消する - 趣味で計算流砂水理)に続いて、河床波の遷移と有効掃流力と無次元有効掃流力の関係を示した岸・黒木による関係式(水理委員会移動床流れの抵抗と河床形状研究小委員会 : 移動床流れにおける河床形態と粗度)を修正します。



岸・黒木による実験式とは?

ご存知のとおり、河床にはdune(砂堆)のような河床波が生じており、出水時には掃流力の増加に伴い河床波の形状がdune~flat-bed(平坦河床)~anti-dune(反砂堆)のように遷移します。

その遷移に伴う流水抵抗の変化を評価したEngelundによる実験式(Hydraulic Resistance of Alluvial Streams)を基に修正を加えたものが、岸・黒木の実験式です。

元式の使用上の課題

実験式なので各係数の打ち切り誤差の影響等により不具合でることはよくあります。

岸・黒木の実験式は以下のように定義されます。

\begin{align}
\rm{dune1} : \tau_*^{\prime} &= 0.21\tau_*^{1/2} \\
\rm{dune2} : \tau_*^{\prime} &= 1.49(R/d)^{-1/4}\tau_* \\
\rm{flat\text{-}bed} : \tau_*^{\prime} &=  \tau_* \\
\rm{antidune} : \tau_*^{\prime} &= 0.264(R/d)^{1/5}\tau_*^{1/2} \\
\rm{transition1} : \tau_*^{\prime} &= 6.5 \times 10^7(R/d)^{-5/2} \tau_*^{11/2}
\end{align}

dune領域とtransition領域、flat-bed領域とantidune領域の区分は次式となる。

\begin{align}
\rm{dune \text{ -- } transition} : \tau_* &= 0.02(R/d)^{1/2} \\
\rm{flat\text{-}bed \text{ -- } antidune} : \tau_* &= 0.07(R/d)^{2/5} 
\end{align}

この式をそのまま図化すると下図のとおりで変曲点付近を拡大すると次のように不連続になっています。なお、R/dは100としています。

import numpy as np
import pandas as pd
import holoviews as hv

hv.extension('bokeh', logo=True) 
def kishikuroki1(taus, Rbyd):
    A1, A2 = 0.02, 0.07 
    
    if taus < A1*(Rbyd)**0.5:
        tausd = 0.21*taus**0.5
    elif taus < A2*(Rbyd)**0.4:
        tausd = 1.49*(Rbyd)**(-0.25)*taus
    else:
        tausd = 0.264*(Rbyd)**0.2*taus**0.5
        
    return tausd

def kishikuroki2(taus, Rbyd):
    A1, A2 = 0.02, 0.07 
    
    if taus < A1*(Rbyd)**0.5:
        tausd = 0.21*taus**0.5
        tausd = tausd if tausd < taus else taus
    elif taus < A2*(Rbyd)**0.4:
        tausd = 6.5*10**7*(Rbyd)**(-2.5)*taus**5.5 
        tausd = tausd if tausd < taus else taus
    else:
        tausd = 0.264*(Rbyd)**0.2*taus**0.5
        
    return tausd

Rbyd = 100

taus = np.logspace(-1.31, 0.7, 10000)

tausd1 = [ kishikuroki1(tt, Rbyd) for tt in taus ]
tausd2 = [ kishikuroki2(tt, Rbyd) for tt in taus ]

g1 = hv.Curve((tausd1, taus)).options(color='r', line_width=2, show_grid=True, line_dash='dashed')
g2 = hv.Curve((tausd2, taus)).options(color='r', line_width=2, show_grid=True)

H1 = hv.HLine(0.02*(Rbyd)**0.5).options(color='k', line_width=1)
H2 = hv.HLine(0.07*(Rbyd)**0.4).options(color='k', line_width=1)

mak = hv.Scatter([[0.094, 0.02*(Rbyd)**0.5], [0.441, 0.07*(Rbyd)**0.4]]).options(size=20, line_color='k', color=None)
gall0 = (H1*H2*g1*g2).options(logx=True, logy=True, width=300, height=250)
gall = gall0*mak.options(logx=True, logy=True, width=300, height=250)

ex1 = gall0.redim.range(x=(0.090,0.100),y=(0.195,0.205)).options(title='拡大図1')
ex2 = gall0.redim.range(x=(0.435,0.445),y=(0.435,0.445)).options(title='拡大図2')

gout = (gall + ex1 + ex2).cols(1).options(shared_axes=False)

out = hv.save(gout,'kishikurokiorg.html')
del out

gout

岸・黒木による実験式の修正

元論文をみると、dune領域とtransition領域、flat-bed領域とantidune領域の区分は、Garde-Rajuの研究成果を参照しているとのことなので、岸・黒木が提案した式のうち、dune2、transition1、anti-duneの式を修正します。

単純に接続する2式の交点を計算し直すだけで、各係数を再定義します。

上式を係数に置き換えて記載すると次のとおりとなる。

\begin{align}
\rm{dune1} : \tau_*^{\prime} &= C_1\tau_*^{1/2} \\
\rm{dune2} : \tau_*^{\prime} &= C_2(R/d)^{-1/4}\tau_* \\
\rm{flat\text{-}bed} : \tau_*^{\prime} &=  \tau_* \\
\rm{antidune} : \tau_*^{\prime} &= C_3(R/d)^{1/5}\tau_*^{1/2} \\
\rm{transition1} : \tau_*^{\prime} &= C_4 (R/d)^{-5/2}\tau_*^{11/2}
% C_1 = 0.21 ,C_2 &= 1.49 ,C_3 = 0.264 ,C_4 = 6.5\times 10^7
\end{align}

また、領域区分の式は次のとおりとなる。

\begin{align}
\rm{dune \text{ -- } transition} : \tau_* &= A_1(R/d)^{1/2} \\
\rm{flat\text{-}bed \text{ -- } antidune} : \tau_* &= A_2(R/d)^{2/5} \\
A_1 &= 0.02, A_2=0.07
\end{align}

係数C_2,C_3,C_4を計算し直すと次のとおりである。

\begin{align}
  C_1 &= 0.21 \\
  C_2 &= C_1/\sqrt{A_1} \\
  C_3 &= \sqrt{A_2} \\
  C_4 &= C_1/{A_1}^5 
\end{align}

係数を比較すると下表のとおりとなる。

original modify
C1 0.21 0.21
C2 1.49 1.4849242
C3 0.264 0.2645751
C4 65000000 65625000

上記の同様のグラフを作成すると下のようになり不連続部が解消されていることがわかる。

A1, A2 = 0.02, 0.07 
C1 = 0.21
C2 = C1/A1**0.5
C3 = A2**0.5
C4 = C1/A1**5
def kishikuroki1(taus, Rbyd):
    
    if taus < A1*(Rbyd)**0.5:
        tausd = C1*taus**0.5
    elif taus < A2*(Rbyd)**0.4:
        tausd = C2*(Rbyd)**(-0.25)*taus
    else:
        tausd = C3*(Rbyd)**0.2*taus**0.5
        
    return tausd

def kishikuroki2(taus, Rbyd):
    
    if taus < A1*(Rbyd)**0.5:
        tausd = C1*taus**0.5
        tausd = tausd if tausd < taus else taus
    elif taus < A2*(Rbyd)**0.4:
        tausd = C4*(Rbyd)**(-2.5)*taus**5.5 
        tausd = tausd if tausd < taus else taus
    else:
        tausd = C3*(Rbyd)**0.2*taus**0.5
        
    return tausd

Rbyd = 100

taus = np.logspace(-1.31, 0.7, 10000)

tausd1 = [ kishikuroki1(tt, Rbyd) for tt in taus ]
tausd2 = [ kishikuroki2(tt, Rbyd) for tt in taus ]

g1 = hv.Curve((tausd1, taus)).options(color='r', line_width=2, show_grid=True, line_dash='dashed')
g2 = hv.Curve((tausd2, taus)).options(color='r', line_width=2, show_grid=True)

H1 = hv.HLine(A1*(Rbyd)**0.5).options(color='k', line_width=1)
H2 = hv.HLine(A2*(Rbyd)**0.4).options(color='k', line_width=1)

mak = hv.Scatter([[0.094, A1*(Rbyd)**0.5], [0.441, A2*(Rbyd)**0.4]]).options(size=20, line_color='k', color=None)
gall0 = (H1*H2*g1*g2).options(logx=True, logy=True, width=300, height=250)
gall = gall0*mak.options(logx=True, logy=True, width=300, height=250)

ex1 = gall0.redim.range(x=(0.090,0.100),y=(0.195,0.205)).options(title='拡大図1')
ex2 = gall0.redim.range(x=(0.435,0.445),y=(0.435,0.445)).options(title='拡大図2')

gout = (gall + ex1 + ex2).cols(1).options(shared_axes=False)

out = hv.save(gout,'kishikurokiorg_mod.html')
del out

gout

これで安心して数値計算等に使えます。

Github

nbviewer.jupyter.org