自然河道断面の浅水流近似一次元流モデル - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

ようやくまとめることができたので記事にします。

f:id:SedimentHydraulics:20210724155015p:plain

この記事のポイント

実河川の一次元計算モデルの構築：ソースコードあり
断面諸元の計算をクラス化するとコーディングは簡単
pythonは科学技術計算に使える

支配方程式のソースコード
実河川の計算例
Environment
Github

支配方程式のソースコード

浅水流近似した一次元流れの式（Saint-Venentの式)を自然河道断面（一般断面とも呼ばれる）に対応できるように河道横断面積分した次の式を基本とします。

$\begin{align} & \frac{\partial Q}{\partial t} + \frac{\partial }{\partial x}\left(\frac{\beta Q^2}{A}\right)+gA\frac{\partial H}{\partial x}+gAi_e = 0 \\ & \dfrac{\partial A}{\partial t}+\dfrac{\partial Q}{\partial x} = 0 \end{align}$

以下の条件よりソースコードを作成しています。

スキームはいつもどおりです。（過去2年間で10河川以上このモデルで計算しているので多分安定しているはずです。）
一次元流れなので逆流はないと想定して分岐を最小化しています。
欠点は上流端の水位（河積）の取り扱いです。横断面形状が極端に不連続な場合は不安定になりますので仮想断面を追加して対応しています。（今回は不要）
断面内諸元の計算は、別途クラスを作成してそちらで対応しています。気が向けば解説を書きますが、慣れている人なら容易に理解できます。

def UnSteadyflow(sections, A, Q, H, Abound, Qbound, dt):
    g = float(9.8)
    imax = len(A)
    Anew, Qnew, Hnew = np.zeros(imax), np.zeros(imax), np.zeros(imax)
    ie = np.zeros(imax)
    Beta = np.zeros(imax)
    
# continuous equation
    for i in range(1, imax-1) : 
        dx = 0.5*(sections[i-1].distance - sections[i+1].distance)
        Anew[i] = A[i] - dt * ( Q[i] - Q[i-1] ) / dx
        
    Anew[imax-1] = Abound
    Anew[0] = Anew[1]
#     Anew[0] = (Anew[1] - A[1]) + A[0]
    
    for i in range(imax) : 
        s = sections[i]
        Hnew[i], _, _ = s.A2HBS(Anew[i], H[i])
        arr = s.calIeAlphaBetaRcUsubABS(Q[i], H[i])
        ie[i] = arr[0]
        Beta[i] = arr[2]
    
# moumentum equation
    for i in range(1, imax-1): 
        ic, im, ip = i, i-1, i+1
        dxp = sections[ic].distance - sections[ip].distance
        dxm = sections[im].distance - sections[ic].distance
        dxc = 0.5*(sections[im].distance - sections[ip].distance)
        
        Cr1 = 0.5*( Q[ic]/A[ic] + Q[ip]/A[ip] )*dt/dxp
        Cr2 = 0.5*( Q[ic]/A[ic] + Q[im]/A[im] )*dt/dxm
        dHdx1 = ( Hnew[ip] - Hnew[ic] ) / dxp
        dHdx2 = ( Hnew[ic] - Hnew[im] ) / dxm
        dHdx = (float(1.0) - Cr1) * dHdx1 + Cr2 * dHdx2
        
        Qnew[ic] = Q[ic] - dt * ( Beta[ic]*Q[ic]**2/A[ic] - Beta[im]*Q[im]**2/A[im] ) / dxc \
                         - dt * g * Anew[ic] * dHdx \
                         - dt * g * A[ic] * ie[ic] 
        
    Qnew[imax-1] = Qnew[imax-2]
    Qnew[0] = Qbound
        
    return Anew, Qnew, Hnew

全てのコードはこちら

実河川の計算例

対象河川

対象河川は8909090001です。

対象範囲は下図のとおりで観測データの存在状況より、3.000～17.600としています。

検討対象出水

水文水質データベースより、流量観測所Y地点の1996年以降で最大流量を記録した2012年7月12日の出水を対象としました。

下流端はO地点の観測水位を与え、上流端は計算範囲上流端の少し下流に位置するY地点の観測流量を与える。

上流端流量および下流端水位は下図のとおりです。

計算条件

マニングの粗度係数は、低水路0.025、高水敷は0.035を与える。
区間距離は河道中心線に沿った距離より算出した。
境界条件は上記のとおり

計算

全コードは少し長くなるので主要な部分のみ解説します。

計算は、初期条件計算と本計算の2つに分類できます。

初期条件計算

まず、初期時点の流量、水位を条件とした不等流計算により水位を設定し、それを初期条件とした非定常計算による助走計算を20時間実施し、その結果を初期値として与えた。

# nonuniform flow
Qt = np.full(len(chD2U), Qup[0], dtype=np.float64)
Huni = model.NonUniformflow(chD2U, Qt, Hdown[0])
Auni = np.array( [chD2U[i].H2ABS(hh)[0] for i, hh in enumerate(Huni)] )
A0 = Auni[0]

# unsteady flow : 20hr
Q = np.full_like(Auni, Qup[0])
A, H = Auni[::-1], Huni[::-1]
for n in range(1, int(3600*20/dt)):
    A, Q, H = model.UnSteadyflow(chU2D, A, Q, H, A0, Qup[0], dt)

本計算

上記の初期条件、前述の境界条件より計算する。

nmax = len(chU2D)
ntmax = len(Qup)
Amat, Qmat, Hmat = np.zeros((ntmax, nmax)), np.zeros((ntmax, nmax)), np.zeros((ntmax, nmax))

Amat[0], Qmat[0], Hmat[0] = A, Q, H

for n in range(1, ntmax):
    A0, _, _ = chU2D[-1].H2ABS(Hdown[n])
    A, Q, H = model.UnSteadyflow(chU2D, A, Q, H, A0, Qup[n], dt)
    Amat[n], Qmat[n], Hmat[n] = A, Q, H

Export

出力はxarrayを用いてNetCDF形式のファイルを作成する。

gdfsectsorted = gdfsect.sort_values('distancefromDB', ascending=False)

ds = xr.Dataset( {'A': (['t','x'], Amat)
                ,'Q': (['t','x'], Qmat)
                ,'H': (['t','x'], Hmat)
                ,'zbmin': (['x'], [c.zbmin() for c in chU2D] )
                 }
                , coords={'x':gdfsectsorted['distancefromDB'].values
                          , 't':dfhydroip.index.values} 
                ,attrs = {'name':gdfsectsorted['name'].values.tolist()}
               )

dout = ds.to_netcdf(r'calout.nc')
del dout