第 24 回数値流体力学シンポジウム
E10-4

3 次元流体・構造・地形変化連成数値計算モデルに基づく遡上津波による陸上構造物周辺の局所洗掘現象に関する数値シミュレーション
Numerical Simulation of Local Scouring around Inland Structure due to Run-up Tsunami
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 two-way coupled fluid-structure-sediment interaction model (FSSM) is applied to predict local
scouring around an inland square structure due to a run-up tsunami. The FSSM consists of a main solver and threemodules. The main solver is an extended Navier-Stokes solver (XNS) to compute incompressible viscous multi-phaseflow including seepage air-water flow in porous media. The first module is a volume of fluid module (VFM) to track
air-water interface motion, the second module is an immersed boundary module (IBM) to track structure motion, andthe third module is a sediment transport module (STM) to track seabed profile evolution and suspended sediment motion. Compared with experimental data measured in hydraulic experiments in terms of water surface elevation,
pore-water pressure, and final land profiles, numerical results demonstrate the computational capability of the FSSM to
predict reasonably well the experimental data including a scour hole around the offshore edge of the structure.

（Immersed Boundary；IB）法(11)に基づく埋め込み境界モジュール（IB Module；IBM），浮遊砂の巻き上げ，移流拡散，沈降の各輸送過程を考慮して浮遊砂濃度を解析し，掃流砂と浮遊砂の輸送に伴う地形変化を追跡する底質輸送モジュール（Sediment Transport
Module；STM）
が，
XNSとの相互作用を考慮できるようにTwo-Way
カップリング手法により組み込まれている．本章では，XNS，VFM，IBM，STM の説明を行うとともに，XNS，VFM，IBM，STM 間のカップリング手法を述べる．

１．緒言
2004 年12 月のスマトラ島沖地震や2010 年2 月のチリ地震に伴う大津波によって陸上構造物の周囲に大規模な洗掘が生じた事例が報告されている．この種の現象に関しては，これまでに護岸の法先や背後に生じる洗掘(1)(2)(3)(4)や透水性斜面上に設置した円柱周辺の洗掘(5)に対する検討が行われている．また，Nakamura ら(6)(7)は遡上津波による陸上構造物周辺の洗掘発生機構を水理実験と数値計算の両面から明らかにしており，その結果に基づいて構築した地形変化予測モデルにより洗掘の再現計算を行っている．しかし，津波の作用に伴う地形変化の計算は行っているものの，その地形変化が津波に与える影響については無視できると仮定して地形変化から津波へのフードバックは行っていないことから，流体と地形変化の相互作用を考慮していない点に課題を残している．
本研究では，遡上津波による矩形構造物周辺の洗掘に関する水理実験(6)を対象に，流体，構造，地形変化の 3 者間の相互作用を考慮した 3 次元流体・構造・地形変化連成数値計算モデル（ThreeDimensional Two-Way Coupled Fluid-Structure-Sediment Interaction(8)
Model；以下 FSSM と略す）
を適用して再現計算を行う．
そして，
水理実験との比較によりその妥当性を検証するとともに，相互作用を考慮していない Nakamura ら(7)による計算結果との比較を行い，遡上津波による地形変化の再現性について検討を行う．

２．
１ 拡張Navier-Stokesソルバー（XNS），
VOF モジュール（VFM），
埋め込み境界モジュール（IBM）
流体の運動は，VFM で求めた気液界面の位置と IBM で求めた構造物の運動の影響を考慮した XNS で，気液界面の追跡は VFM
で，構造物の運動は IBM で取り扱う．以下，XNS，VFM，IBMの支配方程式と数値計算スキームの説明を行う．
２．１．１ 支配方程式
気液界面の位置と構造物の運動の影響を考慮した流体の運動を取り扱う XNS と，気液界面の追跡を取り扱う VFM の支配方程式をまず説明する．VOF 法において VOF 関数 F （0  F  1 ；F  0 ：空気；0  F  1 ：気液界面； F  1 ：水）が各計算格子
での流体率を表しているように，
体積力型 IB 法(11)では構造物を各計算格子での体積率 Fob （ 0  Fob  1 ； Fob  0 ：構造物外部；
0  Fob  1：構造物表面；Fob  1：構造物内部）で表す．そして，

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

Fig. 1 Typical computational domain of the FSSM
1

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4
ひずみ速度テンソル Dij の大きさであり， Lmij と LRij はそれぞれ構造物の内部では構造物の速度になるように， Fob に基づいて相互作用力ベクトル f i ob を作用させる．ここで，実流速ベクトルをob
ob
vif ，構造物の速度ベクトルを viob （  uiob  eijk ob
j rk ； ui ：構ob
造物の並進速度ベクトル；i ：構造物の角速度ベクトル； riob ：構造物の重心からの相対位置ベクトル； eijk ：置換記号）とした
Tとき，体積平均速度ベクトルを vi  u v w  1  Fob  vif 
Fob viob と定義する．また，体積空隙率と面積空隙率が等しく，その空隙率 m の時間変化（ m t  0 ； t ：時間）が十分に小さいと仮定すると，拡張した連続式，NS 方程式，VOF 関数 F の移流
方程式はそれぞれ次式で与えられる．
  mv j 

 q*

x j

vi

m  C 1  m  t 
A



Lma
ij  vi v j  v i v j

Qi  vi q* vi CD1ˆ 1  m 
2md 50

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

(11)



 





 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

(14)

(15)

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

(4)

vi v j v j  
H ij  v i v j  vi v j

(13)

(3)

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

md 502(10)

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

CB 

T

12CD 2 ˆ 1  m  
Lij  vi v j  vi v jとすると，係数 CB ， CS はそれぞれ

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

(9)

(12)

(2)

F ˆ
xi 



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



x j

f i s  と表される．ただし，上付きの－はグリッドフィルタを表す．こ とグリッドフィルタ幅  の比を 
こで，テストフィルタ幅 （    ；上付きの～：テストフィルタ）としたとき，Germano
Identity Lij ， H ij ，R ij ， Mij をそれぞれ  mvi v j F   mv j F 

 Fq*
t
x j

 

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

(1)

m p
m
 mg i   f i s  Ri  f i ob ˆ xi
ˆ

m



(8)

I ijobduiob
 f i ex  f i re
dt

d ob
j
dt

 eijk rj  f kex  f kre 

(17)

(18)

(5)

(6)

と表される．また， ija は Horiuti(14)より
Ra ija  Lma
ij  CB Lij  CS D Dij(7)と与えられる．ただし， Lmij は修正 Leonard 応力テンソル， D は
2

Fig. 2 Artificial damping zone in a typical computational domain
Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4
1 次精度 Euler 陰解法に導入に伴う B n はそれぞれと書かれる．ここで，Xiao ら(20)に倣って構造物の表面に作用する
せん断力の影響が十分に小さいと仮定すると， f i ex は p
f i ex   Fob   ob gi All cell xi

A0ni  (19)

  m n vin v nj x jm n s n CD1 1  m  n n nfi vi v j v j  
ˆ n
2d50
n



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



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

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

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

(22)

n

n

Bn



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

２．１．２．３ 埋め込み境界モジュール（IBM）
IBM では，式(17)と式(18)に 1 次精度 Euler 陽解法を適用し，構造物の並進速度ベクトル uiob と角速度ベクトル iob をそれぞれ
t n 1/ 2 ex n
 fi  fi re n 
m ob

(30)

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

uiob n 1  uiob n n 1n  ob

 ob
j
j

(23)

から求める．そして，構造物の重心位置ベクトル xiob と回転角ベクトル iob を，2 次精度 Crank-Nicolson 法に基づいて，それぞれ次式により更新する．

obi

ここで，上述したように構造物の内部では vi が v になるように強制することから，相互作用力ベクトル f i ob n はFobn  viob n  vip1 n２．１．２．２ VOF モジュール（VFM）
VFM では，式(3)で表される VOF 関数 F の移流方程式に前述したように MARS(10)を適用した．なお，隣接する格子の境界における F の不連続性をできるかぎり押さえながら気液界面の輸送を行うために，1 タイムステップでの界面輸送を全方向に一度に計算する Unsplit Time Integration を時間方向に多段階で行う
Fractional Step Method(25)を採用した．

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

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

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

(20)

vip 2  vip1 (27)

t n 1 2  t n  3 2  t n 1 2 nBn  1 

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

 m n p n
t n 1 2
 mn gi   
vip1  vin  n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

 

(26)

1 
 n
 2mnˆ n Dijn   x  mn ija n   Qin  mn ijnv njˆ x jj(24)

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

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

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

(33)

f i ob n n 1/ 2tmnm  C A 1  m n  ˆnn

vip 2  1  Fobn  vip1  Fobn viob n

Bn

なお，求められた構造物の位置に基づく Fob の更新には，Yuki ら(26)の双曲線正接関数を用いる手法とは異なり，各計算格子を各方向に数分割ずつした小格子を考え，その分割した小格子の中央か(25)

なお，3 次精度 Adams-Bashforth 法の導入に伴う A0ni ， A1ni ， A2ni と

3Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4ら構造物を構成する面までの垂線の長さを基に Fob を算出することで構造物の総体積の増減を抑えた．以上をまとめると，IBM での計算の手順は次のようになる．
1. f i ex （式(19)）と f i re を求める．
2. 式(30)より uiob を，式(31)より iob を求める．
3. 式(32)より xiob を，式(33)より iob を求める．
4. 求めた xiob と iob より， Fob を更新する．

z   0.9

(34)



sin 2  tan 2 






2
s



cos  sin  


s


(40)

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

(41)

1
FD   CD1 w d 502 vr28(42)

FD cos 1  W sin  cos      dW cos (43)

FD sin 1  W sin  sin    (44)となる．また， vb ， Cv f v f ， vr の関係より，

２．２．２．１ 掃流砂量式
ここでは，単位幅，単位時間当たりの i 方向の掃流砂量 qi が砂粒子の平均輸送速度の i 方向成分 vbi に比例すると考えた
Engelund・Fredsøe(28)による掃流砂量式

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

vr cos 1  vb  Cv f v f cos (46)と書かれる．ここで，式(43)から式(46)に対して未知数は vb ，vr ，
 ，1 の 4 つであることから，vb とその方向  ，すなわち vbiが求められ，最終的に式(35)より掃流砂量 qi が求められる．

(35)を採用した．ただし， pEF は掃流形態で移動する砂粒子の割合であり，次式で与えられる．

２．２．２．２ 斜面崩壊計算
Roulund ら(27)は，斜面の傾斜角が砂粒子の安息角 r より微小角度 r だけ上回ったときに斜面が崩壊し始め， r より微小角度 r
だけ下回るまで崩壊し続けることを水理実験で確認している．この現象を基に，Roulund ら(27)は以下の手法を提案している．
斜面の傾斜角  が  r   r を上回ったとき，その全砂粒子（ pEF  1 ）が静水中（ v f  0 ，  1  180 ）で一定速度 vb で

(36)

ここで，  d は砂粒子の動摩擦係数（  tan  d ； d ：動摩擦角），
 * は Shields 数， *c は斜面での限界 Shields 数であり， * は

v 2f


s

と書かれる．ここで，Fig. 3 に示した砂粒子が各タイムステップ
間で一定速度 vb で移動すると仮定すると，砂粒子の移動方向およびそれと直交方向の運動方程式はそれぞれ

２．２．２ 掃流砂輸送モデル
掃流砂輸送モデルは，Roulund ら(27)を参考に構築した掃流砂量
式と斜面崩壊計算から構成される．

* 


s

により評価した(29)．ここで， *0 は水平面での限界 Shields 数， s
は砂粒子の静止摩擦係数（  tan  s ；  s ：静止摩擦角）である．
続いて，式(35)の vbi を求めるために，Fig. 3 に示す掃流形態で移動する砂粒子 1 個を考える．ここで，vb は vbi の大きさ，Cv f v fは砂粒子の位置での流速（ Cv f ：摩擦速度の係数）， vr は vb に対する砂粒子の位置での相対流速，  は斜面下向きに対する v f の方向， は斜面の傾斜角， は vb に対する v f の方向，1 は vb
に対する vr の方向である．また，W は砂粒子の水中重量， FD は砂粒子に作用する抗力を表し，それぞれ

ここで， z s は基準面からの底面の高さ， qx と q y はそれぞれ x ，
y 軸方向の掃流砂量， pN は底質の巻き上げ関数， qzbS は底面での z 軸方向の沈降による浮遊砂輸送フラックスである．

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


s

 *c   *c 0  cos  1 

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

1
qi   d 50 pEF vbi
6s

また，式(36)の *c は

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

z s
1  qx q y pN  qzbS   0

t 1  m  xy k  k exp   k 6   4.535  k  2000 (39)(37)

 s  1 gd50で定義される．ただし， v f は摩擦速度， s は砂粒子の比重（   s  w ；  s ：砂粒子の密度）であり， v f は

vsurf
vf

z

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



 2





(38)

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

Fig. 3 Force balance on a moving sand particle
in bed-load motion on a sloping bed
4

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4れる捨石や消波ブロックを想定している．
底面の境界条件として，Nielsen ら(32)に倣って，底面での z 軸方向の拡散による浮遊砂輸送フラックス qzbD が底質の巻き上げ関数
pN に等しいとする条件を用い，その pN として

斜面下向き（  180 ）に崩壊すると仮定する．式(43)から式(46)に，式(41)，式(42)， v f  0 ，    1  180 を代入して整理すると，斜面崩壊に伴う砂粒子の平均輸送速度 vbslide は

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

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

qislide 


3

d 50

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

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

(49)

C n 1  C n 

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

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

z sn 1  z sn 

6.

(51)

C
xi

(52)

qiS  mwsi C(53)

qiD  m s

(56)

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

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

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

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

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

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

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

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

(55)

(48)

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

m

if  *   *cを用いた(33)．ここで， C p は底質の巻き上げに関する無次元パラメータである．

また，斜面崩壊に伴う底質の質量保存式は，斜面崩壊時には浮遊
砂の巻き上げや沈降が生じないと仮定すると，式(34)に 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.8
C p 
 w0.2   *c t h 1/ 2  qxslide h q yy
1  m n  x

slide hz sh 1  z sh 

iii.
(54)





(58)

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

ここで， d* は砂粒子の無次元中央粒径（  sgd 502  w2 ）である．
なお，ここでは透過性材料として内部への浮遊砂の侵入が考えら

２．３ XNS・VFM・IBM・STM 間のカップリング手法
FSSM では，XNS，VFM，STM，IBM の接続に Fig. 4 に示す
Two-Way カップリング手法を採用した．具体的には，FSSM での計算の手順は次のようになる．
1. XNS を実行する．
2. XNS で得られた流速場を用いて VFM を実行する．
3. 流速場に加えて，VFM で得られた気液界面の位置を用いてSTM を実行する．
4. 流速場と気液界面の位置に加えて，STM で得られた地形変化を用いて IBM を実行する．
5. VFM，STM，IBM で得られた値をフィードバックして XNSを実行する．
6. 手順 2 から 5 を計算終了時刻まで繰り返す．

Fig. 4 Two-way coupling procedure

5

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4
過水平床表面，奥側の側面に Slip 条件を，沖側境界に Sommerfeld
放射条件を，手前側の側面と岸側境界に勾配ゼロの条件を，上面に圧力一定の条件を課した．気液界面を追跡する VFM では，VOF
関数 F の境界条件として，全境界に勾配ゼロの条件を用いた．
STM では，まず掃流砂の境界条件として，砂地盤の沖側端と岸側端からの掃流砂の供給がゼロの条件を課した．
浮遊砂については，
砂地盤の表面に式(55)による境界条件を用い，計算領域の沖側端と岸側端からの浮遊砂の供給はゼロとした．また，式(38)による砂地盤の表面近傍における z 軸方向の格
摩擦速度 v f の計算には，
子幅 2 mm の 1/2 として， z surf  1 mm での流速 vsurf を用いた．
なお，本論では移動する構造物は取り扱わないことから，構造物の運動を取り扱う IBM は使用しなかった．
Table 1 と 2 に，XNS と STM で用いたパラメータをそれぞれ示
す．XNS の LES に関して，グリッドフィルタとテストフィルタ
は Horiuti(14)に倣ってトップハットフィルタ（ボックスフィルタ）を採用し，グリッドスケールのフィルタ幅  は各位置における各方向の格子幅とした．ここで，フィルタ幅比 は Garmano ら(34)に倣って設定した．その他，付加質量係数 C A ，非線形抵抗力係
数 CD1 ，線形抵抗力係数 CD 2 は水谷ら(13)の実験結果に基づいて，減衰係数 C は中村ら(35)による試行計算の結果に基づいて設定した．Table 2 については，摩擦速度の係数 Cv f ，水平面での限界
Shields 数 *c 0 ，静止摩擦係数  s ，動摩擦係数  d ，van Driest 型減衰係数 A ，砂粒子の水中安息角 r ，斜面崩壊開始基準 r ，斜面
崩壊完了基準 r は Roulund ら(27)に基づいて，底質巻き上げ係数
C p は Nielsen(33)に基づいて設定した．

３．計算条件
Nakamura ら(6)の水理実験をモデル化した Fig. 5 に示す計算領域
を用いて解析を行った．同図に示すように，計算領域の岸側端に
高さ 0.33 m の不透過水平床と，その沖側に 1/10 勾配の不透過一様斜面を設定し，その水平床上に不透過護岸（高さ 0.15 m，天端
沖側法面勾配 1/0.2），
砂地盤（高さ 0.15 m，
長さ 1.00 m），
幅 0.03 m，
不透過水平床（高さ 0.15 m，長さ 1.00 m）を沖から順に設定した．
また，構造物は幅 B の矩形とし，護岸背後から 0.28 m の位置に固定した．なお，Nakamura ら(6)の水理実験では y  0 を中心に対称な実験装置が用いられていることから，津波の遡上とそれに伴う地形変化も y  0 に対してほぼ対象であると考え，Fig. 5 に示すように y  0 の片側のみを対象とした．
計算格子は，構造物の沖側隅角部の周辺，かつ砂地盤の表面近傍に対して 10×10×2 mm 角の等間隔格子を採用し，計算負荷の軽減のためにその範囲から離れるにしたがって x ， z 軸方向に格子幅を徐々に広げた不等間隔格子を用いた．流体の解析を行う
XNS では，境界条件として，底面，不透過斜面，護岸表面，不透４．計算結果および考察
本稿では，津波をモデル化した入射波として，Nakamura ら(6)と同様に押し波のみの長周期波 1 波を考える．また，Nakamuraら(6)はその入射波を用いた実験を計 93 ケース行っているが，ここでは沖波波高 H o  64 mm，造波板を押し出す時間 T  6.0 s，静水深 ho  0.455 m，構造物の幅 B  0.14 m，砂粒子の中央粒径
d 50  0.2 mm のケースを対象に再現計算を行った．なお，砂地盤を構成する砂粒子の密度  s は 2.65×103 kg/m3，空隙率 m は 0.40とし，試行計算の結果より浮遊砂の乱流拡散係数  s を 1.0×
10-3 m2/s，等価砂粗度 k s を 0.5d 50 とした．
Fig. 6 に遡上津波の変形の様子を示す．ここで，左図は実験結果を，右図は計算結果を表す．また，右図に示す計算結果では，50.0 kg/m3 ごとの浮遊砂濃度の等値面も示した．
実験結果を示した
Fig. 6 左図より，護岸を越流した津波が砂地盤上を進行し（Fig. 6(a)，(b)），構造物の沖側面で上方に打ち上げられる様子（Fig. 6(c)）が，確認できる．
その後，打ち上がった波が沖向きに砕波し（Fig. 6(d)）
それが後続の越流波の上に乗りつつ徐々に沖側に伝播していく様子（Fig. 6(e)）も確認できる．そして最終的に，Fig. 6(f)に示すよ
うに水位が低下し，津波が治まっていくことが分かる．同様の現象は，同図右図に示す計算結果からも認められることから，本数値モデル FSSM の妥当性が定性的に確かめられた．
Fig. 7 に護岸前面の水位変動 の比較を，Fig. 8 に砂地盤内部の間隙水圧変動 p の比較を示す．ここで，W2 から W5，P1 から P8の位置は Nakamura ら(6)を参照されたい．Fig. 7 より，W5 で水位の上昇を若干過小評価しており，またその後の引き波に伴う水位の低下をすべての場所で過小評価する傾向があることが分かる．
一方，Fig. 8 より，特に P3 で間隙水圧の上昇を若干過小評価していることが分かる．しかし，Fig. 7 の水位変動，Fig. 8 の間隙水圧変動ともに，実験値と良く一致していると判断できることから，
FSSM の妥当性が定量的にも確かめられた．

Fig. 5 Computational domain
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
Table 2 Parameters of the STM
10.0
摩擦速度の係数 Cv f
底質巻き上げ係数 C p
0.00033
水平面での限界 Shields 数 *c 0 0.05
静止摩擦係数  s0.63
動摩擦係数  d
0.51
Kármán 定数  
0.4
Van Driest 型減衰係数 A
25.0
砂粒子の水中安息角 r
32.0º
斜面崩壊開始基準 r
2.0º
斜面崩壊完了基準 r
2.0º
6

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4

(a)

(b)

(c)

(d)

(e)

(f)

(g)
Fig. 6 Snapshots of the run-up tsunami (left: experimental data(6); right: numerical results)
7

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4(a) Experimental data(b) Numerical result computed using the present model (FSSM)
Fig. 7 Comparison of water surface elevation (c) Numerical result computed in Nakamura et al.(7)Fig. 9 Comparison of final scour depth z sf
洗掘が生じているのみであり，再現性はあまり高くないことが分かる．これは，洗掘の大きさに対して格子が十分小さくなかった
ために，この付近に生じる局所的な現象を十分に解像できなかったためであると考えられる．ただし，Fig. 9(b)では構造物の沖側隅角部から側面にかけて洗掘が連続して生じており，またその側面での洗掘が構造物の比較的近くに生じていることから，Fig. 9(a)に示す実験結果とのある程度の傾向の一致が認められることが分かる．しかし，Fig. 9(c)では構造物から離れた y  0.12 m 付近にも洗掘が生じてしまっており，洗掘の傾向を適切に評価できていないことが分かる．続いて， y 軸に平行な帯状の侵食が生じる護岸のすぐ岸側（0  x  0.15 m）については，Fig. 9(b)からは若干過大評価しているものの同様の侵食パターンが確認できる一方で，Fig. 9(c)ではわずかに侵食が生じているのみであり，再現できていないことが分かる．また，構造物の岸側については，全体的に過大評価しているものの，陰になっている構造物のすぐ岸側で堆積が生じ，陰になっていないその横側で侵食が生じるという傾向がFig. 9(a)と(b)で一致していることが確認できる．しかし，Fig. 9(c)では構造物の岸側隅角部の近傍にも洗掘が生じてしまっており，傾向があまり一致していないことが分かる．以上より，構造物の沖側隅角部の洗掘の再現性に関しては課題を残しているものの，流体と地形変化の相互作用を考慮することによって再現性が全体的に向上することが確認できた．

Fig. 8 Comparison of pore-water pressure p
Fig. 9 に最終洗掘深 z sf の比較を示す．
ここで，正の z sf は侵食を，
負の z sf は堆積を表す．また，同図では，流体と地形変化の相互作用を考慮していない Nakamura ら(7)による計算結果も同時に示した．まず，最も大きい地形変化が現れる構造物の沖側隅角部については，Fig. 9(b)と(c)に示すように，いずれの計算結果もわずかに５．結論
本研究では，遡上津波による陸上構造物周辺の洗掘現象に関す
る水理実験(6)を対象に，中村ら(8)が開発した流体，構造，地形変化の 3 者間の相互作用を解析できる 3 次元流体・構造・地形変化連8Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-4
リーNo. 12 (2001), 457 p.
(17) 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.
(18) Eric Cruz, 横木裕宗, 磯部雅彦, 渡辺 晃, “非線形波動方程式に対する無反射境界条件について,” 海岸工学論文集, 40(1993), pp. 46-50.
(19) 岩田隆一, 梶島岳夫, 竹内伸太郎, “Immersed-boundary andvolume-of-fluid カップリング法による粒子列の間を通過する
上昇気泡の解析,” 混掃流, 23, 5 (2010), pp. 531-538.
(20) Xiao, F., Yabe, T., Ito, T., and Tajima, M., “An algorithm for simulating solid objects suspended in stratified flow,” Comp. Phys.
Communications, 102 (1997), pp. 147-160.
(21) 本多和彦, 富田孝史, 西村大司, 坂口 章, “多数の津波漂流物を解析する数値モデルの開発,” 海洋開発論文集, 25 (2009),
pp. 39-44.
(22) Amsden, A. A. and Harlow, F. H., “A simplified MAC technique forincompressible fluid flow calculation,” J. Comp. Phys., 6 (1970), pp.
322-325.
(23) 梶島岳夫, “乱流の数値シミュレーション,” 養賢堂 (1999),
225 p.
(24) 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.
(25) 秋山 守, 有冨正憲, “新しい気液二相流数値解析－多次元流
動解析－,” コロナ社 (2002), 261 p.
(26) 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.
(27) 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.
(28) Engelund, F. and Fredsøe, J., “A sediment transport model forstraight alluvial channels, Nordic Hydrology, 7 (1976), pp. 293-306.
(29) Brooks, H. N., “Discussion on “Boundary shear stresses in curved
trapezoidal channels” by Ippen, A. T. and Drinker, P. A.,” J. Hydraulic Div., ASCE, 89, 3 (1963), pp. 327-333.
(30) 榊山 勉, 伊藤 裕, 小林正典, 板井雅之, 小林英次, “傾斜堤を透過する浮遊砂の濃度評価への数値計算モデルの適用,”
海岸工学論文集, 51 (2004), pp. 456-460.
(31) 中山哲嚴, 牧野弘幸, 新井雅之, 小林 学, 佐藤勝弘, “導水中に含まれる浮遊砂量を低減するための海水導入工の構造に関する研究,” 海岸工学論文集, 54 (2007), pp. 661-615.
(32) 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.
(33) Nielsen, P., “Coastal Bottom Boundary Layers and SedimentTransport,” Advanced Series on Ocean Eng. Vol. 4, World Scientific(1992), 340 p.
(34) 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.
(35) 中村友昭, 水谷法美, 藤間功司, “遡上津波の変形と陸上構造物に作用する津波力に関する 3 次元数値シミュレーション,”
ながれ, 29, 2 (2010), pp. 107-117.

成数値計算モデルを適用して再現計算を行った．その結果，遡上津波の変形過程の様子を実験結果と比較することで定性的な妥当性をまず確認するとともに，水位変動および間隙水圧変動の比較からモデルの妥当性を定量的にも明らかにした．また，津波作用後の最終的な地形を対象に実験結果との比較を行い，相互作用を考慮しなかった場合と比べて地形変化の再現性が向上することを明らかにした．ただし，構造物の沖側隅角部の洗掘が十分に再現できていないなどの課題を残していることから，Nakamura ら(6)(7)
が指摘しているように地盤内部の応力変動の影響を考慮するなど，
さらなる検討を行っていく所存である．
参考文献(1) 西村仁嗣, 堀川清司, “遡上津波の戻り流れによる護岸法先の洗掘,” 海岸工学講演会論文集, 26 (1979), pp. 210-214．
(2) 宇多高明, 小俣 篤, 横山揚久, “津波遡上に関する実験報告
書－構造物の効果を中心として－,” 土木研究所資料, 2486
(1987), 122 p.
(3) 野口賢二, 佐藤愼司, 田中茂信, “津波遡上による護岸越波および前面洗掘の大規模模型実験,” 海岸工学論文集, 44 (1997),
pp. 296-300.
(4) 後藤仁志, 酒井哲郎, 林 稔, 織田晃治, 五十里洋行, “遡上津波の戻り流れによる護岸法先洗掘のグリッドレス解析,” 海
岸工学論文集, 49 (2002), pp. 46-50.
(5) Tonkin, S., Yeh, H., Kato, F., and Sato, S., “Tsunami scour around acylinder,” J. Fluid Mech., 496 (2003), pp. 165-192.
(6) Nakamura, T., Kuramitsu, Y., and Mizutani, N., “Tsunami scour
around a square structure,” Coastal Eng. J., JSCE, 50, 2 (2008), pp.
209-246.
(7) Nakamura, T., Mizutani, N., and Yim, S. C., “A three-dimensional
coupled fluid-sediment interaction model with bed-load/suspendedload transport for scour analysis around a fixed structure, J. Offshore
Mech. and Arctic Eng., ASME, 131, 3 (2009), 031104, 9 p.
(8) 中村友昭, Solomon C. Yim, 水谷法美, “三次元流体・構造・地形変化連成数値計算モデルの開発と没水物体周辺の洗掘現象への適用”, 土木学会論文集 B（海岸工学）, 57 (2010), pp.
406-410.
(9) 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.
(10) 功刀資彰, “自由界面を含む多相流の直接数値解析法,” 機械学会論文集 B 編, 63, 609 (1997), pp. 1576-1584.
(11) Kajishima, T., Takiguchi, S., Hamasaki, H., and Miyake, Y., “Turbulence structure of particle-laden flow in a vertical plane channel
due to vortex shedding,” JSME Int. J. Ser. B, Fluids and ThermalEng., 44, 4 (2001), pp. 526-535.
(12) Brackbill, J. U., Kothe, B. D., and Zemach, C., “A continuum method for modeling surface tension,” J. Comp. Phys., 100 (1992), pp.
335-354.
(13) 水谷法美 , 前田健一郎 , Ayman M. Mostafa, William G.
McDougal, “透水性構造物の抵抗係数の評価と波・潜水透水性
構造物の非線形相互作用の数値解析,” 海岸工学論文集, 43(1996), pp. 131-135.
(14) Horiuti, K., “A new dynamic two-parameter mixed model forlarge-eddy simulation,” Phys. Fluids, 9, 11 (1997), pp. 3443-3464.
(15) 岩田好一朗, 川崎浩司, 金 度三, “水中構造物による砕波の数値解析,” 海岸工学論文集, 42 (1995), pp. 781-785.
(16) 財団法人沿岸開発技術研究センター , “ 数値波動水路（CADMAS-SURF）の研究・開発,” 沿岸開発技術ライブラ

9

Copyright © 2010 by JSFM