噴石模型の落下試験の動画解析
VIDEO MOTION ANALYSIS OF CINDER MODELS THROUGH FALL EXPERIMENT
劉 美智１） 丸山 敬２） 佐々木 寛介３） 井上 実４） 井口 正人５） 藤田 英輔６）
Meizhi LIU1), Takashi MARUYAMA2), Kansuke SASAKI3), Minoru INOUE4), Masato IGUCHI5) and Eisuke FUJITA6)
ABSTRACT
Falling properties of cinder models, non-simple shaped objects, were measured and three-dimensional falling trajectories were obtained by using video movies. Aerodynamic characteristics of cinder models were examined by the trajectory. The falling models imitated the configuration of cinders collected at the site of the experiment, Sakurajima volcano. Two types, one with sharp edges and one without sharp edges, were selected as representatives and a sphere was selected as reference model. The results showed that the drag coefficient varied with the falling velocity becoming smaller as the velocity increased. When the models fell with terminal velocity, the drag coefficients of the cinder models without and with sharp edges were about 0.5 and 0.6, respectively. The drag coefficient of the sphere model was about 0.3. Key Words: Video motion analysis, Cinder models, Fall experiment
１．はじめに 2014 年 9 月に噴火した御嶽山の噴火の際には 63 名もの死者・行方不明者を出す結果となり，被害者の多くは飛 来した火山礫や火山岩塊による損傷死であった 1)。さらに，噴石による飛散物はこれら人的被害だけでなく建物被害 も引き起こしている。噴石の落下による衝突の被害から免れるために，火口近くには避難所が設置されているものの， 噴石の飛散範囲や衝撃力（重さと速度）など不明な点が多いため，噴石の落下衝突に対する耐衝撃設計が十分なさ れているかどうか，明らかではないのが現状である。衝突時の衝撃力推定を行うためには，飛散物の衝突速度を特定 することが基本となり，飛散性状を精度よく予測する必要がある。火山の噴石に関する研究で代表的なものは，「火山 防災マップ作成指針」2)や「火山噴火予知連絡会火山活動評価検討会報告書」3)がある。前者の指針では，噴石の軌 道を計算する際に，形状を球と想定し，一様な空気抵抗として扱っている。後者の報告書では，噴石の軌道を計算す １） 京都大学大学院 工学研究科建築学専攻 大学院生 （〒611-0011 宇治市五ヶ庄京都大学防災研究所） ２） 京都大学防災研究所 気象・水象災害研究部門 教授 （〒611-0011 宇治市五ヶ庄京都大学防災研究所） ３） 京都大学防災研究所 気象・水象災害研究部門 特定准教授 （〒611-0011 宇治市五ヶ庄京都大学防災研 究所） ４） 一般財団法人 日本気象協会 担当部長 （〒170-6055 東京都豊島区東池袋 3-1-1） ５） 京都大学防災研究所 火山活動研究センター 教授（〒611-0011 宇治市五ヶ庄京都大学防災研究所） ６） 国立研究開発法人 防災科学技術研究所 研究部門長・総括主任研究員 （〒300-0006 茨城県つくば市天 王台3-1）
る際に，空気抵抗は考慮していない。噴石の飛散性状を精度よく予測するためにはまず，いまだ不明な部分が多い 噴石の空力特性を明らかにすることが必要である。 これまで，空力特性を明らかにする研究は風洞実験を用いて行われたものがほとんどである。立川，福山4）は平板 について2 次元的な空力特性を求めた。丸山ら 5)は，6 分力天秤を用いて正方形平板の 3 次元的な静的空力特性 を求めた。Richards ら6)は6 分力天秤を用いて長方形平板と棒状の物体の静的空力特性を求めた。このように，風洞 実験により求められる物体の空力特性は，風洞内に飛散物を一定の姿勢で固定し，一様流を与えて計測されること がほとんどである。しかし，実際の飛散物には周囲の気流性状の非定常性および自身の運動による姿勢の変化が生 じるため，そのような条件下での飛散物の動的空力特性を知ることが重要になる。これに関して，松居ら 7)は加速度セ ンサーおよび角速度センサーが内蔵された平板状模型を用意し，ドーム内の無風状態下で模型を自由落下させ，飛 散物の動的空力特性を実験的に測定した。しかし，これら既往研究の飛散物の多くは単純な形状を持つ物体が多く， 噴石のような不整形物体による実験は，ほとんどないのが現状である。 本研究では，噴火時の噴石の飛散運動を精度よく予測するための空力特性を明らかにすることを目的として，不整 形な形状を持つ噴石の自然風中における落下実験を行った。この落下実験で得られる空力特性が，将来的に，同じ 形状の模型を使用した風洞実験から得られる空力特性との差異を知ることも本実験の目的である。実験は桜島の麓 でドローンにより噴石模型を上空に持ち上げて落下させて行った。落下運動の測定結果を解析することにより得られ た噴石のもつ空力特性を用いて，噴石の飛散や衝撃力のより精度の高い推定を行い，噴石に対する安全な建物の 設計，避難行動計画の作成，自治体等の災害対応などに貢献することが期待される。 ２．実験概要 3 種類の噴石模型（球体模型も含む）をドローンを用いて高度約 150m まで持ち上げ（図 1），落下させることにより， 自然風環境下での噴石の落下性状を測定した。落下の様子は複数のビデオカメラで撮影し，3 次元的な落下運動の 解析を行った。画像解析で得た落下軌跡は，模型の形を点としてとらえているので，本研究は回転運動を考慮してい ない。 ２．１ 実験実施場所 落下実験は，鹿児島県桜島黒髪地獄河原で行った（図2）。実験を行った場所は，昭和噴火（1946 年）の際に溶岩 が流れた場所であり，樹木がほとんど生えておらず視界を遮るものがない。そのため，動画撮影としては適切な実験 地点が確保できた。 模型投下地点 図2．実験実施領域 図1．ドローンによる模型吊り上げの様子
表1． 模型の概要 ケース 番号 模型種類 及び番号 質量 m[kg] 見付面積A[m2] case.1 A1 0.41 0.05 case.2 B1 0.42 0.05 case.3 C1 0.39 0.07 case.4 A2 0.41 0.05 case.5 B2 0.46 0.05 case.6 C2 0.47 0.07 ２．２ 噴石模型 実験現場で採取された噴石から，代表的な形を選択した。表面に細かな凹凸があるものの，全体としては比較的 滑らかな形状を持ったもの（図3 模型 A）と，表面に鋭角な角をもったもの（図 3 模型 B）に大別できた。実験では，そ れらを模倣し，図3 に示す 20cm-30cm 程度の大きさの模型として発泡ウレタンで作成した。さらに，空力特性の実験 データが多い球体を比較検証用に採用し，直径30cm の発泡スチロール製の球を実験模型として用いた。 表1 に模型の質量，見付面積，および，対応するケース番号を示す。見付面積 A は模型を回転させて得られた画 像から求めた平均的値を用いた。 ２．３ ビデオカメラによる撮影 図4 に示すように，基準点 O をほぼ正三角形に配置された A, B, C 地点の中心に設け，その周囲 4 ヶ所 A, B, C, N 地点に 4K ビデオカメラを配置した。ビデオカメラ（Panasonic 社製 HC-VX985M）の画素数は 3840 pixel×2160 pixel，フレームレートは 29.97fps である。 基準点O には高さ 1m にマーカーとなる白い球体を設置し，その近くをドローンの発着場所 P とし，そこから上空に 模型を吊り上げ，ほぼ真上から落下させた。カメラの画角は噴石模型の落下軌道が画像内に入るように調整した。 表 2 に動画解析に使用した画像の撮影地点を「〇」で示す。解析に際して，適切な画像が得られなかった場合が あり，4 つの地点全ての画像を用いることができないことがあった。その理由としては，カメラの自動焦点が鳥や虫，雲 の動きを捉えてしまい，模型に焦点が合わず画像がぼやけている。逆光のためコントラストが下がり，特に地面近くで 背景と模型の区別がつかず，模型の位置が判別できない。などが挙げられる。また，case.2 の模型 B1 および case.3 の模型 C1 では，模型の落下軌道が風に流されて大きく外れ，撮影動画の外に出てしまったため，地面近くのデータ が大きく欠如した。
以下には上記のような不都合が少なく，解析に条件の良い測定結果（case.1 の模型 A1, case.5 の模型 B2，case.6
の模型C2）を選んで解析結果を示す。 ケース 番号 地点Ａ 地点Ｂ 地点Ｎ 地点Ｃ case.1 〇 〇 〇 - case.2 - - 〇 〇 case.3 〇 - - 〇 case.4 - 〇 〇 - case.5 〇 〇 〇 - case.6 〇 〇 〇 - 模型C 模型A 模型B 30c m 図3．噴石模型の種類 基準点O C 地点 N 地点 B 地点 一辺約200m A 地点 図4．撮影地点 N 表2. 動画解析に使用した画像の撮影地点
図5．平均水平風速 UV 図6．乱れ強さ Iz ２．４ 風速プロフィール 噴石模型を落下させる際の，周囲の気流性状をドップラーライダー（三菱電機社製DIABREZZA）と超音波風速計 （カイジョー製 SAT900）を使用して観測した。ドップラーライダーは上空にレーザ光を照射して上空の風向風速を計 測する装置であり，40m 以上の高度のデータが得られた。地面付近の風向風速は超音波風速計による観測値を用い てデータを補った。 図5，図 6 に噴石模型（A1，B2，C2）投下時刻における，前後 20 秒の平均水平風速 UV と乱れ強さ Iz のプロフィ ルを示す。今回の実験の範囲は高度0m-150m であり，実験時の落下の範囲内において周囲の水平風速 UV 最大 値は，模型A1，模型 B2，模型 C2 の落下に際して，それぞれ約 1.0m/s，2.4m/s，1.7m/s とほぼ静穏な状況であった。 噴石模型落下直後の高度150m 近くにおける乱れ強さに際しては，模型 A1 の乱れ強さが，模型 B2 模型 C2 より大 きい結果となった。 ２．５ 模型投下に用いたドローン 図 7 に噴石模型を吊り上げ，投下に使用したドローン（ルーチェサーチ社製 SPIDER CS-6）を示す。ドローンには 模型投下装置と位置情報記録用スマートフォンが搭載されている。 図8 はドローンの飛行記録から鉛直方向の変位を示した。噴石模型の落下高度が約 147ｍと，予定していた 150m 高さに近い値を確認できた。 図7．落下実験に用いたドローン 図8．ドローンの飛行記録（鉛直方向の変位） 0 50 100 150 200 0 1 2 3 4 5 A1 B2 C2 h [m ] UV [m/s] 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 A1 B2 C2 h [m ] Iz 0 50 100 150 200 0 2000 4000 6000 8000 高度 [m ] データの数 147m
３． ビデオカメラによる落下軌跡の解析 噴石模型の落下の様子を複数のビデオカメラを用いて撮影し，ビデオ映像を 解析することにより，噴石模型の軌跡を求めた。動画解析には 3 次元運動解析 ソフト（DITECT 社製 Dipp-Motion V）を用いて，噴石模型の落下運動の 3 次元 軌 跡 を 求 め た 。 画 像 解 析 の ア ル ゴ リ ズ ム と し て ，DLT （ Direct Linear Transformation）法 8)を用いた。これは，複数台のカメラと空間座標が既知である 基準点を用いて映像内の3 次元空間を推定する方法であり，解法上，計測空間 内に6 点以上の基準点を配置する必要がある。 今回の実験では，模型落下の直前にドローンをビデオ画像中で移動させて 暫く静止させ，飛行記録からその座標を求めて，画像解析に必要となる基準点 を定めた。図9 に A 地点から撮影したビデオ画像内の基準点の例を示す。 画像解析用の基準点は，発着地点P から真上にドローンを 150m まで上昇さ せ暫く静止させた後，東西南北，4 方向に水平移動させ，そこで上下にドローン を移動させて記録している。上下移動に際しては高度50m，100m，150m の 3 高度で静止して合計 15 の基準点を記 録した後，再び発着地点P に戻る順序で基準点を記録した。図 10 にドローン飛行記録の立体図を示す。図 11 は地 点A，地点 B，地点 N から撮影した画像である。図 10 および図 11 中にある赤丸，白抜き三角，白抜き丸の意味は同 じである。赤丸は基準点位置を示している。さらに，ドローン着地点 P 近くに設置したマーカー（図中白抜き三角）と， 落下後に計測した模型の着地点（図中白抜き丸）を地面付近の基準点として加えた。得られた基準点全てを解析に 取り込むと，上空のデータの重みが増し，求めた落下軌道の精度が悪くなる場合があったので，図11 のように上空の 基準点はできるだけ各基準点が離れるように6 つ使用し，地上の基準点を含めて合計で 8 つの基準点を用いて解析 した。 ４．解析結果 ４．１ 画像解析による３次元軌跡
画像解析により得られた噴石模型A1・case.1，模型 B2・case.5，模型 C2・case.6 の 3 次元軌跡の結果を図 12 に示
す。図の座標はドローン発着地点P を原点としている。xyz 座標は右手系となっており，x 軸の正方向は東，y 軸の正 方向は北，z 軸の方向は高さを表している。どの模型も地面近くまでの軌跡が捉えられている。画像解析から得られた 図10．ドローンデータによる基準点位置 図11．動画解析に使用した基準点 （左図からA 地点，B 地点，N 地点） 図9．A 地点からみた基準点 y [m]_north x [m]_east h [m]
(a) 模型 A1 (b) 模型 B2 (c) 模型 C2
図12. 画像解析から得られた噴石の落下軌跡
軌跡最低点の水平面地点 (xanalyzed, yanalyzed)，および模型落下地点の実測値 (xmeasured, ymeasured) を表 3 に示す。軌跡 が着地時点まで得られなかったのは，地形の段差や前後の草などで遮られたために，ビデオ画像では着地までの映 像が映っていないためである。一方，模型落下地点の実測値は，噴石模型が地面に到達後，落下地点に行き GPS により求めた。したがって，落下した際に模型が跳ねたり転がったりして，落下地点から移動した後の位置である。 表3 より，模型 A1 の落下地点が (xmeasured, ymeasured) = (-4.26m,-10.44m) であるのに対し，画像解析で得られた落 下地点は (xanalyzed,yanalyzed) = (-4.12m,-9.85m) となっており，14cm-59cm の誤差を生じた。模型 B2 では 1.19m-1.65m， 模型C2 では 0.08m-1.93m の違いがあった。誤差は最大約 2m の結果となった。 ４．２ 変位，速度，加速度 図13~15 に噴石模型(A1, B2, C2)の変位，速度，加速度を示す。各図の(a), (b), (c)は変位，速度，加速度を示して いる。速度ベクトル，加速度ベクトルは画像解析から求まった変位の時刻歴から差分で求めた。また，速度および加 速度の導出の際に1 秒の移動平均を取っている。 模型A1 に関しては図 13(a)より，x 軸と y 軸方向にほとんど移動していないことがわかり，模型はほぼ真下に落下 したことが分かる。模型B2 に関しては図 14(a)より，落下開始位置が水平面座標で発着地点 P から（-3.4m, -3.4m）離 れた位置から落下し始め，落下地点が（-4.6m,-12.2m）となった。模型 C2 に関しては図 15(a)より，落下開始位置（-2.6m, -2.6m）から落下地点（-9.0m,-14.8m）となった。 各図の(b)より，鉛直方向の速度に関しては，どの噴石模型も約 3 秒後から終端速度へ安定し始めている（模型 A1 に関しては顕著に表れていないが，約3 秒以降から速度の変化がそれ以前よりも小さくなっている）。3 秒以降の鉛直 方向速度の平均を取ると，噴石模型 A1，B2 に関しては約-15m/s となった。球状模型 C2 では，約-19m/s となり，噴 石模型A1，B2 よりも早い速度で落下していることが分かる。 各図の(c)より，模型落下開始時刻から約 3 秒の範囲では，落下速度の変化が大きくなるとともに，徐々に模型が受 表3. 模型落下地点と画像解析から得られた模型落下地点
model measured data (xmeasured, ymeasured) [m] analyzed data (xanalyzed, yanalyzed) [m] A1 -4.26 -10.44 -4.12 -9.85 B2 -6.25 -11.00 -4.60 -12.19 C2 -9.09 -12.89 -9.01 -14.82
(a) 変位 [m] (b) 速度 [m/s] (c) 加速度 [m/s2] 図13. 噴石模型 A1 の変位，速度，加速度 (a) 変位 [m] (b) 速度 [m/s] (c) 加速度 [m/s2] 図14. 噴石模型 B2 の変位，速度，加速度 (a) 変位 [m] (b) 速度 [m/s] (c) 加速度 [m/s2] 図15. 噴石模型 C2 の変位，速度，加速度 ける空気力が増え，重力加速度から鉛直方向の加速度が小さくなって0 に近づいていくのが分かる。模型 A1 と模型 B2 では，約 3 秒以降で加速度の値が一定になっているのが分かる。球状模型 C2 では，約 5 秒まで全方向で加速 度が他の2 つの模型に比べて大きく変化している。 ４．３ 抗力係数 図16 の(a), (b), (c) にそれぞれの噴石模型(A1, B2, C2) の抗力係数 Cd と，風向と風力方向の偏差𝜃𝜃を示す。 図16(a)，16(b)における抗力に関しては，噴石模型 A1，B2 どちらも落下と共に速度が増すにしたがって抗力係数 の値が小さくなっている。どちらの模型も落下後 3 秒以降で，抗力係数の減り方が緩やかになっていることが分かる。 その後，速度がほぼ一定に達した際の抗力係数の値は，模型A1 では約 0.5，模型 B2 では約 0.6 となった。 図16(c)における球体模型 C2 の抗力係数は，噴石模型よりも早い時刻（落下後約 2 秒）で，ほぼ一定に達してお り，その値は約0.3 となった。球の抗力係数に関してはレイノルズ数との関係が得られている。動粘性係数を乾燥空気 20℃の場合（1.5×10-5m2/s）として，本実験の球体落下速度からレイノルズ数を求めると約 4.0×105 以内となった。レ イノルズ数最大値での抗力係数は約0.1 である。今回行った球体での実験結果は 0.3 と，これまで得られている 0.1 -100 -50 0 50 100 150 200 0 2 4 6 8 10 12 x y z dis pla ce m en t[m ] time[sec] -20 -10 0 10 20 0 2 4 6 8 10 12 x y z ve lo city [m /s ] time[sec] -20 -10 0 10 20 0 2 4 6 8 10 12 x y z ac cel er at io n[ m 2/s] time[sec] -100 -50 0 50 100 150 200 0 2 4 6 8 10 12 x y z di sp lacem en t[m ] time[sec] -20 -10 0 10 20 0 2 4 6 8 10 12 x y z ve lo city [m /s] time[sec] -20 -10 0 10 20 0 2 4 6 8 10 12 x y z ac cel er at io n[ m 2/s] time[sec] -100 -50 0 50 100 150 200 0 2 4 6 8 10 12 x y z di sp lacem en t[m ] time[sec] -20 -10 0 10 20 0 2 4 6 8 10 12 x y z ve lo city [m /s] time[sec] -20 -10 0 10 20 0 2 4 6 8 10 12 x y z ac cel er at io n[ m 2/s] time[sec]
(a) 模型 A1 (b) 模型 B2 (c) 模型 C2 図16. 抗力係数 Cd および風向と風力方向の偏差𝜃𝜃 [°]の時間変化 より大きいが，解析で得られたレイノルズ数は臨界レイノルズ数の近くになっている。臨界レイノルズ数近くでは，風速 の乱れや物体表面の粗さが大きく影響し，レイノルズ数の変化による抗力係数の変化も大きいことを考えると，今回得 られた球体に対する抗力係数の値の結果は妥当な範囲に入っていると考えられる。 噴石模型A1 および B2 の落下速度からレイノルズ数を計算すると，模型 A1 は約 2.9×105以内となり，模型B2 は 約 2.7×105以内となる。本研究は，空気の抗力が支配的となる，2mm-64mm の大きさで定義される「火山礫」を想定 している。粒径64mm の実際の噴石を想定すると，速度が 68m/s 以下で飛んでくることを対象としている。 風向と模型に加わる風力の方向の偏差𝜃𝜃を見ると，模型 A1 の場合，落ち始めて１秒後には約 10°であるが，2 秒後 には約1°と小さくなり，その後 6°前後の範囲で交互に変化している。このように，鋭角な角のない模型 A1 の場合は， ほぼ真下に落下しつつも，風力を受ける向きは風向方向と少しずれて変化している。鋭角な角をもつ模型B2 の場合 は，落下後6 秒近くまでは，偏差は約 7°と一定の値を保っている。鋭角な角のある場合は，風の剥離点が固定されや すくなり，模型周囲の流れ場の形状が変化しにくくなって，風力の方向も安定しやすくなると考えられる。一方，球体 C2 の場合は，落下後約 2 秒の間偏差は約 30°となり，他の模型よりも大きな値となっている。その後約 4 秒にかけて 風向偏差は小さくなっており，模型周囲の流れ場の形が大きく変わっていることがわかる。その後の偏差は約 4°から 8°の範囲になり，物体が受ける風力の方向は風向に近い向きとなった。 ５．まとめ 本研究では，噴石模型の落下性状をビデオカメラで撮影し，得られた動画を用いて模型の落下の3 次元運動を解 析し，得られた落下の軌跡から噴石のような不整形物体の持つ空力特性を求めた。噴石模型の形状は，実験を行っ た現地で採取された噴石を模し，鋭角な角のあるものとないものの 2 種類を代表として選んだ。実験から得られた抗 力係数は落下にともなう速度の増加とともに小さくなり，終端速度に達した際の値は，鋭角な角のない噴石模型で約 0.5，鋭角な角のある噴石模型で約 0.6 となった。また，比較として用いた球体の抗力係数は約 0.3 となった。 謝辞 本研究は，京都大学防災研究所令和元年度一般共同研究 30G-10，および，次世代火山研究・人材育成総合プ ロジェクト・課題 D：火山災害対策技術の開発・サブテーマ 2：「リアルタイムの火山灰ハザード評価手法の開発」の補 助を受けて行われたものである。実験にあたっては小宮拓海氏，瀧下恒星氏，米田格氏，山崎友也氏，加茂正人氏， 土井こずえ氏の協力を得た。ドローン操作にあたっては，志村智也氏，依田隆志氏，二宮裕之氏の協力を得た。ここ に記し，謝意を表します。 参考文献 1) 及川輝樹, 山岡耕春, 吉本充宏, 中田節也, 竹下欣宏, 前野深, 石塚吉浩, 小森次郎, 嶋野岳人, 中野俊，御 嶽山2014 年噴火，火山, 第 60 巻, 第 3 号，pp.411-415, 2015 0 1 2 3 4 0 10 20 30 40 0 2 4 6 8 10 12 Cd theta C d [-] thet a [d egr ee ] time [sec] 0 1 2 3 4 0 10 20 30 40 0 2 4 6 8 10 12 Cd theta C d [-] the ta [ degr ee] time [sec] 0 1 2 3 4 0 10 20 30 40 0 2 4 6 8 10 Cd theta C d [-] the ta [ degr ee] time [sec]
2) 内閣府防災担当，消防庁，国土交通省水管理，国土保全局砂防部，気象庁，火山防災マップ作成指針，2013 3) 気象庁，火山噴火予知連絡会 火山活動評価検討会 報告書，2014 4) 立川正夫，福山雅弘，台風時の飛散物の軌跡と速度に関する研究 -その１平板の空力特性と運動方程式-，日本 建築学会論文報告集，第302 号，pp.1-11, 1981 5) 丸山敬，河合宏允，樋本圭佑，6 分力天秤を用いた正方形平板の空気力特性の測定，日本建築学会大会学術 講演梗概集，pp.177-178，2005
6) Peter J. Richards, Nathan Williams, Brent Laing, Matthew McCarty, Michael Pond, Numerical calculation of the three-dimensional motion of wind-borne debris, Journal of Wind Engineering and Industrial Aerodynamics, Vol.96, pp.2188-2202, 2008
7) 松居健人，丸山敬，西村宏昭，野田博，自立型計測装置を用いた飛散物模型の空力特性の直接計測の試み， 第25 回風工学シンポジウム論文集，pp.169-174，2018
8) Robert Shapiro, Direct Linear Transformation Method for Three-Dimensional Cinematography, Research Quarterly, Vol. 49, No. 2, pp.197-205, 1978
