cip1d_uns_v_2_2_1

cip1d_uns_v_2_2_1.lzhをダウンロードして解凍すると、「cip1d_uns_v_2_2_1.xls」ができます。 この「cip1d_uns_v_2_2_1.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法の基本式と差分化について「計算条件」シートに書いていますので、参考資料とともに参照して下さい。

• 計算結果の出力時間間隔
• 計算終了時間
• 計算刻み時間
• 最小水深
• 繰り返し計算の打ち切り誤差
• 繰り返し計算回数
• 繰り返し計算の緩和係数
• 初期水深をデータで与える : True or False
• 最下流端の水深をデータで与える : True or False

• 距離
• 断面間距離
• 河床高
• 粗度係数
• 川幅
• 初期水位
• 初期水深 : 初期水深をデータで与える場合に必要となる。

• データ数
• 経過時間 : 必ず0から始まる。
• 流量

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

• 入力チェックシート
• 結果シート : 出力間隔毎の計算結果
• 水深結果シート
• 水位結果シート
• 流速結果シート
• 流量結果シート

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

Option Base 0
Option Explicit

'計算キャンセルのフラグ
Public windowflag As Boolean

'初期水深の与えた方のフラグ
Public hflag As Boolean         'ture:シートのデータ      false:等流計算

'初期水深の与えた方のフラグ
Public hinflag As Boolean       'ture:シートのデータ      false:等流計算

'最大配列数の宣言
Public Const im As Long = 300   '河道データ最大数
Public Const ih As Long = 100   '流量データ最大数

'変数の宣言
Public nx As Long               '断面データ数
Public dt As Double             '計算刻み時間
Public hmin As Double           '最小水深
Public time As Double           '累積時間
Public h0 As Double             '初期条件（等流計算時）
Public qp As Double             '流量
Public hdw As Double            '
Public errmax As Double         '打ち切り誤差
Public err As Double            '誤差
Public alh As Double            '緩和係数
Public etime As Double          '計算終了時間
Public np As Double             '流量データ数
Public hdown As Double          '
Public hin As Double            '

'配列の宣言
Public yu(im) As Double
Public yun(im) As Double
Public gux(im) As Double
Public qc(im) As Double
Public z(im) As Double
Public z_up(im) As Double
Public x(im) As Double
Public x_up(im) As Double
Public b(im) As Double
Public b_up(im) As Double
Public h(im) As Double
Public hn(im) As Double
Public hs(im) As Double
Public hsn(im) As Double
Public z0(im) As Double
Public z0_up(im) As Double
Public cfx(im) As Double

Public mn(im) As Double         'マニングの粗度係数
Public dx(im) As Double         '断面間距離
Public dx2(im) As Double        '断面間距離の2乗

'流量関連の配列宣言
Public thyd(ih) As Double       '流量データの時間
Public qhyd(ih) As Double       '流量データ
Public hhyd(ih) As Double       '最下流端の水深

'定数の宣言
Private Const g As Double = 9.8
Private Const pi As Double = 3.1415


'-------------------------------------------------------------------
' CIP法による1次元自由水面流れ計算
' 定常計算
'
' 北大工学部　清水先生のプログラムに加筆
' 参考URL http://ws3-er.eng.hokudai.ac.jp/yasu/
'
' 本プログラムの機能:
' 任意の川幅
' 任意の断面間距離
' 任意の初期水深条件
'
' 制約条件:
' エクセルのシートの制約により
' 計算結果が255列までしか出力できません。
' それ以上はエラーになります。
'
'
' 合同会社TYS
' 技術開発部
' 11/11/2011 ver 2.1.0 公開版
' 11/15/2011 ver 2.1.1 最下流端の水深を計算で与えた場合のバグ修正
'                      hinflagの追加
'                      dslope変数の追加
'-------------------------------------------------------------------
Public Sub cip1d_b()
    
    Dim i As Long
    Dim nq As Long
    Dim sst As Double
    Dim qinput As Double
    Dim slope As Double

    Dim otime As Double
    Dim dtm As Double
    Dim lmax As Long
    
    Dim itout As Long
    Dim itt As Long
    
    Dim icount As Long
    Dim lcount As Long
    Dim ncount As Long
    Dim mcount As Long

    'フォームの表示
    UserForm1.Show vbModeless
    UserForm1.Repaint
    windowflag = False
    
    '計算結果シートのデータ削除
    Call utlModule.SheetDataDel
    
    '流量データの読み込み
    Dim ws_thq As Worksheet
    Set ws_thq = Sheets("流量データ")
    nq = ws_thq.Range("_nq") - 1           '流量データ数
    For i = 0 To nq
        thyd(i) = ws_thq.Range("_time").Offset(i + 1, 0)    '時間
        qhyd(i) = ws_thq.Range("_q").Offset(i + 1, 0)       '流量
        hhyd(i) = ws_thq.Range("_h").Offset(i + 1, 0)       '水深
    Next i
    
    '計算条件の読み込み
    Dim ws_cond As Worksheet
    Set ws_cond = Sheets("計算条件")
    otime = ws_cond.Range("_otime")         '計算出力時間
    dtm = ws_cond.Range("_dtm")             '計算刻み時間
    hmin = ws_cond.Range("_hmin")           '最小水深
    errmax = ws_cond.Range("_errmax")       '打ち切り誤差
    lmax = ws_cond.Range("_lmax")           '繰り返し回数
    alh = ws_cond.Range("_alh")             '緩和係数
    hflag = ws_cond.Range("_hflag")
    hinflag = ws_cond.Range("_hinflag")
    
    hdown = 0
    etime = thyd(nq)                       '計算終了時間
    qp = qhyd(0)                            '初期流量
    
    dt = dtm
    itout = Int(otime / dt)
    time = 0#
    itt = itout
        
    'フォームの初期計算時刻表示
    DoEvents
    UserForm1.Label2.Visible = True
    UserForm1.Label3.Visible = True
    UserForm1.Label4.Visible = True
    UserForm1.Label2.Caption = 0
    UserForm1.Label4.Caption = etime
    UserForm1.Repaint
    
    Call initl(slope)
    
    '入力データのチェック
    Dim ws_out As Worksheet
    Set ws_out = Sheets("入力チェック")
    
    ws_out.Cells(1, 1) = "nx - 1:"
    ws_out.Cells(1, 2) = nx - 1
    ws_out.Cells(2, 1) = "h0:"
    If hflag = False Then
        ws_out.Cells(2, 2) = h0
    Else
        ws_out.Cells(2, 2) = "データシート"
    End If
    
    ws_out.Cells(3, 1) = "i"
    ws_out.Cells(3, 2) = "x_up"
    ws_out.Cells(3, 3) = "z_up"
    ws_out.Cells(3, 4) = "b_up"
    
    For i = 1 To nx - 1
        ws_out.Cells(i + 3, 1) = i
        ws_out.Cells(i + 3, 2) = x_up(i)
        ws_out.Cells(i + 3, 3) = z_up(i)
        ws_out.Cells(i + 3, 4) = b_up(i)
    Next i
        
    
    icount = 0
    ncount = 0
    mcount = 3
    
    '--------------------------------------------------
    '　計算開始
    '--------------------------------------------------
mstart:
    icount = icount + 1
    Application.StatusBar = "Calculate Step = " & CStr(time) & " / " & CStr(etime)

    '流量データの作成
    For i = 1 To nq
        If (time >= thyd(i - 1)) And (time <= thyd(i)) Then
            sst = (time - thyd(i - 1)) / (thyd(i) - thyd(i - 1))
            qinput = qhyd(i - 1) + (qhyd(i) - qhyd(i - 1)) * sst
            If hinflag = False Then
                hdown = z(nx) + (qinput * mn(nx) / b(nx) / Sqr(slope)) ^ (3# / 5#)
            Else
                hdown = z(nx) + hhyd(i - 1) + (hhyd(i) - hhyd(i - 1)) * sst
            End If
        End If
    Next i
    
    qp = qinput
    If hinflag = False Then
        hin = (qp * mn(1) / b(1) / Sqr(slope)) ^ (3# / 5#)
    Else
        hin = h(1)
    End If


    'フォームの計算時刻表示
    DoEvents
    UserForm1.Label2.Caption = Format(time, "####.######")
    UserForm1.Repaint
    
    If windowflag = True Then GoTo mend:
    
    '
    Call hcal(errmax, err, lmax, lcount, alh, hdown)
    Call bound_h(hn, hsn, hdown)
    Call bound_u(yun)
    Call newgrd_u(yun, yu, gux)
    
    '
    If (itt = itout) Then
        itt = 0
        Call output(ncount)
        Call outputm(mcount)
        Debug.Print time, lcount, err
    End If
    itt = itt + 1
    
    '
    Call dcip1d(yun, gux)
    Call bound_u(yun)
    Call shift(yun, yu)
    
    '
    Call hshift(hn, h, hsn, hs)
    Call qccal(qc, yu, hs)
    
    If (time > etime) Then GoTo mend:
    time = time + dt
    GoTo mstart:
    
mend:
    '--------------------------------------------------
    '　計算終了
    '--------------------------------------------------

    Unload UserForm1
     
    windowflag = False
    
    Close

    Application.StatusBar = False
    
End Sub


'--------------------------------------------
'河道データの読み込みと初期条件の設定
'
'--------------------------------------------
Private Sub initl(dslope As Double)

    Dim i As Long
    Dim hs_up As Double
    Dim hc As Double
    Dim slope As Double
    
    Dim ws_bz As Worksheet
    Set ws_bz = Sheets("河道データ")
    
    nx = ws_bz.Cells(1, "B")
    
    dx(0) = 0
    For i = 1 To nx
        b(i) = ws_bz.Range("_b").Offset(i, 0)       '川幅
        z(i) = ws_bz.Range("_z").Offset(i, 0)       '河床高
        z0(i) = ws_bz.Range("_z0").Offset(i, 0)     '初期河床高
        dx(i) = ws_bz.Range("_dx").Offset(i, 0)     '断面間距離
        mn(i) = ws_bz.Range("_mn").Offset(i, 0)     'マニングの粗度係数
        h(i) = ws_bz.Range("_hinit").Offset(i, 0)   '初期水深
        x(i) = ws_bz.Range("_tl").Offset(i, 0)      '累積距離
        
    Next i
    
    dslope = (z(nx - 1) - z(nx)) / dx(nx)
    
    '++++++++++++++++++++++++++++++++++++++++++++++++++
    If hflag = False Then
        '水深の初期条件を等流水深で決める場合
        'slope = (z(1) - z(2)) / dx(2)
        'h0 = (qp * mn(1) / b(1) / Sqr(slope)) ^ (3# / 5#)
        'hc = ((qp / b(1)) ^ 2 / g) ^ (1# / 3#)
        
        'For i = 1 To nx
        '    h(i) = z0(i) + h0
        '    hs(i) = h(i) - z(i)
        '    If (hs(i) < hmin) Then
        '        hs(i) = hmin
        '        h(i) = hs(i) + z(i)
        '    End If
        '    hn(i) = h(i)
        '    hsn(i) = hs(i)
        'Next i
        
        Call uniformflow(qp)
        
        For i = 1 To nx
            h(i) = z0(i) + h(i)
            hs(i) = h(i) - z(i)
            hn(i) = h(i)
            hsn(i) = hs(i)
        Next i
        '++++++++++++++++++++++++++++++++++++++++++++++++++
    Else
        '++++++++++++++++++++++++++++++++++++++++++++++++++
        '水深の初期条件を読み込む場合
        For i = 1 To nx
            h(i) = z0(i) + h(i)
            hs(i) = h(i) - z(i)
            hn(i) = h(i)
            hsn(i) = hs(i)
        Next i
    End If
    '++++++++++++++++++++++++++++++++++++++++++++++++++
    
    dx(0) = dx(1)
    
    For i = 1 To nx - 1
        x_up(i) = (x(i) + x(i + 1)) * 0.5
        z_up(i) = (z(i) + z(i + 1)) * 0.5
        z0_up(i) = (z(i) + z(i + 1)) * 0.5
        b_up(i) = (b(i) + b(i + 1)) * 0.5
    Next i
    
    For i = 1 To nx - 1
        hs_up = (hs(i) + hs(i + 1)) * 0.5
        yu(i) = qp / (hs_up * b_up(i))
        qc(i) = hs_up * yu(i) * b_up(i)
    Next i
    yu(0) = yu(1)
    yu(nx) = yu(nx - 1)
    
    For i = 0 To nx
        yun(i) = yu(i)
        gux(i) = 0#
        dx2(i) = dx(i) * dx(i)
    Next i

End Sub


'--------------------------------------------
' 水深の計算
'
'--------------------------------------------
Private Sub hcal(errmax As Double, err As Double, lmax As Long, _
                 lcount As Long, alh As Double, hdown As Double)
    
    Dim i As Long
    Dim k As Long
    Dim qu(im) As Double
    Dim div As Double
    Dim hsta As Double
    Dim serr As Double
    
    For k = 1 To lmax
        Call cntfl(yun, hsn, cfx)
        Call rhs
        For i = 1 To nx - 1
            qu(i) = yun(i) * (hsn(i) + hsn(i + 1)) * 0.5 * b_up(i)
        Next i
        qu(1) = qp      '最上流端で一定流量
        qu(0) = qp
        
        err = 0#
        For i = 2 To nx - 1
            div = (qu(i - 1) - qu(i)) / b(i)
            hsta = h(i) + div / dx(i) * dt
            serr = Abs(hsta - hn(i))
            err = err + serr
            hn(i) = hsta * alh + hn(i) * (1# - alh)
            hsn(i) = hn(i) - z(i)
        Next i
        Call bound_h(hn, hsn, hdown)
        If (err < errmax) Then GoTo h200:
    Next k
h200:
    lcount = k
    
End Sub


'--------------------------------------------
' 非移流項の計算
'
'--------------------------------------------
Private Sub cntfl(u() As Double, hs() As Double, cfx() As Double)
    
    Dim i As Long
    Dim hs_up As Double
    Dim dbdx As Double
    
    For i = 1 To nx - 1
        hs_up = (hs(i) + hs(i + 1)) * 0.5
        dbdx = (b(i + 1) - b(i)) / dx(i)
        cfx(i) = -mn(i) ^ 2 * g * u(i) * Abs(u(i)) / hs_up ^ (4# / 3#) - g * hs_up / (2# * b_up(i)) * dbdx
    Next i
    cfx(0) = cfx(1)
    cfx(nx) = cfx(nx - 1)
    
End Sub


'--------------------------------------------
' 移流項の計算
'
'--------------------------------------------
Private Sub rhs()
    
    Dim i As Long
    Dim dhdx As Double
    
    For i = 1 To nx - 1
        dhdx = (hn(i + 1) - hn(i)) / dx(i)
        yun(i) = yu(i) + (cfx(i) - g * dhdx) * dt
    Next i
    Call bound_u(yun)
    
End Sub


'--------------------------------------------
' 水位の境界条件
'
'--------------------------------------------
Private Sub bound_h(h() As Double, hs() As Double, hdown As Double)
    
    hs(1) = hs(2)
    h(1) = z(1) + hs(1)
    hs(nx) = hs(nx - 1)     'nx-1と同じ値を使う
    h(nx) = z(nx) + hs(nx)  '
    'h(nx) = hdown          '等流の計算値を使う
    'hs(nx) = hdown - z(nx) '

End Sub


'--------------------------------------------
' 流速の境界条件
'
'--------------------------------------------
Private Sub bound_u(u() As Double)
    
    'u(0) = qp / (b(1) * hin)
    
End Sub


'--------------------------------------------
' 微分係数の計算
'
'--------------------------------------------
Private Sub newgrd_u(yn() As Double, y() As Double, gx() As Double)
    Dim i As Long
    
    For i = 1 To nx - 1
        'gux(i) = gux(i) + (yun(i + 1) - yun(i - 1) - yu(i + 1) + yu(i - 1)) * 0.5 / dx
        gx(i) = gx(i) + (yn(i + 1) - yn(i - 1) - y(i + 1) + y(i - 1)) * 0.5 / dx(i)
    Next i
    
End Sub


'--------------------------------------------
' CIP法計算のメイン
'
'--------------------------------------------
Private Sub dcip1d(f() As Double, gx() As Double)
    
    Dim u(im) As Double
    Dim fn(im) As Double
    Dim gxn(im) As Double
    
    Dim i As Long
    Dim hs_up As Double
    Dim dx3 As Double
    Dim xx As Double
    Dim isn As Long
    Dim im1 As Long
    Dim a1 As Double
    Dim b1 As Double
    
    For i = 1 To nx - 1
        hs_up = (hs(i) + hs(i + 1)) * 0.5
        u(i) = yu(i)
    Next i
    
    For i = 1 To nx - 1
        dx3 = dx2(i) * dx(i)
        xx = -u(i) * dt
        'isn = sign(1#, u(i))
        isn = Sgn(u(i))
        If isn < 0 Then
            isn = -1
        Else
            isn = 1
        End If
        
        im1 = i - isn
       a1 = ((gx(im1) + gx(i)) * dx(i) * isn - 2# * (f(i) - f(im1))) / (dx3 * isn)
       b1 = (3# * (f(im1) - f(i)) + (gx(im1) + 2# * gx(i)) * dx(i) * isn) / dx2(i)
       fn(i) = ((a1 * xx + b1) * xx + gx(i)) * xx + f(i)
       gxn(i) = (3# * a1 * xx + 2# * b1) * xx + gx(i)
    Next i
    
    For i = 1 To nx - 1
        f(i) = fn(i)
        gx(i) = gxn(i) - (gxn(i) * (u(i + 1) - u(i - 1))) * 0.5 * dt / dx(i)
    Next i
    
End Sub


'--------------------------------------------
' 計算結果の出力その１
'
'--------------------------------------------
Private Sub output(ncount As Long)
    
    Dim i As Long
    Dim h_up(im), hs_up(im)
    
    Dim ws_out As Worksheet
    Set ws_out = Sheets("結果")
    
    For i = 1 To nx - 1
        h_up(i) = (hn(i) + hn(i + 1)) * 0.5
    Next i
    For i = 1 To nx - 1
        hs_up(i) = h_up(i) - z_up(i)
    Next i
    
    ncount = ncount + 1
    ws_out.Cells(ncount, 1) = "time:"
    ws_out.Cells(ncount, 2) = time
    ws_out.Cells(ncount, 3) = "qp:"
    ws_out.Cells(ncount, 4) = qp
    
    ncount = ncount + 1
    ws_out.Cells(ncount, 1) = "i"
    ws_out.Cells(ncount, 2) = "累積距離"
    ws_out.Cells(ncount, 3) = "水深"
    ws_out.Cells(ncount, 4) = "流速"
    ws_out.Cells(ncount, 5) = "流量"
    ws_out.Cells(ncount, 6) = "川幅"
    ws_out.Cells(ncount, 7) = "河床高"
    ws_out.Cells(ncount, 8) = "水位"
    For i = 1 To nx - 1
        ncount = ncount + 1
        ws_out.Cells(ncount, 1) = i
        ws_out.Cells(ncount, 2) = x(i)
        ws_out.Cells(ncount, 3) = hs_up(i)
        ws_out.Cells(ncount, 4) = yun(i)
        ws_out.Cells(ncount, 5) = qc(i)
        ws_out.Cells(ncount, 6) = b_up(i)
        ws_out.Cells(ncount, 7) = z_up(i)
        ws_out.Cells(ncount, 8) = z_up(i) + hs_up(i)
    Next i
    
End Sub


'--------------------------------------------
' 計算結果の出力その２
'
'--------------------------------------------
Public Sub outputm(mcount As Long)
    
    Dim i As Long
    Dim n As Long
    Dim h_up(im), hs_up(im)
    
    Dim ws_h_up As Worksheet
    Set ws_h_up = Sheets("水深結果")
    Dim ws_yun As Worksheet
    Set ws_yun = Sheets("流速結果")
    Dim ws_qc As Worksheet
    Set ws_qc = Sheets("流量結果")
    Dim ws_i As Worksheet
    Set ws_i = Sheets("水位結果")
    
    For i = 1 To nx - 1
        h_up(i) = (hn(i) + hn(i + 1)) * 0.5
    Next i
    For i = 1 To nx - 1
        hs_up(i) = h_up(i) - z_up(i)
    Next i
    
    'mcount = mcount + 1
    ws_h_up.Cells(1, mcount) = time
    ws_yun.Cells(1, mcount) = time
    ws_qc.Cells(1, mcount) = time
    ws_i.Cells(1, mcount) = time

    n = 0
    For i = 1 To nx - 1
        n = n + 1
        
        ws_h_up.Cells(i + 1, 1) = i
        ws_h_up.Cells(i + 1, 2) = x(i)
        ws_h_up.Cells(i + 1, mcount) = hs_up(i)
        
        ws_yun.Cells(i + 1, 1) = i
        ws_yun.Cells(i + 1, 2) = x(i)
        ws_yun.Cells(i + 1, mcount) = yun(i)
        
        ws_qc.Cells(i + 1, 1) = i
        ws_qc.Cells(i + 1, 2) = x(i)
        ws_qc.Cells(i + 1, mcount) = qc(i)
    
        ws_i.Cells(i + 1, 1) = i
        ws_i.Cells(i + 1, 2) = x(i)
        ws_i.Cells(i + 1, mcount) = hs_up(i) + z0(i)
    Next i
    mcount = mcount + 1
    
End Sub


'--------------------------------------------
' 流速の新しいステップへの入れ替え
'
'--------------------------------------------
Private Sub shift(yn() As Double, y() As Double)
    
    Dim i As Long
    
    For i = 1 To nx - 1
        y(i) = yn(i)
    Next i
    
End Sub


'--------------------------------------------
' 水深の新しいステップへの入れ替え
'
'--------------------------------------------
Private Sub hshift(hn() As Double, h() As Double, hsn() As Double, hs() As Double)
    
    Dim i As Long
    
    For i = 1 To nx - 1
        h(i) = hn(i)
        hs(i) = h(i) - z(i)
        If (hs(i) < hmin) Then
            hs(i) = hmin
            h(i) = z(i) + hmin
        End If
        hsn(i) = hs(i)
    Next i
    
End Sub


'--------------------------------------------
' 流量の計算
'
'--------------------------------------------
Private Sub qccal(qc() As Double, u() As Double, hs() As Double)
    
    Dim i As Long
    
    For i = 1 To nx - 1
        If (i = 1) Then
            qc(i) = qp
        Else
            qc(i) = u(i) * (hs(i) + hs(i + 1)) * 0.5 * b_up(i)
        End If
    Next i
    
End Sub


'************************************************
' 等流計算
' 初期水位として使用する
'
'************************************************
Private Sub uniformflow(qq As Double)
    Dim i As Long
    Dim zs0x As Double
    Dim ib As Double
    
    zs0x = (z0(nx) - z0(1)) / x(nx)     '計算対象河道の平均河床勾配
    
    '等流水深
    For i = 1 To nx - 1
        ib = (z0(i + 1) - z0(i)) / dx(i)
        If ib < 0 Then
             h(i) = (mn(i) * qq / b(i) / (Sqr(Abs(ib)))) ^ (3# / 5#)            '河床勾配が正の場合
        Else
             h(i) = (mn(i) * qq / b(i) / (Sqr(Abs(zs0x)))) ^ (3# / 5#)          '河床勾配が負の場合
        End If
    Next i

End Sub