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.lzhをダウンロードして解凍すると、「mac1d_uns_v_1_0_0.xls」ができます。 この「mac1d_uns_v_1_0_0.xls」の使用方法を説明します。 ダウンロードについては、水理水文ソフトウェアのHPを参照してください。
本マクロは、各条件シートに条件を入力するだけで、Maccormack法による1次元不定流計算が行えます。
「計算条件」シートに計算パラメータを入力します。また、Maccormack法の基本式と差分化について「計算条件」シートに書いていますので、参考資料とともに参照して下さい。
計算条件シート
計算結果は、以下の6つのシートに書き出されます。 グラフについては、使用者各自に作成して下さい。
結果シート
最終結果シート
水深結果シート
水位結果シート
流速結果シート
流量結果シート
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) '定常計算を行い各段面の初期値を作成する '流量データのはじめの1ステップを定常計算とする。 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 '-------------------------------------------- ' 計算結果の出力その1 ' '-------------------------------------------- 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 '-------------------------------------------- ' 計算結果の出力その2 ' '-------------------------------------------- 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