円形水槽に作用する地震時動液圧

随時編集予定です。

2012年1月21日付で記事にした「技術計算関数-ベッセル関数・Γ（ガンマ）関数」。
足かけ3年越しで再開し終わりまでもう少しのところになったので記録に残します。
この作業については岡山理科大学　西村先生に助言を頂き，又友人のＦさんとＶＢＡ書き換え作業を共にしたこと。
ソースコードを公開されたＷＡＮｔａｒｏ様三者の方々に深く感謝いたします。

　前回の記事の一部。

　　　【地震時,水の入った水槽の側面・底面に作用する動水圧を求める
　　　　速度ポテンシャル法による式に用います。
　　　　まだ，奥が深くてこの式による答えが出ていません。
　　　　ネット検索でＶＢによるソースプログラムがありました。
　　　　式とＶＢの抜粋です。
　　　　詳しくは ”ホームページ　ＷＡＮｔａｒｏ雑記帳”をご覧下さい。】

1：ＶＢソースコードからＶＢＡコードに書き換えることで分かったこと。
　　　　○ｚ：水槽底盤を0とし上向きを正とする座標→この引数・戻り値ｚはｐ（θ，ｚ）を求める式のみ使用。
　　　　　ベッセル関数Ｉｖ（ｚ）式に使用するｚ：λi・ｒ/H
　　　　　Γ関数　　　　Γ（ｚ）式　同上　　ｚ：ｍ+1，ν+ｍ+1
　　　
　　　　○Math.を削除

　　　　○引数FIはなぜかVBAでは使用できない。

　　　　○エクセル関数をそのままVBAで使用できない関数は定義しておく必要がある。
　　　　　代わりにWorksheetFunctionプロパティを使用も可（使用できない関数もある）。
　　　　　例：π（PI())　ｅ（ネイピア数）　Power　Sqrt　sinh

　　　　○Returnに続いて一行で構文が書けるのはVBnet，VBAではエラーとなる。　　
　　　　　　
　　　　○Pow:POWER関数　　　CDbl:CDBL関数　　FIはfloor関数でなく上記式Fi値


2：速度ポテンシャル法に基づく簡略式等
　　　　　

　　　
3：ソースコード
　　　　　ここまでの変数宣言等は省略
　　　　　　NN = 30
　　　　　　dz = 0.5
　　　　　　 z = -dz
　　　　　　Do
　　　　　　　z = z + dz 'サージタンク底面からの距離
　　　　　　　　sum = 0.0
　　　　　　　　　For i = 0 To NN
　　　　　　　　　　lambda = (2.0 * CDbl(i) + 1.0) / 2.0 * Math.PI
　　　　　　　　　　　sum = sum + Math.Pow(-1.0, CDbl(i)) / lambda * FI(r / H, lambda) * Math.Cos(lambda * z / H)
　　　　　　　　　Next i
　　　　　　　 '地震時動水圧強度(速度ポテンシャル)
　　　　　　　　pp = gamma * kh * r * Math.Cos(theta) * sum


　　　　　　　Private Function FI(ByVal rbh As Double, ByVal lambda As Double) As Double
　　　　　　　　Dim frbh As Double
　　　　　　　　frbh = 2.0 / rbh * FMBF(1, lambda * rbh) / (lambda * FMBF(0, lambda * rbh) - FMBF(1, lambda * rbh) / rbh)
　　　　　　　　Return frbh
　　　　　　　　End Function

　　　　　　　Private Function FMBF(ByVal nu As Integer, ByVal z As Double) As Double
　　　　　　　'第1種変形Bessel関数
　　　　　　　　Dim sum As Double
　　　　　　　　Dim i As Integer
　　　　　　　　sum = 0.0
　　　　　　　　　For i = 0 To NN
　　　　　　　　　　sum = sum + 1.0 / GAMMAF(CDbl(i + 1)) / GAMMAF(CDbl(nu + i + 1)) * Math.Pow(0.5 * z, CDbl(2 * i))
　　　　　　　　　Next i
　　　　　　　　Return Math.Pow(0.5 * z, CDbl(nu)) * sum
　　　　　　　　End Function

　　　　　　　　Private Function GAMMAF(ByVal z As Double) As Double
　　　 'ガンマ関数
　　　　　　　　Dim a As Double
　　　　　　　　　a = z / Math.E * Math.Sqrt(z * Math.Sinh(1.0 / z) + 1.0 / 810.0 / Math.Pow(z, 6.0))
　　　　　　　　Return Math.Sqrt(2.0 * Math.PI / z) * Math.Pow(a, z)
　　　　　　　　End Function

4：書き換えVBAコード（Fi値,Iν（ベッセル関数），Γ（ガンマ関数部分）
　　　　　　　

5：同計算結果（水深3ｍ～水面まで記載省略）
　　　　　　　

6：結果の検証（水道用プレストレスコンクリートタンク設計施工指針・解説による式に基づいてＥＸＣＥＬ計算表を作成しました。）

　　　　　　　ループ回数の違いによる誤差があるものの小数第2位までは一致しています。
　　　　　　　グラフはエクセルグラフ機能使用後ＪＰＧ変換回転貼付け。