第 26 回数値流体力学シンポジウム
講演番号 E07-1

3 次元流体・構造・地形変化連成数値計算モデルを用いた海上を浮遊するコンテナの風による漂流挙動に関する数値シミュレーション
Numerical Simulation of Wind-Induced Drift Behavior of Shipping Container Floating on
Water Surface using Three-Dimensional Coupled Fluid-Structure-Sediment Interaction Model○ 中村友昭, 名古屋大学高等研究院, 名古屋市千種区不老町, E-mail: tnakamura@nagoya-u.jp
水谷法美, 名古屋大学社会基盤工学専攻, 名古屋市千種区不老町, E-mail: mizutani@civil.nagoya-u.ac.jp
Tomoaki NAKAMURA, Inst. for Advanced Res., Nagoya Univ., Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan
Norimi MIZUTANI, Dept. of Civil Eng., Nagoya Univ., Furo-cho, Chikusa-ku, Nagoya 464-8603, JapanA three-dimensional coupled fluid-structure-sediment interaction model (FSSM) is applied to simulate wind-induced
drift behavior of a shipping container floating on the water surface. The FSSM consists of an extended Navier-Stokes
solver for incompressible viscous multi-phase flow, a volume-of-fluid module for air-water interface tracking, animmersed-boundary module for movable structure motion, and a sediment transport module for suspended sedimentmotion and seabed profile evolution. Numerical results show that the drift of the container under wind action can becomputed using the FSSM regardless of the weight of the shipping container and the wind speed, and demonstrate thecomputational capability of the FSSM to simulate wind-induced drift behavior of floating containers.

の考察を述べ，最後に本研究の結論を第 5 章にまとめる．

１．緒言2011 年東北地方太平洋沖地震に伴い発生した大津波により，仙台港に蔵置されていた約4,400個のコンテナのうち約2,000個が海上へ流出した．それらのうち，約 340 個は仙台港内で沈没し，約200 個は石巻市から岩沼市にかけて打ち上げられたことが確認されている(1)(2)．このように，沈没したり打ち上げられたりしたコン
テナは，最初のうちは津波の余波により漂流し，その後は風とそれに伴う風波により運ばれたと推測される．
一方，巨大地震の発生時には，港湾は物資輸送の拠点や救援活動の基地として期待されている．しかし，被災地の支援のために入出港する船舶にとって海上を漂流するコンテナは航行の障害となることから，流出したコンテナの早期の発見と回収が課題である．その際に，コンテナの位置の正確な予測ができれば，コンテ
ナのより効率的な回収が可能になると考えられるものの，そのためには流出したコンテナの津波による漂流だけではなく，風とそれに伴う風波による漂流と，その後の海底への沈没や海岸への打ち上げまでを包括的に取り扱える予測手法の確立が不可欠である．
それらのうち，津波が治まってからのコンテナの漂流挙動に対して支配的と考えられる風とそれに伴う風波の影響については，コンテナに作用する力として風と流れによる抗力を考慮した数値計算モデルの構築と，現地でのコンテナの漂流事例および漂流実験との比較による妥当性の検討が Daniel ら(3)によって行われている．
しかし，比較対象とした現地データには複雑に変化する気象や海流の影響が含まれていることから，熊谷・鈴木(4)と熊谷・根木(5)は理想的な条件下での検討を行うために風洞水槽を用いた水理実
験を行っているものの，風とそれに伴う風波によるコンテナの漂流挙動の解明とそれに基づく予測手法の確立には到っていない．
本研究では，水谷ら(6)による水理実験の風洞水槽をモデル化した計算領域に，中村ら(7)と同様のスケールでモデル化したコンテ
ナの模型を浮かべ，そのコンテナに風を作用させたときの漂流挙
動を流体・構造・地形変化の 3 者間の相互作用を解析できる 3 次
元流体・構造・地形変化連成数値計算モデル（Three-Dimensional
Coupled Fluid-Structure-Sediment Interaction Model；以下，FSSM と略す）(8)(9)により解析する．そして，風によるコンテナの漂流現象に対する FSSM の適用性を検討する．以下，第 2 章では FSSM の詳細を，第 3 章では計算条件の詳細を，第 4 章では計算結果とそ

２．3 次元流体・構造・地形変化連成数値計算モデル
3 次元流体・構造・地形変化連成数値計算モデル（FSSM）(8)(9)の計算領域を Fig. 1 に示す．FSSM は，メインソルバーと 3 つのモジュールから構成されている．メインソルバーは，透過性材料の間隙内部の流体を含む非圧縮性粘性流体からなる全気液相に適用できるように一般化し，さらに可動構造物や地形変化を取り扱えるように拡張した連続式と Navier-Stokes（NS）方程式に基づく
Large-Eddy Simulation（LES）である（Extended NS Solver；以下，
XNS と略す）．
その XNS には，気液界面を追跡する Volume of Fluid（VOF）法(10)を改良した Multi-Interface Advection and Reconstruction
Solver（MARS）(11)に基づく VOF モジュール（VOF Module；以下，
VFM と略す），可動構造物を取り扱う体積力型埋め込み境界（Immersed Boundary；IB）法(12)に基づく IB モジュール（IBModule；以下，IBM と略す）
，渦と浸透滲出流の両者の効果(9)を導入した上で，浮遊砂の巻き上げ・移流拡散・沈降の各輸送過程を考慮して浮遊砂の解析を行い，掃流砂と浮遊砂の輸送に伴う地形変化の計算を行う底質輸送モジュール（Sediment Transport
Module；以下，STM と略す）が，XNS との相互作用を考慮できるように Two-Way カップリング手法により組み込まれている．本章では，XNS，VFM，IBM，STM，カップリング手法の説明を行うとともに，パラメータの設定値について述べる．
２．
１ 拡張Navier-Stokesソルバー（XNS），
VOF モジュール（VFM），
埋め込み境界モジュール（IBM）
流体の運動は XNS で，気液界面の追跡は VFM で，可動構造物の運動は IBM で取り扱う．以下，XNS，VFM，IBM の支配方程式と数値計算スキームの説明を行う．
２．１．１ 支配方程式
VFM で求めた気液界面の位置，IBM で求めた可動構造物の運動，
STMで求めた地形変化を考慮して流体の運動を計算するXNSと，気液界面を追跡する VFM の支配方程式をまず説明する．
LES では，Grid Scale（GS）と呼ばれる大きなスケールの渦を直接計算し，普遍性が期待できる Subgrid Scale（SGS）と呼ばれ1

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1
1  F   a ；  w ：水の粘性係数； a ：空気の粘性係数）， C A は
透過性材料の付加質量係数(13)，Ri は透過性材料による線形・非線
形抵抗力ベクトル(13)， f i s は CSF（Continuum Surface Force）モデル (14) に 基 づ く 表 面 張 力ベ クト ル ，  ij は DTM （ Dynamic
Two-Parameter Mixed Model）(15)に基づく乱流応力テンソル，Dij はひずみ速度テンソルの GS 成分（   vi x j   v j xi ）， q* は単位時間当たりの造波ソース強度(16)，Qi は造波ソースベクトル，i
T
は減衰領域での減衰関数ベクトル（  0 0   w ；  ：減衰関数）であり，上付きの a はテンソルの非等方成分を表す．ただし，式(6)の導出に際して空隙率の空間勾配は考慮しているものの，空隙率が空間的に変化する場所で，式(6)右辺第（ m xi  0 ）1 項の圧力勾配項と第 2 項の重力加速度項が静水状態で釣り合わない現象が生じてしまうことから，CADMAS-SURF(17)に倣って圧力勾配項については m xi  0 とした．また，式(6)において  を
表面張力係数， を気液界面の曲率， avg を気液界面での流体の密度（    w   a  2 ），CD 2 ，CD1 をそれぞれ線形，非線形抵抗
力係数，ˆ を流体の動粘性係数（  ˆ ˆ  F w  1  F  a ； w ：水の動粘性係数； a ：空気の動粘性係数）， d50 を透過性材料の中央粒径とすると， f i s ， Ri ， Qi はそれぞれ

Fig. 1 Typical computational domain of the FSSM.
る小さなスケールの渦をモデル化する．実流速ベクトルを vif ，上付きの－をグリッドフィルタとすると，
vif  vif  vif 
T

(1)

v
w  と SGS 成分 vi  に分解でのように GS 成分 vi  uきる．また，VOF 法において VOF 関数 F （ 0  F  1 ； F  0 ：空気；0  F  1 ：気液界面； F  1 ：水）が各計算格子での流体率を表しているように，体積力型 IB 法(12)では可動構造物を各計
算格子での体積率 Fob （ 0  Fob  1 ； Fob  0 ：可動構造物外部；
0  Fob  1：可動構造物表面；Fob  1：可動構造物内部）で表す．
そして，可動構造物の内部では可動構造物の速度になるように，
Fob に基づいて相互作用力ベクトル f i ob を作用させる．ここで，ob
ob
可動構造物の速度ベクトルを viob （  uiob  eijk ob
j rk ； ui ：可動ob
構造物の並進速度ベクトル；i ：可動構造物の角速度ベクトル；
riob ：可動構造物の重心からの相対位置ベクトル；eijk ：置換記号）
Tとしたとき，体積平均速度ベクトル vi （  u v w ）を
f

f

f

f

f

vi  vi  vi  1  Fob  vif  Fob viob

Ri  

vi  1  Fob  vif (4)  q*x j

m  C 1  m  vt iA mvi v j









 

(12)





(13)

 とグリッドフィルタ幅と表される．ここで，テストフィルタ幅  の比を  （    ；上付きの～：テストフィルタ）とするとき，Germano Identity Lij ， H ij ， R ij ， Mij をそれぞれ(5)Lij  vi v j  vi v j

(14)

 
H ij  vi v j  vi v j

(15)





(6)

 





 
R ij  vi  vi v j  v j  vi  vi v j  v j

(16)

 
Mij   2 D Dij  D Dij

(17)



F  mv j F

 Fq*
t
x j

(11)

LRij  vi  vi v j  v j  vi  vi v j  v j
1 
2mˆ Dij  m ija   Qi  miˆ x j
x j

m

(10)

Lmij  vi v j  vi v j

x j



v j 
2  
 mˆ

x j 
3 xi ただし， Lmij と LRij はそれぞれ



m p
m

 mgi   f i s  Ri  f i ob ˆ xi
ˆ


(9)

Ra
 ija  Lma
ij  CB Lij  CS D Dij

である．また，体積空隙率と面積空隙率が等しく，その空隙率 mの時間変化（ m t  0 ；t ：時間）が十分に小さいと仮定する．
以上より，直接計算できる GS 成分 vi を vi から分離するために，拡張した連続式，NS 方程式，VOF 関数 F の移流方程式にグリッドフィルタを作用させると，それぞれ次式が得られる．

 mv j

12CD 2 ˆ 1  m C ˆ 1  m 
vi  D1vi v j v j
2md50
md 502と表される．また，修正 Leonard 応力テンソルを Lmij ，ひずみ速度テンソルの GS 成分 Dij の大きさを D とすると， ija は Horiuti(15)より次式で与えられる．

と定義する．ここで， vi と vi はそれぞれ vi の GS 成分（ Tu v w ）と SGS 成分であり， vi と vi はそれぞれ



(3)

(8)

Qi  vi q* (2)vi  1  Fob  vif  Fob viob

F ˆ
xi  avg

f i s  とすると，係数 CB ， CS はそれぞれ
CB 

(7)

ここで， p は圧力の GS 成分， xi は位置ベクトル（ TT x y z  ）， gi は重力加速度ベクトル（  0 0  g  ； g ：重力加速度）， ̂ は流体の密度（  F  w  1  F   a ；  w ：水の密度；  a ：空気の密度）， ̂ は流体の粘性係数（  F  w 

 L - H R M M   L - H  M R M

CS  

a
ij

a
ij

a
ij

kl

kl

a
ij

a
ij

ij

a
kl

kl

R ijaR ijaM kl M kl  R ijaMijR kla M kl

 L - H  M R R   L - H R R M
a
ij

a
ij

ij

a
kl

a
kl

a
ij

a
ij

a
ij

R ijaR ijaM klM kl  R ijaM ijR kla M kl

a
kl

kl

(18)

(19)

と表される．ただし， CS が負になることで計算が不安定になら2

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1（Simplified Marker and Cell）法(23)を採用した．ただし，式(6)右辺
第 1 項の圧力勾配項と第 2 項の重力加速度項については梶島(24)を参考に 1 次精度 Euler 陽解法を，式(9)右辺第 1 項の線形抵抗力項については安定化のために 1 次精度 Euler 陰解法を，その他の項については 3 次精度 Adams-Bashforth 法を適用した．すなわち，
タイムステップ数を上付きの n ，第 n ステップから第  n  1 ス
テップまでの時間の刻み幅を t n 1/ 2 とすると，予測段階の第 1 ス
テップ，予測段階の第 2 ステップ，修正段階はそれぞれないように，Horiuti(15)に倣って CS を Clipping した．また，減衰関
数  は Hinatsu(18)と Cruz ら(19)を参考に


C

g

   h0  x 
 z

sin 
,1 in artificial
 min 
2   h    damping zones (20)otherwiseとした．ただし，C は減衰係数，h は静水深， は水位変動，
は減衰領域の長さ， x は減衰領域の開始位置からの距離， z は
．なお，上述したように構造物の底面からの距離である（Fig. 2）
内部では vi が viob となるように強制することから，岩田ら(20)が述
べているように構造物の内部での ̂ やˆ は任意であり，計算結果に与える影響は十分に小さいことを確認している．
続いて，可動構造物の運動を取り扱う IBM の支配方程式を説明する．可動構造物に作用する外力ベクトルを f i ex ，底面との接触による反力ベクトルを f i re ，可動構造物の質量を mob ，可動構
造物の重心回りの慣性テンソルを I ijob とすると，可動構造物の並進速度ベクトル uiob に関する運動量保存則と角速度ベクトル iob
に関する角運動量保存則はそれぞれ

duiobm f i ex  fi redtob

I ijobd  objdt

 eijk rj  f kex  f kre 

vi

p1

vi

vi

p2

n 1

 vi

 vi

p2

(24)



t n 1 2
m n ob nfiBnnm  C A 1  m  ˆ n

(25)



t n 1 2
m n  n 1 2
m n  C A 1  m n  ˆ n xi

(26)

p1

(21)

n

p1

Bn

p2

と書ける．ただし， vi と vi は vi の予測値， n 1 2 は圧力増分n 1n（  p  p ）であり，  n 1 2 は式(26)の両辺に  xi を作用させた次に示す Poisson 方程式から求められる．

(22)


 
mn
m n  n 1 2
 n
Bn    nnxi  m  C A 1  m  ˆ
xi


p2
n
 m vixi  q*n 1t n1 2と書かれる．ここで，Xiao ら(21)に倣って可動構造物の表面に作用するせん断力の影響が十分に小さいと仮定すると， f i ex は p ob g   xi y j zk
f ex   Fob  xi , j ,k n m n  p n
t n 1 2
 vi  n
 n
 m n gi   
n 
ˆ
x


m
C
m


1
 
i
A

2
n 1 2 n 1 2
t


t

n
n

A1i 
A2ni   B n
 A0 i 

 
2
6

 

(23)

と表される．ただし，ob は可動構造物の密度，xi ，y j ，zk
はそれぞれ x，y，z 軸方向の格子幅である．また，可動構造物が底面と接触した場合には，本多ら(22)に倣って，それぞれの接触点での相対速度の反力方向成分が接触後にゼロになると仮定して
f i re を求めた．ただし，例えば可動構造物の自由度が 6 のとき，可動構造物を構成している平面のうちの 1 つの面に対して，接触点が 4 点以上あるか，または同一直線上の接触点が 3 点以上ある
場合には， f i re が一意に定まらず求められないことから，4 点目
以降の接触点，または同一直線上の場合には 3 点目以降の接触点については，それらの接触点での f i re が底面へのめり込み量に比
例すると仮定して計算を行った．

(27)





ここで，上述したように構造物の内部では vi が viob になるように強制することから，相互作用力ベクトル f i ob n は
fi



ob n



Fobn viob n  vi

p1



(28)

t n1/ 2mn
m n  C A 1  m n  ˆ nBnとなる．したがって，式(25)と式(28)より，予測段階の第 2 ステッ
プは最終的に次のように表される．

vip2

 1  Fobn  vi

p1

 Fobn viob n(29)

なお，3 次精度 Adams-Bashforth 法の導入に伴う A0ni ， A1ni ， A2ni と
1 次精度 Euler 陰解法に導入に伴う B n はそれぞれ

A 

２．１．２ 数値計算スキーム
２．１．２．１ 拡張 Navier-Stokes ソルバー（XNS）
XNS では，式(5)と式(6)で表される vi と p の連成計算に SMACn

 m n vi v j

n
0i

n

  m f  C 1  m  v v v  
n

n

ˆ

x j

n

n

i

2d 50





n

D1

sn

i

j

n

j

(30)

n
1 

2m nˆ n Dij  n mn ijan   Qin  mn inˆ x jx j

A1ni  t n  3 2  t n  3 2  2t n 1 2  A0ni   t n  3 2  t n 1 2  A0ni12  t n 1 2  A0ni 2
2



 t

n 3 2t n 1 2  t n  3 2  t n 1 2 

A2ni  2 t n  3 2 A0ni   t n  3 2  t n 1 2  A0ni1  t n 1 2 A0ni 22

tn 3 2

Fig. 2 Artificial damping zone in a typical computational domain.

3 (32)

 t n  1 2   t n  3 2  t n  1 2 

12CD 2ˆ 1  m 
t n 1 2nn
d 502
m  C A 1  m nBn  1 

(31)

n

(33)

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1
明を行うとともに，数値計算スキームを説明する．

と表される．一方，空間方向の離散化には，式(30)右辺第 1 項の対流項については 5 次精度 MUSCL（Monotone Upstream-Centered
Schemes for Conservation Laws）(25)を，他の項については 2 次精度
中央差分法を採用した．また，Poisson 方程式の解法には MICCG（Modified Incomplete Cholesky Conjugate Gradient）法を適用した．
以上をまとめると，XNS での計算の手順は次のようになる．
p1
1. 式(24)より vi を求める．
p2
2. 式(29)より vi を求める．
n 1
3. 式 (27) よ り  n 1 2 を 求 め ， 求 め た  n 1 2 よ り pn（  p   n 1 2 ）を更新する．
n 1
4. 式(26)より vi を求める．

２．２．１ 底質の質量保存式
底質の質量保存式は次式で表される．


1  qx q yzs pN  qzbS   0

t 1  m  xy

ここで， zs は基準面からの底質表面の高さ， qx と q y はそれぞれ
x ， y 軸方向の掃流砂量， pN は底質の巻き上げ関数， qzbS は底質表面での z 軸方向の沈降による浮遊砂輸送フラックスである．
２．２．２ 掃流砂輸送モデル
掃流砂輸送モデルでは，漂砂現象に与える渦と浸透滲出流の影響を考慮するために，掃流砂量 qi に対して渦の効果が，限界Shields 数 *c と掃流砂量 qi に対して浸透滲出流による層流抵抗と乱流抵抗の効果が組み込まれている(9)．
底質粒子として，非粘着性で均一の粒子径の球体を考える．また，底面流速および浸透滲出流速 w（滲出を正とする）による慣性力と，限界 Shields 数を上回る状況では比較的大きい底面流速が生じていると考えられることから，底面流速による層流抵抗の影
響は無視する．このとき，底質粒子 1 個に作用する力は
1. 底質粒子の水中重量 W

２．１．２．２ VOF モジュール（VFM）
VFM では，式(7)で表される VOF 関数 F の移流方程式に前述したように MARS(11)を適用した．なお，隣接する格子の境界における F の不連続性をできるかぎり押さえながら気液界面の輸送を行うために，1 タイムステップでの界面輸送を全方向に一度に計算する Unsplit Time Integration を時間方向に多段階で行う
Fractional Step Method(26)を採用した．
２．１．２．３ 埋め込み境界モジュール（IBM）
IBM では，式(21)と式(22)に 1 次精度 Euler 陽解法を適用し，可動構造物の並進速度ベクトル uiob と角速度ベクトル iob をそれぞれ次式から求める．

uiob n 1  uiob n n 1/ 2

t
mob

f

ex n
i

 fi re n t n 1/ 2n 1n  ob ob eijk rjn  f kex n  f kre n 
 ob
j
j
I ij

1
W     s   w  gd503
6

iob n 1  iob n 

(39)

(34)

2. 底面流速による乱流抵抗力 FD と揚力 FL
(40)

(35)

1
FD   CD1 w d 502 vr2  CFD vr2
81
FL   CL  w d502 vr2  k LCFD vr2
8(41)

そして，可動構造物の重心位置ベクトル xiob と回転角ベクトルiob を，2 次精度 Crank-Nicolson 法に基づいて，それぞれ

xiob n 1  xiob n (38)

3. 浸透滲出流速 w による作用力 Fw

t n 1/ 2 ob nui  uiob n 1 
2

(36)

t n 1/ 2 ob ni  iob n 1 
2

(37)

1
Fw   wCD1 w d 502 w w  3 wCD 2  w wd 50 w
8(42)から構成される．ここで，vr は底質粒子の位置での vb に対する相対底面流速（vb：底質粒子の平均輸送速度），CL は揚力係数，kL
は CL と CD1 の比（= CL / CD1）
，  s は底質粒子の密度であり，により更新する．なお，求められた可動構造物の位置に基づく Fobの更新には，Yuki ら(27)の双曲線正接関数を用いる手法とは異なり，各計算格子を各方向に数分割ずつした小格子を考え，その分割した小格子の中央から可動構造物を構成する面までの垂線の長さを基に Fob を算出することで可動構造物の総体積の増減を抑えた．
以上をまとめると，IBM での計算の手順は次のようになる．
1. f i ex （式(23)）と f i re を求める．
2. 式(34)より uiob を，式(35)より iob を求める．
3. 式(36)より xiob を，式(37)よりiob を求める．
4. 求めた xiob とiob より， Fob を更新する．

1
CFD   CD1  wd 502
8(43)である．また，式(40)に示した FD のは遮蔽係数であり，式(41)に示した FL のと等しいと仮定している．さらに，その類推から底質内部にある粒子によってもある程度遮蔽されることで 1 個の底質粒子のみが存在すると考えた場合よりも抵抗力が減少すると考え，式(42)に示した Fw にも遮蔽係数w を導入している．
２．２．２．１ 限界 Shields 数
浸透滲出流速による作用力 Fw が底質粒子の水中重量 W を下回っているとき，移動限界時の底質粒子には Fig. 3 に示す外力が作用している．ここで，摩擦速度を vf，その係数を Cvf，斜面下向きに対する vf の方向を，vb に対する vf の方向を，vb に対する vrの方向を1 としたとき，同図中の FDc，FLc，vfc，vbc，vrc，c，c，
1c は移動限界時の FD，FL，vf，vb，vr，，，1 を表す．また，
は斜面の傾斜角，s は底質粒子の最大静止摩擦係数（  tan  s ； s ：最大静止摩擦角）である．このとき，移動限界時より vbc = 0
つまり vrc = Cvf vfc かつc = 1c であり，vbc と同方向および直交方向

２．２ 底質輸送モジュール（STM）
STM は，
掃流砂と浮遊砂の輸送に伴う地形変化の計算を行う底質の質量保存式，Roulund ら(28)を参考に構築した地形変化の計算に必要な掃流砂量を求める掃流砂輸送モデル，地形変化の計算に必要な底質の巻き上げと沈降による浮遊砂輸送フラックスを求めるとともに，浮遊砂の巻き上げ・移流拡散・沈降の各輸送過程を考慮して浮遊砂の濃度分布を計算する浮遊砂輸送モデル，Roulundら(28)を参考に構築した地形変化の計算後に底質斜面の崩壊を計算する斜面崩壊モデルから構成されている．以下，それぞれの説

4

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1の力のつり合いより次式が得られる．
vrc2 

 cos   sin  cos     W  Fsc cos  c   s kL  CFD

c



w

* 
(44)
vsurf

(45)

vf

 *c  s  1 gd 50

v 2fc 0
 *c 0 s  1 gd50

2
rc 0

v
Cvf2  s  1 gd 50

z   0.9


2
rc
2
rc 0

v
v

 sW

1   s kL  CFD

(46)

(47)

(48)

(49)


s


s

(54)

vi  vi  CSGS vi  vi  CSGS (vi  vi )ここで， pEF は掃流形態で移動する底質粒子の割合，vbi は vb の i
方向成分であり，底質粒子の動摩擦係数を d （  tan  d ； d ：動摩擦角），Shields 数を * とすると， pEF は(55)と書ける．また，Multi-Level Filtered Model(31)に基づいてより高次の Transfer Field を用いると，3 次と 5 次の近似値 vi はそれぞれvi  vi  CSGS {2(vi  vi )  (vi  vi )}(50)

vi  vi  CSGS {4(vi  vi )  6(vi  vi )と与えられる，ただし， d   s である．また， * はFig. 3 Critical force on a sediment particle
on a sloping bed.

sとなる．ここで，左辺 vi は SGS 成分のうち比較的大きなスケー
ルの変動を表している．また，Transfer Field と呼ばれる右辺(vi  vi ) は，GS 成分 vi と GS 成分から小スケールの渦を消去した
成分 vi の差であることから，GS 成分のうち比較的小さなスケー
ルの変動を表しており，XNS で直接計算可能な GS 成分 vi から求められる．したがって，式(54)は SGS 成分のうち比較的大きなス
ケールの変動と，GS 成分のうち比較的小さなスケールの変動が相似であることを表していることから，SGS 成分のうち比較的大きなスケールの変動を vi ，すなわち Transfer Field で近似できるとし，また SGS 成分の残りの変動については Transfer Field に比例すると仮定し，その比例定数を CSGS とする．以上より，GS 成分 viに SGS 成分 vi の影響を考慮した vi の近似値 vi は２．２．２．２ 掃流砂量
単位幅，単位時間当たりの i 方向の掃流砂量 qi は，Engelund・Fredsøe(29)による次式から求められる．

if  *   *c0

pEF   6
   *   *c  if  *   *c
 d


s

vi  vi  vi

以上より， *c は Fw  W のとき式(44)から式(47)より，Fw > W のとき式(48)より求められる．

1
qi   d 50 pEF vbi
6 k  k exp   k 6   4.535  k  2000 (53)

漂砂現象に与える渦の影響を考慮しないとき，式(52)の vsurf はXNS で求めた vi から計算される．一方，漂砂は本来 vi にグリッドフィルタを施した GS 成分 vi ではなく，あらゆるスケールの渦を含む vi の変化に対応して生じると考えられることから，直接解像できない SGS 成分 vi の影響も無視できないと考えられる．
Scale Similarity Model(30)によれば，式(2)の両辺にグリッドフィルタを作用させて整理すると，

一方，Fw が W を上回っているとき，底質粒子に静止摩擦力が作用しなくなることから，次式が得られる．

 *c
0
 *c 0(52)

ただし，vsurf は底面から微小距離 z  （  z surf v f  w ）だけ上方での接線方向流速の大きさ，   は Kármán 定数， A は van Driest 型
減衰係数であり， ks を粗度 Reynolds 数（  k s v f  w ； k s ：等価砂粗度）とすると， z  は次式で与えられる．

と表される．ただし，vrc2 0 は w = 0 かつ = 0 での vrc2 であり，式(44)に Fw =  = c = 1c = 0 を代入した次式から求められる．
vrc2 0 

0

ただし，    s すなわち  r   r   s である．なお，  r と  r は後述する．ここで，未知数は vrc2 とc であることから，繰り返し計算により vrc2 とc が求められる．したがって，浸透滲出流を考慮しないとき（w = 0）の水平面（ = 0）での限界 Shields 数を *c 0 ，底質粒子の比重を s（   s  w ）とすると， *c は

vrc2
C  s  1 gd 50z

 2  1 1 

22
 1  4 2  z   z   1  exp   z   z   A   dz  （ cは cと同符号かつ  c   c ）

2
vf

(51)

 s  1 gd50により定義され，その v f は次式から求められる．

W  Fw  sin  sin  c c  arctan 

2


C
v
W
F
sin
cos



w
c
 FD rc 

v 2fc

v 2f

(56)

(57)

4(vi  vi )  (vi  vi )}(b) Fw > W(a) Fw  W
Fig. 4 Force on a sediment particle in bed-load motion on a sloping bed.
5

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1となる．したがって，GS 成分 vi の代わりに式(55)から式(57)のいずれかから求めた vi を用いて式(52)の vsurf を算出することで，
Shields 数 * に対して渦の影響の考慮を行っている．
続いて，式(49)の vbi を求めるために，Fig. 4 に示す掃流形態で移動する底質粒子を考える．まず Fig. 4(a)に示した Fw が W を下回っているとき，vb と同方向および直交方向の力のつり合いより
FD cos 1  W  Fw  sin  cos     d W  Fw  cos   FL (58)

FD sin 1  W  Fw  sin  sin    (59)

qiS  mwsiCと表される．ただし，  s は浮遊砂の乱流拡散係数， wsi は底質粒
T子の沈降速度ベクトル（   0 0  ws  ； ws ：底質粒子の沈降速度）であり， ws は次に示す Rubey の式により与えられる．
 2 3636 
ws  sgd 50 



 3 d
d* 
*


(60)

vr cos 1  vb  Cv f v f cos 

(61)

0

1.5
pN     *   *c   s  10.6 g 0.6 d 500.8C p w0.2   *c である．ここで，式(58)から式(61)に対して未知数は vb，vr，，1 の 4 つであることから，繰り返し計算により vb とその方向，すなわち vbi が求められる．したがって，前述した Fw  W での *cとあわせて，最終的に式(49)と式(50)より qi が求められる．
一方，Fw が W を上回っているとき，底質粒子に動摩擦力が作用しなくなることから，Fig. 4(b)に示すように底質粒子に作用する力は FD と(Fw – W) sin となり，次式が成立する．
1     Fw  W  sin 
CFD

(63)

となる．また， Cv f v f ， vb ， vr の幾何学的関係より，

(64)

2


2  vr cos   Cvf v f 1   arccos 
 1
2
2
2
2
 vr sin    vr cos   Cvf v f （はと逆符号）(65)

vbslide 

qislide q  mviC

qiD   m s

C
xi


3

d 50

1 d 50 s  1 g  sin   d cos  
3 CD1

slide
z s
1  qxslide q y 



0
t 1  m  x
y 

(66)

(72)

(73)

(74)

により表されることから，最終的に式(73)と式(74)より斜面崩壊に伴う地形変化が計算できる．

ここで，C は浮遊砂濃度，qiC ，qiD ，qiS はそれぞれ移流，拡散，沈降による浮遊砂輸送フラックスであり，それぞれCi

4 d 50 s  1 g  sin   d cos  
3 CD1

また，斜面崩壊に伴う底質の質量保存式は，斜面崩壊時には浮遊
砂の巻き上げや沈降が生じないと仮定すると，式(38)に pN 
qzbS  0 を代入した

２．２．３ 浮遊砂輸送モデル
浮遊砂輸送モデルでは，榊山ら(32)や中山ら(33)を参考に，透過性材料の間隙内部にも適用できるように一般化した次式に示す浮遊砂の質量保存式を支配方程式として用いている．
C qCj  q Dj  q Sj   0
t x j(71)

ここで，    d すなわち r   r   d である．したがって，上式と pEF  1 を式(49)に代入すると，最終的に斜面崩壊に伴う掃流砂
量 qislide は次のように求められる．

となる．したがって，上式より vb とその方向が求められることから，前述した Fw > W での *c とあわせて，最終的に式(49)と式(50)より qi が求められる．
以上より，掃流砂量 qi は vbi と前述した *c から求められ，そのvbi は Fw  W のとき式(58)から式(61)より，Fw > W のとき式(62)から式(65)より求められる．

mif  *   *c

２．２．４ 斜面崩壊モデル
Roulund ら(28)は，斜面の傾斜角  が底質粒子の安息角 r より微小角度 r だけ上回ったときに斜面が崩壊し始め， r より微小角度  r だけ下回るまで崩壊し続けることを水理実験で確認している．この現象を基に，Roulund ら(28)は以下の手法を提案している．
斜面の傾斜角  が  r   r を上回ったとき，その底質粒子（ pEF  1 ）が静水中（ v f  0 ，  1  180 ）で一定速度 vb で
斜面下向き（   180 ）に崩壊すると仮定する．また，斜面崩壊時における FL と Fw の影響は十分小さいと仮定する．このとき，式(58)から式(61)に，式(39)，式(40)， FL  Fw  v f  0 ，   
1  180 を代入して整理すると，斜面崩壊に伴う底質粒子の平均輸送速度 vbslide は次式で与えられる．

(62)

vb  Cvf v f cos   vr cos 1

if  *   *cを用いた(35)．ここで， C p は底質の巻き上げに関する無次元パラメータである．なお， * を介して渦の影響が， *c を介して浸透滲出流の影響が考慮されている．

このとき，底面接線方向の力のつり合いより
vr 

(70)

ここで，d* は底質粒子の無次元中央粒径（  sgd502  w2 ）である．
なお，ここでは透過性材料として内部への浮遊砂の侵入が考えられる捨石や消波ブロックを想定している．
底質表面の境界条件として，Nielsen ら(34)に倣って，底質表面での z 軸方向の拡散による浮遊砂輸送フラックス qzbD が底質の巻き上げ関数 pN に等しいとする条件を用い，その pN としてが得られる．また， Cv f v f ， vb ， vr の幾何学的関係より，
vr sin 1  Cv f v f sin (69)

２．２．５ 数値計算スキーム
STM での計算の手順は次の通りである．
1. 式(67)より qiC を，式(68)より qiD を，式(69)より qiS を求める．
なお，式(67)と式(69)の離散化には安定性を考慮して 1 次精度風上差分法を，式(68)の離散化には 2 次精度中央差分法を
適用した．

(67)
(68)
6

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1

2. 求めた qiS より qzbS を，式(71)より pN を求める．
3. 求めた qzbS と pN を底質表面での境界条件として，式(66)を1 次精度 Euler 陽解法により離散化した次式より C を求める．
C n 1  C n t n 1/ 2  qCj n  q Dj n  q Sj n 
m n x j(75)

4. 式(49)より qi を求める．なお，式(58)から式(61)の計算には，
Roulund ら(28)に倣って Newton-Raphson 法を適用した．
5. 式(38)を 1 次精度 Euler 陽解法により離散化した次式より zsを更新する．
zsn 1  zsn 6.

nt n 1/ 2  qxn q y pNn  q zbS n n y1  m  x

(76)Fig. 5 Two-way coupling procedure.

 が r   を上回ったとき，以下の斜面崩壊計算を行う．
なお，ここでは Roulund ら(28)に倣って斜面崩壊が瞬間的に起こると仮定し，手順 1 から 5 までとは独立したタイムス
テップ数 h と時間の刻み幅 t h 1/ 2 を使って計算を進める．
i. 式(73)より斜面崩壊に伴う qislide を求める．
ii. 式(74)を 1 次精度 Euler 陽解法により離散化した次式よりzs を更新する．
rt h 1/ 2  qxslide h q yy
1  m n  x

slide hzsh 1  zsh Table 1 Parameters of the XNS.
フィルタ幅比2.0
重力加速度 g9.81 m/s2
水の密度  w
9.97×102 kg/m3
空気の密度  a
1.18 kg/m3
水の動粘性係数 w
8.93×10-7 m2/s
空気の動粘性係数 a1.54×10-5 m2/s
表面張力係数
7.20×10-2 N/m
付加質量係数 C A-0.04
非線形抵抗力係数 CD1 0.45
線形抵抗力係数 CD 225.0
減衰係数 C3.0(77)iii. 斜面の傾斜角  がすべての場所で r   r を下回るまで，
手順 i から ii を繰り返す．

Table 2 Parameters of the STM.
摩擦速度の係数 Cv f
10.0
底質巻き上げ係数 C p
0.00033
揚力係数と非線形抵抗力係数の比 kL 0.85
底面流速の遮蔽係数
0.4
浸透滲出流速の遮蔽係数w
0.4
水平面での限界 Shields 数 *c 00.05
最大静止摩擦係数  s0.63
動摩擦係数 d0.51
0.4
Kármán 定数 
25.0
Van Driest 型減衰係数 A
底質粒子の水中安息角 r
32.0º
斜面崩壊開始基準 r
2.0º
斜面崩壊完了基準 r
2.0º

２．３ XNS・VFM・IBM・STM 間のカップリング手法
FSSM では，XNS，VFM，STM，IBM の接続に Fig. 5 に示す
Two-Way カップリング手法を採用した．具体的には，FSSM での計算の手順は次のようになる．
1. XNS を実行する．
2. XNS で得られた流速場を用いて VFM を実行する．
3. 流速場に加えて，VFM で得られた気液界面の位置を用いてSTM を実行する．
4. 流速場と気液界面の位置に加えて，STM で得られた地形変化を用いて IBM を実行する．
5. VFM，STM，IBM で得られた値をフィードバックして XNSを実行する．
6. 手順 2 から 5 を計算終了時刻まで繰り返す．
２．４ パラメータ
Table 1 と 2 に，XNS と STM で用いたパラメータをそれぞれ示
す．Table 1 に示した XNS に関して，LES のグリッドフィルタとテストフィルタには Horiuti(15)に倣ってトップハットフィルタ（ボックスフィルタ）を採用し，グリッドスケールのフィルタ幅  は各位置における各方向の格子幅と等しいとした．また，フィルタ
幅比  は Garmano ら(36)に倣って設定した．その他，付加質量係数
C A ，非線形抵抗力係数 CD1 ，線形抵抗力係数 CD 2 は水谷ら(13)の実験結果に基づいて，減衰係数 C は中村ら(37)による試行計算の結果に基づいて設定した．Table 2 に示した STM については，摩擦速度の係数 Cvf ，水平面での限界 Shields 数 *c 0 ，最大静止摩擦係数  s ，動摩擦係数 d ，van Driest 型減衰係数 A ，底質粒子の水中安息角  r ，斜面崩壊開始基準  r ，斜面崩壊完了基準  r はRoulund ら(28)に基づいて，底質巻き上げ係数 C p は Nielsen(35)に基づいて，揚力係数と非線形抵抗力係数の比 kL と底面流速の遮蔽係
数は中川・辻本(38)に基づいて，浸透滲出流速の遮蔽係数w は中村・水谷(9)による試行計算の結果に基づいて設定した．

３．計算条件
風洞水槽を用いた水谷ら(6)による水理実験を基に計算領域を決定し，また津波によるコンテナの漂流特性の検討を行った中村ら(7)による水理実験を基にコンテナのモデルを決定した．
Fig. 6 に計算領域の概略図を示す．同図に示すように，水谷ら(6)と同様に，計算領域内に反流を分離するための隔壁床を設け，その上部の静水深は 0.45 m，さらにその上部の通風部の高さは 0.25m とした．ただし，その長さは計算負荷の軽減のために水谷ら(6)よりも短い 3.20 m とし，
その風上側と風下側にそれぞれ 1.00 m の減衰領域を設けた．さらに，その風上側には 1.00 m の助走区間を設け，梶島(24)を参考に助走区間内に境界層の発達を促すための障害物（長さ 40 mm，高さ 40 mm）を 3 つ設置した．一方，コンテ
ナは，中村ら(7)と同様に，国際規格に基づく 40 ft コンテナを Froudeの相似則に基づいて 1/75 でモデル化した 163 mm×32 mm×35
mm とし，その質量は貨物半載時と満載時を想定したそれぞれ
37.8 g と 77.6 g とした．なお，中村ら(7)の水理実験ではコンテナの7

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1

Fig. 6 Computational domain.
Fig. 7 Motion of the shipping container in the still water.
重心がコンテナのほぼ中央にくるように模型を製作しているものの，ここでは簡単のために安定な浮体としてのコンテナを対象とするために重心を初期の浮心と一致させた．また，コンテナの長軸を風向と同方向に設置した場合でも，漂流開始とともに風向と直交方向になるまで回転することが確認されている(4)ことから，
ここではコンテナの長軸を風向と直交方向とした場合のみを対象とした．それに伴い，計算領域の奥行きは，コンテナの長軸方向の長さ 163 mm の 3 倍に相当する 0.489 m とした．
計算格子は，y 軸方向に対してコンテナの長軸方向の長さ 163
mm を 1/16 した 10.1875 mm の一様な格子幅y を用いた．また，減衰領域と助走区間を除く範囲の静水面の近傍（0  x  3.20 m，0.04  z  0.04 m）に対して，x 軸方向の格子幅x はコンテナの短
軸方向の長さ 32 mm を 1/8 した 4 mm，z 軸方向の格子幅z はコ
ンテナの高さ 35 mm を 1/7 した 5 mm とした．また，それ以外の範囲については計算負荷のさらなる軽減のために格子幅を徐々に広げ，助走区間ではx = 20 mm，隔壁床の近くではz = 20 mm とした．流速と圧力の境界条件として，底面，上面，側面などの不透過壁の表面に対して No Slip 条件を，下流側境界に対して圧力一定の条件を課し，上流側境界に対しては，前述したように障害in
物で境界層の発達を促すことから，ここでは uairの一様風速を与えた．また，VOF 関数の境界条件には，全境界に勾配ゼロの条件を用いた．なお，本研究では地形変化は取り扱わないことから，
FSSM のうち底質輸送モジュール STM は使用しなかった．
コンテナの初期位置は，前述したようにコンテナの長軸を風向と直交方向，すなわち y 軸方向と一致させるとともに，その重心が減衰領域と助走区間を除く計算領域の中央である x = 1.60 m，y= 0.00 m となるように設定した．また，コンテナの喫水は，静水状態での動揺ができる限り抑えられるように，半載のとき 7.20mm，満載のとき 14.90 mm とした．すなわち，コンテナの重心と浮心の初期の位置が，半載のとき z = - 3.60 mm，満載のとき z = 7.45 mm となるように設定した．この詳細は次章で述べる．一方，コンテナに作用させる風は，計算開始 10 秒後以降に高さ 0.25 mの通風部での断面平均風速が平常時を想定した低風速の 3.0 m/s
inまたは 4.0 m/s となるように uairを設定した．ただし，3 次元計算のために計算負荷が高く，風とそれに伴う風波が十分に安定するまでの長時間の計算は現時点では計算所要時間の観点から現実的ではないこと，また減衰領域と助走区間を除く計算領域が前述したように 3.20 m しかなく，長時間の計算を行うとコンテナが下流側の減衰領域に達してしまうことから，以下では風作用初期における風によるコンテナの漂流挙動に着目して検討を行う．

４．計算結果及び考察
静水状態でのコンテナの動揺を Fig. 7 に示す．
ここで，xxc ，xzc ，c y はそれぞれコンテナの重心位置の水平変位，鉛直変位，y 軸
回りの回転角である．また，正の xxc は風下側への変位を，正のxzc は上向きへの変位を，正の  yc は時計回りの回転を表す．同
図より，計算開始直後に xzc に振動が認められるものの，1 秒後過ぎには小さくなっていることが確認できる．このとき，半載の場合には 0.0025 mm ほど浮き上がっており，満載の場合には 0.018
mm ほど沈み込んでいるものの，コンテナの高さ 35 mm と比較して十分に小さいことが分かる．また， xzc と  yc には変動が認められない．したがって，静水状態でのコンテナの動揺が風による漂流現象に与える影響は十分に小さいと考えられる．
in
コンテナが風下へ最も大きく漂流したケース（半載， uair
= 4.0m/s）におけるコンテナの漂流挙動を Fig. 8 に例示する．また，Fig.
8 と同じ風の条件下でコンテナを満載とした場合の結果を Fig. 9に示す．Fig. 8 と Fig. 9 より，時間の経過とともにコンテナが徐々に風下へ流されていることが確認できる．また，風の作用とコン
テナの漂流の影響を受け，コンテナ周辺の水面が波立ち始めていることが分かる．さらに，Fig. 8 と Fig. 9 を比較すると，半載のコンテナの方が軽く，また風を受ける面積が広いことから，同じ時
刻でも満載の場合と比較して半載の方が漂流した距離が長いこと
inが分かる．図示しないが，コンテナの質量が等しいとき， uairの増加に伴ってもコンテナの漂流距離が伸びることを確認している．
以上より，流体と構造の相互作用を解析できる FSSM により，風とそれに伴うコンテナの漂流を計算できることが示唆される．
５．結論
本研究では，風洞水槽(6)をモデル化した計算領域に 40 ft コンテ
ナのモデル(7)を浮かべ，そのコンテナに風を作用させたときの漂流挙動を，流体・構造・地形変化の 3 者間の相互作用を解析でき
る 3 次元流体・構造・地形変化連成数値計算モデル FSSM(8)(9)を用いて解析した．FSSM は，流体の運動を計算する拡張 NS ソルバ
ー，気液界面の位置を追跡する VOF モジュール，可動構造物の運動を計算する IB モジュール，
浮遊砂濃度と地形変化を計算する
底質輸送モジュールから構成されており，それらの接続には流体・構造・地形変化の相互作用が考慮できるように Two-Way カッ
プリング手法が採用されている．FSSM による解析の結果，コンテナの質量の減少やコンテナに作用する風速の増加とともにコンテナの漂流距離が伸びることを確認した．また，FSSM を用いる
ことで風によるコンテナの漂流現象を計算できることを示した．

8

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1

in= 4.0 m/s.
Fig. 8 Drift of the half-loaded shipping container for uairin
= 4.0 m/s.
Fig. 9 Drift of the full-loaded shipping container for uair

今後は，風によるコンテナの漂流特性の解明に向け，FSSM を用いてさらなる検討を行っていく所存である．
最後に，本研究は公益財団法人 日本科学協会 笹川科学研究助
成による補助を受けたことを付記し，謝意を表する．
参考文献(1) 朝日新聞社, “海水浴場のコンテナ、撤去始まる 津波で仙台港から流出,” http://www.asahi.com/special/10005/TKY20110823
0237.html, 参照 2011-08-28.
(2) 読売新聞社, “漂着コンテナ 撤去本格化,” http://www.yomiuri.
co.jp/e-japan/miyagi/news/20110826-OYT8T00048.htm, 参 照
2011-08-28.
(3) Daniel, P., Jan, G., Cabioc’h, F., Landau, Y., and Loiseau, E., “Driftmodeling of cargo containers,” Spill Science & Technology Bulletin,
7, 5-6 (2002), pp. 279-288.
(4) 熊谷兼太郎, 鈴木 武, “風洞水槽におけるコンテナ模型の流動に関する基礎的実験,” 海洋開発論文集, 26 (2010), pp.
951-956.

9(5) 熊谷兼太郎, 根木貴史, “風洞水槽におけるコンテナ模型の流動に与える風と流れの効果に関する実験,” 土木学会論文集
B3（海洋開発）, 67, 2 (2011), pp. I_154-I_159.
(6) 水谷夏樹, 鍛治允啓, 宮島昌弘, “砕波を伴わない風波の自己調整機能に対する気流の剥離の役割について,” 土木学会論文集 B2（海岸工学）, 67, 2 (2011), pp. I_056-I_060.
(7) 中村友昭, 若松靖弘, 水谷法美, “反射波の影響を考慮した遡上津波によるエプロン上のコンテナの漂流特性に関する研究,” 土木学会論文集 B2 （海岸工学）, 67, 2 (2011), pp.
I_276-I_280.
(8) 中村友昭, 小竹康夫, 水谷法美, “3 次元流体・構造・地形変化連成数値計算モデルを用いた施工時における上部斜面堤の波浪動揺に関する数値シミュレーション,” 第 25 回数値流体
力学シンポジウム, E11-5 (2011), 10 p.
(9) 中村友昭, 水谷法美, “渦と浸透滲出流の影響を考慮した漂砂
計算手法と遡上津波による陸上構造物周辺の洗掘現象への適用に関する研究,” 土木学会論文集 B3（海洋開発）, 68, 1(2012), pp. 12-23.
Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E07-1
schemes for solving the compressible Euler and Navier-Stokes
equations,” Computers Fluids, 22, 2-3 (1993), pp. 259-270.
(26) 秋山 守, 有冨正憲, “新しい気液二相流数値解析－多次元流
動解析－,” コロナ社 (2002), 261 p.
(27) Yuki, Y., Takeuchi, S., and Kajishima, T., “Efficient immersedboundary method for strong interaction problem of arbitrary shape
object with the self-induced flow,” J. Fluid Sci. Tech., JSME, 2, 1
(2007), pp. 1-11.
(28) Roulund, A., Sumer, B. M., Fredsøe, J., and Michelsen, J.,“Numerical and experimental investigation of flow and scour
around a circular pile,” J. Fluid Mech., 534 (2005), pp. 351-401.
(29) Engelund, F. and Fredsøe, J., “A sediment transport model forstraight alluvial channels, Nordic Hydrology, 7 (1976), pp. 293-306.
(30) Bardina, J., “Improved turbulence models based on large eddy
simulation of homogeneous, incompressible turbulent flows,” Ph.D.
dissertation, Stanford University (1983).
(31) Horiuti, K., “Transformation properties of dynamic sub-grid-scalemodels in a frame of reference undergoing rotation,” J. Turbulence,
7, 16 (2006), 27 p.
(32) 榊山 勉, 伊藤 裕, 小林正典, 板井雅之, 小林英次, “傾斜堤を透過する浮遊砂の濃度評価への数値計算モデルの適用,”
海岸工学論文集, 51 (2004), pp. 456-460.
(33) 中山哲嚴, 牧野弘幸, 新井雅之, 小林 学, 佐藤勝弘, “導水中に含まれる浮遊砂量を低減するための海水導入工の構造に関する研究,” 海岸工学論文集, 54 (2007), pp. 661-615.
(34) Nielsen, P., Svendsen, I. A., and Staub, C., “Onshore-offshore
sediment movement on a beach,” Proc., Int. Conf. on Coastal Eng.,
ASCE, 16 (1978), pp. 1475-1492.
(35) Nielsen, P., “Coastal Bottom Boundary Layers and SedimentTransport,” Advanced Series on Ocean Eng. Vol. 4, World Scientific(1992), 340 p.
(36) Germano, M., Poimelli, U., Moin, P., and Cabot, W., “A dynamic
subgrid-scale eddy viscosity model,” Phys. Fluids A, 3, 7 (1991), pp.
1760-1795.
(37) 中村友昭, 水谷法美, 藤間功司, “遡上津波の変形と陸上構造物に作用する津波力に関する 3 次元数値シミュレーション,”
ながれ, 29, 2 (2010), pp. 107-117.
(38) 中川博次, 辻本哲郎, “ 移動床流れの水理,” 技報堂出版(1986), 310 p.

(10) Hirt, C. W. and Nichols, B. D., “Volume of fluid (VOF) method for
dynamics of free boundaries,” J. Comp. Phys., 39 (1981), pp.
201-225.
(11) 功刀資彰, “自由界面を含む多相流の直接数値解析法,” 機械学会論文集 B 編, 63, 609 (1997), pp. 1576-1584.
(12) Kajishima, T., Takiguchi, S., Hamasaki, H., and Miyake, Y.,“Turbulence structure of particle-laden flow in a vertical planechannel due to vortex shedding,” JSME Int. J. Ser. B, Fluids andThermal Eng., 44, 4 (2001), pp. 526-535.
(13) 水谷法美 , 前田健一郎 , Ayman M. Mostafa, William G.
McDougal, “透水性構造物の抵抗係数の評価と波・潜水透水性
構造物の非線形相互作用の数値解析,” 海岸工学論文集, 43(1996), pp. 131-135.
(14) Brackbill, J. U., Kothe, B. D., and Zemach, C., “A continuum
method for modeling surface tension,” J. Comp. Phys., 100 (1992),
pp. 335-354.
(15) Horiuti, K., “A new dynamic two-parameter mixed model forlarge-eddy simulation,” Phys. Fluids, 9, 11 (1997), pp. 3443-3464.
(16) 岩田好一朗, 川崎浩司, 金 度三, “水中構造物による砕波の数値解析,” 海岸工学論文集, 42 (1995), pp. 781-785.
(17) 財団法人沿岸開発技術研究センター , “ 数値波動水路（CADMAS-SURF）の研究・開発,” 沿岸開発技術ライブラリーNo. 12 (2001), 457 p.
(18) Hinatsu, M., “Numerical simulation of unsteady viscous nonlinear
waves using moving grid system fitted on a free surface,” J. KansaiSociety of Naval Architects, 217 (1992), pp. 1-11.
(19) Eric Cruz, 横木裕宗, 磯部雅彦, 渡辺 晃, “非線形波動方程式に対する無反射境界条件について,” 海岸工学論文集, 40(1993), pp. 46-50.
(20) 岩田隆一, 梶島岳夫, 竹内伸太郎, “Immersed-boundary andvolume-of-fluid カップリング法による粒子列の間を通過する
上昇気泡の解析,” 混掃流, 23, 5 (2010), pp. 531-538.
(21) Xiao, F., Yabe, T., Ito, T., and Tajima, M., “An algorithm forsimulating solid objects suspended in stratified flow,” Comp. Phys.
Communications, 102 (1997), pp. 147-160.
(22) 本多和彦, 富田孝史, 西村大司, 坂口 章, “多数の津波漂流物を解析する数値モデルの開発,” 海洋開発論文集, 25 (2009),
pp. 39-44.
(23) Amsden, A. A. and Harlow, F. H., “A simplified MAC technique forincompressible fluid flow calculation,” J. Comp. Phys., 6 (1970), pp.
322-325.
(24) 梶島岳夫, “乱流の数値シミュレーション,” 養賢堂 (1999),
225 p.
(25) Yamamoto, S. and Daiguji, H., “Higher-order-accurate upwind

10

Copyright © 2012 by JSFM