岩垣による限界掃流力の式の不連続部を解消する - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

河床変動の数値計算では、よく使用する岩垣による限界掃流力の式について不具合があったので若干修正しました。なんかこういうのって土木っぽい。。。

モチベーション

岩垣の式は実験式のため、図化すると関数の接続部で不連続になります。普通に使う分には特に問題もないですが、収束計算等に使うとエラーが生じる場合があるので修正しておきます。

拡大図

岩垣の式の係数の修正

岩垣の式は次式のとおりである。

$\begin{align} R_* \ge 671.0 \: : \:\: & u_{*c}^2 = 0.05G_wgd\\ 162.7 \le R_* \le 671.0 \: : \:\: & u_{*c}^2 = (0.01505G_wg)^{\dfrac{25}{22}}\nu^{\dfrac{3}{11}}d^{\dfrac{31}{22}} \\ 54.2 \le R_* \le 162.7 \: : \:\: & u_{*c}^2 = 0.034G_wgd \\ 2.14 \le R_* \le 54.2 \: : \:\: & u_{*c}^2 = (0.1235G_wg)^{\dfrac{25}{32}}\nu^{\dfrac{7}{16}}d^{\dfrac{11}{32}}\\ R_* \le 2.14 \: : \:\: & u_{*c}^2 = 0.14G_wgd \end{align}$

ここに、 $G_w=\sigma/\rho-1$ 、 $R_*=\dfrac{G_w^{1/2}g^{1/2}d^{3/2}}{\nu}$ である。（ $R_*$ は無次元粒径かな？)

岩垣の式はshields図表の近似式となるが、shields図表のx軸:限界状態の粒子レイノルズ数? $\dfrac{u_{*c}d}{\nu}$ とy軸: 無次元限界掃流力 $\tau_{*c}=\frac{u_{*c}^2}{G_wgd}$ を使って整理すると次のとおりになる。

$\begin{align} R_* \ge 671.0 \: : \:\: & \tau_{*c} = 0.05\\ 162.7 \le R_* \le 671.0 \: : \:\: & \tau_{*c} =\alpha_1 \left(\dfrac{G_wgd^3}{\nu^2}\right)^{\beta_1} = \alpha_1 \left[ \left(\dfrac{u_{*c}d}{\nu}\right)^2 \tau_{*c}^{-1}\right]^{\beta_1} \\ 54.2 \le R_* \le 162.7 \: : \:\: & \tau_{*c} = 0.034 \\ 2.14 \le R_* \le 54.2 \: : \:\: & \tau_{*c} =\alpha_2 \left(\dfrac{G_wgd^3}{\nu^2}\right)^{\beta_2}=\alpha_2 \left[ \left(\dfrac{u_{*c}d}{\nu}\right)^2 \tau_{*c}^{-1}\right]^{\beta_2} \\ R_* \le 2.14 \: : \:\: & \tau_{*c} = 0.14 \end{align}$

ここに、 $\alpha_1=0.01505^{\dfrac{25}{22}},\beta_1=\dfrac{3}{22},\alpha_2=0.1235^{\dfrac{25}{32}},\beta_2=-\dfrac{7}{32}$ である。また、 $R_*=\dfrac{G_w^{1/2}g^{1/2}d^{3/2}}{\nu}=\left(\dfrac{u_{*c}d}{\nu}\right) \tau_{*c}^{-0.5}$ となる。

元論文に詳しく書いてないですが、基本的には実験値を直線で近似したもので、上式の1,3,5番目は粒子レイノルズ数によらず一定、2,4番目は両対数軸上で1次関数で近似したものとなっている。

今回の課題である不連続部は $\alpha,\beta$ の近似精度に依存すると仮定してこれらの見直しを行う。

参考として、2番目の式を展開する。

$x = \log_{10} \left(\dfrac{u_{*c}d}{\nu}\right)$ 、 $y = \log_{10} \tau_{*c}$ として式変形すると以下となる。

$\begin{align} y = \dfrac{2\beta_1}{1+\beta_1} x + \dfrac{1}{1+\beta_1}\log_{10} \alpha_1 \end{align}$

$y_1=0.034,y_2=0.05, x_1=162.7\sqrt{y_1}, x_2=671.0\sqrt{y_2}$ として、両対数軸上で直線近似により $\alpha_1,\beta_1$ を計算する。

$\begin{align} \beta_1 &= \dfrac{A}{2-A} \\ \alpha_1 &= 10^{(1+\beta_1)B} \end{align}$

ここに、 $A$ ：近似直線の傾き、 $B$ ：切片である。

以上より得られた係数を列挙すると下表のとおりである。

	元論文	修正値
$\alpha_1$	0.008492171377400728	0.008502572284837527
$\beta_1$	0.13636363636363635	0.13609748853308548
$\alpha_2$	0.19514831719540923	0.19535327890603318
$\beta_2$	-0.21875	-0.2189567763631981

冒頭と同様に図化すると以下のとおり不連続部が解消されている。大局的には何も変わらないが数値計算上はこちらを使うほうが良い。

拡大図

修正した岩垣の式のコード

pythonのコードは以下のとおり

オリジナル

def iwagaki_org(d, g=9.8, Gw=1.65, nu=0.000001):
    #d:diameter[m]
    t = Gw*g*d**3/nu**2
    Rs = np.sqrt(t)
    
    R1, R2, R3, R4 = 671.0, 162.7, 54.2, 2.14
    Y1, Y2, Y3 = 0.05, 0.034, 0.14
    alpha1, beta1 = 0.008492171377400728, 0.13636363636363635
    alpha2, beta2 = 0.19514831719540923, -0.21875
    
    if Rs >= R1:
        taus = Y1
    elif Rs >= R2:
        taus = alpha1*t**beta1
    elif Rs >= R3:
        taus = Y2
    elif Rs >= R4:
        taus = alpha2*t**beta2
    else:
        taus = Y3
        
    return taus, Rs*np.sqrt(taus)

修正版

def iwagaki_mod(d, g=9.8, Gw=1.65, nu=0.000001):
    
    def calcoef(R1, R2, Y1, Y2):
        X1, X2 = R1*np.sqrt(Y1), R2*np.sqrt(Y2)
        A = (np.log10(Y2) -  np.log10(Y1))/ (np.log10(X2) -  np.log10(X1))
        B = np.log10(Y1)- A*np.log10(X1)
        beta = A/(2-A)
        alpha = 10**((1+beta)*B)
        return alpha, beta
    
    #d:diameter[m]
    t = Gw*g*d**3/nu**2
    Rs = np.sqrt(t)
    
    R1, R2, R3, R4 = 671.0, 162.7, 54.2, 2.14
    Y1, Y2, Y3 = 0.05, 0.034, 0.14
    
    alpha1, beta1 = calcoef(R2,R1,Y2,Y1)
    alpha2, beta2 = calcoef(R4,R3,Y3,Y2)
    
    if Rs >= R1:
        taus = Y1
    elif Rs >= R2:
        taus = alpha1*t**beta1
    elif Rs >= R3:
        taus = Y2
    elif Rs >= R4:
        taus = alpha2*t**beta2
    else:
        taus = Y3
        
    return taus, Rs*np.sqrt(taus)

参考：図化のコード

import numpy as np
import holoviews as hv
hv.extension('bokeh')

dd = np.logspace(-5, -2, 10000)
arr = np.array([iwagaki_org(d) for d in dd])
arr2 = np.array([iwagaki_mod(d) for d in dd])

fig = hv.Curve((arr[:,1],arr[:,0]), label='org').options(logx=True, logy=True, width=500, show_grid=True).redim.range(y=(0.03,0.3))

fig2 = hv.Curve((arr2[:,1],arr2[:,0]), label='mod').options(logx=True, logy=True, width=500, show_grid=True).redim.range(y=(0.03,0.3))

g = fig*fig2

out = hv.save(fig,'iwagaki_org.html')
out = hv.save(g,'iwagaki_mod.html')
del out

グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

はてブロへのhtmlグラフの埋め込みは以下を参照

computational-sediment-hyd.hatenablog.jp