修士論文要旨（2011年度）
浅水長波流れ解析のための Discontinuous Galerkin 有限要素法に関する研究
Studies on the Discontinuous Galerkin FEM for Shallow Water Flow Analysis
土木工学専攻 33号 牧野 優作 Yusaku MAKINO
1. はじめに
河川の氾濫，高潮，津波などの数値シミュレーション には，浅水長波方程式が広く用いられている．そのため 浅水長波方程式を高精度に解くことは水害による被害を 予測する上で重要である．偏微分方程式の代表的な離散 化解析手法には，差分法，有限体積法，有限要素法などが 一般に用いられているが，浅水長波流れ問題は，複雑な地 形を計算対象領域とするため，任意形状の適合性に優れ た有限要素法や有限体積法などは有効な手法と考えられ る．これまで，有限要素法においては，重み関数と近似関 数の関数値が連続となるCG(Continuous Galerkin)法 に基づく有限要素法が一般的であったが，近年，要素間 で関数値の不連続を許容し，要素境界において局所的な Fluxの収支を考慮して解析を行う，DG(Discontinuous Galerkin)法に基づく有限要素法1) が注目され，特に不 連続性のある問題に対して有効に用いられている．
そこで本研究は，浅水長波流れの高精度な数値解析手 法の構築を行うことを目的とし，浅水長波方程式の離 散化手法にDG法を適用する．なお，要素形状として は三角形で，補間関数としては1次関数(三角形1次要 素)を用いる．陸域における移動境界手法としては固定 メッシュを用いるEuler的手法を用いる．数値解析例と して段波問題，ダムブレイク問題を取り上げ，厳密解お よび従来のCG法(SUPG法に基づく安定化有限要素 法)2) との結果の比較によりDG法の有効性について検 討する．
2. 数値解析手法
2.1 支配方程式支配方程式には，以下に示す浅水長波方程式を用いる．
∂U
∂t +∇ ·F(U) =S(U) (1) ここで，Uは保存変数，F(U)は流束関数，S(U)はソー ス項であり,以下のように定義される．
U= [ h, uh, vh ]T, (2)
F(U) = [ F1(U), F2(U) ]T, (3)
F1(U) =
uh u2h+12gh2
uvh
,F2(U) =
vh uvh v2h+12gh2
図– 1 座標系
㑆㑐ᢙ ⷐ⚛ౝㇱ㗔ၞ
ⷐ⚛Ⴚ⇇
図– 2 1次要素モデル(1次元)
S(U) =
0
−gh∂z∂x
−gh∂z∂y
(4)
ここで，hは全水深，u，vはx，y軸方向の断面平均流 速，gは重力加速度であり，zは基準面からの高さであ る（図-1参照）．
2.2 DG法による空間方向の離散化1)
DG法は，1次要素を用いる場合，要素間で関数値の C0不連続を許容しており，要素境界においてFluxの 収支を考慮して解析を行う手法である．図-2にDG法 の1次要素モデルを示す．式(1)に対してDG法を適用 すると，以下の弱形式が得られる．
Z
e Ω
Wh· ∂Uh
∂t dΩ− Z
e
Ω∇Wh·F(Uh)dΩ +
Z
∂Ωe
Wh·( ˆF(Uh)·n)dΩ = Z
e Ω
Wh·S(Uh)dΩ (5) ここで，Wh= [whh, whqx, whqy]Tは要素ごとに定義され る不連続な重み関数である．なお，補間関数としては一 次関数を用いる．また，Ωe は要素内部領域の集合，∂Ωe は要素境界の集合，n= [nx, ny]Tは∂Ωeの外向き法線 ベクトル，式(5)の左辺第3項はフラックス項であり，
次節で説明を行う．
2.3 数値フラックス
本論文では数値フラックスとして，次式で定義される Local Lax-Friedrichs Flux3) を用いる．
F(Uˆ h)·n
=1
2[F(Uh−) +F(Uh+)−amax¡
Uh+−Uh−¢ ] (6) ここで，amaxは以下のように定義される．
amax= max(|u+nx+v+ny−c+|,
|u−nx+v−ny−c−|,
|u+nx+v+ny+c+|,
|u−nx+v−ny+c−|) (7)
c+=p
gh+, c−=p gh−
また，f±は，以下のように定義される．
f±|∂Ωe= lim
ϵ→0f(x±ϵn) (8) 2.4 時間方向の離散化
式(5)を三角形1次要素を用いて離散化すると，以下 の時間に関する連立常微分方程式が得られる．
dU
dt =L(U) (9)
そして，上式に対して，3次精度陽的Runge-Kutta法 を適用すると，以下の離散化式が得られる．
U(1) =Un+ ∆tL(Un) (10)
U(2) = 3 4Un+1
4
³
U(1)+ ∆tL(U(1))´
(11)
U(n+1)=1
3Un+2 3
³
U(2)+ ∆tL(U(2))´
(12)
2.5 Slope Limiter処理
上記で示した離散化手法により解を求めると，不連続 面においてオーバーシュートやアンダーシュートが生じ る．このような数値不安定性を回避するために，Slope Limiter処理4) を導入する．この処理をRunge-Kutta 法の各ステップ後に行う．
2.6 移動境界手法5)
津波遡上や洪水氾濫などの現象を解く場合，移動する 水際境界を考慮する必要がある．陸域の移動境界手法 として，移動メッシュを用いるLagrange的手法と固定 メッシュを用いるEuler的手法に大別される．本研究で は，Euler的手法を用いる．
本手法では，陸域において微小な水深を設ける．これ より，計算領域全体で浅水長波方程式を解き，得られた 水深の大きさから陸水判定を行う．要素の平均水深が微
ಣℂ೨ ಣℂᓟ
図– 3 水際線の処理（ケース2）
ಣℂ೨ ಣℂᓟ
図– 4 水際線の処理（ケース3）
小水深以下の場合，その要素は陸域とする．また，以下 に示す処理を行うことにより水域の要素については水深 を微小水深以上にする．この処理を施すことにより水際 線の移動を滑らかに実現する．
要素内の水深の分布は次の3ケースに大別される．
• ケース1
すべての節点の水深が微小水深より大きい．
• ケース2
要素の平均水深が微小水深より小さい．(図-3) ここで，h0は微小水深，¯hは着目する要素の平均 水深，qiは着目する要素のi番目の節点における 処理前の流量，qˆiは着目する要素のi番目の節点 における処理後の流量である．
• ケース3
要素の平均水深が微小水深より大きい，かつ節点 の一つ以上の水深が微小水深より小さい(図-4) ここで，∆hは微小水深と微小水深以下の節点の 水深の差である．
次に，ケース別の処理法を示す．
• ケース1の場合
処理前の水深，流量をそのまま用いる．なお，こ の要素は水域要素とする．
• ケース2の場合
水深はすべての節点で要素の平均水深¯hとなり，
流量はすべて0となる．なお，この要素は陸域要 素とする．
• ケース3の場合
微小水深以下の節点 (i = 2) の水深を微小水深 h0にする．この時の増加量∆hを微小水深以上 の節点(i= 1,3)の水深から均等に差し引く．流 量は，微小水深以下の節点(i= 2)の流量を微小 水深以上の節(i= 1,3)点の流量に均等に割り振 る．なお，この要素は水域要素とする．
上記の処理をRunge-Kutta法の各ステップ後に行う．
3. 数値解析例
本手法の計算精度と安定性を検討するために，段波 問題，ダムブレイク問題を取り上げ，厳密解，CG法 による計算結果との比較を行う．なお，CG法の空間 方向の離散化にはSUPG法に基づく安定化有限要素法
2)，時間方向の離散化にはDG法と同様に3次精度陽的 Runge-Kuuta法を用い，要素としてはDG法，CG法 ともに三角形一次要素を用いた．
3.1 段波問題
解析モデルを図-5に示す．水路奥行き幅は1.0[m]と した．x，y方向分割幅はともに0.1[m]，微小時間増分 量は0.001[s]とした．図-6は計算に用いたメッシュの 拡大図である．
解析結果として，図-7，図-8はそれぞれ1秒後の水 面形状，流速分布についてDG法のSlope Limiter処理 の有無による違いを示したものである．図より，Slope Limiter処理を行うことにより，オーバーシュート，アン ダーシュートを抑えていることがわかる．また，図-9， 図-10は同じく1 秒後の水面形状，流速分布について Slope Limiter処理を行ったDG法とSUPG法に基づ く安定化有限要素法を適用したCG法による解析結果の 比較を示したものである．図より，DG法による結果は CG法に比べて減衰がなく，厳密解とよい一致を示して いることがわかる．
3.2 ダムブレイク問題
解析モデルを図-11に示す．水路奥行き幅は1.0[m]
と し た ．x 方 向 分 割 幅 は 0.05[m]，y 方 向 分 割 幅 は 0.1[m]，微 小 時 間 増 分 量 は 1.0×10−4[s]，微 小 水 深 は1.0×10−5[m] と し た ．な お ，CG 法 の 微 小 水 深 は 1.0×10−2[m]とした．
解析結果として，図-12，図-13はそれぞれ1秒後の 水面形状，流速分布についてDG法とSUPG法に基づ く安定化有限要素法を適用したCG法による解析結果 の比較を示したものである．図より，DG法による結果 はCG法よりも微小水深をより小さく取れるため，厳密 解とよい一致を示していることがわかる．
1.0[m]
5.0[m]
10.0[m]
0.4[m]
図– 5 解析モデルおよび初期条件
図– 6 メッシュ拡大図
図– 7 水面形状
4. おわりに
本研究では，浅水長波流れの高精度な数値解析手法の 構築を目的とし，浅水長波方程式の離散化手法にDG法 を適用し，その有効性を検討した．数値解析例として，
段波問題とダムブレイク問題を取り上げ，以下の結論を 得た．
• Slope Limiter処理を行うことにより，不連続面 でのオーバーシュート，アンダーシュートを抑え られることが確認された．
• DG法は，従来のCG法に比べて減衰が少なく，厳 密解とよい一致を示していることが確認された．
図– 8 流速分布
図– 9 水面形状
図– 10 流速分布
• 移動境界を有する問題においてDG法は，厳密解 とよい一致を示していることが確認された．
今後の課題として，より複雑な問題への適用などが挙 げられる．
参考文献
1) D.Schwanenberg, M.Harms : Discontinuous Galerkin finite-element method for transcritical two-dimensional shallow water flows,Journal of Hydraulic Engineering., Vol.130, No.5, pp.412-421, 2004.
2) S.Takase, K.Kashiyama, S.Tanaka and T.E.Tezduyar： Space-time SUPG formulations of the shallow water equations, International Journal for Numerical Meth-
[m]
=O?
5.0[m]
図– 11 解析モデルおよび初期条件
図– 12 水面形状
図– 13 流速分布
ods in Fluids, Vol.64, pp1379-1394, 2010.
3) G.Kesserwani, R.Ghostine, J.Vazquez, A.Ghenaim, R.Mose：Riemann solvers with Runge-Kutta discontin- uous Galerkin schemes for the 1D shallow water equa- tions,Journal of Hydraulic Engineering, Vol.134, No.2, pp.243-255, 2008.
4) L.J.Durlofsky, B.Engquist and S.Osher：Traiangle based adaptive stencils for the solution of hyperbolic conservation lows, Journal of Computational Physics, 98, pp.64-73, 1992.
5) S.Bunya, E.J.Kubatko, J.J.Westerink and C.Dawson： A wetting and drying treatmentfor the Runge-Kutta discontinuous Galerkin solution to the shallow water equations, Computer Methods in Applied Mechanics and Engineering, Vol.198, pp1548-1562, 2009.
