Hydro Disaster Software Index
Wiki管理
合同会社TYS
技術開発部
TEL 042-812-5920
水理水文ソフトウェア
GIS関連のソフトウェア
Fortranで水理計算
Excel VBAで水理計算
TYS-RBAモデル
合同会社TYS
— mail to TYS Support
累計：
本日：
昨日：
mac1d_uns_v_1_0_0
mac1d_uns_v_1_0_0.lzhをダウンロードして解凍すると、「mac1d_uns_v_1_0_0.xls」ができます。 この「mac1d_uns_v_1_0_0.xls」の使用方法を説明します。 ダウンロードについては、水理水文ソフトウェアのHPを参照してください。
マクロの説明
本マクロは、各条件シートに条件を入力するだけで、Maccormack法による1次元不定流計算が行えます。
計算条件とデータの入力
計算条件
「計算条件」シートに計算パラメータを入力します。また、Maccormack法の基本式と差分化について「計算条件」シートに書いていますので、参考資料とともに参照して下さい。
計算条件シート
- Δt : 計算刻み時間
- 計算時間
- 出力間隔
- EPS(乱流粘性係数)
- 人工粘性係数 : 1.0から5.0の値経験的に決定する
- 初期水深をデータで与えるか否か : True or Flase
- 断面数
再度計算実行した場合は、計算結果シート（複数）の値がすべて上書きされます。
河道データ
河道シート
- 距離
- 断面間距離
- 河床高
- 粗度係数
- 川幅
- 初期水位
- 初期水深 : 初期水深をデータで与える場合に必要となる。
流量データ
流量シート
- データ数
- 経過時間 : 必ず0から始まる。
- 流量
計算結果
計算結果は、以下の６つのシートに書き出されます。 グラフについては、使用者各自に作成して下さい。
- 結果シート : 出力間隔毎の計算結果
- 最終結果シート
- 水深結果シート
- 水位結果シート
- 流速結果シート
- 流量結果シート
結果シート
最終結果シート
水深結果シート
水位結果シート
流速結果シート
流量結果シート
VBAソースファイル
Maccromack法に1次元不定流計算の計算部分のソースコードを公開します。このマクロについては、ソースファイルはすべて公開しています。
Option Explicit
Option Base 1
'計算キャンセルのフラグ
Public windowflag As Boolean
'最大配列数の宣言
Public Const ih As Long = 100 '流量データ最大数
'ベクトル式
Type vector2d
x As Double
y As Double
End Type
'変数の宣言
Private ied As Long '外部境界を含む断面分割数 +2
Private ndata As Long '外部境界を含む断面分割数 +2
Private dt As Double '計算刻み時間
Private dtdx As Double
Private q As Double '流量
Private eps As Double '乱流粘性係数
Private kv As Double '人工粘性係数
Private istep As Long '繰り返し計算回数
Private otime As Double '出力時間間隔
Private nq As Long
Private qinput As Double
Public etime As Double
Public time As Double
'配列の宣言
Private z0() As Double '河床高
Private hinit() As Double '水深
Private level() As Double '水位
Private s0x() As Double '河床勾配ib
Private totalx() As Double '累積距離
Private manning() As Double 'マニングの粗度係数
Private b0() As Double '川幅
Private dx() As Double '断面間距離
'流量関連の配列宣言
Public thyd(ih) As Double '流量データの時間
Public qhyd(ih) As Double '流量データ
Public hhyd(ih) As Double '最下流端の水深
'ベクトルの定義
Private U() As vector2d
Private preU() As vector2d
Private E() As vector2d
Private preE() As vector2d
Private C() As vector2d
Private preC() As vector2d
Private V() As vector2d
Private preV() As vector2d
Private Const ist As Long = 2 '計算開始メッシュ
Private Const g As Double = 9.8 '重力加速度
'************************************************
' MacCormak法による1次元非定常流の計算
'
' 人工粘性項としては、拡散法を適用
'
' 合同会社TYS
' 技術開発部
' 11/08/2011
'
' 参考図書 エクセル河川数値解析入門 山海堂
' 河村哲也 著 ISBN4-381-01696-3
'
'************************************************
Public Sub macRiver1d()
Dim i As Long
Dim j As Long
Dim k As Long
Dim n As Long
Dim itout As Long
Dim itt As Long
Dim ncount As Long
Dim mcount As Long
Dim sst As Double
'初期水深の与えた方のフラグ
Dim hflag As Boolean 'ture:シートのデータ false:等流計算
'計算条件の読み込み
Dim ws_jk As Worksheet
Set ws_jk = Sheets("計算条件")
dt = ws_jk.Range("_dt")
etime = ws_jk.Range("_istep")
otime = ws_jk.Range("_otime")
'q = ws_jk.Range("_q")
eps = ws_jk.Range("_eps")
kv = ws_jk.Range("_kv")
hflag = ws_jk.Range("_hflag")
itout = Int(otime / dt)
itt = itout
ied = ist + ws_jk.Range("_xdata")
ndata = ied + 1
'配列の再宣言
ReDim z0(ndata) As Double
ReDim hinit(ndata) As Double
ReDim level(ndata) As Double
ReDim s0x(ndata) As Double
ReDim totalx(ndata) As Double
ReDim manning(ndata) As Double
ReDim b0(ndata) As Double
ReDim dx(ndata) As Double
ReDim U(ndata) As vector2d 'U(x,y)
ReDim preU(ndata) As vector2d
ReDim E(ndata) As vector2d 'E(x,y)
ReDim preE(ndata) As vector2d
ReDim C(ndata) As vector2d 'C(x,y)
ReDim preC(ndata) As vector2d
ReDim V(ndata) As vector2d 'V(x,y)
ReDim preV(ndata) As vector2d
'配列の初期化
For i = 1 To ndata
z0(i) = 0#
hinit(i) = 0#
level(i) = 0#
s0x(i) = 0#
totalx(i) = 0#
manning(i) = 0#
b0(i) = 0#
dx(i) = 0#
Next i
'計算結果シート内のデータ削除
Call SheetDataDel
'計算データの読み込み
Call readData
'流量データの読み込み
Dim ws_thq As Worksheet
Set ws_thq = Sheets("流量データ")
nq = ws_thq.Range("_nq") '流量データ数
For i = 1 To nq
thyd(i) = ws_thq.Range("_time").Offset(i, 0) '時間
qhyd(i) = ws_thq.Range("_q").Offset(i, 0) '流量
hhyd(i) = ws_thq.Range("_h").Offset(i, 0) '水深
'Debug.Print i, thyd(i), qhyd(i)
Next i
'河床勾配の計算
Call setS0x
'++++++++++++++++++++++++++++++++++++++++++++++++++++++
'初期水深を等流あるいは不等流計算で求める場合
'サブルーチンを用意しいますが、チェックしていませんので
'利用者各自でテストして下さい。
If hflag = False Then
Call Uniformflow
'Call non_Uniformflow
End If
'++++++++++++++++++++++++++++++++++++++++++++++++++++++
'初期条件の設定
Call setUInit
Call setEInit
Call setCInit
time = 0
ncount = 0
mcount = 3
'初期状態の出力
'Call output(ncount)
'Call outputm(mcount)
'定常計算を行い各段面の初期値を作成する
'流量データのはじめの１ステップを定常計算とする。
For n = 1 To 50000
Call setDT
Call predict
Call correct
'最上流端で同じqを与える
U(ist).y = qhyd(1)
preU(ist).y = qhyd(1)
U(1).y = qhyd(1)
preU(1).y = qhyd(1)
'Debug.Print "n="; n
Next n
'--------------------------------------------------
' 非定常計算開始
'--------------------------------------------------
mstart:
Application.StatusBar = "Calculate Step = " & CStr(time) & " / " & CStr(etime)
'流量データの作成
For i = 2 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
'ステップで流量を増減させる場合
qinput = qhyd(i - 1)
End If
Next i
If windowflag = True Then GoTo mend:
Call setDT
Call predict
Call correct
'計算途中の出力
If (itt = itout) Then
itt = 0
Call output(ncount)
Call outputm(mcount)
Debug.Print time
End If
itt = itt + 1
'定常計算の場合、最上流端で同じqを与える
U(ist).y = qinput
preU(ist).y = qinput
U(1).y = qinput
preU(1).y = qinput
If (time > etime) Then GoTo mend:
time = time + dt
GoTo mstart:
mend:
'--------------------------------------------------
' 計算終了
'--------------------------------------------------
'計算結果出力（最終の計算結果）
Dim ws_o As Worksheet
Set ws_o = Sheets("最終結果")
For i = 1 To ied
ws_o.Range("_X").Offset(i) = totalx(i)
ws_o.Range("_I").Offset(i) = i
ws_o.Range("_OB").Offset(i) = z0(i)
ws_o.Range("_OSL").Offset(i) = hinit(i) + z0(i)
ws_o.Range("_OL").Offset(i) = U(i).x / b0(i) + z0(i)
ws_o.Range("_suishin").Offset(i) = U(i).x / b0(i)
Next
Close
Application.StatusBar = False
Set ws_o = Nothing
Set ws_jk = Nothing
End Sub
Private Sub set_C()
Dim i As Long
Dim h As Double
Dim uu As Double
Dim ie As Double '摩擦勾配
Dim wk As Double
Dim qq As Double
For i = ist To ied - 1
C(i).x = 0
h = U(i).x / b0(i)
qq = U(i).y
uu = qq / (b0(i) * h)
If h > 0 Then
ie = manning(i) ^ 2 * uu * Abs(uu) / h ^ (4# / 3#)
Else
ie = 0
End If
'If h > 0 Then
' ie = manning(i) ^ 2 * qq ^ 2 / h ^ (10# / 3#)
'Else
' ie = 0
'End If
wk = eps * (U(i - 1).y - 2 * U(i).y + U(i + 1).y) / dx(i) ^ 2
C(i).y = g * b0(i) * h * (s0x(i) - ie) + wk
Next i
C(ist - 1) = C(ist)
C(ied) = C(ied - 1)
End Sub
'************************************************
' 予測子段階で使う人工粘性係数の計算
'
'
'************************************************
Private Sub set_V()
Dim i As Long
Dim kx As Double
Dim h As Double
Dim uu As Double
Dim ustar As Double
Dim ief As Double
Dim hdx As Double
Dim hux As Double
For i = ist To ied - 1
h = U(i).x / b0(i)
uu = U(i).y / (h * b0(i))
'摩擦速度の計算
'If (h > 0) Then
' ief = ((h * manning(i)) ^ 2#) / (h ^ (4# / 3#)) 'エネルギー勾配
' ustar = Sqr(g * h * ief)
'Else
' ustar = 0#
'End If
If (h > 0) Then
ustar = Sqr(g) * manning(i) * uu / (h ^ (1# / 6#))
Else
ustar = 0#
End If
hdx = U(i - 1).x / b0(i - 1)
hux = U(i + 1).x / b0(i + 1)
kx = kv * ustar * h / dx(i)
V(i).x = kx * (hdx - 2 * h + hux)
V(i).y = kx * (U(i - 1).y - 2 * U(i).y + U(i + 1).y)
Next i
V(ist - 1) = V(ist)
V(ied) = V(ied - 1)
End Sub
Private Sub set_E()
Dim i As Long
Dim h As Double
Dim uu As Double
For i = ist To ied - 1
h = U(i).x / b0(i)
uu = U(i).y / (h * b0(i))
E(i).x = U(i).y
E(i).y = 0.5 * g * b0(i) * h ^ 2 + uu ^ 2 * (b0(i) * h)
Next i
E(ist - 1) = E(ist)
E(ied) = E(ied - 1)
End Sub
Private Sub set_preC()
Dim i As Long
Dim h As Double
Dim uu As Double
Dim qq As Double
Dim ie As Double 'エネルギー勾配
Dim wk As Double
For i = ist To ied - 1
preC(i).x = 0
h = preU(i).x / b0(i)
qq = preU(i).y
uu = qq / (h * b0(i))
If h > 0 Then
ie = manning(i) ^ 2 * uu * Abs(uu) / h ^ (4# / 3#)
Else
ie = 0
End If
'If h > 0 Then
' ie = manning(i) ^ 2 * qq ^ 2 / h ^ (10# / 3#)
'Else
' ie = 0
'End If
wk = eps * (preU(i - 1).y - 2 * preU(i).y + preU(i + 1).y) / dx(i) ^ 2
preC(i).y = g * b0(i) * h * (s0x(i) - ie) + wk
Next i
preC(ist - 1) = preC(ist)
preC(ied) = preC(ied - 1)
End Sub
'************************************************
' 修正子段階で使う人工粘性係数の計算
'
'
'************************************************
Private Sub set_preV()
Dim i As Long
Dim kx As Double
Dim h As Double
Dim uu As Double
Dim ustar As Double
Dim ief As Double
Dim hdx As Double
Dim hux As Double
For i = ist To ied - 1
h = preU(i).x / b0(i)
uu = preU(i).y / (h * b0(i))
'摩擦速度の計算
'If (h > 0) Then
' ief = ((h * manning(i)) ^ 2#) / (h ^ (4# / 3#)) 'エネルギー勾配
' ustar = Sqr(g * h * ief)
'Else
' ustar = 0#
'End If
If (h > 0) Then
ustar = Sqr(g) * manning(i) * uu / (h ^ (1# / 6#))
Else
ustar = 0#
End If
hdx = preU(i - 1).x / b0(i - 1)
hux = preU(i + 1).x / b0(i + 1)
kx = kv * ustar * h / dx(i)
preV(i).x = kx * (hdx - 2 * h + hux)
preV(i).y = kx * (preU(i - 1).y - 2 * preU(i).y + preU(i + 1).y)
Next i
preV(ist - 1) = preV(ist)
preV(ied) = preV(ied - 1)
End Sub
Private Sub set_preE()
Dim i As Long
Dim h As Double, uu As Double
For i = ist To ied - 1
h = preU(i).x / b0(i)
uu = preU(i).y / (h * b0(i))
preE(i).x = preU(i).y
preE(i).y = 0.5 * g * b0(i) * h ^ 2 + uu ^ 2 * (b0(i) * h)
Next i
preE(ist - 1) = preE(ist)
preE(ied) = preE(ied - 1)
End Sub
'************************************************
' 修正子段階の計算
'
'
'************************************************
Private Sub correct()
Dim i As Long
Call set_preE
Call set_preV
Call set_preC
For i = ist To ied - 1
dtdx = dt / dx(i)
U(i).x = (U(i).x + preU(i).x) _
- dtdx * (preE(i + 1).x - preE(i).x + preV(i + 1).x - preV(i).x) _
+ dt * preC(i).x
U(i).y = (U(i).y + preU(i).y) _
- dtdx * (preE(i + 1).y - preE(i).y + preV(i + 1).y - preV(i).y) _
+ dt * preC(i).y
U(i).x = U(i).x / 2#
U(i).y = U(i).y / 2#
Next i
U(ist - 1) = U(ist)
U(ied) = U(ied - 1)
End Sub
'************************************************
' 予測子段階の計算
'
'
'************************************************
Private Sub predict()
Dim i As Long
Call set_E
Call set_V
Call set_C
For i = ist To ied - 1
dtdx = dt / dx(i)
preU(i).x = U(i).x _
- dtdx * (E(i).x - E(i - 1).x - V(i).x + V(i - 1).x) _
+ dt * C(i).x
preU(i).y = U(i).y _
- dtdx * (E(i).y - E(i - 1).y - V(i).y + V(i - 1).y) _
+ dt * C(i).y
Next i
preU(ist - 1) = preU(ist)
preU(ied) = preU(ied - 1)
End Sub
'************************************************
' CFL条件による計算刻み時間の設定
'
'
'************************************************
Private Sub setDT()
Dim i As Long
Dim uu As Double
Dim h As Double
Dim da As Double
Dim ddx As Double
For i = ist To ied - 1
h = U(i).x / b0(i)
uu = U(i).y / (h * b0(i))
If (h > 0#) Then
da = dx(i) / (uu + Sqr(g * h) + 2 * eps / dx(i))
If dt > da Then
dt = da
End If
End If
Next i
End Sub
Private Sub setCInit()
Dim i As Long
Dim h As Double
Dim uu As Double
Dim qq As Double
Dim ie As Double 'エネルギー勾配
For i = 1 To ied
C(i).x = 0
h = hinit(i)
qq = qhyd(1)
uu = qhyd(1) / (h * b0(i))
If h > 0 Then
ie = manning(i) ^ 2 * uu * Abs(uu) / h ^ (4# / 3#)
Else
ie = 0
End If
'If h > 0 Then
' ie = manning(i) ^ 2 * qq ^ 2 / h ^ (10# / 3#)
'Else
' ie = 0
'End If
C(i).y = g * b0(i) * h * (s0x(i) - ie)
preC(i) = C(i)
Next i
End Sub
Private Sub setEInit()
Dim i As Long
Dim h As Double
Dim uu As Double
For i = 1 To ied
h = hinit(i)
uu = qhyd(1) / (h * b0(i))
E(i).x = qhyd(1)
E(i).y = (g * b0(i) * h ^ 2) / 2# + uu ^ 2 * (b0(i) * h)
preE(i) = E(i)
Next i
End Sub
Private Sub setUInit()
Dim i As Long
For i = 1 To ied
U(i).x = b0(i) * hinit(i)
U(i).y = qhyd(1)
preU(i) = U(i)
Next i
End Sub
Private Sub setS0x()
Dim i As Long
For i = ist To ied - 1
s0x(i) = (z0(i - 1) - z0(i + 1)) / (dx(i - 1) + dx(i + 1))
Next i
s0x(ist - 1) = s0x(ist)
s0x(ied) = s0x(ied - 1)
End Sub
Private Sub readData()
Dim i As Long
Dim tl As Double
Dim ws_data As Worksheet
Set ws_data = Sheets("河道データ")
For i = 1 To ied
totalx(i) = ws_data.Cells(i + 2, "B") '累積距離
dx(i) = ws_data.Cells(i + 2, "C") '断面区間距離
z0(i) = ws_data.Cells(i + 2, "D") '河床高
manning(i) = ws_data.Cells(i + 2, "E") 'マニングの粗度係数
b0(i) = ws_data.Cells(i + 2, "F") '川幅
level(i) = ws_data.Cells(i + 2, "G") '水位
hinit(i) = level(i) - z0(i) '水深
Next i
End Sub
'************************************************
' 等流計算
' 初期水位として使用する
'
'************************************************
Private Sub Uniformflow()
Dim i As Long
Dim qq As Double
Dim zs0x As Double
qq = qhyd(1)
zs0x = (z0(1) - z0(ied)) / totalx(ied)
'等流水深
For i = 1 To ied
If s0x(i) > 0 Then
hinit(i) = (manning(i) * qq / b0(i) / (Sqr(s0x(i)))) ^ (3# / 5#) '河床勾配が正の場合
Else
hinit(i) = (manning(i) * qq / b0(i) / (Sqr(zs0x))) ^ (3# / 5#) '河床勾配が負の場合
End If
level(i) = z0(i) + hinit(i)
Next i
End Sub
'************************************************
' 不等流計算(ルンゲ・クッタ法)
' 初期水位として使用する
'
'************************************************
Private Sub non_Uniformflow()
Dim i As Long
Dim qq As Double
Dim aa As Double
Dim bb As Double
Dim cc As Double
Dim fh As Double
Dim fdh As Double
Dim dh As Double
Dim h00 As Double
Dim zang As Double
Dim hs(1000) As Double
qq = q
'不等流水深
zang = Abs(z0(2) - z0(1)) / dx(1)
hinit(1) = (manning(1) * qq / b0(1) / (Sqr(zang))) ^ (3# / 5#) 'h00
hs(1) = hinit(1)
For i = 2 To ied
cc = -z0(i) + z0(i - 1) - hs(i - 1) - qq ^ 2 / (2# * g * b0(i) ^ 2 * hs(i - 1) ^ 2) _
- dx(i) * manning(i) ^ 2 * qq ^ 2 / (2# * b0(i) ^ 2 * hs(i - 1) ^ (10# / 3#))
hs(i) = hs(i - 1)
g120: aa = qq ^ 2 / (2# * g * b0(i) ^ 2)
bb = -dx(i) * manning(i) ^ 2 * qq ^ 2 / (2# * b0(i) ^ 2)
fh = hs(i) + aa / hs(i) ^ 2 + bb / hs(i) ^ (10# / 3#) + cc
fdh = 1# - 2# * aa / hs(i) ^ 3 - 10# / 3# * bb / hs(i) ^ (13# / 3#)
dh = -fh / fdh
If Abs(dh) < 0.00001 Then GoTo g135:
hs(i) = hs(i) + dh
If hs(i) <= 0 Then
hs(i) = hs(i) - dh / 2#
End If
GoTo g120:
g135: hinit(i) = hs(i)
'Debug.Print i, hinit(i)
level(i) = z0(i) + hinit(i)
Next i
End Sub
'--------------------------------------------
' 計算結果の出力その１
'
'--------------------------------------------
Private Sub output(ncount As Long)
Dim i As Long
Dim ws_out As Worksheet
Set ws_out = Sheets("結果")
ncount = ncount + 1
ws_out.Cells(ncount, 1) = "time:"
ws_out.Cells(ncount, 2) = time
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 ied
ncount = ncount + 1
ws_out.Cells(ncount, 1) = i
ws_out.Cells(ncount, 2) = totalx(i)
ws_out.Cells(ncount, 3) = U(i).x / b0(i)
ws_out.Cells(ncount, 4) = U(i).y / U(i).x
ws_out.Cells(ncount, 5) = U(i).y
ws_out.Cells(ncount, 6) = b0(i)
ws_out.Cells(ncount, 7) = z0(i)
ws_out.Cells(ncount, 8) = z0(i) + U(i).x / b0(i)
Next i
End Sub
'--------------------------------------------
' 計算結果の出力その２
'
'--------------------------------------------
Public Sub outputm(mcount As Long)
Dim i As Long
Dim n As Long
Dim ws_h As Worksheet
Set ws_h = Sheets("水深結果")
Dim ws_i As Worksheet
Set ws_i = Sheets("水位結果")
Dim ws_u As Worksheet
Set ws_u = Sheets("流速結果")
Dim ws_q As Worksheet
Set ws_q = Sheets("流量結果")
'mcount = mcount + 1
ws_h.Cells(1, mcount) = Int(time + dt) & "step"
ws_i.Cells(1, mcount) = Int(time + dt) & "step"
ws_u.Cells(1, mcount) = Int(time + dt) & "step"
ws_q.Cells(1, mcount) = Int(time + dt) & "step" 'time & "step"
n = 0
For i = 1 To ied
n = n + 1
ws_h.Cells(i + 1, 1) = Int(i) & "断面"
ws_h.Cells(i + 1, 2) = totalx(i)
ws_h.Cells(i + 1, mcount) = U(i).x / b0(i)
ws_i.Cells(i + 1, 1) = Int(i) & "断面"
ws_i.Cells(i + 1, 2) = totalx(i)
ws_i.Cells(i + 1, mcount) = U(i).x / b0(i) + z0(i)
ws_u.Cells(i + 1, 1) = Int(i) & "断面"
ws_u.Cells(i + 1, 2) = totalx(i)
ws_u.Cells(i + 1, mcount) = U(i).y / U(i).x
ws_q.Cells(i + 1, 1) = Int(i) & "断面" 'i & "断面"
ws_q.Cells(i + 1, 2) = totalx(i)
ws_q.Cells(i + 1, mcount) = U(i).y
Next i
mcount = mcount + 1
End Sub
特に明示されていない限り、本Wikiの内容は次のライセンスに従います： CC Attribution-Share Alike 3.0 Unported
