流砂量計算ソフトの開発：鶴巻有一郎（2005.8.1：母の誕生日を祝して）
まえがき
地上での土砂の流出は、雨水による侵食、崩壊、地すべりがある。河川内の土砂輸送現象を見ると、流水力によって流送される砂礫の現象と、河床堆積層の全体が急勾配のため安定性を欠き全体が流動する土石流の現象がある。
ここで扱うものは、流水力によって流送される砂礫の現象の一部分である。水理学では、流水力によって流送される砂礫を、三種類に分類しており、掃流砂、浮遊砂、ウォシュロードである。「掃流砂」は、河床の砂礫を流水力によって、押し、河床上を転動、飛び跳ねする河床近傍に見られる比較的大きな径の砂礫の現象を言う。「浮遊砂」は河床の砂礫が浮遊し流下するものの、流速が遅くなると沈殿し始める砂礫を言う。「ウォシュロード」は流速が０になってもなかなか沈殿しないシルトである。一般に、流砂量計算公式と言うと、「掃流砂」、「浮遊砂」に関するものが多く、ここで扱うものはこの二つである。
河床材料には混合粒径を取扱い、流砂量の算定に当たっては、岩垣公式、佐藤・吉川・芦田の式、Brownの式、 Einsteinの式、LaneとKalinskeの式等を用いる。これら各式を定点の安定断面形における年間流況の下の年間総流出土砂量を算出する。その結果を見ると、各式による値は、10倍から100倍程度異なり、驚かされる。また、Einsteinの理論の難解さに興味惹かれるが、丁寧に解説したものは見当たらず、全式の誘導を行うには至っていない。このため、理論の説明は行わない。プログラム作成は、将来、河床変動解析プログラム開発に繋げる配慮と豊富なコメントを記した。
演算結果を再整理した例
ところで、現在、HPの検索で流砂量計算ソフトを検索すると実際にその物ずばりは、見つからない。最も、検索機構で、本HPを「水理学」や「水工水理学」の中から探そうとしても、皆無であるから、「流砂量計算ソフト」は存在するかもしれないが、あったとしても、極小数なのであろう。
このため、本HPをご覧になられ、かつ、流砂量に興味を持ち、勉強中の方はラッキーなのである。確かに、需要の少なさと、その厄介さから、特定の会社・研究所の独占物なのかもしれない、と思うと益々、公開したくなってきました。さらに、「水理学」の検索で現われるものは、書籍の宣伝、大学講座の案内、有名先生の経歴、、、と無意味なものが多い中、世に隠れた充実したHPを目指したいものである。本HPが良かったと思われる方は、口コミの宣伝や、リンクしていただければ幸いである。
流砂量計算ソフトと河床変動解析ソフトは、若かりし頃、Fortranで作りましたが、ここに、VBAで新たに作ったものを公開します。ソフトは、末尾に掲載していますので、各自でVBAを立ち上げ、貼り付けてくだされば、ウィルスの心配も不要である。
なお、流砂量公式は、次の書籍を参考にしました。
石原藤次郎・本間仁編 応用水理学中I，丸善株式会社pp5-36
水理学演習下巻荒木正夫、椿東一郎共著、森北出版
Ⅰ.1 入力データの整理作業
1.1 各種公式図表の入力
水理学演習下巻荒木正夫、椿東一郎共著、森北出版；pp197.図10.9から
をx=τc/τ0 ,y=f(τc/τ0)として読み取り近似式を求めると
0>x>0.5に対し
y = -0.3853716 x4 - 1.1386687 x3 - 0.5997683 x2 + 0.0306731 x + 1.0000094
0.5>x>0.903974に対し
y = -8.965847 x4 + 29.424557 x3 - 33.835262 x2 + 14.770071 x - 1.345754
0.903974>xに対し
y = -9.515844 x4 + 38.126980 x3 - 55.739295 x2 + 34.688015 x - 7.512319
となる。この関係式はサブルーチン「 DOKEN(tauc, tau0, mani, s,
qBD, ua) '掃流砂量公式：佐藤・吉川・芦田の式 」に組み込まれている。
アインシュタイン公式については、
石原藤次郎・本間仁編 応用水理学中I，丸善株式会社pp5-36を参照
のd65U*/11.6νとｘ欄の数値をsheet7のセル（B56，C71）に記入してください。サブルーチン 「Einxx(d65, myu, ua, x, Rx) 'アインシュタインの浮遊砂量の係数ｘ,Xの算定」 で使用されます。
のLOG(ｘ)とLOG(Y)欄の数値をsheet7のセル（D113，E139）に記入してください。サブルーチン 「EinY(d65, ua, RY) 'RY=Y＝揚圧力補正係数をX=d65U*/11.6υの関係図から求める」 で使用されます。
のLOG(ｘ)とLOG(Y)欄の数値をsheet7のセル（F169，G194）に記入してください。サブルーチン 「EinGUZAI(Rx, d, guzai) 'アインシュタインの浮遊砂量の係数guzai=ξの算定」 で使用されます。
のLOG(x)とLOG(y)欄の数値をsheet7のセル（F208，G222）に記入してください。サブルーチン 「EinFaiPai(FaiGYakuA, PaiA) 'アインシュタインの１/Ψ*からΦ*を求める」 で使用されます。
のLOG(x)とLOG(y)欄の数値をsheet7のセル（F244，G267）に記入してください。サブルーチン 「EinRaRwa(Ra, Rwa) 'R*からRwaを求める」 で使用されます。
また、アインシュタインの浮遊砂量の濃度分布を算定するI1，I2は、図表によらず、直接にシンプソン公式を用いて下記の式を積分している。
参考書籍での上式のZ文字は大文字のZを、ｙは小文字のｚと記され、見誤り易いことから、ここでは、ｙを用いて表示した。
積分範囲a/hから１の間の分割数はa/hからa/h間隔に20個、次の間隔は10 a/hを20個、次の間隔は100 a/hを20個、次は残りの1までの間を20等分した合計８０分割で算定している。各分割過程の関数値をsheet6に表示している。このサブルーチンは「Integral102(aBh, h, Rz, I1, I2)
'アインシュタインの浮遊砂量isqs算定のための積分値算定」で行われている。
Lane & Kalinskeのw0/U*とPの関連も図表読み取りではなく、直接、100分割の数値積分を行っている。Sheet6にその結果を示す。このサブルーチンは「IntegralKalin(w0, uab, v, P) 'Lane & Kalinskeによる浮遊砂量isqs算定のための積分値Pの算定」で行われている。関係式は、次の通りである。
ここにC0の単位はgr/m3、Fは沈降速度w0なる砂が河床材料に占める割合（％）である。
前節の読み取り図表を用いているため、その再現性とプログラムのチェックを兼ねて、アインシュタイン混合粒径における、浮遊砂量と掃流砂量を「石原藤次郎・本間仁編 応用水理学中I，丸善株式会社pp24-36」の例題を行うこととする。本ソフト使用の計算結果と例題の答えを比較すると下図の如くである。なお、例題の答えを黄色セル、本ソフトの値の比較セルを青セルで表す。
以上の如くであり、ほぼ満足する結果である。なお、上記の演算結果はsheet5に表示される。また、ここでの例題ではd65＝0.35mmであるため、サブルーチン「EinRbd(v, Ie, x, d65, R, ua, myu, Rx, Rbd) 'アインシュタインのU*に係わるRb'の算定」の冒頭で「'd65 = 0.035 '応用水理学中I：石原・本間編pp.34-36の例題をチェック計算するときは’をはずせば良い」を挿入している。しかし、最大篩目径に0.9mmを入力した場合のd65は0.30858mmと算出され、例題の径にはならないが、この算定値はサブルーチン「BedMaterial(kkn, s, da, dm, d65, myu) '河床材料の粒度分布から平均粒径、65%粒径を求める」で求められる最確値である。
目次↑
既に、本HPで公開済みの不等流計算ソフトではXY座標を用いた河川断面形が採用されている。しかしながら、ここでは、平均的な流砂量を求めていることから、変形断面形ではなく、台形断面形に変換し、その河床幅に発生する流砂量を扱うこととする。
断面形の変換に当たっては、XY座標の変形断面形で求められた面積、径深、水面幅、標高値を採用するももの、その値を正しく再現する流量、マニングの粗度係数、河川勾配をまず定める必要がある。
流量は、その河川の年平均洪水流量を、かつ、そのときのマニングの粗度係数を採用する。それは、流砂量の支配的な流量として、ここでは採用する。さらに、水位には、等流状態の水位を算定するが、その水深には、径深を採用することとする。すなわち、流砂量には、平均的支配流量時の平均的な台形断面形を作成し採用するのである。
具体的には、次のように行う。
① 公開済み不等流計算ソフトのXY座標を用いた河川断面形を整理した結果を表すsheet5をコピーして、本ソフトのsheet２に貼り付ける。
② 支配流量Q、 マニングの粗度係数ｎ、河川勾配I＝sin(θ)を定める。
③ 等流水位を算定するために
の右辺項を算定し、sheet２セル（E14,M14）の中からこの値を超える値の列を探した後、その列の15行と16行の係数を用いてｈを求める。これが、等流水深h0である。
④ h0に対する径深Rを、その列の12行と13行の係数を用いて定め、さらに
を算定してAを定める。
⑤ 17行の「hに対応するA3/Bの値」の行の、その列とその左側列の値と、9行目のｈの関連からh0に対応するA3/Bを算定する。
⑥ このA3/Bの値と④で求めたAからBを算定する。このBは水面幅である。
⑦ Aを台形の面積とすると
から河床幅Bbが定められる。
⑧ XY座標断面形の河床標高Ｚ(x,y)にh0を加えたのが等流水位である。
以上が断面形の整理である。
この台形断面を他の流量に使用した場合の等流水深とその場合の断面積は、次のように求める。
⑨ まず，上式の右辺項を求める。これをCと置く、
⑩ 径深Rと水深hがほぼ等しい河川断面として扱う。このため上式は次のように変換される。
⑪ 支配流量時の等流水深h0と水面幅Bs0の関係とh,Bsの関連は次のようである。
⑫ 上式を満足するh,Rを算定する。添え字１は第一近似値を、２は第二近似値を表す。
なお、sheet２への貼り付け例は下図のようである。
目次↑
流砂量の算定に当たっては、岩垣公式の限界掃流力を算出するサブルーチン「IWAGAKI da, uac2, tauc」をもうけているが、その結果を特に出力していない。掃流砂量には佐藤・吉川・芦田の式を採用し、サブルーチン「DOKEN tauc, tau0, mani, s, qBD, ua」で演算を行っている。掃流砂量と浮遊砂量を含む式は、Brownの式と Einsteinの式を扱う。Brownの式は、サブルーチン「BROWN ua, s, dm, qBB, Fai」に示す。Einsteinの式は、サブルーチン「Einstein kkn, d65, myu, ua, v, Ie, R, qTE, EqTh, h, TqsLK」に示す。また、浮遊砂量のみの式は、LaneとKalinskeの式である。これは、Einsteinのサブルーチンの終端で演算されている。すなわち、サブルーチン「LaneKalinske w0, uab, v, ib, h, qsLK 'Lane & Kalinskeの浮遊砂量」に示されている。
入力はsheet1の黄色のセルである。内容と位置は次の通り。
流量分割総数＝セル(C3)； 断面内分割数(天端から最深部間の分割数)＝セル(H2)；河川断面形見直しの流量（xY断面形を台形断面に変換）＝セル(E4)；流量（m3/s) ＝セル(A８，A行)；継続時間(hr) ＝セル(B８，B行)；粗度係数＝セル(J8)；河床標高（EL.ｍ）＝セル(L8)；河床勾配＝セル(Q8)；篩の目(mm) ＝セル(T7，AA7)、なお、最大値は必ず数値を入力すること；△ｐ通過百分率（％）＝セル(U８，AA８)
その他、sheet2,sheet7への入力はすでに述べた通りである。
Sheet1への入力例
演算結果を再整理した例（緑はこのために追加、黄色と結果を接続表示）
演算結果は、定点の安定断面形における年間流況の下の年間総流出土砂量である。その結果を見ると、各式による値は、10倍から100倍程度異なる。そこで、さらに、他の面からの検討を加え、対象河川に比較的適合する流砂量公式を探す必要がある。まず流域の年間侵食高を検討する。わが国の貯水池容量は、比較的小さく、ダムへの堆砂量は、掃流砂量と浮遊砂量であり、ウォシュロードはほとんど堆積しない。その堆積量から、年間の流出土砂量が把握される。まず、そのようなデータと比較検討する。
さらに、河床は変動するものであり、かつ、粒径分布構成も流量に応じて変動（アーマリング）する。すなわち、河道区間の上流端には、上流で生産される土砂量が流入する。そのような、河床変動計算を何年間にも亘って、行い始めて、その河川の特性が把握される。次回は、本プログラムをそこまで発展させたいものである。
流砂量計算プログラム：３９９KB RiverBedLoad.xls
プログラム・ソフトのダウンロードは都合により移動しましたここをクリックしてください。
流砂量計算プログラム
（06.5.15）鶴巻有一郎
演算実行は「RiverBedFluctuation」をクリックする。ソフトの内容は以下の通りである。
