土木学会論文集 B Vol.64 No.2, , C C 3 Key Words : free-su

3

次元自由水面流中の接触を伴う

任意形状物体運動に対する数値解法

牛島 省

1

_{・福谷 彰}

2

_{・牧野 統師}

3 1_{正会員 京都大学大学院准教授 社会基盤工学専攻（〒 615-8540 京都市西京区京都大学桂 C クラスタ）} E-mail: [email protected] 2_{正会員 大阪府（〒 540-8570 大阪市中央区大手前 2 丁目）} 3_{学生員 京都大学大学院 社会基盤工学専攻 修士課程（〒 615-8540 京都市西京区京都大学桂 C クラスタ）} 3次元自由水面流中において，接触を伴う複雑形状物体の運動を扱う数値解法を考察した．この解法では，物 体を含む自由水面流れを固気液多相場として扱うことにより，物体に作用する流体力と物体が流体に及ぼす運 動量が考慮される．流体中の任意形状物体は剛体として扱われ，四面体要素により表現される．慣性テンソル 等の物体の特性量や，流体との力学的な相互作用の評価には四面体要素を利用し，物体間の衝突判定には表面 近傍に配置された接触判定球を用いる．この解法を一様流中に置かれた円柱周りの流れ，直方体形状の物体の 水中落下過程，四脚ブロックの水中での積み重なり，また波動流れによる衝突を伴う立方体ブロックの輸送問題 に適用した．これらの問題に対して，実験結果との比較により，提案された数値解法の有効性が確認された．

Key Words : free-surface flow, fluid-solid interaction, collision detection, arbitrarily-shaped body, rigid body motion, 3D MICS

1.

緒言

自由水面流れの影響を受ける物体の運動を把握する ことは，河川や海岸等における防災や環境問題，また 各種の工学的な問題を解明する上で重要な課題となっ ている．具体的な問題としては，津波や洪水時に発生 する漂流物や巨石等の輸送，浮体構造物の動揺，また 自然石で作られた水理構造物や護岸・消波ブロックの 流れに対する安定性などがあげられる．これらの問題 では，一般に物体形状は複雑であり，さらに運動に伴 い物体間の衝突が生ずる場合もある．したがって，現 象を適切に評価するには，自由水面を有する 3 次元流 体中において，衝突を伴う複雑形状物体の運動を適切 に扱う手法が必要とされる．また，物体に作用する流 体力を求める際には，抗力係数等を用いる簡易手法で はなく，自由水面流れと物体との力学的な相互作用が 適切に考慮された手法を用いることが望まれる． これらの現象を数値的に予測する手法として，これ までに流体と物体の連成作用が考慮された 2 次元解法 が検討されている1), 2), 3), 4)_{．また，著者らは，3 次} 元自由水面流中で接触を伴う球体の運動を扱う多相場 の解法 (MICS)5)_{を提案した．しかし，既往の解法で} は，3 次元流体中における複雑形状物体間の衝突までを 含めた問題を扱うには到っていない．このため，本研 究では MICS を発展させて，自由水面流れと任意形状 物体の相互作用のみならず，物体間の衝突までを扱え る解法とし，その有効性を検討する． 接触を伴う問題を扱う固体モデルとして，衝突判定 を容易に行うために，球体を連結して物体形状を模擬 する方法 (以下，球体連結モデル) が提案されている6)_． しかし，この方法では形状表現に限界があり，慣性テン ソルの設定や，物体に作用する流体力の評価が適切に行 えない場合が多い．一方，より正確な方法として，多面 体として物体形状を表現し，ボクセル7)_{や，単純形状の} 凸多面体群を用いる衝突判定法が検討されている．後者 では，物体の幾何形状を包含する境界ボリュームを階層 化して，対象物を単純な凸多面体に分解し，Lin-Canny のアルゴリズム8)_{や V クリップアルゴリズム，またよ} り効率的とされる GJK 法9)_{などを用いて，凸多面体間} の衝突判定が行われる．このような衝突判定方法は，物 体頂点部分の接触評価などが正確である反面，一般に数 値処理が複雑で，計算時間を要するという問題がある． 自由水面流れと物体運動の両者を扱う場合には，物 体形状を正確に表現することが可能で，しかも衝突判 定の計算負荷がなるべく小さい固体モデルが望ましい． このため，本研究では，任意形状物体を四面体要素の 集合として表し，物体間の衝突判定には物体表面付近 に配置された接触判定球を利用する固体モデルである， T 型モデル10)を用いる．この固体モデルでは，衝突判 定以外の処理には四面体要素を利用するので，形状表 現や流体との力学的な相互作用の評価を，球体連結モ デルよりも精度良く行うことができる．さらに，凸多 面体モデルで用いられる衝突判定法と比較して，処理 が簡単であり，計算負荷が小さいという利点がある．

本研究では，最初に既報5)_{の計算手法に対して，多} 相場中の固相の取り扱いと，固相に対する流体力の評 価法を改良した手法を示す．次に，四面体要素で表現さ れた任意形状物体運動の計算方法，流体と物体との力 学的な相互作用の取り扱い，そして物体間の接触力の 評価方法を示す．得られた解法を，固定円柱周りの流 れ，直方体あるいは複雑形状物体の水中落下過程，さ らに波動流れによる衝突を伴う物体輸送に適用し，実 験結果との比較を通じて解法の有効性を検討する．

2.

多相場の数値解法

既報5)_{では固相の密度を考慮して多相場全体の計算} が行われたが，本研究では物体を含む自由水面流れ場 に対して，最初に気相と液相のみから構成される多相 場の基礎式を解き，次に得られた結果に物体密度を考 慮して物体に作用する流体力を求める解法とした．物 体は剛体として扱われ，得られた流体力を用いて剛体 運動計算が行われる．計算された物体の運動量は多相 場に反映される． (1) 多相場の基礎式と予測段階の計算 気相と液相から構成される場の支配方程式は，多相 場の基礎式5)_{を用いて次のように記述される．} ∂ρf ∂t + ∂ ∂xj (ρfuj) = 0 (1) ∂uj ∂xj = 0 (2) ∂ui ∂t + ∂ ∂xj (uiuj) = fi− 1 ρf ∂p ∂xi + 1 ρf ∂ ∂xj [ ∂ ∂xj (µui) + ∂ ∂xi (µuj) ] (3) ここで，t は時間，xiは直交座標系の座標成分，uiは 気相と液相の質量平均速度であり，p，ρf，µ は順に気 相と液相の体積平均圧力，密度，粘性係数である．ま た，fiは外力加速度の xi成分である． 式 (1) は質量保存則より導かれ，保存形スキームを用 いて数値解を求めることにより，気液相の密度 ρf，す なわち自由水面形状が定められる．式 (3) をコロケート 格子上で離散化し，MAC 系解法11)_{に基づく予測段階，} 圧力計算段階，修正段階の演算を順に適用する．予測 段階では，セル中心における流速の推定値 u∗_i を陰的解 法である C-ISMAC 法12)_{により求める．時間方向に離} 散化された基礎式は次のように表される． u∗_i − un i ∆t = fi− 1 ρf ∂pn ∂xi − α ∂ ∂xj (u∗_iun_j)− (1 − α) ∂ ∂xj (un_iun_j) + β ρf ∂ ∂xj [ ∂ ∂xj (µu∗i) + ∂ ∂xi (µu∗j) ] +1− β ρf ∂ ∂xj [ ∂ ∂xj (µun_i) + ∂ ∂xi (µun_j) ] (4) ここに，α と β は，0 以上 1 以下の値を取るパラメータ であり，計算安定性を考慮して適切な値を設定すること により，計算時間が短縮される12)_{．C-ISMAC 法では，} 式 (4) に次の関係を用いる． u∗_i = un_i + ˜ui (5) 式 (5) を式 (4) へ代入すれば，次式が得られる． ˜ ui ∆t+ α ∂ ∂xj (˜uiunj)− β ρf ∂ ∂xj [ ∂ ∂xj (µ˜ui) + ∂ ∂xi (µ˜uj) ] = fi− 1 ρf ∂pn ∂xi − ∂ ∂xj (uniunj) + 1 ρf ∂ ∂xj [ ∂ ∂xj (µun_i) + ∂ ∂xi (µun_j) ] (6) 式 (6) 左辺に 1 次の空間離散法を適用し，右辺に高次 精度の保存形スキームを用いる．この操作により，比 較的簡単な離散化式から高精度の数値解 ˜uiを求めるこ とが可能となる12)_{．得られた ˜u} iを式 (5) に代入するこ とにより，セル中心における u∗_i が得られる．そして， CBP スキーム11)_{で生ずる圧力振動を除去するため，次} 式により u∗_i から圧力勾配項を取り除く． ˆ ui= u∗i + 1 ρf ∂pn ∂xi ∆t (7) 次に，適当な空間内挿関数 fbを用いて，セル中心の ˆui をセル境界上の値 fb(ˆui) とし，その位置で次式のよう に圧力勾配項を考慮する． ub,i= fb(ˆui)− 1 ρf ∂pn ∂xi ∆t (8) 以上が予測段階で行われる演算である． (2) 圧力計算段階と修正段階の計算 次に，圧力計算段階では，式 (8) と非圧縮条件である 式 (2) から圧力場を求める．式 (8) に基づき，n + 1 ス テップのセル境界における流速成分 un+1 b,i が p n+1_を用 いて次式で与えられるとする． un+1_b,i = fb(ˆui)− 1 ρf ∂pn+1 ∂xi ∆t (9) 式 (8) と式 (9) の差を取れば，次式が得られる． un+1_b,i − ub,i=− 1 ρf ∂ϕ ∂xi ∆t (10) ここに，ϕ = pn+1_{− p}n_{である．式 (10) を式 (2) に用い} て，C-HSMAC 法13)_{の基礎式となる以下の 3 式を導く．} ∂ ∂xi ( 1 ρ ∂ψ ∂xi ) = D k ∆t (11)

pk+1= pk+ ψ (12) uk+1_b,i = ukb,i− ∆t ρ ∂ψ ∂xi (13) ここに，上添字 k は C-HSMAC 法の反復ステップ数を 表し，Dk_{= ∂u}k b,i/∂xiである．与えられたしきい値 ϵϕ に対して，すべての計算セルにおいて_|Dk_{| < ϵ} ϕとな るまで上記の 3 式を用いる反復演算を行う．本研究で は，圧力の初期値を静水圧 p0としている．このため， 反復演算終了後の圧力は，次式のように表される． pn+1= p0+ ϕ (14) 以上のようにして pn+1_{が得られれば，これより n + 1} ステップの流速場が定められる．この演算が修正段階 に相当する． (3) 物体部分を含む多相場の基礎式 流体計算セル内に物体部分が含まれる場合には，式 (3) に対応する次の運動方程式が成り立つとする． ρb Dui Dt =−ρbgδ3,i− ∂p ∂xi + ∂ ∂xj [ ∂ ∂xj (µui) + ∂ ∂xi (µuj) ] (15) ここに，g は重力加速度，x3は鉛直上方に向かう直交 座標系の成分，δi,jはクロネッカーのデルタである．ま た，ρbはセル内に占める物体部分の質量が考慮された 多相場の密度であり，四面体サブセル法14)_{により定め} られる．四面体サブセル法は，流体計算セルを細分割し てサブセルを生成し，四面体内に含まれるサブセルの 個数から，セル内の物体体積を求める手法である．式 (14) を用いれば，式 (15) は次のように表される． Dui Dt =− ρb− ρf ρb gδ3,i− 1 ρb ∂ϕ ∂xi + 1 ρb ∂ ∂xj [ ∂ ∂xj (µui) + ∂ ∂xi (µuj) ] (16) 式 (16) 右辺第 1 項は浮力，第 2 項は流動により生ずる 圧力，第 3 項は粘性力に相当する．本解法では式 (16) は直接計算せず，以下の 3. (3) で示すように，物体に作 用する流体力として式 (16) 右辺の各項が用いられる．

3.

任意形状物体モデルと数値計算法

(1) 四面体要素分割と基本姿勢の慣性テンソル 本計算手法では，流体中に含まれる任意形状の物体 は剛体として扱われ，四面体要素の集合として表現さ れる．図–1 の例のように，最初に CAD ソフトウエア を用いて物体形状を定め，次にその出力データを四面 体格子生成ソフトウエアの入力データとして，物体を 四面体要素に分割し，四面体群の頂点座標を設定する． 物体と流体との相互作用を適切に扱うには，物体全体の (a) CADソフトウエアにより作成された物体形状 (b)四面体要素により表された物体形状 図–1 任意形状物体の形状表現 スケールに対して流体計算セルが十分細かいことが必 要である．一方，物体を構成する四面体要素に対して は，流体計算セルから生ずる制約条件はなく，物体形 状を適切に表現できる四面体要素分割を行えばよい． 物体 k の重心点座標は，それを構成する四面体群の 頂点座標から定められる．この四面体群の頂点座標を 平行移動し，物体 k の重心点が原点となるようにした 状態を物体 k の基本姿勢とする．基本姿勢における物 体 k の慣性テンソル I0,kは，物体重心点を基準とする 四面体 s の慣性テンソル I0,k,sの総和として，次のよ うに求められる． I0,k= Nk ∑ s=1 I0,k,s (17) ここで，Nkは物体 k を構成する四面体数である． (2) 剛体運動の計算方法 剛体の運動は，重心点の並進運動と重心点のまわり の回転運動に分けて扱う．図–2 に概略的に示すように， 最初に物体重心点を原点とする基本姿勢からの回転運 動を行い，次に得られた姿勢の物体を所定の位置まで 並進運動させる．この方法では，基本姿勢における慣 性テンソルのみを回転運動の計算で使用するので，各 時刻における慣性テンソルの計算が不要となる．

図–2 剛体運動計算の概略 重心点を原点とする場合に，物体 k の回転運動によ る姿勢を求める方法を以下に示す．剛体の回転運動は， 次の Euler の運動方程式に従う． ˙ Lk = Nk (18) ここで，上付きのドットは時間微分を表す．Nkは物体 k に作用するトルクであり，Lkは角運動量ベクトルで ある．ある時刻 t の Lkは，基本姿勢に対する角運動量 ベクトル L0,kを用いて次のように表される． Lk= RkL0,k (19) ここに，Rkは基本姿勢から時刻 t の姿勢への変換を表 す回転行列である．物体の基本姿勢に対する角速度ベ クトルを ω0,k とすると，L0,k = I0,kω0,k と表される ことから，式 (18) は次式のように表される． Rk ( ω0,k× I0,kω0,k+ I0,kω˙0,k ) = Nk (20) 式 (20) より，角加速度 ˙ω0,kは，次式から求められる． ˙ ω0,k= I−1_0,k ( R−1_k Nk− ω0,k× I0,kω0,k ) (21) また，基本姿勢に対する角速度ベクトル ω0,kを時刻 t の姿勢に対する角速度ベクトル ωkとするには，次の回 転操作を行えばよい． ωk= Rkω0,k (22) 得られた ωkを用いて物体の姿勢が定められる．なお， 式 (21) および式 (22) 右辺等に含まれる回転演算を行う 際には，回転行列を乗ずるのではなく，四元数を利用 する．四元数を用いる演算は，計算の精度や効率の面 で有利であることが示されている15)_． 一方，物体重心点の並進運動は，次式に従う． Mk˙vk= Fk (23) ここで，Mkと vkは物体 k の質量と速度，Fkは外力 ベクトルである．式 (23) を時間積分し，n + 1 ステッ プの物体 k の重心点を表す位置ベクトル xc,kを求める． そして，回転運動の計算で得られた姿勢の物体 k の重 心点を点 xc,kまで並進させることにより，n + 1 ステッ プの物体位置と姿勢が定められる．

C

F

Ck

object-k

r

Ck

x

c, k

図–3 流体計算セルと物体を構成する四面体要素の関係 (3) 物体に作用する流体力 物体に作用する流体力は，式 (16) 右辺から定められ る．図–3 に，流体計算セルと物体を構成する四面体要 素の概略的な関係を示す．ある流体計算セル C 内にお いて，物体 k に作用する xi方向の流体力成分を FCki とする．この FCkiは，セル体積 VC，四面体サブセル 法により求められたセル中に占める物体 k の体積割合 αk，物体密度 σkを用いて次式から求められる． FCki= αkσkVC [ −ρb− ρf ρb gδ3,i− 1 ρb ∂ϕ ∂xi +1 ρb ∂ ∂xj { ∂ ∂xj (µui) + ∂ ∂xi (µuj) }] (24) 上式より流体力ベクトル FCkが定められる．物体 k を 含む流体計算セル C に対する FCkの総和が物体 k に 作用する流体力となり，式 (23) の Fkに用いられる． また，計算セル C において，物体 k に作用する流体 力モーメントを NCkとする．NCkは，式 (24) により 求められた流体力ベクトル FCkと，物体重心点 xc,kか らセル中心点に向かうベクトル rCkの外積として，次 式により計算される． NCk= rCk× FCk (25) 流体力と同様に，関係するセルに対する NCkの総和が 物体 k に作用する流体力モーメントとなり，式 (21) の Nkに考慮される． (4) 物体運動量を多相場で考慮する方法 上記の剛体運動計算で得られた物体の運動量を用い て，多相場の流速を定める．流体計算セルと物体の関 係を図–4 に概略的に示す．物体の並進および角速度ベ クトル vkと ωkから，ある流体計算セル C の位置にお ける物体の速度ベクトル vCkを次式により求める． vCk= vk+ ωk× rCk (26) この vCkと，セル内の気液相部分の速度 uf を質量平 均して，セル C の多相場の流速ベクトルとする．セル

C

v

object-k

r

Ck k

ω

_k

x

c, k

v

Ck 図–4 物体運動量を多相場に考慮する方法 内には，複数の物体が含まれる場合があるので，セル 内に存在するすべての物体部分の質量 Mbと運動量ベ クトル mb，また気体および液体部分の質量 Mfを次式 から定める． Mb = ∑ k σkαkVC (27) mb= ∑ k σkαk∆VCvCk (28) Mf = ρf ( 1−∑ k αk ) VC (29) 以上の物理量を用いて，流体計算セル C における多相 場の流速ベクトル uC は，次式の質量平均速度ベクト ルとして計算される． uC = Mfuf+ mb Mf+ Mb (30) (5) 物体の接触判定 本研究の解法では，物体を四面体要素により表現す るため，物体表面は多面体として表わされる．多面体 表示された物体間の衝突判定法がこれまでに提案され ているが9)_{，数値処理が複雑であるものが多い．この} ため，本研究では四面体要素の重心点に中心を有する 接触球を配置し，これを利用して DEM16)_{と同様に衝} 突判定と接触力の評価を行う．このように，物体運動 計算や流体力の評価を行う際には物体を構成する四面 体要素を用い，一方衝突判定には接触判定球を用いる 固体モデルは T 型モデルといわれる10)_． 図–5 に接触判定球の配置例を示す．従来の研究では， この接触球に相当する球体要素のみを用いて物体を表 現する方法がいくつか提案されているが，本解法のよ うに，接触判定球は衝突判定のみに用いることとすれ ば，物体の慣性テンソルや質量などの特性量や，流体 との力学的な相互作用は，より精度良く評価される． 接触判定球は物体間の衝突判定のみに用いられるの で，その配置や形状は任意に設定可能である．このた め，本研究では図–6 に概略的に示すように，物体表面 図–5 接触判定球の設定例

d

n

n

v

r

v

rt

v

rn

j

i

object − m

object − k

x

c, k

x

c, m

r

k, i

r

m , j 図–6 接触判定球の配置の概略 に相当する面を有する四面体要素のみに接触判定球を 配置することとし，中心を四面体重心点と一致させ，半 径は球体体積が四面体体積と等しくなるように定めた． なお，接触判定球を四面体内接球としたり，あるいは 四面体要素と独立して物体表面近傍に密に配置するこ とにより，球面が物体外部に突出したり，球体間に隙間 が生ずるなどの問題を避けることも可能である． 接触判定球間に作用する接触力は，法線および接線方 向のバネとダッシュポット，さらに接線方向の摩擦スラ イダから構成されるモデルにより評価される．図–6 に 示す物体 k に属する接触判定球 i と物体 m の接触判定 球 j の相対速度ベクトル vrは，各物体の並進速度と， 各物体の重心点回りの回転運動を考慮して次式のよう に求められる． vr= vm− vk+ ωm× rm,j− ωk× rk,i (31) 図–6 に示すように，式 (31) の rm,jと rk,iは，各物体 の重心点から 2 つの接触判定球の接触点に向かうベク トルを表す．この vrを用いて，法線方向の接触力 fn は次式から定められる． f_n=−Kndnn + Dnvrn (32)

ここに，Knと Dnは法線方向のバネ定数および粘性係 数であり，既報17)_{と同様にして定める．d} nは図–6 に 示す接触判定球の重なり深さで，球体の中心間距離と球 体径から求められる．また，n は接触球 i の中心から接 触球 j の中心へ向かう単位ベクトル，vrnは (vr· n) n で与えられる相対速度ベクトルである． 同様に，接線方向の接触力 ftは次式から求められる． f_t= Ktdt− Dtvrt (33) この Ktと Dtは接線方向のバネ定数および粘性係数 であり，vrtは接線方向の相対速度ベクトルで，vrt = vr− vrnで与えられる．また，dtは変位ベクトルで， 接触時間にわたる積分として次式から定められる． dt= ∫ t2 t1 vrt dt (34) 接線方向接触力 f_tが法線方向接触力 f_nと摩擦係数 µ から求められる摩擦力を越えるときには滑りが生ずる． 滑りが生じた場合，あるいは 2 つの接触判定球が離れ た場合には，変位ベクトル dtはゼロクリアされる．そ れらを除く接触時刻 t1から t2にわたり，式 (34) の数 値積分が計算され，dtが求められる． 以上のようにして得られた接触力 f_nと f_tは，物体 重心点の並進運動の外力と，回転運動のモーメントと して，式 (23) および式 (21) に考慮される．また，境界 面と物体との接触判定にも上記の接触判定球が用いら れる．DEM の計算法に基づくため，同一物体中の複数 の接触判定球が同時に接触している場合でも，すべて の接触力を考慮することができる． なお，物体の運動方程式を数値積分する際には，時 間増分 ∆tbを十分小さくしないと正しい接触判定が行 えなくなる．本報では，∆tbの基準を，接触判定球の最 小径 dminと相対速度の大きさに基づく次式の接触判定 球クーラン数 Cdsから定めた． Cds= |v rn|∆tb dmin (35) 計算中はこの Cdsをモニタリングし，流体計算の時間 増分 ∆t よりも，物体運動計算の ∆tbを小さく設定す る可変時間ステップ法17)_{を利用した．}

4.

数値解法の適用性の検討

(1) 固定円柱回りの流れ 四面体要素を用いて表現された物体周辺の流体計算 が適切に行われることを確認するため，一様流中に固 定された円柱回りの 3 次元流れに計算手法を適用する． 一様流中の単一円柱回りの流れは，円柱直径 D と一様 流速 U を用いて定義されるレイノルズ数 Red= U D/ν の値に応じて変化する．既往の実験結果では，Red ≈ 40 を越えると，2 次元的な非定常後流渦が発生し，さ らに Red が 150 から 200 を越えると，後流は 3 次元非 定常流れに遷移することが知られている18)_{．このよう} な後流渦の特性と，自由水面のない流体中の円柱に対 する抗力係数 Cd が，物体を四面体要素の集合として 表現する本解法で適切に求められることを確認する． 計算では円柱径 D を 0.1 m とし，計算領域の流下方 向長さと幅をそれぞれ 1.2 m および 0.8 m，円柱の軸 方向高さを 0.5 m として，各方向の計算セル数を 120， 80 および 50 とした．円柱は，898 個の四面体要素によ り表現される．流入流速を固定して，動粘性係数を変 化させることにより各 Redにおける計算を行った． 図–7 に，計算で得られた円柱背後の流れの渦度分布 を示す．Red= 30 では後流は一対の定常渦となり，円 柱軸方向に一様な 2 次元的な分布を示す．これに対し て，Red = 100 では円柱背後に非定常な交互渦が形成 される．この非定常渦は，円柱軸方向にほぼ同一の挙 (a) Red = 30 (b) Red= 100 (c) Red= 300 図–7 円柱背後に生ずる渦度の等高線

動を示す 2 次元非定常渦であるが，図–7 (c) に示すよ うに，Red= 300 では軸方向にも一様でない 3 次元後流 渦となる．計算で得られた上記の円柱背後の流れは，実 験で確認されている後流渦の挙動とよく一致している． 図–8 に抗力係数 Cd と Redの関係を示す．計算結果 は，実験結果とほぼ一致し，本解法により物体に作用 する抗力が適切に計算されることが確認された． (2) 単一直方体の水中落下 次に，運動する物体と流れの相互作用に関する検証 を行うため，単一の直方体形状の物体が水中を自由落 下する過程を対象とする実験と数値計算を行った． 物体の初期姿勢と座標系を図–9 に示す．座標系は物 体初期姿勢の重心点を原点として空間に固定され，鉛 直上方を z 軸とする．物体の各辺の長さ a_{× b × c は，} 80× 38 × 15 mm であり，比重は 1.36 である．この 物体を水が満たされたダクト内で自由落下させる実験 を行った．ダクト底面は，x× y 方向長さが 180 × 90 mm であり，水深は 0.8 m である．初期の物体重心点 はダクト底面から 0.77 m の高さに設定し，図–9 に示 す状態を初期姿勢として落下させた． 図–8 抗力係数とレイノルズ数の関係

y

x

z

g

a

b

c

図–9 物体初期姿勢と座標系 実験では，物体の初期姿勢の設定誤差等により，自 由落下の過程には多少のばらつきが生じたが，典型的 な運動のパターンとして，周期的な x 軸回りの回転運 動と，y 方向の大きな変動を伴う状態が観察された．こ の y 方向の変動は落下とともに発達し，ダクト下方で は物体は側面に衝突した．この運動は，物体背後に交 互に形成される後流渦と y 軸方向の物体運動が連成す るために生ずると考えられる．なお，x 方向には有意な 変動は認められなかった． 図–10 落下過程の物体の姿勢(左は実験結果，右は計算結果． いずれも時間間隔は0.1秒) 図–11 物体落下過程の計算結果(各図の時間間隔は0.25秒)

数値計算では，直方体を 6 個の四面体要素に分割した． 流体計算セル数は，x× y × z の各方向に 36 × 18 × 160 とし，サブセル数は 5_{× 5 × 5 とした．} ダクト側面の x 軸方向から見た実験結果と計算結果 を図–10 に示す．これらの結果は，落下過程の各時刻 の物体の姿勢を重ね合わせて示したものである．また， 各時刻の流況と物体の姿勢の計算結果を図–11 に示す． これらの結果に示されるように，計算により得られた 物体の落下運動は，1 秒弱の周期を持つ x 軸回りの回転 運動と y 方向の変動を伴うものとなり，ダクト下方で はその運動が発達して側面に衝突する．この落下過程 は，実験結果とよく一致している．計算では，図–11 に 示されるように，物体運動と連成して交互に発生する 後流渦が再現されており，流れと物体運動の相互作用 が適切に扱われていると考えられる． 図–12 に物体重心点の z 座標の時系列を示す．物体 の平均的な沈降速度は約 0.2 m/s であり，物体が壁面 に衝突するまでの各時刻の重心点位置は，実験と計算 でよく一致している． 図–12 物体重心点高さzと時間tの関係 図–13 V (計算結果)およびVm(実験結果)と時間tの関係 図–13 に計算で得られた y 方向の物体速度 V の時系 列と，実験で計測された同方向の最大速度 Vmを示す． 実験では，物体初期姿勢の設定誤差等により y 方向の 運動の発達過程には多少の相違が見られたため，5 ケー スの実験から得られた全ての Vmを図–13 にプロット している．計算による V の時系列の極大あるいは極小 値を，実験で得られた Vmと比較すると，後者の方が 絶対値はやや大きいが，概ね良い一致を示している． (3) 消波ブロック模型の水中投入 水中に投入された消波ブロック模型が積み重なる過 程の計算を行い，複雑形状物体が流体中で相互に接触 する運動に対する解法の適用性を検討する． 実験に用いた消波ブロック模型は，モルタル製の 4 脚 テトラポッドであり，高さは約 40 mm，比重は約 2.1 である．このブロックを矩形容器内の水中に投入する． 矩形容器の底面形状は 90× 180 mm，高さは 190 mm であり，内部に水を満たして静水時の水深を 100 mm とした．容器側壁面にブロックの 3 脚を接触させ，ブ ロックの重心点を底面から 0.15 m の高さとした状態を 初期姿勢とする．この初期状態から，約 1 秒間隔でブ ロックを 10 個水中へ投入する実験を行った． 図–14 に，実験で得られた全ブロック投入後の堆積 状況を示す．水中に投入されたブロックは，後から落下 してくるブロックと衝突し，脚部を複雑に組み合わせ ながら堆積する．実験により得られた典型的なブロッ クの積み上がり状態は，図–14 に示されるようになり， ブロックの堆積角度は約 45 度となった． 計算では，格子生成ソフトウエアにより，1 個の消波 ブロックを 416 個の四面体要素に分割した．ブロック の接触判定球の摩擦係数と反発係数は，それぞれ 0.7 お よび 0.5 とし，法線方向および接線方向のバネ定数 Kn と Ktは，いずれも 1.0× 105(N/m) と設定した．水と 空気が占める矩形容器内の領域に対して，18× 36 × 38 の計算セルを設定し，実験と同様の初期条件から，時 図–14 ブロックの堆積状況(実験結果の写真)

(a) t = 0 (s) (b) t = 2 (s) (c) t = 4 (s) (d) t = 6 (s) (e) t = 8 (s) (f) t = 10 (s) 図–15 水中投入されたブロックの計算結果 間増分 ∆t を 1.0× 10−2_{秒として非定常計算を行った．} 図–15 に計算結果を示す．ブロックを順に投下して いくと，ブロック脚部が互いに組み合わさり，容器端 部付近に堆積し始める．最終的には，図–15 (f) に示 すように，実験と同様の堆積角度でブロックが積み重 なる状況が再現された．以上の結果から，本計算手法 により，複雑形状物体の接触を伴う運動が適切に扱え ることが示された． (4) 波動流れによる物体輸送 波動流れによる物体輸送の過程を対象とする実験を 行い，その結果を用いて計算手法の適用性を確認する． 実験では，造波により自由水面流れを発生させて，立 方体ブロックが，他のブロックや境界面との接触を伴 いながら移動する状況を把握した． 実験に用いた水槽を図–16 に示す．水槽左端 (下流) 側には，造波板があり，上流側には高さ 0.1 m のボック スが固定されている．水槽幅 W は 0.19 m であり，造 波板からボックスまでの距離 L1は 0.7 m，またボック スの x 方向長さ L2は 0.38 m である．造波板は PC 制 御のスライダにより上流側へ平行移動した後，停止す る条件とした．

c

D

D

_c

W

h

_o

box

L

₁

L

₂ b2

c

b1

c

b3

c

wave generator b4

c

A B B A x z y x 図–16 造波水槽(Cb1= 15.5 mm，Cb2= 13 mm，Cb3= 16 mm，Cb4= 18 mm) ボックスの上面には同一形状の 2 つの半円柱 A, B が 固定されている．半円柱の直径 Dcは 123 mm，高さは 45 mm である．ボックス上面には，一辺 30 mm，比 重 0.5 の 8 個の立方体ブロックが置かれ，それらが自 由水面流れにより輸送される．初期の水深 h0を 0.1 m とし，ボックス下流端から 0.1 m 下流側の位置の最大 波高が約 28 mm となる造波条件とした． 数値計算では，底面から高さ 0.25 m までの空気部分 を含む領域に対して，108× 38 × 50 の計算セルを設定 した．時間増分 ∆t は 1.0× 10−2 _{秒とし，実験と同じ} 初期状態から時間発展的な計算を行った．立方体ブロッ クと半円柱は，それぞれ 168 個および 490 個の四面体 要素により表現される． 図–17 に，造波開始後の各時刻における計算結果を 示す．ボックス上面の自由水面流れにより，ブロック群 は移動を開始する．半円柱 A に一部のブロックが衝突 し，その周辺で 2 個のブロックが停止する．他のブロッ クは，半円柱 A に衝突した流れにより，半円柱 B 側へ 向きを変えて輸送され，いくつかのブロックは半円柱 B に衝突して停止する．また，それらより後方にある ブロックは，前面のブロックに衝突して停止する． 図–18 に，3 回の実験により得られたすべてのブロッ ク中心点位置と計算結果を合わせて示す．図中の直線 で結ばれた丸印は，計算により得られた 0.1 秒毎のブ ロック重心点位置を表している．また，正方形は，時 刻 2.0 秒におけるブロック停止時の姿勢を表す．一方， 図–18 の大きい丸印は，実験により得られたブロック 停止時の重心点位置を表す．図中には，3 回の実験で得 られた 24 個の重心点位置が示されている．実験の初期 条件の設定誤差や流体の乱れ等により，重心点位置に は試行毎に多少のばらつきが見られるが，平均的には ある領域に集中して分布している．図–18 より，計算 により得られたブロックの停止位置は，実験結果が分 布する領域と，ほぼ一致していることが確認できる．

(a) t = 1.0 (s) (b) t = 1.1 (s) (c) t = 1.2 (s) (d) t = 1.3 (s) (e) t = 1.4 (s) (f) t = 1.5 (s) 図–17 波動流れによるブロック輸送の計算結果 0.000 0.095 0.190 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 y (m) x (m) 図–18 実験結果と計算結果の比較(大きい●印は3回分の実 験結果，他は計算結果)

5.

結言

本研究では，3 次元自由水面流れにおいて，接触を 伴う任意形状物体の運動を扱う解法を検討した．この 解法では，多相場の流体計算と剛体運動計算が行われ， 流体と物体との力学的な相互作用が考慮される．この ため，流体力の評価には抗力係数等の経験定数が不要 となる．また，任意形状物体は四面体要素により表現 され，物体間の衝突判定には物体表面付近に配置され た接触判定球を利用する．その結果，凸多面体モデル で用いられる衝突判定法と比較して計算負荷が小さく， しかも慣性テンソルなどの物体の特性量や流体との力 学的な相互作用は，球体要素のみを用いる手法と比較 してより精度良く評価される． 本解法を一様流中に置かれた円柱周りの流れ，直方 体形状の物体の水中落下過程，四脚ブロックの水中で の積み重なり，また波動流れによる立方体ブロックの 輸送問題に適用した．これらの広範な問題に対して解 法の適用性を検討した結果，計算結果は実験結果を概 ねよく再現することが確認された． 参考文献 1) 細山田得三,津田朗宏：直交矩形格子を用いた移動物体周 りの流れの数値計算,応用力学論文集, Vol. 5, pp. 699– 707, 2002. 2) 田中聖三,樫山和男：バックグラウンドメッシュに基づ くメッシュ再構築手法を用いた自由表面を有する流体– 構造連成問題のためのALE有限要素法,応用力学論文 集, Vol. 8, pp. 295–302, 2005. 3) 陸田秀美,常山鉄平,土井康明：重合ソロバン格子法によ る固体と自由表面の相互作用に関する数値解析,応用力 学論文集, Vol. 8, pp. 251–258, 2005.

4) Doig, R., Okazawa, S. and Fujikubo, M.: Solid-fluid interaction analysis by using a multi-material Eule-rian finite element method, J. Applied Mechanics, Vol. 9, pp. 151–159, 2006.

5) 牛島省,山田修三,藤岡奨,禰津家久：３次元自由水面流

れによる物体輸送の数値解法(3D MICS)の提案と適用

性の検討, 土木学会論文集, Vol. 810/II-74, pp. 79–89, 2006.

6) 熊谷兼多郎,小田勝也,藤井直樹：津波によるコンテナの 漂流挙動シミュレーションモデルの適用性,海岸工学論 文集, Vol. 53, pp. 241–245, 2006. 7) 鈴木克幸,久保田純,大坪英臣：ボクセルベース衝突判定 アルゴリズムを用いた剛体運動シミュレーション,応用 力学論文集, Vol. 6, pp. 131–139, 2003.

8) Lin, M. and Canny, J.: A fast algorithm for incre-mental distance calculation, Proc. IEEE Int. Conf.

on Robotics and Automaion, pp. 1008–1014, 1991.

9) Glbert, E. G., Johnson, D. W. and Keerthi, S. S.: A fast procedure for computing the distance between complex objects in three-dimensional space, IEEE

Journal of Robotics and Automation, Vol. 4, No. 2,

pp. 193–564, 1988. 10) 牛島省,福谷彰,牧野統師,禰津家久：3次元流体中を運 動する接触と変形を考慮した任意形状固体モデルの数値 解法,応用力学論文集, Vol. 10, pp. 139–146, 2007. 11) 牛島省,竹村雅樹,禰津家久：コロケート格子配置を用い たMAC系解法の計算スキームに関する考察,土木学会 論文集, No. 719/II-61, pp. 11–19, 2002. 12) 牛島省,禰津家久：陰解法を用いたコロケート格子によ

PREDICTION METHOD FOR MOVEMENTS WITH COLLISIONS

OF ARBITRARILY-SHAPED OBJECTS IN 3D FREE-SURFACE FLOWS

Satoru USHIJIMA, Akira FUKUTANI and Osashi MAKINO

This paper deals with a computational method for the movements of arbitrarily-shaped solid objects with collisions in three-dimensional free-surface flows. Since the mechanical interactions between objects and free-surface flows are adequately taken into account in the numerical method, no empirical coefficients, like resistance coefficients, are needed. An arbitrarily-shaped object is represented by tetrahedron elements, which are used to calculate its motion and estimate interactions with fluids, while contact-spheres are utilized for collision detections. The prediction method was applied to various experimental results and its validity has been demonstrated.

る高次精度の流体解析手法の提案,土木学会論文集, No.

719/II-61, pp. 21–30, 2002.

13) 牛島省,奥山洋平：非圧縮性流体計算におけるC-HSMAC

法とSOLA法の収束特性,土木学会論文集, No. 747/II-65, pp. 197–202, 2003. 14) 牛島省,牧野統師,禰津家久：四面体サブセル法を用いる 市街地に流入する氾濫流の3次元数値計算,水工学論文 集, Vol. 51, pp. 787–792, 2007. 15) 牛島省,福谷彰,藤岡奨,禰津家久：3次元流体中を移動す る任意形状物体の数値解析手法,水工学論文集, Vol. 51, pp. 847–852, 2007.

16) Cundall, P. A. and Strack, O. D. L.: A discrete nu-merical model for granular assemblies, Geotechnique, Vol. 29, No. 1, pp. 47–65, 1979.

17) 牛島省,竹村雅樹,山田修三,禰津家久：非圧縮性流体解

析に基づく粒子₋流体混合系の計算法(MICS)の提案,

土木学会論文集, No. 740/II-64, pp. 121–130, 2003. 18) Chen, S. S.: Flow-Induced Vibration of Circular

Cylindrical Structures, Hemisphere, 1988.