Hydro Disaster Software Index
Wiki管理
合同会社TYS
技術開発部
TEL 042-812-5920
水理水文ソフトウェア
GIS関連のソフトウェア
Fortranで水理計算
Excel VBAで水理計算
TYS-RBAモデル
合同会社TYS
— mail to TYS Support
累計: 本日: 昨日:
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
計算結果は、以下の14のシートに書き出されます。 グラフについては、使用者各自に作成して下さい。
入力チェックシート
結果シート
水深結果シート
水位結果シート
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