高次補間に基づく
DG有限要素法による浅水長波流れ解析手法の構築
Development of a numerical method for shallow water flow based on higher order discontinuous Galerkin finite element method15N3100020D 花澤 広貴 Hiroki HANAZAWA
Key W ords:DG-F EM, Shallow water f low, Runge Kutta method, Slope limiter
1.
はじめに
近年, 地震に伴う津波, 台風や豪雨に伴う高潮, 洪水 などの水害が深刻化してきている. これらの被害を予測 する手段として, 数値シミュレーションが広く用いられ ており, その精度の向上が求められている. 津波, 高潮, 洪水などの現象を表す数値モデルは,双曲型偏微分方程 式に分類される. 双曲型問題は不連続的な解を有しやす く,これを精度よく表現するためには比較的細かなメッ シュが必要とされ, 計算負荷がかかる. 双曲型方程式 に対する有効な解析手法の一つとして, Discontinuous Galerkin有限要素法1)(以下DG法)が注目されている. DG法は要素間の不連続性を許容する手法であり,保存 則を満足するとともに, 高次要素の適用が容易という特 徴を有する. そのため, 比較的粗いメッシュを用いても 不連続な現象をより精度よく表現することが可能な手法 であると考えられる.
そこで本研究では, 高精度な浅水長波流れ解析手法の 構築を目的とし, 浅水長波方程式に対して高次補間(2 次)に基づくDG法の適用を行った. 本手法の妥当性と 有効性について検討するため, 数値解析例として段波問 題と孤立波伝播問題を取り上げた.
2.
数値解析手法
(1) 支配方程式支配方程式には,以下に示す保存型2次元浅水長波方 程式を用いる.
∂U
∂t +∂F(U)
∂x +∂G(U)
∂y =S(U) (1) ここで, U, F(U), S(U)はそれぞれ未知変数, 流束関 数, ソース項に関するベクトルであり, 以下のように表 される.
U=
" H uH vH
#
,F(U) =
uH
u2H+12gH2 uvH
(2)
G(U) = 2 4
vH uvH v2H+12gH2
3 5,S(U) =
" 0
−gH(S0x+Sf x)
−gH(S0y+Sf y)
# (3)
S0x= ∂z
∂x, S0y= ∂z
∂y (4)
Sf x=gn2u|u|
H13 , Sf y= gn2v|v|
H13 (5) ここで, H は全水深, u,vはそれぞれx,y方向の断面平 均流速, gは重力加速度である. また,S0x, S0yはそれぞ れx,y方向の底面勾配,Sf x,Sf yはそれぞれx,y方向の 底面せん断応力,nはマニングの粗度係数である. (2) DG法による空間方向の離散化
支配方程式に対して, 空間方向の離散化としてDG法 を適用する. 重み付き残差法を適用し,第二項,第三項に 対して部分積分を行い, Gaussの発散定理を用いると弱 形式は以下のように表される.
Z
Ω
W∂U
∂t dΩ− Z
Ω
³ F∂W
∂x +G∂W
∂y
´ dΩ +
Z
Γ
W(Fnx+Gny)dΓ = Z
Ω
WSdΩ (6) ここで,Wは要素毎に定義される不連続な重み関数,nx, ny は要素境界における外向き単位法線ベクトルの成分 である. 各要素内の補間は, 図-1に示す, 直交性を有す
るDubinerの基底関数2)と自由度の線形結合で表され
る. 式(6)において, 各項の積分には数値積分法の一つ であるLegendre-Gauss法を適用する. 図-2,図-3に示 すような一般化座標系への変換を行い, 要素内部及び要 素境界における積分点と重みを用いて積分計算を行う.
図– 1 階層型直交基底関数2)
2016年度 中央大学理工学部都市環境学科修士論文発表会要旨(2017年2月)
図– 2 1次要素の積分点(一般化座標)
図– 3 2次要素の積分点(一般化座標)
(3) 数値フラックス
DG法では要素境界で物理フラックスを複数有する ため, 数値フラックスを用いて近似を行う. 式(6)左辺 第三項のフラックスと外向き単位法線ベクトルの内積 を数値フラックスに置き換え, 本論文では以下に示す Local-Lax Friedrich flux3) を用いる.
Eˆ =1 2
hEn(U+) +En(U−)i
−amax
2 (U+−U−)(7) ここで,
En =Fnx+Gny (8) であり, amaxは以下のように表される.
amax=maxh
|λ1(U+)|,|λ2(U+)|,|λ3(U+)|,
|λ1(U−)|,|λ2(U−)|,|λ3(U−)|i (9)
λ1=unx+vny−p
gH (10)
λ2=unx+vny (11)
λ3=unx+vny+p
gH (12)
図– 4 要素境界値の定義
また, U+,U−は要素境界における値であり,図-4に示 すように, 要素境界であらかじめ設定された法線ベクト ルの方向の値をU+ とし, もう一方の値をU− と定義 する.
(4) 時間方向の離散化
以上を整理すると,式(6)は以下のような有限要素方 程式となる.
M ˙U=AfF(ΦnU) +AgG(ΦnU)
+CS(ΦnU)−B ˆE (13) ここで,M,Af,Ag,C,Bはそれぞれ質量行列,フラッ クスF, Gに関する行列, ソース項に関する行列, 要素 境界に関する行列, また,Φnは各積分点における基底関 数の行列である. 式(13)を全要素で重ね合わせる際, 要 素境界では数値フラックスを積分点で評価する. このと き,任意境界辺の単位法線ベクトルが着目した要素に対 して内向きの場合は, 積分点の順序及び法線ベクトルの 向きを正しく考慮するためにB˜m=−BmIaとする. こ こで,mはローカルな辺番号である. また, Iaは辺上積 分点数分の正方行列であり, 単位行列を反転させたもの である. 以上より, 式(13)は以下のような時間に関する 常微分方程式として表される.
dU
dt =L(U) (14)
式(14)に対して,時間方向の離散化として式(15)〜(19) に示す陽的Runge-Kutta法(2/3公式)を適用する. 1 次要素を用いる場合は2段階2次, 2次要素を用いる場 合は3段階3次と,空間精度に対して次数を1つ上げて 離散化を行う.
2段階2次Runge-Kutta法 Un+
2
3 =Un+2
3∆tU˙n (15) Un+1=Un+∆t
4 ( ˙Un+ 3 ˙¯Un+23) (16)
2016年度 中央大学理工学部都市環境学科修士論文発表会要旨(2017年2月)
3段階3次Runge-Kutta法 U¯n+12 =Un+∆t
2
U˙n (17)
U¯n+1=Un+ ∆t(−U˙n+ 2 ˙¯Un+12) (18)
Un+1=Un+∆t
6 ( ˙Un+ 4 ˙¯Un+12 + ˙¯Un+1) (19) (5) slope limiter処理手法
段波などの不連続面を有する問題を解いた際に発生す るオーバーシュートやアンダーシュートを抑制するため に, 本研究ではslope limiter処理手法5) 6)を適用する. この処理をRunge-Kutta法の各段階で行っていく.
3.
数値解析例
(1) 段波問題
本手法の妥当性を確認するために, 段波問題を取り 上げる. 図-5 に示すような初期水位を与え, 境界条 件として壁面でslip条件を与える. メッシュ分割は,
∆x = 0.25m, ∆y = 2.5m とし, 微小時間増分量は 0.001secとした.
図-6, 図-7に1秒後の水面形状及び流速分布を示す. これらより, DG法の1次要素及び2次要素の結果は厳 密解と良い一致を示しており, 本手法の妥当性が確認さ れた. また,図-8,図-9より, slope limiter処理を施すこ とによって, 不連続面でのオーバーシュート, アンダー シュートを抑制できていることが確認できる.
図– 5 解析モデル(段波問題)
図– 6 1秒後の水面形状(limiter処理なし)
図– 7 1秒後の流速分布(limiter処理なし)
図– 8 1秒後の水面形状(limiter処理あり)
図– 9 1秒後の流速分布(limiter処理あり)
(2) 孤立波伝播問題
高次要素を用いることの有効性を確認するために, 図-10に示すStreetの水理実験モデル7)を用いた孤立 波の伝播問題を取り上げる. 式(20)に示す孤立波及び 式(21), (22)に示す流量を初期条件として与え, ∆x=
∆y = 0.1mとし, 2次要素の要素数は1次要素の半分 としている. 壁面はslip条件とし, 微小時間増分量は 0.001secとした. なお,本例題ではlimiter処理は行って いない.
ζ=η Ã
sech r3H
4h3(x−ct)
!2
(20)
qx=c(h+ζ) (21) c=p
gh (22)
2016年度 中央大学理工学部都市環境学科修士論文発表会要旨(2017年2月)
図– 10 解析モデル(孤立波伝播問題)
図– 11 1次要素20分割, 2次要素10分割時の空間波形
図– 12 1次要素30分割, 2次要素15分割時の空間波形
図– 13 1次要素60分割, 2次要素30分割時の空間波形
また,領域分割数を20, 30, 60の3パターンで変更し1 次要素と2次要素の空間精度の比較を行った. なお, そ れぞれのパターンで2次要素の要素数は1次要素の要素
数の半分としている. 参照解としてはSUPG法に基づ く安定有限要素法8)を用いた.
図-11〜図-13に領域分割幅を変更した際の空間波形 を示す. これらより, DG法の結果は参照解と概ね良い 一致を示していることがわかる. また, 2次要素を用い た結果は1次要素に比べて少ない要素数を用いた場合で も, 1次要素と概ね同等の精度が出ていることが確認で きる.
4.
おわりに
本研究では，高精度な浅水長波流れ解析手法の構築を 目的として，浅水長波方程式に対して高次補間(2次)に 基づくDG有限要素法の適用を行った．本手法の妥当 性と有効性を確認するため, 段波問題と孤立波伝播問題 を取り上げ以下の結論を得た．
• 段波問題において, DG法の1次要素及び2次要 素の結果は厳密解と良い一致を示し,その妥当性 を確認した.
• 孤立波伝播問題において, 2 次要素を用いた結 果は, 1次要素に比べて少ない要素数を用いた場 合でも1次要素と同等の精度が出ることを確認 した.
今後の課題として, 遡上問題を取り扱うための移動境 界条件処理手法の導入及び実地形問題への適用が挙げら れる.
参考文献
1) B.Cockburn, G.E.Karniadakis, and C.W. Shu : The Development of Discontinuous Galerkin Methods, T heory Comput. Appl., pp.3-50, 1999.
2) 牧野優作,田中聖三,桜庭雅明,樫山和男: Discontinuous Galerkin有限要素法の浅水長波流れ解析への適用,第61 回理論応用力学講演会, OS13-01, 2012.
3) M.Dubiner : Spectral Method on Triangles and Other Domains,J. Sci. Comput., 6, pp.345-390, 1991.
4) B.Cockburn, C.W. Shu : The Runge-Kutta discontin- uous Galerkin method for conservation laws V: multi- dimensional systems,J. Comput. P hys.141, pp.199- 224, 1998.
5) L.J.Durlofsky, B.Engquist and Osher : Triangle based adaptive stencils for the solution of hyperbolic conser- vation lows,J. Compt. P hys., 98, pp.64-73, 1992.
6) L.Krividonova, J.Xin, J.-F.Remacle, N.Chevaugeon and J.E.Flacherty : Shock detection and limiting with discontinuous Galerkin methods for hyperbolic con- servation laws,Appl. N umer. M ath., Vol.48, Issues.3- 4, pp.323-338, 2004.
7) R.L.Street, S.J.Burgers and P.W.Whitford : The Be- havior of Solitary Waves on a Stepped Slope, Dept.
Civ. Eng. Stanf ord U niv. T ech. Rep., No.93, 1968.
8) 利根川大輔,樫山和男: 安定化有限要素法による非線形分 散波理論に基づいた津波遡上解析手法の構築研究,応用力 学論文集, Vol.12, pp.127-134, 2009.
2016年度 中央大学理工学部都市環境学科修士論文発表会要旨(2017年2月)
