第 25 回数値流体力学シンポジウム
E11-5

3 次元流体・構造・地形変化連成数値計算モデルを用いた
施工時における上部斜面堤の波浪動揺に関する数値シミュレーション
Numerical Simulation of Wave-Induced Motion of Floating Sloping-Top Caissonusing Three-Dimensional Coupled Fluid-Structure-Sediment Interaction Model○ 中村友昭, 名古屋大学高等研究院, 名古屋市千種区不老町, E-mail: tnakamura@nagoya-u.jp
小竹康夫, 東洋建設株式会社鳴尾研究所, 西宮市鳴尾浜 1-25-1, E-mail: kotake-yasuo@toyo-const.co.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
Yasuo KOTAKE, Naruo Res. Inst., Toyo Construction Co., Ltd., 1-25-1, Naruohama, Nishinomiya 663-8142, 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 is applied to simulate wave-induced motion ofa sloping-top caisson floating above a rubble mound. As a result, numerical results show that the motion of the floating
caisson is affected by the wave period and the wave height, and suggest that the floating caisson can collide with the
top of the mound for large wave period and large wave height. In addition, the model is a useful tool in developing
effective countermeasures to reduce the motion of the floating caisson.

（Immersed Boundary；IB）法(6)に基づく埋め込み境界モジュール（IB Module；以下，IBM と略す）
，渦(2)と浸透滲出流(3)の効果を導入した上で，浮遊砂の巻き上げ，移流拡散，沈降の各輸送過程を考慮して浮遊砂の解析を行い，掃流砂と浮遊砂の輸送に伴う地形変化を計算する底質輸送モジュール（Sediment Transport Module；
以下，STM と略す）が，XNS との相互作用を考慮できるように
Two-Way カップリング手法により組み込まれている．本章では，
XNS，VFM，IBM，STM，カップリング手法の説明を行うとともに，パラメータの設定値について述べる．

１．緒言
上部斜面堤は，上部斜面に作用する波力により堤体の安定性の向上を図れることから，断面を小型化することができ，工費縮減に効果的な構造形式である．その一方で，上部斜面堤が係留状態にある据付時に波浪が作用すると，断面を小型化していることも影響して大きな動揺が生じる可能性があり，ケーソン据付前のマ
ウンドの天端に接触し，損傷する危険性がある．本研究では，一般的な断面形状の上部斜面堤を取り扱い，簡単のために係留索が設置されていない堤体の動揺特性を，流体，構造，地形変化の 3
者間の相互作用を解析可能な 3 次元流体・構造・地形変化連成数
値計算モデル（Three-Dimensional Coupled Fluid-Structure-SedimentInteraction Model；以下，FSSM と略す）(1)(2)(3)により検討を行う．

２．
１ 拡張Navier-Stokesソルバー（XNS），
VOF モジュール（VFM），
埋め込み境界モジュール（IBM）
流体の運動は XNS で，気液界面の追跡は VFM で，可動構造物の運動は IBM で取り扱う．以下，XNS，VFM，IBM の支配方程式と数値計算スキームの説明を行う．

２．3 次元流体・構造・地形変化連成数値計算モデル
3 次元流体・構造・地形変化連成数値計算モデル（FSSM）(1) (2)(3)の計算対象領域を Fig. 1 に示す．FSSM は，メインソルバーと 3つのモジュールから構成されている．メインソルバーは，透過性材料の間隙内部の流体を含む非圧縮性粘性流体からなる全気液相に適用できるように一般化し，さらに可動構造物の運動や地形変化を取り扱えるように拡張した連続式と Navier-Stokes（NS）方程式に基づく Large-Eddy Simulation（LES）である（Extended NS
Solver；以下，XNS と略す）
．その XNS には，気液界面を追跡す
る Volume of Fluid（VOF）法(4)を改良した Multi-Interface Advection(5)に基づくVOF モジュール（VOF
and Reconstruction Solver（MARS）
Module；VFM），可動構造物を取り扱う体積力型埋め込み境界

２．１．１ 支配方程式
気液界面の位置，可動構造物の運動，地形変化を考慮した流体の運動を取り扱う XNS と，気液界面の追跡を取り扱う VFM の支配方程式をまず説明する．
LES では，Grid Scale（GS）と呼ばれる大きなスケールの渦を直接計算し，普遍性が期待できる Subgrid Scale（SGS）と呼ばれる小さなスケールの渦をモデル化する．実流速ベクトルを vif ，上付きの－をグリッドフィルタとすると，
vif  vif  vif T(1)のように GS 成分 vi  uvw  と SGS 成分 vi  に分解で
きる．また，VOF 法において VOF 関数 F （0  F  1 ； F  0 ：空気；0  F  1 ：気液界面； F  1 ：水）が各計算格子での流体率を表しているように，体積力型 IB 法(6)では可動構造物を各計算
格子での体積率 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 ：可動構造物の角速度ベクトル；fFig. 1 Typical computational domain of the FSSM.
1

f

f

f

f

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5

riob ：可動構造物の重心からの相対位置ベクトル；eijk ：置換記号）
Tとしたとき，体積平均速度ベクトル vi （  u v w ）を
vi  vi  vi  1  Fob  vif  Fob viob

(2)

(3)

vi  1  Fob  vif (4)Ra ija  Lma
ij  CB Lij  CS D Dij

 q

*

Lmij  vi v j  vi v j mvi v jm  C A 1  m  vti  x
j

m





Ri  

12CD 2 ˆ 1  m 
md

2
50

vi 

CD1ˆ 1  m 
2md 50(14) 
H ij  vi v j  vi v j

(15)

 





(16)

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

(17)

CB 

(7)

F ˆ
xi 

(13)

とすると，係数 CB ， CS はそれぞれ

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

f i s   
R ij  vi  vi v j  v j  vi  vi v j  v jF  mv j F

 Fq*
t
x j




Lij  vi v j  vi v j





(6)

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

 とグリッドフィルタ幅と表される．ここで，テストフィルタ幅   ；上付きの～：テストフィルタ）としたと
 の比を  （  き，Germano Identity Lij ， H ij ， R ij ， Mij をそれぞれ



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





(12)

LRij  vi  vi v j  v j  vi  vi v j  v j

(5)

x j

(11)

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

 mv j(10)と表される．また，修正 Leonard 応力テンソルを Lmij ，ひずみ速度テンソルの GS 成分 Dij の大きさを D とすると， ija は Horiuti(9)より次式で与えられる．

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

v 
2  
 mˆ j 
3 xi x j 

Qi  vi q*  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 R M kl M kl  R MijR M kl
a
ij

a
ij

a
ij

a
kl

a
kl

kl

(18)

(19)

と表される．ただし， CS が負になることで計算が不安定にならないように，Horiuti(9)に倣って CS を Clipping した．また，減衰関
数  は Hinatsu(12)と Cruz ら(13)を参考に
Cg   h0  x 
 z

sin 
,1 in artificial
 min 


2
h




 damping zones (20)otherwiseとした．ただし，C は減衰係数，h は静水深， は水位変動，
は減衰領域の長さ， x は減衰領域の開始位置からの距離， z は
．なお，上述したように構造物の底面からの距離である（Fig. 2）
内部では vi が viob となるように強制することから，岩田ら(14)が述べているように構造物の内部での ̂ やˆ は任意であり，計算結果に与える影響は十分に小さいことを確認している．

(8)vi v j v jFig. 2 Artificial damping zone in a typical computational domain.

(9)2

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5
続いて，可動構造物の運動を取り扱う IBM の支配方程式を説明する．可動構造物に作用する外力ベクトルを f i ex ，底面との接触による反力ベクトルを f i re ，可動構造物の質量を m ob ，可動構
造物の重心回りの慣性テンソルを I ijob とすると，可動構造物の並進速度ベクトル uiob に関する運動量保存則と角速度ベクトル iob
に関する角運動量保存則はそれぞれ
duiobm f i ex  f i re
dt
ob

d ob
j
ob

I ij

dt


 
mn
m n  n 1 2
 n
Bn    
n
n

xi
xi  m  C A 1  m  ˆp2n
 m vixi  q*n 1
t n 1 2

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

 eijk rj  f kex  f kre 

fi

(22)

と書かれる．ここで，Xiao ら(15)に倣って可動構造物の表面に作用するせん断力の影響が十分に小さいと仮定すると， f i ex は p
f ex   Fob    ob g   xi y j zk
i , j ,k
 x


vi

(23)


A 

vi

n 1

 vi

 vi

p2

(25)



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

(26)

p1

n 1/ 2



(28)

t
m
m n  C A 1  m n  ˆ nnBn

 1  Fobn  vi

p1

 Fobn viob n

(29)

n

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

n

ˆ n

x j

i





n

D1

sn

n

i

2d 50

j

n

j

(30)

n
1 

2m nˆ n Dij  n mn ija n   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 

Bn

2

t

n 3 2 (32)


12CD 2ˆ 1  m 
t n 1 2
d 502
m  C A 1  m n n(31)

t n 1 2  t n  3 2  t n 1 2 nBn  1 n(33)と表される．一方，空間方向の離散化には，式(30)右辺第 1 項の対流項については 5 次精度 MUSCL（Monotone Upstream-Centered
Schemes for Conservation Laws）(19)を，他の項については 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 を求める．

(24)

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

n

p2





p1

p1

A2ni  2 t n  3 2 A0ni   t n  3 2  t n 1 2  A0ni1  t n 1 2 A0ni 2 n m n  p n
t n 1 2
 vi  n
 m n gi   m  C A 1  m n   ˆ n xin 1 2 2 An  B n
t n 1 2 n  t
A1i 
  A0ni 
2i 

 
2
6

  
p2n

 m n vi v jn
0i

２．１．２．１ 拡張 Navier-Stokes ソルバー（XNS）
XNS では，式(5)と式(6)で表される vi と p の連成計算に SMAC
（Simplified Marker and Cell）法(17)を採用した．ただし，式(6)右辺
第 1 項の圧力勾配項と第 2 項の重力加速度項については梶島(18)を参考に 1 次精度 Euler 陽解法を，式(9)右辺第 1 項の線形抵抗力項については安定化のために 1 次精度 Euler 陰解法を，その他の項については 3 次精度 Adams-Bashforth 法を適用した．すなわち，
タイムステップ数を上付きの n ，第 n ステップから第  n  1 ス
テップまでの時間の刻み幅を t n 1/ 2 とすると，予測段階の第 1 ス
テップ，予測段階の第 2 ステップ，修正段階はそれぞれvi

Fobn viob n  vi

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

２．１．２ 数値計算スキーム

p1



ob n

となる．したがって，式(25)と式(28)より，予測段階の第 2 ステッ
プは最終的に次のように表される．

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

vi

(27)





２．１．２．２ VOF モジュール（VFM）
VFM では，式(7)で表される VOF 関数 F の移流方程式に前述したように MARS(5)を適用した．なお，隣接する格子の境界における F の不連続性をできるかぎり押さえながら気液界面の輸送を行うために，1 タイムステップでの界面輸送を全方向に一度に計算する Unsplit Time Integration を時間方向に多段階で行う
Fractional Step Method(20)を採用した．

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

3

２．１．２．３ 埋め込み境界モジュール（IBM）
IBM では，式(21)と式(22)に 1 次精度 Euler 陽解法を適用し，可動構造物の並進速度ベクトル uiob と角速度ベクトル iob をそれぞCopyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5
れ次式から求める．
t n 1/ 2 ex n
 fi  fi re n 
m ob

(34)

t n 1/ 2
eijk rjn  f kex n  f kre n 
I ijob(35)

uiob n 1  uiob n 
n 1
n
 ob
  objj

２．２．２．１ 掃流砂量式
漂砂現象に与える渦と浸透滲出流の影響を考慮するために，掃流砂量 qi に対して渦の効果が組み込まれており(2)，浅野・米永(23)と同様に限界 Shields 数 *c と掃流砂量 qi に対して浸透滲出流による層流抵抗と乱流抵抗の効果が組み込まれている(3)．ここで，底質粒子は非粘着性で均一の粒子径の球体とし，底面流速および浸透滲出流速 w（滲出を正とする）による慣性力は無視する．このとき，底質粒子に作用する外力は以下から構成される．
1. 底質粒子の水中重量 W

そして，可動構造物の重心位置ベクトル xiob と回転角ベクトル
 iob を，2 次精度 Crank-Nicolson 法に基づいて，それぞれ
xiob n 1  xiob n  iob n 1  iob n 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
W     s   w  gd 503
6
ここで，  s は底質粒子の密度である．
2. 底面流速による乱流抵抗力 FD と揚力 FL

により更新する．なお，求められた可動構造物の位置に基づく Fobの更新には，Yuki ら(21)の双曲線正接関数を用いる手法とは異なり，各計算格子を各方向に数分割ずつした小格子を考え，その分割した小格子の中央から可動構造物を構成する面までの垂線の長さを基に Fob を算出することで可動構造物の総体積の増減を抑えた．
以上をまとめると，IBM での計算の手順は次のようになる．
1. f i ex （式(23)）と f i re を求める．
2. 式(34)より uiob を，式(35)より iob を求める．
3. 式(36)より xiob を，式(37)より iob を求める．
4. 求めた xiob と iob より， Fob を更新する．

(40)

1
FL   CL  w d 502 vr2  k LCFD vr2
8

(41)

1
Fw   wCD1 w d 502 w w  3 wCD 2  w w d 50 w
8(42)

ただし，上式右辺第 1 項は乱流抵抗力，右辺第 2 項は層流抵抗力であり，両者で遮蔽係数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 と同方向および直交方向の力のつり合いより次式が得られる．

(38)

ここで， z s は基準面からの底面の高さ， qx と q y はそれぞれ x ，
y 軸方向の掃流砂量， pN は底質の巻き上げ関数， qzbS は底面での z 軸方向の沈降による浮遊砂輸送フラックスである．
２．２．２ 掃流砂輸送モデル
掃流砂輸送モデルは，Roulund ら(22)を参考に構築した以下の掃流砂量式と斜面崩壊計算から構成される．

Fig. 3 Critical force on a sand particle
on a sloping bed.

1
FD   CD1 w d 502 vr2  CFD vr2
8

ここで，vr は底質粒子の位置での vb に対する相対底面流速，CL は揚力係数，kL は CL
（vb：底質粒子の平均輸送速度）と CD1 の比（= CL / CD1）であり，FD と FL で遮蔽係数が等しいとした．また，限界 Shields 数を上回る状況では比較的大きい底面流速が生じていると考えられることから，底面流速による層流抵抗の影響は十分小さいと仮定した．
3. 浸透滲出流速 w による作用力 Fw

２．２ 底質輸送モジュール（STM）
STM は，地形変化を計算する底質の質量保存式，地形変化の計算に必要な掃流砂量を求めるとともに，地形変化の計算後に斜面の崩壊計算を行う掃流砂輸送モデル，地形変化の計算に必要な底質の巻き上げと沈降による浮遊砂輸送フラックスを求め，浮遊砂の解析を行う浮遊砂輸送モデルから構成される．以下，それぞれの説明を行うとともに，数値計算スキームを説明する．

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

(39)

(b) Fw > W(a) Fw  W
Fig. 4 Force on a moving sand particle in bed-load motion on a sloping bed.
4

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5

vrc2 

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

c



w





W  Fw  sin  sin  c c  arctan 2 CFD vrc  W  Fw  sin  cos  c z   0.9

(43)

(44)

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

 sW

(45)

(46)

1  s kL  CFD

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

FD cos 1  W  Fw  sin  cos    

２．２．２．１．２ 掃流砂量
単位幅，単位時間当たりの i 方向の掃流砂量 qi は，Engelund・Fredsøe(24)による次式から求められる．

ここで， pEF は掃流形態で移動する底質粒子の割合，vbi は vb の i
方向成分であり，底質粒子の動摩擦係数を  d （  tan  d ； d ：動摩擦角），Shields 数を * とすると， pEF は

(49)

v 2f

vf

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

vr sin 1  Cv f v f sin (57)

vr cos 1  vb  Cv f v f cos (58)

1    (59)

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

vr 

z

 2 1 1 
022
 1  4 2  z   z   1  exp   z   z   A   dz 

 



(55)

(50)

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

vsurf  d W  Fw  cos   FL 

である．ここで，式(55)から式(58)に対して未知数は vb，vr，，1 の 4 つであることから，繰り返し計算により vb とその方向，すなわち vbi が求められる．したがって，前述した Fw  W での *cとあわせて，最終的に式(48)と式(49)より qi が求められる．
一方，Fw が W を上回っているとき，底質粒子に動摩擦力が作用しなくなることから，Fig. 4(b)に示すように底質粒子に作用する力は FD と(Fw – W) sin となり，次式が成立する．

と与えられる，ただし， d   s である．また， * は

* 

(54)

が得られる．また， vb ， Cv f v f ， vr の関係より，(48)

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

(53)

と書ける．
そこで，
漂砂現象に与える渦の影響を考慮するときは，上式で求めた vi から式(51)の vsurf を計算する．
続いて，式(48)の vbi を求めるために，Fig. 4 に示す掃流形態で移動する底質粒子を考える．まず Fig. 4(a)に示した Fw が W を下回っているとき，vb と同方向および直交方向の力のつり合いより(47)

1qi   d 50 pEF vbi
6svi  vi  CSGS vi  vi  CSGS (vi  vi )一方，Fw が W を上回っているとき，底質粒子に静止摩擦力が作用しなくなることから，次式が得られる．

 *c
0
 *c 0sより，左辺 vi は SGS 成分のうち比較的大きなスケールの変動を表しており，Transfer Field と呼ばれる右辺 (vi  vi ) は，GS 成分から小スケールの渦を消去した成分 vi を GS 成分 vi から引いていることから，GS 成分のうち比較的小さなスケールの変動を表している．ここで，式(53)右辺の Transfer Field は XNS で求めた GS
成分 vi から直接計算できることから，SGS 成分 vi が vi ，すなわ
ち Transfer Field で近似できるとし，またその影響の程度を表す無次元パラメータを CSGS とすると，LES で直接計算できる GS 成分
vi に SGS 成分 vi の影響を考慮した vi の近似値 vi はと表される．ただし，vrc2 0 は w = 0 かつ = 0 での vrc2 であり，式(43)に Fw =  = c = 1c = 0 を代入した次式から求められる．

vrc2 0 svi  vi  vi

2
rc

v
v
2
C
s
s
gd

1
 *c 
 50  vf   1 gd50  vrc2
2
v fc 0vrc2 0vrc2 0 *c 02
 s  1 gd50 Cvf  s  1 gd50

s

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

 k  k exp   k 6   4.535  k  2000 (52)



(51)

 Fw  W  sin 
CFD

(60)

となる．また，Fig. 4(b)より，

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

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

(61)

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

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5となる．したがって，上式より vb とその方向が求められることから，前述した Fw > W での *c とあわせて，最終的に式(48)と式(49)より qi が求められる．
以上より，掃流砂量 qi は vbi と前述した *c から求められ，そのvbi は Fw  W のとき式(55)から式(58)より，Fw > W のとき式(59)から式(62)より求められる．

 2 3636 
ws  sgd 50 



 3 d
d* 
*


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

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

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

vbslide 


3

d 50

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

(63)

(64)
C n 1  C n (65)z sn 1  z sn 

(66)

を支配方程式として用いた．ここで，C は浮遊砂濃度，qiC ，qiD ，
qiS はそれぞれ移流，拡散，沈降による浮遊砂輸送フラックスであり，それぞれ次式で表される．

(68)

qiS  mwsi C

(69)

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

(73)

 が r   r を上回ったとき，以下の斜面崩壊計算を行う．

t h 1/ 2  qxslide h q yy
1  m n  x

slide hz sh 1  z sh 





(74)

iii. 斜面の傾斜角  がすべての場所で r   r を下回るまで，
手順 i から ii を繰り返す．

(67)

C
xi

qiD  m s(72)

なお，ここでは Roulund ら(22)に倣って斜面崩壊が瞬間的に起こると仮定し，手順 1 から 5 までとは独立したタイムス
テップ数 h と時間の刻み幅 t h 1/ 2 を使って計算を進める．
i. 式(64)より斜面崩壊に伴う qislide を求める．
ii. 式(65)を 1 次精度 Euler 陽解法により離散化した次式より
z s を更新する．

２．２．３ 浮遊砂輸送モデル
榊山ら(26)や中山ら(27)を参考に，透過性材料の間隙内部にも適用できるように一般化した浮遊砂の質量保存式qiC  mvi C

t n 1/ 2  qCj n  q Dj n  q Sj n 
m n x j

4. 式(48)より qi を求める．なお，式(55)から式(58)の計算には，
Roulund ら(22)に倣って Newton-Raphson 法を適用した．
5. 式(38)を 1 次精度 Euler 陽解法により離散化した次式より z sを更新する．

6.

C qCj  q Dj  q Sj   0
t x j(71)

２．２．４ 数値計算スキーム
STM での計算の手順は次の通りである．
1. 式(67)より qiC を，式(68)より qiD を，式(69)より qiS を求める．
なお，式(67)と式(69)の離散化には安定性を考慮して 1 次精度風上差分法を，式(68)の離散化には 2 次精度中央差分法を適用した．
2. 求めた qiS より qzbS を，式(71)より pN を求める．
3. 求めた qzbS と pN を底面での境界条件として，式(66)を 1 次
精度 Euler 陽解法により離散化した次式より C を求める．

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

m

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

また，斜面崩壊に伴う底質の質量保存式は，斜面崩壊時には浮遊
砂の巻き上げや沈降が生じないと仮定すると，式(38)に pN 
qzbS  0 を代入した
slide
1  qxslide q y 
z s



0
t 1  m  xy 

if  *   *c

0

1.5
pN     *   *c   s  10.6 g 0.6 d 500.8C p w0.2   *c 

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

qislide (70)

ただし，  s は浮遊砂の乱流拡散係数， wsi は底質粒子の沈降速度
T
ベクトル（   0 0  ws  ；ws ：底質粒子の沈降速度）であり，
ws は次に示す Rubey の式により与えられる．

6

２．３ XNS・VFM・IBM・STM 間のカップリング手法
FSSM では，XNS，VFM，STM，IBM の接続に Fig. 5 に示す
Two-Way カップリング手法を採用した．具体的には，FSSM での計算の手順は次のようになる．
1. XNS を実行する．
2. XNS で得られた流速場を用いて VFM を実行する．
3. 流速場に加えて，VFM で得られた気液界面の位置を用いてSTM を実行する．
4. 流速場と気液界面の位置に加えて，STM で得られた地形変
Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5
 r ，斜面崩壊開始基準 r ，斜面崩壊完了基準 r は Roulund ら(22)に基づいて，底質巻き上げ係数 C p は Nielsen(29)に基づいて，揚力係数と非線形抵抗力係数の比 kL と底面流速の遮蔽係数は中川・辻本(32)に基づいて，浸透滲出流速の遮蔽係数w は中村・水谷(3)による試行計算の結果に基づいて設定した．

化を用いて IBM を実行する．
5. VFM，STM，IBM で得られた値をフィードバックして XNSを実行する．
6. 手順 2 から 5 を計算終了時刻まで繰り返す．
２．４ パラメータ
Table 1 と 2 に，XNS と STM で用いたパラメータをそれぞれ示
す．Table 1 に示した XNS の LES に関して，グリッドフィルタとテストフィルタは Horiuti(9)に倣ってトップハットフィルタ（ボックスフィルタ）を採用し，グリッドスケールのフィルタ幅  は各位置における各方向の格子幅と等しいとした．また，フィルタ幅
その他，
付加質量係数 C A ，比  は Garmano ら(30)に倣って設定した．
非線形抵抗力係数 CD1 ，線形抵抗力係数 CD 2 は水谷ら(8)の実験結果に基づいて，減衰係数 C は中村ら(31)による試行計算の結果に基づいて設定した．Table 2 に示した STM については，摩擦速度
水平面での限界 Shields 数 *c 0 ，
最大静止摩擦係数  s ，の係数 Cvf ，動摩擦係数  d ，van Driest 型減衰係数 A ，底質粒子の水中安息角３．計算条件
本研究では，Fig. 6 に示すように，簡単のために係留索が設置されていない一般的な断面形状の上部斜面堤を対象に解析を行った．ここで，FSSM は上述したように 3 次元モデルであるが，後述するように細かい計算格子を用いたことから，計算負荷の軽減のために，ここでは単位奥行きの断面 2 次元計算を行った．
Fig. 6 に示すように，約 500 kg の捨石を想定した中央粒径 d50 =0.356 m からなるマウンドを設定し，その上にコンクリートを想定した密度s = 2.3×103 kg/m3 の上部斜面堤を浮かべた．また，同図の左端（x = – 46.90 m）に造波ソースを設定するとともに，同図に示す領域の沖側と岸側には，後述する入射波の 2 波長分以上に相当する約 150 m の減衰領域を確保した．座標軸は，Fig. 6 に示

Fig. 6 Initial condition of the sloping-top caisson floating above the rubble mound.

Fig. 5 Two-way coupling procedure.
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(a) Horizontal displacement x

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
動摩擦係数  d
0.51
0.4
Kármán 定数  
25.0
Van Driest 型減衰係数 A
底質粒子の水中安息角 r
32.0º
斜面崩壊開始基準 r
2.0º
斜面崩壊完了基準 r2.0º(b) Vertical displacement z

7(c) Rotation angle 
Fig. 7 Motion of the floating sloping-top caisson in the still water.
Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5
入射波は規則波とし，海上工事で施工可能な風波を想定して入射波周期 T を 3.0，5.0，7.0 s の 3 種類とした．また，入射波高 H
は施工作業の中止基準を参考に 1.0 m とし，比較のために T = 7.0 s
では H = 1.8 m のケースも行った．静水深 h は HWL を想定した
15.5 m とした．波の作用時間は 60 s 間とし，その間の堤体の動揺特性を検討した．

すように，水平岸向き座標を x 軸，鉛直上向き座標を z 軸とした．
なお，マウンドの捨石は動かないものとし，堤体とマウンドの接触は考慮しなかった．また，マウンドの空隙率 m は 0.40 とした．
計算格子は，x 軸方向に対して，堤体が移動する可能性がある
範囲については 50 mm の等間隔格子を採用し，それ以外の範囲については計算負荷のさらなる軽減のために不等間隔格子を用いた．
また，z 軸方向に対しては 100 mm の等間隔格子を用いた．流速と圧力の境界条件として，底面に Slip 条件を，上面に圧力一定の条
また，件を，沖側境界と岸側境界に Sommerfeld 放射条件を課した．
VOF 関数の境界条件には，全境界に勾配ゼロの条件を用いた．なお，本論では地形変化は取り扱わないことから，FSSM のうち底質輸送モジュール STM は使用しなかった．
ここでは，堤体の動揺が問題となる状況を対象とするために，静水状態で堤体がマウンドの天端から 0.50 m 浮いた据付直前の状況を取り扱った．ただし，実施工では堤体内部を隔壁で仕切った各区画の水位を調整し，動揺を最小限に抑えながら沈設作業を行うと考えられること，また静水状態で堤体が大きく動揺してしまうと，波作用時の動揺特性の把握が困難になることから，静水状態で堤体の動揺ができる限り抑えられるように堤体内部の各区画の水量を 0.5 mm 単位で調節した．その結果，堤体内部の沖側および中央の区画と岸側の区画の水位差を 321.0 mm としたときに，堤体の動揺が最も抑えられることが分かった．そのときの堤体の重心位置の水平変位x，鉛直変位z，回転角の時間変化をFig. 7 に示す．ここで，t は時刻であり，正のx は岸向きの変位を，正のz は上向きの変位を，正のは時計回りの回転を表す．Fig. 7より，計算開始 60 s 後にx = 7 mm ほど沖側に移動するとともに，
z = 8 mm ほど浮き上がっているものの，堤体の移動は非常に小さく抑えられていることが確認できる．そこで，区画間の水位差を 321.0 mm としたときを対象に以下解析を行った．

４．計算結果及び考察
Fig. 8 に H = 1.0 m，T = 3.0 s での，Fig. 9 に最も大きな動揺が生じた H = 1.8 m，T = 7.0 s での堤体の動揺の様子を例示する．Fig. 8より，H = 1.0 m，T = 3.0 s のケースでは堤体内に水塊が入りそうなほどの打ち上げが斜面上で生じているものの，堤体の動揺は小さいことが確認できる．一方，Fig. 9 に示した H = 1.8 m，T = 7.0 sのケースでは，波谷が堤体の前面を通過した後に堤体が反時計回りに回転し，堤体の上部が沖側に大きく傾くことが分かる．その後，堤体の前面に波峰が作用すると，堤体が時計回りに回転し，堤体の上部が大きく岸側に傾くことが確認でき，Fig. 8 と比較して大きな動揺が生じていることが分かる．
Fig. 10 に堤体の重心位置の水平変位x，鉛直変位z，回転角の時間変化を示す．同図より，x，z，の振幅は周期が長くなるほど，また波高が大きくなるほど増加する傾向が確認できる．
したがって，z との振幅がともに大きな H = 1.8 m，T = 7.0 s のケースでは堤体がマウンドの天端に接触する危険性が高いことが示唆される．また，Fig. 10(a)より，堤体の岸沖方向の漂流挙動は周期によって異なっており，周期が短い T = 3.0 s や T = 5.0 s では
岸側に，周期が長い T = 7.0 s では沖側に漂流していることが分かる．その一方で，z は周期や波高によらず若干沈み込んだ状態を中心に上下運動しており，は沖側や岸側に偏ることなく初期状
態（ = 0 rad）を中心に回転運動していることが確認できる．

(a)

(b)

(c)

(d)

(e)

(f)

(g)
Fig. 8 Wave-induced motion of the sloping-top caisson for T = 3.0 s, H = 1.0 m, h = 15.5 m: (a) t = 57.0 s (t / T = 9.0),(b) t = 57.5 s, (c) t = 58.0 s, (d) t = 58.5 s, (e) t = 59.0 s, (f) t = 59.5 s, and (g) t = 60.0 s (t / T = 10.0).
8
Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
Fig. 9 Wave-induced motion of the sloping-top caisson for T = 7.0 s, H = 1.8 m, h = 15.5 m: (a) t = 49.0 s (t / T = 7.0),(b) t = 50.0 s, (c) t = 51.0 s, (d) t = 52.0 s, (e) t = 53.0 s, (f) t = 54.0 s, (g) t = 55.0 s, and (h) t = 56.0 s (t / T = 8.0).

以上より，堤体の動揺量は周期や波高の増大とともに大きくなり，マウンドの天端に接触する危険性が高くなることを確認するとともに，据付時の堤体の動揺特性を検討するための本数値計算
モデル FSSM の有用性を示した．
５．結論
本論では，一般的な断面形状の上部斜面堤を取り扱い，簡単のために係留索が設置されていない状況を対象に，堤体の動揺特性
を流体，構造，地形変化の 3 者間の相互作用を解析可能な 3 次元
流体・構造・地形変化連成数値計算モデルにより検討を行った．
本数値計算モデルは，流体の運動を解析するメインソルバー，気液界面を追跡するモジュール，可動構造物の運動を解析するモジ
ュール，浮遊砂と地形変化を解析するモジュールから構成されており，それらの接続には流体，構造，地形変化の相互作用が考慮本数できるように Two-Way カップリング手法が採用されている．
値計算モデルによる検討の結果，堤体の動揺特性は周期や波高によって異なり，本研究で対象とした海上工事で施工可能な波浪条
件では，周期や波高の増大とともにマウンドの天端に接触する危険性が高くなることを示した．また，据付時に堤体の動揺特性を検討するための本数値計算モデルの有用性を確認した．ただし，より現実的な係留索を設置した場合の堤体の動揺特性は解析していないなど課題を残していることから，本数値計算モデルの改良を含めさらなる検討を行っていく所存である．

(a) Horizontal displacement x(b) Vertical displacement z

参考文献(1) 中村友昭, 水谷法美, “3 次元流体・構造・地形変化連成数値計算モデルに基づく遡上津波による陸上構造物周辺の局所
洗掘現象に関する数値シミュレーション,” 第 24 回数値流体

9(c) Rotation angle 
Fig. 10 Wave-induced motion of the floating sloping-top caisson.
Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-5

pp. 39-44.
(17) Amsden, A. A. and Harlow, F. H., “A simplified MAC technique forincompressible fluid flow calculation,” J. Comp. Phys., 6 (1970), pp.
322-325.
(18) 梶島岳夫, “乱流の数値シミュレーション,” 養賢堂 (1999),
225 p.
(19) Yamamoto, S. and Daiguji, H., “Higher-order-accurate upwind
schemes for solving the compressible Euler and Navier-Stokes
equations,” Computers Fluids, 22, 2-3 (1993), pp. 259-270.
(20) 秋山 守, 有冨正憲, “新しい気液二相流数値解析－多次元流
動解析－,” コロナ社 (2002), 261 p.
(21) 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.
(22) 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.
(23) 浅野敏之, 米永敏之, “前浜砂層内の浸透流を考慮した波打ち帯の漂砂の解析,”海岸工学論文集, 47 (2000), pp. 511-515.
(24) Engelund, F. and Fredsøe, J., “A sediment transport model forstraight alluvial channels, Nordic Hydrology, 7 (1976), pp. 293-306.
(25) Bardina, J., “Improved turbulence models based on large eddy
simulation of homogeneous, incompressible turbulent flows,” Ph.D.
dissertation, Stanford University (1983).
(26) 榊山 勉, 伊藤 裕, 小林正典, 板井雅之, 小林英次, “傾斜堤を透過する浮遊砂の濃度評価への数値計算モデルの適用,”
海岸工学論文集, 51 (2004), pp. 456-460.
(27) 中山哲嚴, 牧野弘幸, 新井雅之, 小林 学, 佐藤勝弘, “導水中に含まれる浮遊砂量を低減するための海水導入工の構造に関する研究,” 海岸工学論文集, 54 (2007), pp. 661-615.
(28) 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.
(29) Nielsen, P., “Coastal Bottom Boundary Layers and SedimentTransport,” Advanced Series on Ocean Eng. Vol. 4, World Scientific(1992), 340 p.
(30) 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.
(31) 中村友昭, 水谷法美, 藤間功司, “遡上津波の変形と陸上構造物に作用する津波力に関する 3 次元数値シミュレーション,”
ながれ, 29, 2 (2010), pp. 107-117.
(32) 中川博次, 辻本哲郎, “ 移動床流れの水理,” 技報堂出版(1986), 310 p.

力学シンポジウム, E10-4 (2010), 9 p.
(2) 中村友昭, 水谷法美, “Large-Eddy Simulation における漂砂計
算への渦の効果の導入と遡上津波による洗掘現象への適用に関する研究,” 土木学会論文集 B3(海洋開発), 67, 4 (2011),
印刷中.
(3) 中村友昭, 水谷法美, “浸透滲出流による層流・乱流抵抗の影響を考慮した漂砂モデルとその適用に関する研究,” 土木学会論文集 B2(海岸工学), 67, 3 (2011), 印刷中.
(4) 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.
(5) 功刀資彰, “自由界面を含む多相流の直接数値解析法,” 機械学会論文集 B 編, 63, 609 (1997), pp. 1576-1584.
(6) 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.
(7) Brackbill, J. U., Kothe, B. D., and Zemach, C., “A continuum
method for modeling surface tension,” J. Comp. Phys., 100 (1992),
pp. 335-354.
(8) 水谷法美 , 前田健一郎 , Ayman M. Mostafa, William G.
McDougal, “透水性構造物の抵抗係数の評価と波・潜水透水性
構造物の非線形相互作用の数値解析,” 海岸工学論文集, 43(1996), pp. 131-135.
(9) Horiuti, K., “A new dynamic two-parameter mixed model forlarge-eddy simulation,” Phys. Fluids, 9, 11 (1997), pp. 3443-3464.
(10) 岩田好一朗, 川崎浩司, 金 度三, “水中構造物による砕波の数値解析,” 海岸工学論文集, 42 (1995), pp. 781-785.
(11) 財団法人沿岸開発技術研究センター , “ 数値波動水路（CADMAS-SURF）の研究・開発,” 沿岸開発技術ライブラリーNo. 12 (2001), 457 p.
(12) 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.
(13) Eric Cruz, 横木裕宗, 磯部雅彦, 渡辺 晃, “非線形波動方程式に対する無反射境界条件について,” 海岸工学論文集, 40(1993), pp. 46-50.
(14) 岩田隆一, 梶島岳夫, 竹内伸太郎, “Immersed-boundary andvolume-of-fluid カップリング法による粒子列の間を通過する
上昇気泡の解析,” 混掃流, 23, 5 (2010), pp. 531-538.
(15) 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.
(16) 本多和彦, 富田孝史, 西村大司, 坂口 章, “多数の津波漂流物を解析する数値モデルの開発,” 海洋開発論文集, 25 (2009),

10

Copyright © 2011 by JSFM