====== cip1d_uns_bed_v_1_0_0 ======
cip1d_uns_bed_v_1_0_0.lzhをダウンロードして解凍すると、「**cip1d_uns_bed_v_1_0_0.xls**」ができます。
この「**cip1d_uns_bed_v_1_0_0.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次元モデル:不定流モデル:image31.jpg?650|計算条件シート}}
;#;
計算条件シート
;#;
{{ :河床変動計算:1次元モデル:不定流モデル:image32.jpg?325 |計算フロー}}
;#;
計算フロー
;#;
=== 流量データ ===
{{ :河床変動計算:1次元モデル:不定流モデル:image33.jpg?650|流量データシート}}
;#;
流量データシート
;#;
=== 河道データ ===
{{ :河床変動計算:1次元モデル:不定流モデル:image34.jpg?650|河道データシート}}
;#;
河道データシート
;#;
==== 計算結果 ====
計算結果は、以下の14のシートに書き出されます。 グラフについては、使用者各自に作成して下さい。
* 計算結果出力シート : 出力間隔毎の計算結果
* 計算河床変動量シート
* 計算河床高シート
* 計算水深シート
* 計算水位シート
* 計算流積シート
* 計算流量シート
* 計算流速シート
* 計算フルード数シート
* 計算平均粒径シート
* 計算流砂量シート
* 計算累積流砂量シート
* 計算粒径別累積流砂量シート : 出力間隔毎の計算結果
* 計算粒径区分割合シート : 出力間隔毎の計算結果
{{ :河床変動計算:1次元モデル:不定流モデル:image35.jpg?650|入力チェックシート}}
;#;
入力チェックシート
;#;
{{ :河床変動計算:1次元モデル:不定流モデル:image36.jpg?650|結果シート}}
;#;
結果シート
;#;
{{ :河床変動計算:1次元モデル:不定流モデル:image37.jpg?650|水深結果シート}}
;#;
水深結果シート
;#;
{{ :河床変動計算:1次元モデル:不定流モデル:image37.jpg?650|水位結果シート}}
;#;
水位結果シート
;#;
===== VBAソースファイル =====
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