不等流計算

独立行政法人 土木研究所 寒地土木研究所(旧名:北海道開発局土木試験所河川研究室)で発行された「現場のための水理学 (単断面における不等流計算)」を参考に抜粋してまとめました。オリジナルは土木試験所月報No.411号(1987.8)~415号(1987.12)まで5回に渡って連載されたものです。

詳細な内容を勉強される方は、現場のための水理学から入手してください。

なお、掲載したExcel VBAについては、弊社にてオリジナルプログラムを修正・加筆しました。

基礎式

任意の2つの断面間で流れの状態を考える。(下図を参照)想定した2つの断面を、各々断面1、断面2とし、各断面の諸量には断面番号の下付き数字をつける。また,図中にhfは2つの断面の総水頭差、すなわち、各々の断面で流体がもつエネルギ-の差を表わして、流体が断面間を流下するときに失なわれるエネルギ-(水頭)である。この損失を摩擦損失水頭と呼び。

河道縦断模式図

以上のことから,任意の2つの断面間では,次式が成り立つことがわかる。

f hgh Vgh V +2= z + +2z + +・(2.3.1)

ここで,z;河床高(m),h;水深(m),V;流速(m/s),hf;摩擦損失水頭(m),g;重力加速度(m/s2)すなわち,上式が不等流を表す式であります。そして,この式は水頭に換算されたエネルギ-が任意の個所で等しいという、“エネルギ-保存の法則“となっている。 また,ここで不等流の場合,河床勾配ib,水面勾配iwおよびエネルギ-勾配ie は,水深や流速が違う分,各々異なるということに注意しておいて下さい。(2.3.1)式が任意断面間でなりたつ不等流の式です。

次に(2.3.1)式を一般化する。 まず,図-2.1 では,対象とする2つの断面間の距離はl(m)でしたが,これを微少長さ⊿x(m)とします。つまり,微少区間を考えることにより,より細かな現象を捉えることを可能とするわけです。 次に,摩擦損失水頭hf ですが,これは図-2.1 を見てもらえば,エネルギ-勾配がie=hf /l,すなわち,l を⊿x とおくと,ie=hf /⊿x となることがわかり,これから, ・・・・・・・・・・・・・・・・・・・・・・(2.2) と表すことができます。 以上のことを念頭において,式(2.1)を書きなおすと, ・・・・・・・・・(2.3) 上式において、右辺から左辺を引いて,各項を⊿x で割ると, ・・・・・・・・・・・・・・・・・・・・・・・・(2.4) ここで,⊿x をさらに微少量として,すなわち,限りなく0 に近い長さとしてdx で表わし,また,z やh,V2/2g の変化量もdz,dh,d(V2/2g)で表わすと式(2.4)は, ・・・・・・・・・・・・・・・・・(2.5) となります。また,水位H はH=h+z であるので, ・・・・・・・・・・・・・・・・・・(2.6) となります。このように表わした式(2.5)および式(2.6)が,不等流における流れの微少変化を表した”微分表示”というものです。 ところで,式(2.5)および式(2.6)にあるエネルギ-勾配は、実際にどのようにして算出すればよいでしょう か。このie は,下記の式(2.7)を用いて算出することができます。 ・・・・・・・・・・(2.7) ここで,V;流速(m/s), n;粗度係数,R;径深(m),ie;エネルギ-の勾配 これが,不等流の平均流速式です。 なお,このような形をもつ平均流速式をマニング型の式といい,今後,頻繁に使われることになります。また,この中にでてくるn を(マニングの)粗度係数といい,河床の粗さや形状などに左右される抵抗を表わすものといわれており,厳密には,きわめて高度な科学的見地より決められねばなりませんが,実用上は,実測値からの逆算や推定などからあらかじめ決められている既知値として扱っていくこととします。さらに,R は径深というもので,流水断面積もしくは流積A(m2)を潤辺(断面において水に潤っている辺の長さ)S(m)で割ったもので, ・・・・・・・・・(2.8) と表わされるものです。この辺の概念は,図-2.2 をみて把握して下さい。一般の大部分の河川のように,水深に比して川幅が広い場合は,径深R は水深h に等しいとみなしてさしつかえありません。すなわち, ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・(2.9) となり,このような断面を 広短形断面 と称しています。 さて,式(2.7)からie を求めると, ・・・・・・・・・・・・・(2.10) また,平均流速V は流量Q (m3/s)を流積A (m2)で割ったものであるので,式(2.10)は, ・・・・・・・・・・・・・(2.11) となり,これを式(2.6)に代入して,V=Q/A などとおくと, ・・・・・・・・・・・(2.12) となります。この式(2.12)は,最終的に不等流を表わす基礎式となり,今後,頻繁に使われるものなので,確実に覚えておいて下さい。 ここで,特に広矩形断面の場合においては,R~h とおくことができました。また,流積A は川幅B(m)と水深h(m)の積で表わされることから,式(2.12)は, (2.13) となり,これが広矩形断面をもつ水路の不等流計算の基礎式となります。

単断面の不等流計算

それでは,いよいよ実際の不等流計算に臨んでみることにしましょう。初めて計算を行う人のために,本節では不等流計算の中で最も初歩的な広矩形単断面の例から考えていきたいと思います。

もう一度,図-2.1 および式(2.13)を見て下さい。対象とする最小の計算区間を図-2.1 のように,上流側断面1,下流側断面2の間と設定します。ここで,式(2.13)を実際の計算に用いることができるように,諸量の差をもって表わすこととします。ここで注意しておきたいのは流量Q は一定としていること,第三項目は断面1での値と断面2での値の平均となっていることです。

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・(2.25)

ゆえに,両辺に⊿x をかけ,左辺に下流側,右辺に上流側の諸量を示す項をもってくると

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・(2.26)

また,水位H を河床高z と水深h の和で表すと,

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・(2.27)

となり,このような表し方を差分表示といって,今後,基礎方程式を実際に計算機を用いて解く場合,広く使われる手法です。

次に,式(2.27)の中でなにが既知量か,なにが未知量かを考えてみましょう。まず,流量Q,河床高z,河幅B,断面間の距離⊿x は既知量として与えられます。また,重力加速度g,粗度係数n も定数として与えられます。結局,水深h が未知量として求めるべきものとなりますが,式が1本しかないのに,上流側の水深h1 と下流側の水深h2 を未知量として,両方同時に求めることはできません。そこで,どちらか一方の水深が既知量として与えられなければならないのですが,ここにいたって,2-2 節で示された流れの性質が問題となってくるのです。すなわち常流の場合,流れの変化は下流から上流に及び,射流ではこれが上流方向に及ばないため,常に流れは上流から下流に変化します。この事実を計算にも適用し,常流の場合下流から上流方向へ,射流の場合上流から下流方向へ計算を進めていくことにします。つまり,常流の場合下流側のh2 が既知量として与えられ,上流側のh1 を求めていけばよいわけです。なお,断面がい くつかある場合は,下流瑞の水深が境界条件として与えられ,上流に向かって断面間の計算を逐次行っていくような方法をとることになります。

例題

河床高および河幅が表-2.7 に示されるような広矩形断面河川に,流量Q=1500m3/s が流下した場合の水面形を求めよ。ただし,粗度係数n=0.025,下流端の水深を2.5m とする。

河道諸元

断面番号下流端からの距離河床高河幅
(m)(m)(m)
1 0 0.0 300
2 500 0.5 320
3 1000 0.9 280
4 1200 0.8 250
5 1800 2.0 300
6 2100 2.3 300
7 2500 3.0 320
8 3000 3.0 350
9 3300 3.5 300
10 3800 4.0 250

常流であるから,下流から上流に向かって計算するので,下流側水深h2 を既知とし,上流側水深h1 を求めることから,式(2.27)を参考として,

・・・(1)

ただし,

・・・(2)

・・・(3)

・・・(4)

ここでは、1階のニュ-トン法を用いることとする。

(1)はf(h1)=0 となるとき,式(2.27)と同一になる。つまり,f(h1)=0 を満たすようなh1 が断面間の関 係を満足する上流側水深の正しい値といえるわけである。

f(h)=0 となるh を求める方法としては種々の方法があるが,今回はニュ-トン法を用いた。以下にその説明を述べる。

ⅰ)h1 になんらかの仮定値を代入し,f(h1)を計算する。 ⅱ)f(h1)> ε のとき,h1=h1+⊿h1 なる更新を行う。 ただし,⊿h1 は上記のような更新を行ったときにf(h1+⊿h1)を0 とするようなものを採用しなけれ ばならない。そこで,f(h1+⊿h1)を一階微分の項までテイラ-展開し,

f(h1+⊿h1)≒f’(h1)⊿h1=0 ・・・・・・・・・・・・・・・・②

これから,

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・③

とすればよい。 本問の場合,f’(h1)は①式を微分することにより,

・・・・・・・・・・・・・④

で求められる。 ⅲ)ⅱ)で更新したh1 をもって,再度f(h1)の計算を行う。 ⅳ)f(h1)が十分零に近づくまで,すなわち| f(h1)|< εとなるまで上記の手順を繰返す。なお,ε は打切り誤差といい,解の精度を勘案して十分零に近い値,例えばmm の精度で解をだすならε =0.001というような値をとるものである。

VBAプログラム

初期条件シート

計算結果シート

プログラムソース

Option Explicit
Option Base 1
Private Const g As Double = 9.8        '重力加速度
Private Const ndim As Integer = 100    '最大断面数

Private Sub pro03()
    '広矩形断面における不等流計算
    'Newton-Raphson法を用いた水位計算
    '微分項--->第1項まで

    Dim d(ndim) As Double
    Dim dl(ndim) As Double
    Dim z(ndim) As Double
    Dim b(ndim) As Double
    Dim sh(ndim) As Double
    Dim hh(ndim) As Double
    Dim sn(ndim) As Double
    
    Dim j As Integer
    Dim ndanmen As Integer
    Dim q As Double
    Dim n As Double
    Dim eps As Double
    
    Dim h0 As Double
    Dim aa As Double
    Dim bb As Double
    Dim cc As Double
    Dim fh As Double
    Dim f1 As Double
    Dim dh As Double
    
    Dim ws_i As Worksheet
    Dim ws_o As Worksheet
    
    Set ws_i = Sheets("初期条件")
    Set ws_o = Sheets("計算結果")

    '計算条件
    ndanmen = ws_i.Range("_ndanmen")                '河道断面数
    q = ws_i.Range("_q")                            '流量 m^3/s
    n = ws_i.Range("_n")                            '粗度係数
    sh(1) = ws_i.Range("_sh_1")                     '第一断面水深
    eps = ws_i.Range("_eps")                        '打ち切り誤差

    '河道データ読み込み
    For j = 0 To ndanmen - 1
        d(j + 1) = Range("_kyori").Offset(j + 1)    '追加距離の読み込み
        z(j + 1) = Range("_z").Offset(j + 1)        '河床高の読み込み
        b(j + 1) = Range("_b").Offset(j + 1)        '河幅の読み込み
        'Debug.Print j + 1, z(j + 1)
    Next j
    
    hh(1) = sh(1) + z(1)                            '第一断面水位
    
    dl(1) = 0
    '区間距離の計算
    For j = 2 To ndanmen
        dl(j) = d(j) - d(j - 1)
    Next j
    
    '計算開始
    For j = 1 To ndanmen
        sn(j) = 0
        If (j = 1) Then GoTo n34:
        h0 = sh(j - 1)
        '係数aa,bb,ccの計算
        aa = q ^ 2 / (2 * g * b(j) ^ 2)
        bb = -(q ^ 2) * n ^ 2 * dl(j) / (2 * b(j) ^ 2)
        cc = z(j) - (z(j - 1) + sh(j - 1) _
          + q ^ 2 / (2 * g * b(j - 1) ^ 2 * sh(j - 1) ^ 2) _
          + q ^ 2 * n ^ 2 * dl(j) / (2 * b(j - 1) ^ 2 * sh(j - 1) ^ (10# / 3#)))
        'f(hh),f'(hh)の計算
n28:    fh = h0 + aa / (h0 ^ 2) + bb / (h0 ^ (10# / 3#)) + cc            ' f(hh)
        f1 = 1# - 2# * aa / (h0 ^ 3#) - 10 * bb / (3# * h0 ^ (13# / 3#)) ' f'(hh)
        dh = -fh / f1                                                    ' δh
        If (Abs(fh) < eps) Then GoTo n33:                                ' 条件判断
        h0 = h0 + dh                                                     ' hh=hh+dh
        sn(j) = sn(j) + 1
        GoTo n28:
n33:    sh(j) = h0
        hh(j) = sh(j) + z(j)
n34:
    Next j
    
    '結果の出力
    ws_o.Cells(1, 1) = "j"
    ws_o.Cells(1, 2) = "追加距離(m)"
    ws_o.Cells(1, 3) = "区間距離(m)"
    ws_o.Cells(1, 4) = "河幅(m)"
    ws_o.Cells(1, 5) = "河床高(m)"
    ws_o.Cells(1, 6) = "水深(m)"
    ws_o.Cells(1, 7) = "水位(m)"
    ws_o.Cells(1, 8) = "収束回数"
    For j = 1 To ndanmen
        ws_o.Cells(1 + j, 1) = j
        ws_o.Cells(1 + j, 2) = d(j)
        ws_o.Cells(1 + j, 3) = dl(j)
        ws_o.Cells(1 + j, 4) = b(j)
        ws_o.Cells(1 + j, 5) = z(j)
        ws_o.Cells(1 + j, 6) = sh(j)
        ws_o.Cells(1 + j, 7) = hh(j)
        ws_o.Cells(1 + j, 8) = sn(j)
    Next j

End Sub


Private Sub CommandButton1_Click()
    'シート上の「計算開始」ボタンをクリックした場合に実行
    Call pro03
End Sub

一般断面の不等流計算

独立行政法人 土木研究所 寒地土木研究所(旧名:北海道開発局土木試験所河川研究室)で発行された「現場のための水理学 (一般断面における不等流計算)」を現場のための水理学から入手してください。 基本的には、単一断面が基本となりますので、FortranやBasicからのVBAへの変更は比較的簡単に行えます。

VBAプログラム

計算条件

計算結果

プログラムソース

Option Explicit
Option Base 1
Private Const ndim As Integer = 100     '断面数
Private Const soku As Integer = 20      '一断面当たりの測点数
Private Const g As Double = 9.8         '重力加速度


Private Sub nftoryu()
    '任意断面での不等流計算プログラム
    'original programed by y.murakami 2005/11/6 寒冷地土木研究所
    'improbed by TYS 2011/10/15
    '
    '計算条件
    '●下流端の水位
    '●上流からの流量
    
    Dim lx(ndim) As Double
    Dim dx(ndim) As Double
    Dim n(soku) As Integer
    Dim hh(soku) As Double
    Dim tn(soku) As Integer
    Dim y(ndim, soku) As Double
    Dim z(ndim, soku) As Double
    Dim sn(ndim, soku) As Double
    
    Dim yy(soku) As Double
    Dim zz(soku) As Double
    Dim nn(soku) As Double
    
    Dim ndanmen As Integer
    Dim h0 As Double
    Dim q As Double
    Dim eps As Double
    Dim kk As Double
    
    Dim hhh As Double
    Dim aa As Double
    Dim bb As Double
    Dim sf As Double
    Dim lf As Double
    Dim f1 As Double
    Dim f2 As Double
    Dim dh As Double
    
    Dim j As Integer
    Dim i As Integer
    Dim nk As Integer
    
    Dim ws_i As Worksheet
    Dim ws_o As Worksheet
    
    Set ws_i = Sheets("初期条件")
    Set ws_o = Sheets("計算結果")

    ndanmen = ws_i.Range("_ndanmen")            '断面数
    '
    h0 = ws_i.Range("_h0")                      '下流端水位(m)
    q = ws_i.Range("_q")                        '流量(m^3/s)
    eps = ws_i.Range("_eps")                    '許容誤差(m)
    kk = ws_i.Range("_kk")                      '緩和係数
    
        
    Dim nrow As Integer
    nrow = 15
    
    For j = 1 To ndanmen
        
        n(j) = ws_i.Cells(nrow + j - 1, 3)              '断面毎の測点数
        lx(j) = ws_i.Cells(nrow + j - 1, 4)             '下流端断面からの累加距離
    
        nrow = nrow + 1
        For i = 1 To n(j)                               '水平方向距離(m)
            y(j, i) = ws_i.Cells(nrow + j - 1, i + 1)
        Next i
        nrow = nrow + 1
        For i = 1 To n(j)                               '河床高さ(m)
            z(j, i) = ws_i.Cells(nrow + j - 1, i + 1)
        Next i
        nrow = nrow + 1
        For i = 1 To n(j)                               'マニングの粗度係数
            sn(j, i) = ws_i.Cells(nrow + j - 1, i + 1)
        Next i
        'nrow = nrow + 1
    Next j
    
    '区間距離
    For j = 2 To ndanmen
        dx(j) = lx(j) - lx(j - 1)
    Next j

    '
    hh(1) = h0
    tn(1) = 0
    For j = 2 To ndanmen
        hhh = hh(j - 1)

        '*********f2(下流側 計算)**********
        aa = q ^ 2 / (2 * g)
        bb = q ^ 2 * dx(j) / 2#
        nk = n(j - 1)
        For i = 1 To nk
            nn(i) = sn(j - 1, i)
            zz(i) = z(j - 1, i)
            yy(i) = y(j - 1, i)
        Next i
        Call water_height(nk, hhh, nn, zz, yy, sf, lf) 'φψの計算
        f2 = hhh + aa * sf / lf ^ 3# + bb * 1# / lf ^ 2

        '*********f1(上流側 計算)**********
        aa = q ^ 2 / (2 * g)
        bb = q ^ 2 * dx(j) / 2#
        hhh = hh(j - 1) + z(j, Int(n(j) / 2#))
        nk = n(j)
        For i = 1 To nk
            nn(i) = sn(j, i)
            zz(i) = z(j, i)
            yy(i) = y(j, i)
        Next i
n30:    Call water_height(nk, hhh, nn, zz, yy, sf, lf) 'φψの計算
        f1 = hhh + aa * sf / lf ^ 3# - bb * 1# / lf ^ 2
        dh = f1 - f2
        If (Abs(dh) < eps) Then
            hh(j) = hhh
        Else
            hhh = hhh - dh * kk
            tn(j) = tn(j) + 1
            GoTo n30:
        End If
        Debug.Print j, hh(j), f1, f2, dh
    Next j

    '結果の出力
    ws_o.Cells(1, 1) = "j"
    ws_o.Cells(1, 2) = "追加距離(m)"
    ws_o.Cells(1, 3) = "区間距離(m)"
    ws_o.Cells(1, 4) = "水位(m)"
    ws_o.Cells(1, 5) = "収束回数"
    For j = 1 To ndanmen
        ws_o.Cells(1 + j, 1) = j
        ws_o.Cells(1 + j, 2) = lx(j)
        ws_o.Cells(1 + j, 3) = dx(j)
        ws_o.Cells(1 + j, 4) = hh(j)
        ws_o.Cells(1 + j, 5) = tn(j)
    Next j
    
End Sub

Private Sub water_height(nk As Integer, hhh As Double, _
                         nn() As Double, zz() As Double, _
                         yy() As Double, sf As Double, lf As Double)
    '
    '水面処理(φψの計算)*******************
    'sf;ψ
    'lf;φ
    
    Dim sh(soku) As Double
    Dim i As Integer
    Dim dy As Double
    Dim sfi As Double
    Dim sli As Double
    
    sf = 0
    lf = 0
    For i = 1 To nk
        sh(i) = hhh - zz(i)
    Next i
    For i = 1 To nk - 1
        dy = yy(i + 1) - yy(i)
        If zz(i) > hhh And zz(i + 1) > hhh Then GoTo w50:
        If zz(i) > hhh And hhh >= zz(i + 1) Then
            sfi = 1# / nn(i) ^ 3# * (sh(i + 1) / 2#) ^ 3# _
                * sh(i + 1) * dy / (zz(i) - zz(i + 1)) / 2#
            sli = 1# / nn(i) * (sh(i + 1) / 2#) ^ (5# / 3#) _
                * sh(i + 1) * dy / (zz(i) - zz(i + 1)) / 2#
        End If
        If zz(i) <= hhh And hhh < zz(i + 1) Then
            sfi = 1# / nn(i) ^ 3# * (sh(i) / 2#) ^ 3# _
                * sh(i) * dy / (zz(i + 1) - zz(i)) / 2#
            sli = 1# / nn(i) * (sh(i) / 2#) ^ (5# / 3#) _
                * sh(i) * dy / (zz(i + 1) - zz(i)) / 2#
        End If
        If zz(i) <= hhh And zz(i + 1) <= hhh Then
            sfi = 1# / nn(i) ^ 3# * ((sh(i) + sh(i + 1)) / 2#) ^ 3# * dy
            sli = 1# / nn(i) * ((sh(i) + sh(i + 1)) / 2#) ^ (5# / 3#) * dy
        End If
        sf = sf + sfi
        lf = lf + sli
w50:
    Next i

End Sub


Private Sub CommandButton1_Click()
    Call nftoryu
End Sub

ダウンロード

 
水理計算/1.河川の流れ/3.不等流計算.txt · 最終更新: 2012/07/30 23:07 by tys
 
特に明示されていない限り、本Wikiの内容は次のライセンスに従います: CC Attribution-Share Alike 3.0 Unported
Recent changes RSS feed Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki