====== 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**」の使用方法を説明します。
ダウンロードについては、[[http://www.disaster-software.net|水理水文ソフトウェア]]の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法の基本式と差分化について「計算条件」シートに書いていますので、参考資料とともに参照して下さい。
{{ :河川不定流計算:1次元不定流:image10.jpg?650|計算条件シート}}
;#;
計算条件シート
;#;
* 計算結果の出力時間間隔
* 計算終了時間
* 計算刻み時間
* 最小水深
* 繰り返し計算の打ち切り誤差
* 繰り返し計算回数
* 繰り返し計算の緩和係数
* 初期水深をデータで与える : True or False
* 最下流端の水深をデータで与える : True or False
\\
{{ :河川不定流計算:1次元不定流:image11.jpg?300 |計算フローチャート}}
;#;
計算フローチャート
;#;
再度計算実行した場合は、計算結果シート(複数)の値がすべて上書きされます。
=== 河道データ ===
{{ :河川不定流計算:1次元不定流:image13.jpg?650|河道シート}}
;#;
河道シート
;#;
\\
* 距離
* 断面間距離
* 河床高
* 粗度係数
* 川幅
* 初期水位
* 初期水深 : 初期水深をデータで与える場合に必要となる。
=== 流量データ ===
{{ :河川不定流計算:1次元不定流:image12.jpg?650|流量シート}}
;#;
流量シート
;#;
\\
* データ数
* 経過時間 : 必ず0から始まる。
* 流量
==== 計算結果 ====
計算結果は、以下の6つのシートに書き出されます。 グラフについては、使用者各自に作成して下さい。
* 入力チェックシート
* 結果シート : 出力間隔毎の計算結果
* 水深結果シート
* 水位結果シート
* 流速結果シート
* 流量結果シート
{{ :河川不定流計算:1次元不定流:image14.jpg?650|入力チェックシート}}
;#;
入力チェックシート
;#;
\\
{{ :河川不定流計算:1次元不定流:image15.jpg?650|結果シート}}
;#;
結果シート
;#;
\\
{{ :河川不定流計算:1次元不定流:image17.jpg?650|水深結果シート}}
;#;
水深結果シート
;#;
\\
{{ :河川不定流計算:1次元不定流:image18.jpg?650|水位結果シート}}
;#;
水位結果シート
;#;
\\
{{ :河川不定流計算:1次元不定流:image19.jpg?650|流速結果シート}}
;#;
流速結果シート
;#;
\\
{{ :河川不定流計算:1次元不定流:image20.jpg?650|流量結果シート}}
;#;
流量結果シート
;#;
\\
===== VBAソースファイル =====
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
'--------------------------------------------
' 計算結果の出力その1
'
'--------------------------------------------
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
'--------------------------------------------
' 計算結果の出力その2
'
'--------------------------------------------
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