cip1d_uns_bed_v_1_0_0

cip1d_uns_bed_v_1_0_0.lzhをダウンロードして解凍すると、「cip1d_uns_bed_v_1_0_0.xls」ができます。 この「cip1d_uns_bed_v_1_0_0.xls」の使用方法を説明します。 ダウンロードについては、水理水文ソフトウェアのHPを参照してください。

本マクロは、各条件シートに条件を入力するだけで、CIP法による1次元不定流河床変動計算が行えます。

参考URL

http://ws3-er.eng.hokudai.ac.jp/yasu/

http://ja.wikipedia.org/wiki/CIP%E6%B3%95

http://homepage.mac.com/galois21/fortran/lesson/index.html

http://www.civil.chuo-u.ac.jp/lab/eisei/member/2008/suirigaku8.pdf

http://ws3-er.eng.hokudai.ac.jp/yasu/pdf_files/Textbook%20for%20master%20course%28Japanese%29.pdf

「計算条件」シートに計算パラメータを入力します。また、CIP法の基本式と差分化について「計算条件」シートに書いていますので、参考資料とともに参照して下さい。

　計算結果は、以下の１４のシートに書き出されます。　グラフについては、使用者各自に作成して下さい。

• 計算結果出力シート : 出力間隔毎の計算結果
• 計算河床変動量シート
• 計算河床高シート
• 計算水深シート
• 計算水位シート
• 計算流積シート
• 計算流量シート
• 計算流速シート
• 計算フルード数シート
• 計算平均粒径シート
• 計算流砂量シート
• 計算累積流砂量シート
• 計算粒径別累積流砂量シート : 出力間隔毎の計算結果
• 計算粒径区分割合シート : 出力間隔毎の計算結果

　CIP法による1次元不定流計算の河床変動の計算部分のソースコードを公開します。このマクロについては、ソースファイルはすべて公開しています。

' ***********************************
' 河床変動量の計算ルーチン
'
' ***********************************
Public Sub calRVA()
    
    Dim i As Integer
    
    '最下流端の流量を定義
    qc(nx) = qc(nx - 1)

    '掃流砂量と浮遊砂量の計算
    Call calQBQF

    Call calConcent
    
    Call calZ
        
    Call calZsib
    
    '累積値の計算
    For i = 1 To nx - 1
        qbsum(i) = qbsum(i) + qb(i)
        qfsum(i) = qfsum(i) + qf(i)
        dqbsum(i) = dqbsum(i) + dqb(i)
        dzdtsum(i) = dzdtsum(i) + dzdt(i) * dt
    Next i

End Sub


' ***********************************
' 河床材料データの読み込み
'
' ***********************************
Public Sub bedinput()

    Dim i As Integer
    Dim k As Integer
    
    '粒径分布は対数正規分布(8区分)
    For k = 1 To ndx
        dd(1) = 0.1
        dd(2) = 0.2
        dd(3) = 0.5
        dd(4) = 0.7
        dd(5) = 1#
        dd(6) = 2#
        dd(7) = 4#
        dd(8) = 10#
    Next k
    
    Const d10_0 As Double = 0.2     'mm
    Const d50_0 As Double = 0.5     'mm
    Const d90_0 As Double = 1#      'mm
    
    For i = 1 To nx
       d10(i) = d10_0
       d50(i) = d50_0
       d90(i) = d90_0
    Next i

    nd = ndx - 1
    For k = 2 To nd
        d(k) = (dd(k) + dd(k - 1)) * 0.5 / 1000#
    Next k
    
    
    '初期粒径分布の設定
    Call initgrain(nx, nd)
    
    '沈降速度の計算
    Call fallvel(nd)
    
    '限界流速の計算(critical shear velocity)
    For k = 2 To nd
        Call usc(d(k), uci(k))
        tsci0(k) = uci(k) ^ 2 / (s * g * d(k))
    Next k

    '初期平均粒径の計算
    For i = 1 To nx
        dm0(i) = 0#
        dm(i) = 0#
        For k = 1 To nd
            dm0(i) = dm0(i) + d(k) * sib(i, k)
            dm(i) = dm(i) + d(k) * sib(i, k)
        Next k
    Next i

    '平均粒径の計算
    Call setDM

End Sub


' ***********************************
'
Public Sub setDM()
    
    Dim i As Integer
    Dim k As Integer
    Dim ucm As Double
    
    '
    '------  dm ------
    '
    For i = 1 To nx
        dm(i) = 0#
        For k = 1 To nd
            dm(i) = dm(i) + d(k) * sib(i, k)
        Next k
    Next i

    '
    ' ------  tau*cm ------
    '
    For i = 1 To nx
        Call usc(dm(i), ucm)
        tscm(i) = ucm ^ 2 / (s * g * dm(i))
    Next i


End Sub



' ***********************************
'
Private Sub usc(d As Double, uc As Double)
    Dim rst As Double
    
    rst = Sqr(s * g) * d ^ (3# / 2#) / snu
    If rst <= 2.14 Then uc = 0.14 * s * g * d
    If rst > 2.14 Then uc = (0.1235 * s * g) ^ (25# / 32#) * snu ^ (7# / 16#) * d ^ (11# / 32#)
    If rst > 54.2 Then uc = 0.034 * s * g * d
    If rst > 162.7 Then uc = (0.01505 * s * g) ^ (25# / 22#) * snu ^ (-3# / 11#) * d ^ (31# / 22#)
    If rst > 671# Then uc = 0.05 * s * g * d
    
    uc = Sqr(uc)

End Sub


' ***********************************
'
Private Sub omega(bs As Double, tsi As Double, ome As Double)
    Dim a1 As Double, a2 As Double, a3 As Double
    Dim ad As Double
    Dim t As Double, zx As Double, px As Double
    Dim er1 As Double, er As Double
    Dim xxx As Double
    
    a1 = 0.4361836
    a2 = -0.1201676
    a3 = 0.937298
    
    If tsi <= 0.0000001 Then
        ome = 0#
        Exit Sub
    End If
    
    ad = bs / tsi - 2#
    If ad >= 0# Then
        xxx = ad * Sqr(2#)
    End If
    If ad < 0# Then
        xxx = -ad * Sqr(2#)
    End If
    
    t = 1# / (1# + 0.33627 * xxx)
    zx = 1# / Sqr(2# * pi) * Exp(-xxx ^ 2# / 2#)
    px = 1# - zx * (a1 * t + a2 * t ^ 2 + a3 * t ^ 3)
    
    er1 = 2# - 2# * px
    If ad >= 0# Then
        er = er1 / 2#
    End If
    If ad < 0# Then
        er = (2# - er1) / 2#
    End If
    
    If (bs = 0# Or er = 0#) Then
        ome = 0#
    Else
        ome = tsi / bs / (2# * Sqr(pi)) * Exp(-ad ^ 2) / er + tsi * 2# / bs - 1#
    End If
    
    'if(ome.lt.1e-10) ome=0.

End Sub


' ***********************************
'
Private Function alf(bet As Double)
    'if bet > 20. then
    '   alf = bet
    'else
       alf = bet / (1# - Exp(-bet))
    'end if

End Function


' ***********************************
'
Private Function ppx(dd As Double, d10 As Double, d50 As Double, d90 As Double)
    Dim s1 As Double, s2 As Double
    Dim x As Double
      
    s1 = 1.282 / (Log10(d90) - Log10(d50))
    s2 = -1.282 / (Log10(d10) - Log10(d50))

    If dd >= d50 Then
        x = (Log10(dd) - Log10(d50)) * s1
    Else
        x = (Log10(dd) - Log10(d50)) * s2
    End If
    
    ppx = pc(x)

End Function


' ***********************************
'
Private Function tfunc(x As Double)
    If x >= 0# Then
        tfunc = 1# / (1# + 0.33267 * x)
    Else
        tfunc = 1# / (1# - 0.33267 * x)
    End If

End Function


' ***********************************
'
Private Function pc(x As Double)

    Dim a1 As Double
    Dim a2 As Double
    Dim a3 As Double
    Dim tx As Double
    
    
    a1 = 0.4361836
    a2 = -0.1201676
    a3 = 0.937298
    tx = tfunc(x)
      
    If x >= 0# Then
       pc = 1# - 1# / Sqr(2# * pi) * _
       Exp(-x ^ 2 / 2#) * (a1 * tx + a2 * tx ^ 2 + a3 * tx ^ 3)
    Else
       pc = 1# / Sqr(2# * pi) * _
       Exp(-x ^ 2 / 2#) * (a1 * tx + a2 * tx ^ 2 + a3 * tx ^ 3)
    End If
    
End Function


' ***********************************
'
Private Function max(a As Double, b As Double)
    max = a
    If max < b Then max = b
End Function


' ***********************************
'
Private Function min(a As Double, b As Double)
    min = a
    If min > b Then min = b
End Function


' ***********************************
'
Private Function Log10(ByVal Val As Double)
    Log10 = Log(Val) / Log(10)
End Function


' ***********************************
' 初期粒径の設定
'
' ***********************************
Private Sub initgrain(nx As Long, nd As Long)
    Dim i As Long
    Dim k As Long
    
    For i = 1 To nx
        p0(i, 1) = 0#
        p0(i, nd) = 1#
        For k = 2 To nd - 1
            p0(i, k) = ppx(dd(k), d10(i), d50(i), d90(i))
        Next k
    Next i
    For i = 1 To nx
        For k = 2 To nd
            sib(i, k) = p0(i, k) - p0(i, k - 1)
            sib0(i, k) = sib(i, k)
        Next k
        For k = 1 To nd
            p(i, k) = p0(i, k)
        Next k
    Next i
End Sub


' ***********************************
' 沈降速度の計算
'
' ***********************************
Private Sub fallvel(nd As Long)
    Dim k As Long
    
    For k = 2 To nd
        If d(k) > 0.001 Then
            w0(k) = 32.8 * Sqr(d(k) * 100#) * 0.01
        Else
            w0(k) = Sqr(2# / 3# * s * g * d(k) + (6# * snu / d(k)) ^ 2) - 6# * snu / d(k)
        End If
    Next k

End Sub


' ***********************************
' 粒径別掃流砂量と浮遊砂量の計算
' 芦田・道上の式
' ***********************************
Private Sub calQBQF()
    
    Dim i As Integer
    Dim k As Integer
    Dim tsci As Double
    Dim usci As Double
    Dim xi As Double
    Dim xxi As Double
    Dim bs As Double
    Dim ome As Double
    
    Dim si(nix) As Double
    Dim a(nix) As Double
    Dim u(nix) As Double
        
    For i = 1 To nx
        a(i) = b(i) * hs(i)
        u(i) = qc(i) / a(i)
        si(i) = (mn(i) * qc(i) / b(i)) ^ 2 / hs(i) ^ 3.3333
        us(i) = Sqr(g * hs(i) * si(i))
        tsm(i) = us(i) ^ 2 / (s * g * dm(i))
    Next i
    
    '
    'qb (bed load) & qf (suspended load)
    '
    For i = 1 To nx
        qb(i) = 0#
        qf(i) = 0#
        For k = 2 To nd
            tsi(i, k) = us(i) ^ 2 / (s * g * d(k))
            xi = (Log10(23#) / Log10(21# * d(k) / dm(i) + 2#)) ^ 2
            tsci = xi * tscm(i)
            If tsi(i, k) <= tsci Then
                qbi(i, k) = 0#
                qfi(i, k) = 0#
            Else
                usci = Sqr(tsci * s * g * d(k))
                qbi(i, k) = sib(i, k) * 17# * Sqr(s * g * d(k) ^ 3) * _
                            tsi(i, k) ^ 1.5 * (1# - tsci / tsi(i, k)) * (1# - usci / us(i))
                If qbi(i, k) < 1E-20 Then qbi(i, k) = 0#
                qb(i) = qb(i) + qbi(i, k)

                xxi = tsci / tsci0(k)
                bs = bs0 * xxi
                Call omega(bs, tsi(i, k), ome)
                qfi(i, k) = sib(i, k) * sk * (als * ome * s * g * d(k) / (rs * us(i)) - w0(k))
                If qfi(i, k) < 1E-20 Then qfi(i, k) = 0#
                qf(i) = qf(i) + qfi(i, k)
            End If
        Next k
    Next i
    
End Sub


'
' ------ cal. of c ------
'
Private Sub calConcent()

    Dim i As Integer
    Dim k As Integer
    
    Dim bet As Double
    Dim alfx As Double
    Dim c(nix, ndx) As Double
    Dim ct(nix) As Double
    Dim last As Integer
    
    last = nx
    
    '濃度の計算
    For k = 2 To nd
    
        bet = 15# * w0(k) / us(last)
        alfx = alf(bet)
        c(last, k) = qfi(last, k) / (w0(k) * alfx)
        cb(last, k) = c(last, k) * alfx
        
        For i = last - 1 To 1 Step -1
            bet = 15# * w0(k) / us(i)
            alfx = alf(bet)
            c(i, k) = (yu(i) * c(i + 1, k) + dx(i) * b(i) * qfi(i, k)) / (yu(i) + dx(i) * b(i) * w0(k) * alfx)
            If c(i, k) <= 0.0000000001 Then
                c(i, k) = 0#
            End If
            cb(i, k) = c(i, k) * alfx
        Next i
        
    Next k
    
    'wcとctの計算
    For i = 1 To nx
        ct(i) = 0#
        wc(i) = 0#
        For k = 2 To nd
            ct(i) = ct(i) + c(i, k)
            wc(i) = wc(i) + w0(k) * cb(i, k)
        Next k
    Next i

End Sub


'
' ------ cal. of z ------
'
Private Sub calZ()
    
    Dim i As Integer
    Dim k As Integer
    Dim dq As Double
    
    For i = 2 To nx
    
        dq = (qb(i - 1) * b(i - 1) - qb(i) * b(i)) / (b(i) * dx(i))
        dzdt(i) = (dq + wc(i) - qf(i)) / (1# - ramda)
        dqb(i) = dzdt(i)
        'dzdt(i) = -1 * (dqb(i) + wc(i) - qf(i)) / (1# - ramda)
        
    Next i
    
    'dqb(1) = (qb(1) * b(1) - 0#) / (b(i) * dx(i))
    dq = (qb(1) * b(1) - qb(1) * b(1)) / (b(1) * dx(1))
    dzdt(1) = (dq + wc(1) - qf(1)) / (1# - ramda)
    dqb(1) = dzdt(1)

End Sub


'
' ------ cal. of z, sib ------
'
Private Sub calZsib()
    
    Dim i As Integer
    Dim k As Integer
    
    Dim dletz As Double
    Dim xi As Double
    Dim dqbi As Double
    Dim dpi As Double
    Dim ax As Double
    
    ax = 1# / dt
    
    For i = 2 To nx
        dletz = dt * dzdt(i)
        z(i) = z(i) + dletz
        
        For k = 2 To nd
            If dletz > 0# Then
                xi = sib(i, k)
            Else
                xi = sib0(i, k)
            End If
                dqbi = (qbi(i, k) * b(i) - qbi(i - 1, k) * b(i - 1)) / (b(i) * dx(i))
            dpi = ax * (dt / (1# - ramda) * (dqbi + w0(k) * cb(i, k) - qfi(i, k)) - xi * dletz)
            sib(i, k) = sib(i, k) + dpi
            If sib(i, k) < 0.0000000001 Then sib(i, k) = 0#
        Next k
        
    Next i

End Sub