長岡工業高等専門学校研究紀要

論第 52 巻（2016）

文

波動水槽内の浅水流再現シミュレーション－ SUPG法による解析結果と水位計測値の比較 －倉橋 貴彦1・野上 雅人2・今井 伸哉3・吉荒 太一4・衛藤 俊彦5
1 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
2 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
3 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
4 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
5 環境都市工学科（Department of Civil Engineering, National Institute of Technology, Nagaoka College）

Numerical Simulation of Shallow Water Flow in Wave Tank
- Comparison of Numerical Results by SUPG method and Observed Water Elevation Takahiko KURAHASHI1, Masato NOGAMI2, Shinya IMAI3,Taichi YOSHIARA4 and Toshihiko ETO5
Abstract
In this study, we present numerical experiments of shallow water flow based on the SUPG method. As
the governing equation, the shallow water equation is employed, and remarks on the effect of thediscretization method and the nonlinear term are carried out. In addition, comparison of the numericalresults and the measurement result based on the experiment using the wave tank is carried out, and theapplicability of the numerical results for simulation of shallow water flow is investigated.
Key Words : finite element method, shallow water flow, SUPG method

１．序論整備が行われている．本論文では，浅水長波方程式の非線形項の有無や，数値計算法の適切性を検証す

一般に，河川の流れや海の流れ等，水波の波長に比べて水深が十分に深い場合，浅水長波流れとしてるために，波動水槽 2)により得られた水位計測値と数値解析結果を比較し，考察を行う．

取り扱い，浅水長波方程式による流れのシミュレーションが行われる．浅水長波方程式は偏微分方程式の分類としては双曲型の偏微分方程式であり，数値的に解く場合は，計算安定性に十分留意する必要が

２．SUPG法およびθ法による浅水長波方程式の離散化ある 1)．シミュレーションの対象としては，海岸の方であれば，津波や高潮，河川の方であれば洪水を対象とした数値予測が行われ，近年では，各自然災
害に対する危険箇所把握のため，ハザードマップの本論文で対象とする浅水長波方程式の運動方程式，連続式を式(1)，(2)にそれぞれ示す．式(1)，(2)はEinsteinの総和規約により表示している．

- 21 -

倉橋 貴彦・野上 雅人・今井 伸哉・吉荒 太一・衛藤 俊彦

ui  u j ui , j  g ,i  0(1)ui  N1ui1  N 2ui 2  N 3ui 3   ui (10)  N11  N 22  N33  T (11)T  hui ,i  0(2)ここに，uiはx,y方向流速，ηは水位変動量，gは重力T加速度，hは基準水深を示す．初期条件および境界

 *  N11*  N 22*  N33*  T   

条件は式(3)，(4)の様に定義する．

ui t 0   uˆiin (3)(4)

空間方向に対する離散化にはSUPG法 3) を用いる．
SUPG法では，重み関数を移流流速の影響も考慮しここに，N1，N2，N3は形状関数を示す．この補間関および式(15)に示す有限要素方程式が得られる．

M ui  u j S j ui  gSi   uk M  ui   uk u j H jk ui   uk g H ki    0
kた形で定義し，式(5)，(6)のように重み付き残差方(14)M  h Si ui   uk M    uk h H ki ui   0

程式を誘導する．

 u  u u u  u u  g d  0 (5)(6)   u    hu d  0*i*i ,kki*ke*,kji, j(13)数を重み付き残差方程式に代入することで，式(14)u x  uˆ x , u y  0 on 1u y  0 on 2e  (12)ui*  N1ui*1  N 2ui*2  N3ui3   ui

(15)

k,ii ,i

ここに，変数の上に付けたバーは要素内平均値を示しており，また式(14)および式(15)における係数行列は式(16)のように書くことができる．

ここに，τは安定化パラメータを示し，以下の式に

M   T dS    T d,i i  S   T d,j j T M  k    , k  dT H jk    , j  , k  dH ki     , k  ,i T dより与えられる．
12e

 22
  2   2 u i        t   he     (7)nuie   uei2(8)i 1(16)

ここに，Δtは時間増分量を示す．また，heは要素長を示し，式(9)で与えられ，式(9)中のAe は要素の面積を示す．

次に，式(14)および式(15)に対してθ法 4)を適用し，時間方向に離散化をする．θ法による離散化に際し，

he  2 Ae(9)式(14)および式(15)における各変数を式(17)～(20)のように示す．

x,y方向流速ui，水位変動量ηまたそれぞれの変数に対する重み関数ui*，η*の補間関数として，三角形一次 要 素 に 対 す る 補 間 関 数 を 適 用 す る （ 式 (10) ～(13)）．
- 22 -ui   ui  ui n 1n

t(17)

波動水槽内の浅水流再現シミュレーション

ui    uin1 1   uin (18)  (19)n 1  nt

図-1

    n1 1    n (20)

式(17)～(20)を式(14)および式(15)に代入することにより，空間方向および時間方向に離散化した有限要素方程式を得る（式(21) および式(22)）．

M   u M u  t u S   u u H u  tg S   u H   M   u M u  t 1   u S   u u H u  tg 1   S   u H  n 1i

kkjjkj

有限要素分割図

表-1 解析条件
時間増分量 Δt[s]0.001

パラメータ θ時間ステップ

1.0
3000

基準水深 h[m]重力加速度 g[m/s2]
x 方向分割数
y 方向分割数

100
9.81
1002n 1i

jkn 1ikkki

kni

jjkj(21)ni

jkni 0kM   u M  kkkin 1  th  S i   u k H ki  uin 1  

 M   u k M  k  n(22)  th 1   Si   u k H ki  uin 0

t=0.0s における水位・流速ベクトル分布（図-2-1）

t=0.1s における水位・流速ベクトル分布（図-2-2）

３．数値実験
３．１

段波の解析

図-1 に示す有限要素メッシュを用いて，段波の解析を行った．解析条件を表-1 に示す．x=0m，10m の境界では，ux=0m/s，y=0m，0.5m の境界ではuy=0m/s の条件を与えた．水位変動量 η の初期条件
としては，x=5m を境に，右側で－1m，左側で＋1mの条件を与えた．解析結果を図-2 に示す．結果より，時間経過とともに，水位変動量の波形が発散することなく安定に計算できていることがわかる．

- 23 -倉橋 貴彦・野上 雅人・今井 伸哉・吉荒 太一・衛藤 俊彦

３．２

波動水槽における水位の計測結果と有限要素解析結果の比較

次に，波動水槽における水位の計測結果と有限要素解析結果の比較を示す．図-3および図-4に，水槽の写真と計測点の位置を示す．また，水槽の左端より 表 -2 に 示 す 条 件 に 基 づ い た 波 を 生 成 し ，No.1~No.3の点において時間経過に対する水位変動量を計測する．また，計測開始は，波が安定した状態から計測を行った．各点における計測結果を図-5に示す．

t=0.2s における水位・流速ベクトル分布（図-2-3）

図-3

No.3

水槽の写真

No.2

No.1
Inflow direction

0.50m

y0.20m0.25mx2.00m

2.00m
4.00m

t=0.3s における水位・流速ベクトル分布: Observation point（図-2-4）図-4観測点の位置表-2 実験条件

各時刻における水位分布の比較
図-2（図-2-5）段波の解析結果- 24 -

基準水深 h[m]0.46

スライダーボリューム
周期 T[s]5.0
5.0

計測間隔[s]
計測時間[s]0.0110

波動水槽内の浅水流再現シミュレーション

図-5 水位変動量の計測結果

図-7

流入境界条件

ここで，実験モデルに対する再現シミュレーショ

図-8にθ=1.0とした完全陰解法による解析結果（2ンを行う．図-6に計算モデルのメッシュ図，表-3に秒おきにおける水位および流速ベクトル分布の鳥瞰解析条件をそれぞれ示す．解析領域は水槽モデルと図，各観測点における水位の経時変化）を示す．ま同じにし，長さ25m水槽を表したモデルとする．モ

た，図-9にθ=0.0としたEuler法の各観測点におけるデルの右端の各節点に，実験のNo.1における計測デ水位の経時変化を示す．（※θ=0.0とした場合の結ータは0.01秒間隔であるため，解析の時間増分量をΔt=0.001[s]に合わせて各計測値間を10等分したもの

果においては，不安定な計算結果となっているため，を流入境界条件として適用する．図-7に流入部にお果は，節点上の結果を線形補間することにより求めける水位変動量ηの境界条件を示す．初期条件は，たものである．結果よりθ=1.0とした完全陰解法の解析領域においてη=0と設定し解析を行う．

結果はθ=0.0とした陽解法の結果に比べ，安定に計鳥瞰図は除外する．）観測点における水位の解析結算が行えていることがわかる．次に，観測点No.2，No.3における水位の計測結果とθ=1.0とした完全陰解法の比較を図-10および図-11に示す．結果として，水位に若干差はあるが，良い一致をした結果が得られていることがわかる．

図-6 メッシュ分割図
表-3 解析条件
時間増分量 Δt[s]パラメータ θ時間ステップ
基準水深 h[m]0.001
完全陰解法：1.0
Euler 法：0.0
完全陰解法：10000
Euler 法：7000
0.46

重力加速度 g[m/s2]
節点数

9.81753

要素数1000

t=0.0s における水位・流速ベクトル分布- 25 -（図-8-1）

倉橋 貴彦・野上 雅人・今井 伸哉・吉荒 太一・衛藤 俊彦

t=2.0s における水位・流速ベクトル分布（図-8-2）t=8.0s における水位・流速ベクトル分布（図-8-5）t=4.0s における水位・流速ベクトル分布（図-8-3）t=10.0s における水位・流速ベクトル分布（図-8-6）

各観測点における水位の経時変化（図-8-7）図-8 解析結果（θ=1.0とした完全陰解法）
t=6.0s における水位・流速ベクトル分布（図-8-4）- 26 -波動水槽内の浅水流再現シミュレーション波の問題に対する数値実験を行い，波動水槽のモデルに対する水位の計測結果と数値計算結果の比較を行った．結果として，段波の問題では，不連続な水面波形を有する初期値問題にも関わらず，発散せず安定に解析を行えることを確認できた．次に，波動水槽モデルに対して，実験において得られた水位の計測結果と数値解析の結果の比較を行った．結果として，θ=1.0とした完全陰解法とθ=0.0としたEuler法の比較を行ったところ，θ=1.0とした完全陰解法では安定に解析を行うことができることを確認できた．
また，実験において得られた計測結果と解析結果の比較では，水位の値に若干の差はあるが良い一致を図-9 解析結果（θ=0.0としたEuler法）示す結果を得ることができた．本検討では，粘性項や摩擦項を取り除いた解析を行っているため，それら各項の影響を考慮するとより良い結果が得られるものと考えられる．

謝辞：本論文を執筆するにあたり，科学研究費補助金（基盤(C)）15K05786の援助を受けた．ここに謝意を表す．また，本論文の解析結果は，九州大学情報基盤研究開発センターの高性能アプリケーション
サーバSR16000を使用し計算を行ったものである．
計算機の使用について，センター関係者の方へ謝意
図-10を表す．

No.2における計測結果と解析結果の比較

参考文献
1)

竹内 則雄, 樫山 和男, 寺田 賢二郎，”計算力学 有限要素法の基礎”，森北出版, 2003.

2)

倉橋 貴彦, 吉荒 太一, 衛藤 俊彦, ”波動水槽による水位計測値を用いた水槽内における流れ場の推定解析－ カルマンフィルタ理論を用いた有限要素解析による流れ場の推定アルゴリズムの紹介 －”, 長岡工業高等専門学校紀要，Vol.51, pp.1-9, 2015.

3)

Tezduyar,T.E., ”Stabilized finite element formulations forincompressible flow computations”, Advanced in AppliedMechanics, Vol.28, pp.1-44, 1991.

4)
図-11

No.3における計測結果と解析結果の比較

川原睦人，”有限要素法流体解析”，日科技連出版社,
1985．
(2016. 8. 1 受付)

４．結論
本研究では， 移流項を考慮した浅水長波方程式
に対してSUPG法およびθ法による離散化を行い，段- 27 -