高次要素を用いる圧力安定化有限要素法の解析と高レイノルズ数流れ問題への応用 (数値解析学の最前線 : 理論・方法・応用)

全文

(1)36 高次要素を用いる圧力安定化有限要素法の解析と高レイノルズ数流れ問題への応用. Analysis of a pressure‐stabilized finite element method with higher‐order elements and application to high‐Reynolds‐number flow problems 内海晋弥. 早稲田大学基幹理工学部,[email protected] Shinya Uchiumi. School of Fundamental Science and Engineering. Waseda University. 1. はじめに本稿では,Navier‐Stokes 問題と Oseen のための高次要素を用いる圧力安定化. 局所線形化流. 速Lagrange‐Galerkin スキームを導入する.Oseen 問題に対するスキームでは,粘性係数依存性に注目した理論的な誤差評価を示す.加えて,ある Navier‐Stokes 流れの数値計算結果を紹介し, 導入したスキームの現実的な問題における有効性を示す. 非圧縮粘性流体の運動を記述する Navier Stokes 方程式と,その移流項を線形化した Oseen 方程. 式に対する数値計算スキームを考える.特に粘性係数が小さい場合 , すなわち,レイノルズ数が高い場合を考察する.このような状況で安定な計算を行うためには,物質微分項の近似方法の選択が重. 要な論点の一つであることが知られている.特性曲線の方法と有限要素法を結合させたLagrange‐. Galerkin 法 (LG 法) はそのような状況に対処する有効な手法のつである [11, 12, 13, 14]. 解くべき連立一次方程式に現れる係数行列は対称であり,ゆえに,共役勾配 (CG) 法や最小残差 (MINRES) 法[1] といった,効率が良い線形ソルバーを用いることができる.LG 法はスキームの -. 定式化と実装との間に乖離があることが問題であったが,最近,局所線形化流速を用いて,数値積分を使うことなく厳密に計算できかつ収束性を証明できるスキームが作成され,この問題が解. 消された [1314]. これらの流体問題の有限要素解析において,流速と圧力を近似する空間の組の選択も重要な論点である.この考察は,より簡単な定常 Stokes 問題においても必要である.純粋な Galerkin スキー. ムでは,その組は inf‐sup 条件を満たすことが要求され,Taylor‐Hood 要素とも呼ばれる P_{k}/P_{k-1} 要素 (k\geq 2) がそれを満たす代表的なものの一つである [2]. ここで, P_{k} は三角形 k 次要素を表す. 一方で,流速と圧力に同じ次数の有限要素空間を用いる P_{k}/P_{k} 要素 (k\geq 1) を使う際には適切な安定化項が必要となる.定常 Stokes 問題を考えるとき, P_{1}/P_{1} 要素に対しては Brezzi‐Pitkäranta. [3] によって安定化項が導入され,Burman [4] によって高次要素へ拡張された.非定常 Stokes 問題に対しては,Burman‐Fernández [5] で [4] と同種の安定化項を用いた対称なスキームが解析された.一方で,運動方程式の残差に基づくタイプの安定化項も研究されてきた [9_{\dot{\ovalbox{\t \smal REJECT} 10] . しかし,スキームの実装や解析が前者のタイプの安定化項よりも,特に非定常問題に対して,複雑である. Navier‐Stokes 問題や Oseen 問題における移流項の制御の問題とは別に,粘性係数依存性はよ. り単純な定常 Stokes 問題にも現れる.その依存性の改善方法の一つは grad‐div 安定化項の付加. である.非定常 Oseen 問題に対して de Frutos ら [7] により,時間方向に差分近似 , 空間方向に Galerk甫近似を用い,grad‐div 項を付加したスキームが考察され,その粘性係数依存性に対する効果が解析された.また,非圧縮性を厳密に満たす Scott‐Vogelius 要素が Navier‐Stokes 問題に.

(2) 37 対して用いられ,有効性が報告されている [6]. 本稿では,Oseen 問題に対する P 扉 P_{k} 圧力安定化 Lagrange‐Galerkin スキームを導入し, Navier‐Stokes 問題と合わせて数値結果を示す.ここでは Burman [4] の安定化項を加えたスキームを考察する.このスキームは非定常 Stokes 問題に対して考察されていた [5] のスキームと同種のものであるが,彼らのスキームとは安定化項のパラメータの取り方が異なる.また,本スキー. ムは Lagrange Galerkin 法の利点である対称性が引き継がれている.さらに [14] と同じく局所線形化流速が用いられており,スキームに忠実な実装が可能である.第2節では,Oseen 問題のた. めのスキームに対して,[15] に基づき,粘性係数依存性に注目した誤差評価を示す.第3節では, 円板回りの Navier‐Stokes 流れの数値計算結果を報告する.そこでは, P 雇 P_{k-1} 要素を用いるスキームと比較して, P 副琉圧力安定化法が良好な結果を示すことが観察できる.. 2. Oseen 問題のための Lagrange‐Galerkin スキームとその誤差評価 (u, p) : \Omega\cross(0, T)arrow \mathbb{R}^{d}\cross \mathbb{R} を未知関数とする Oseen 問題 :. \frac{\partial u}{\partial t}+(w\cdot\nabla)u-\nu\triangle u+\nabla p=f, (x, t)\in\Omega\cross(0, T) \nabla\cdot u=0, (x, t)\in\Omega\cross(0, T). ,. (Os). ,. u=0, (x, t)\in\partial\Omega\cross(0, T). ,. u(\cdot, 0)=u^{0}, x\in\Omega \Omega\subset \mathbb{R}^{d}, d=2,3 は多角形または多面体領域, T>0_{\dot{r}}0<\nu\leq 1 はそれぞれ時. を考える.ここに,. 刻,粘性係数を表す定数, \Omega. w,. f :. \Omega\cros (0, T)ar ow \mathbb{R}^{d_{\ovalbox{\t \small REJECT}} .u^{0}. : \Omegaarrow \mathbb{R}^{d} は与えられた関数である. \partial\Omega は. の境界を表す.. Oseen 問題に対する数値計算法を考える.まず,特性曲線に沿った物質微分の離散化を導入する.. w. を滑らかとする.特性曲線 X(t;x, s) は常微分方程式系. \frac{dX}{dt}(t;x, s)=w(X(t;x, s), t) , t<s, X(s;x, s)=x の解として定義される.これを用いると,物質微分項. ( \frac{\partial}{\partial t}+(w\cdot\nabla) u. を. ( \frac{\partial u}{\partial t}+(w\cdot\nabla)u)(X(t), t)=\frac{d}{dt}u(X(t), t) と書ける.. \triangle t>0. を時間刻みとする.. t^{n}\equiv n\triangle t_{:}.u^{n}(x)\equiv u(x, n\triangle t) とし, f^{n} なども同様に定め. る.流速場 w^{\star_{1} : \Omegaarrow \mathbb{R}^{d} に対して写像 X_{1}(w^{*}) を. X_{1}(w^{*})(x)\equiv x-w^{*}(x)\triangle t で定める.. X_{1}(w(x, t)) は X(t-\triangle t;x, t) の Euler 近似である.このとき. \frac{\partial u^{n} {\partial t}+(w^{n}\cdot\nabla)u^{n}=\frac{u^{n}-u^{n.-1} \circ X_{1}(w^{n-1})}{\triangle t}+0(\triangle t) が成り立つ.ここでする.. \circ. は関数の合成を表す.また,以下では N_{T}\equiv\lfloor T/\triangle t\rfloor を総時間ステップ数と.

(3) 38 次に,有限要素空間を導入する.. \{T_{h}\}_{h\downarrow 0} を \overline{\Omega} の正則な三角形 (四面体) 分割列とし,各要素. K\in T_{h} は少なくとも1つ境界上に無い頂点を持つとする.. w_{h}^{(m)}\subset H^{1}(\Omega). に対して,. 作用素とする.. h. を要素の最大直径とする.整数 m\geq 1. を五上の P_{m} 有限要素空間とし,. k\geq 2 に対して,. \Pi_{h}^{(m)}. を P_{m} 有限要素空間への補間. V_{h}\cross\overline{Q}_{h}:=( W_{h}^{(k)})^{d}\cap H_{0}^{1}(\Omega)^{d}) \cross(W_{h}^{(k-1)}\cap L_{0}^{2}(\Omega)) を流速,圧力に対応する P_{k}/P_{k-1} 有限要素空間とする.これは (一般化された) Taylor‐Hood 要. 素とも呼ばれ,離散化された Stokes 問題の適切性を保証する inf‐sup 条件を満たす [2]. 双一次形式. a, b. を. a(u_{\grave{\ovalbox{\tt\small REJECT}}}v):=\nu(\nabla u, \nabla v) , b(v, q) : =-(\nabla\cdot v, q) で定める.ここで. (\cdot, ). は. L^{2}(\Omega), L^{2}(\Omega)^{d}. または. L^{2}(\Omega)^{d\cross d}. の内積である.. スキーム OsTH. u_{h}^{0}\in 砺を u^{0} の近似とする.次を満たす. \{(u_{h}^{n}, p_{h}^{n})\}_{n=1}^{N_{T} \subset V_{h}\cross\overline{Q}_{h}. を求めよ.. ( \frac{u_{h}^{n}-u_{h}^{n-1}\circ X_{1}(H_{h}^{({\imath})}w^{n.-1})}{\triangle t}, v_{h})+a(u_{h}^{n}, v_{h})+b(v_{h}, p_{h}^{n})=(f^{n}, v_{h}) , \foral v_{h}\in V_{h},. b(u_{h}^{n}, q_{h})=0, \forall q_{h}\in\overline{Q}_{h}.. P_{k}/P た有限要素空間と安定化項を導入する.. k\geq 1 に対して,. V_{h}\cross Q_{h}:=((W_{h}^{(k)})^{d}\cap H_{0}^{1}(\Omega)^{d})\cross(W_{h} ^{(k)}\cap L_{0}^{2}(\Omega)) を P_{k}/P_{k} 有限要素空間とする.この要素は inf‐sup 条件を満たさないので圧力安定化項を必要とする.圧力安定化項 \mathcal{C}_{h} を. \mathcal{C}_{h}(p, q):=\sum_{K\in T_{h} h_{K}^{2k}\sum_{|\alpha|=k}(D^{\alpha} p, D^{\alpha}q)_{K} K の直径であり, (\cdot, )_{K} は L^{2}(K) の内積である.この項は Bunnan [4] によって導入されている. P_{1}/P_{1} 要素に対する Brezzi‐Pitkäranta [3] の安定化項の高次要素へ. で定める.ここに, h_{K} は要素の拡張である.. スキーム OsPstab.. u_{h}^{0}\in 琉を. u^{0} の近似とする.次を満たす. \{(u_{h}^{n}, p_{h}^{n})\}_{n=1}^{N_{T} \subset V_{h}\cross Q_{h}. ( \frac{u_{h}^{n}-u_{h}^{7?.-1}\circ X_{1}(H_{h}^{(\`{I})}w^{71.-1})}{\triangle t}, v_{h})+a(u_{h}^{n}, v_{h})+b(v_{hP_{h}^{n})}=(f^{n}, v_{h}). ,. を求めよ.. \forall v_{h} \in Vń,. b(u_{h}^{n}, q_{h})-\delta_{0}\mathcal{C}_{h}(p_{h}^{n}, q_{h})=0, \forall q_{h}\in Q_{h}.. ここで,. \delta_{0}>0 は安定化パラメータである.. 注意1.. 1. スキームから生じる連立一次方程式の係数行列は対称である.. 2. 非定常 Stokes 問題に対して Burman‐Ferllández [5] では [4] と同種の安定化項を用いた対称なスキームが解析されていた.彼らは安定化項のパラメータ \delta_{0} を. たが,その選択では以下で示す. 3.. k=1. \nu. \nu. に依存して取ってい. に依存しない誤差評価を得ることは困難と思われる.. のとき,スキーム OsPstab はNotsu Tabata [11] により作成と解析が行われている.. 下では,粘性係数に注目した新しい誤差評価を示す..

(4) 39 4. 写像 X_{1}(\cdot) では,元の流速場 w^{n-1} の代わりに局所線形化流速場 H_{h}^{(1)}w^{n-1} が使われている.これの導入により (u_{h}^{n-1_{\circ}}X_{1}(\Pi_{h}^{(1)}w^{n-1}), v_{h}) は厳密に積分することができる [13_{:-}14].. 元の流速場を使った時,厳密な積分値を求めることは困難である.. 本スキームの粘性係数依存性に注目した誤差評価を述べる.. 定理.. u_{h}. をスキーム OsPstab の解とし,(Os) の解 (u, p) は十分滑らかとする.流速場. 滑らかで \partial\Omega で 0 とする.時間刻み \triangle t は 0<\triangle t\leq\triangle t_{0_{\dot{\ovalbox{\t \small REJECT}} }\triangle t_ {0}. 満たすとする.ここで. \alpha. u_{h}^{0}. c. が存在して. \Vert u-u_{h}\Vert_{\ell^{\infty}(L^{2})}, \sqrt{\nu}\Vert ▽ (u-u_{h})\Vert_{\ell^{2}(L^{2})}\leq c(\triangle t+h^{k}+h^{2}). \psi=\{\psi^{n}\}_{n=0}^{N_{T}. が成立する.ここに,. を. は u^{0} の P_{k} 有. h, \triangle t に依存しない正定数. \nu,. は十分. :=(4\alpha_{*}|W|_{C([0,T];W^{1,\infty}(\Omega)^{d})})^{-1}. 。はメッシュの最小角のみに依存する正定数である.. 限要素空間における Lagrange 補間とする.このとき,. w. (^{*}). に対して,. \Vert\psi\Vert_{\ell^{\infty}(L^{2})}\equiv\max\{\Vert\psi^{n}\Vert_{L^{2} (\Omega)};n=0, N_{T}\},. \Vert\psi\Vert_{\el^{2}(L^{2}) \equiv(\trianglet\sum_{n=1}^{N_{T} \Vert\psi^ {n}\Vert_{L^{2}(\Omega)}^{2})^{1/2} である.. 注意2.. 1. 定数. 2. スキーム. c. は厳密解. OsT^{-}H. u, p. に依存する.. に対しても (^{*} ) と同じ収束次数の評価を得ているが,定数 c が. \nu. に依存する.. 3. (^{*}) 右辺の h^{2} の項は,局所線形化流速と元の流速との差の評価から現れる.. 定理の証明は [15] を参照されたい.. 3. 数値結果. 1. 本節では Navier‐Stokes 問題の数値例を示す.Navier‐Stokes 問題は (Os) において流速場. u. それぞれ,スキーム OsTH とスキーム OsPstab において,. w^{n-1} を. u_{h}^{n-1}. を未知. で置き換えることによ. り得られる. スキーム NSTH. 次を満たす. \{(u_{h}^{n}, P_{h}^{\mathfrak{n} )\}_{\gamma 1=1}^{N_{T} \subset V_{h} \cros \overline{Q}_{h}. を求めよ.. ( \frac{u_{h}^{n}-u_{h}^{\tau\iota-1}\circ X_{1}(\Pi_{h}^{(1)}u_{h}^{n-1})} {\triangle t}, v_{h})+a(u_{h}^{n}, v_{h})+b(v_{h},p_{h}^{n})=(f^{n}, v_{h}), \foral v_{h}\in V_{h},. b(u_{h}^{n}, q_{h})=0, \forall q_{h}\in Q_{h}.. スキーム NSPstab. 次を満たす. \{ (u_{h}^{n}, p_{h}^{n})\}_{n=1}^{N_{T} \subset V_{h}\cross Q_{h}. ( \frac{u_{h}^{n}-u_{h}^{n-1}\circ X_{1}(\Pi_{h}^{(1)}u_{h}^{n-1})}{\triangle t}, v_{h})+a(. を求めよ.. 嶋 v_{h})+b(v_{hP_{h}^{n})}=(f^{n}, v_{h}), \forall v_{h}\in V_{h},. b(u_{h;}^{n}q_{h})-\delta_{0}\mathcal{C}_{h}(P_{h}^{n}, q_{h})=0, \forall q_{h} \in Q_{h}.. ここで,. w. に置き換えることにより得られる.対応するスキーム NSTH とスキーム NSPstab は,. \delta_{0}>0 は安定化パラメータである..

(5) 40. 1-\frac{3^{Q} {4}0. u=(^{1}. (-5,-2) -0 \langle 10,-2). 図1: 領域と境界条件.. 図2: メツシュ 1. これまでところ,スキーム NSPstab に対する粘性係数依存性を排する理論的な誤差評価は得ら. れていない.[15] では,創生解から生成される Navier‐Stokes 問題に対してこれらのスキームの誤差が数値的に比較され,スキーム NSPstab の誤差がより小さいことが報告されている.特に粘性係数が小さい時にこの差は顕著である.以下では,円板回りにおける Navier‐Stokes 流れの数. 値実験を行い,両者のスキームの精度を比較する. 領域 \Omega と Dirichlet 境界条件を図1のように与える.Navier‐Stokes 問題において. \nu=10^{-3_{:}}. f=0 とする.. 2つのスキームにおいて,時間刻みを. \triangle t=0.01. とする.FreeFem. ++. [8] を使って2つのメッ. シュを作成した.円板周辺のみを細かく分割したメッシュ 1と,全体的に細かく分割したメッシュ. 2を用いる (図 2_{\ovalbox{\t smalREJC.T} 表1). 両スキームにおいてん =2 とする.すなわち,スキーム NSTH では P_{2}/P_{1} 要素を,スキーム NSPstab では P_{2}/P_{2} 要素を用いる.スキーム NSPstab では安定化パラメータを \delta_{0}=10^{-3} とした.スキームの初期値 u_{h}^{0} は,図1 と同じ境界条件を課した Stokes 問題の解 (図3) とした.なお,この Stokes 問題では粘性係数を \nu=1 と設定している. 図4は t^{n}=20 における鰐の流線と, u_{h1}^{n} (1, ) のプロファイルを示している.メッシュ 1で P_{2}/P_{2} 要素を用いて得られた解は,メッシュ 2で P_{2}/P_{1} 要素を用いて得られた解とほぼ一致している.一方,メッシュ 1で P_{2}/P_{1} 要素を用いて得られた解は,これらとは異なる値をとっている. 全体的に細かいメッシュ 2によって得られた解はより厳密解に近いと考えられるので, P_{2}/P_{2} 要 \cdot. 素を用いて得られた解がより精度が良いと考えられる.. 4. おわりに第2節では,非定常 Oseen 問題に対して P_{h}./P_{k} 要素を用い,安定化項を加えるスキームを導入. し,その粘性係数依存性に注目した誤差評価を述べた.Navier‐Stokes 問題に対するスキームでは理論的な結果を得られていないが,第3節で見た円板回りの流れの数値結果では,. P_{2}/P_{2} 要素を.

(6) 41 41. メッュ1;2の \overline{T}^{\backsla h\ovalbox{\t smal REJ CT}- タ \underline{表1:_{\backslash }} メッシュ 1 メッシュ 2 \ovalbox{\t smalREJ CT}. \sqrt[ovalbx{\t smalREJCT}]{\ovalbx{\t smalREJCT}. \ovalbox{\t\smal REJ CT}. 総要素数. 総自由度 (P2/P1) 総自由度 (P2/P2). 14_{\vee}.188 65.046. 120_{-}.738. 86.568. 160_{\dot{\ovalbox{\t \smal REJECT}} 824. 26_{\vee}.564. x_{1} u_{1\iota 1}(1,\cdot). 図3:. u_{h}^{0}. の流線 (上,左下) と. u_{h1}^{0}(1, \cdot). のプロファイル (右下)..

(7) 42 -\cdot-\cdot-\cdot-P2lP1 , meshl -----P2/P2 , meshl. P2/P1 , mesh2. a^{\sim}. S. -1. 0. 1. 2. -1. 0. 1. 2. x_{1} uh1(1,\cdot) 図圭 u_{h}^{n}(t^{n}=20) の流線 (_{-}\vdash_{-} , 左下) と. 用い安定化項を加えるスキームが,. u_{h\`{I} ^{n}(1, \cdot). のプロファイル (右下).. P_{2}/P_{1} 要素を用いるスキームよりも良好な結果が得られるこ. とが観察された.. 参考文献 [1] M. Benzi, G.H. Golub_{i} and J. Liesen. Numerical solution of saddle point problems. Acta Numerica.. 14:1-137_{:} 2005.. [2] D. Boffi, F. Brezzi. and M. Fortin.. Mixed Finite Element Methods and Applications.. Springer. Berıin Heidelberg. 2013.. [3] F. Brezzi and J. Pitkäranta. On the stabilization of finite clement approximations of the Stokes equations. In W. Hackbusch. editor. Efficient Solutions of Elliptic Systems, pages 11‐19. Vieweg: 1984.. [4] E. Burman. Pressure projection stabilizations for Galerkin approximations of Stokes’ and Darcy’s problem. Numerical Methods for Partial Differential Equations. 24(1):127-143_{\ovalbox{\t \small REJECT}}. 2008.. [5] E. Burman and M. A. Fernández. Galerkin finite element methods with symmetric pressure stabilization for the transient Stokes equations: Stability and convergence analysis. SIAM Journal on Numerical Analysis, 47(1):409-439 , 2009.. [6] M. A. Case, V. J. Ervin. A. Linke, and L. G. Rebholz. A connection between Scott‐Vogelius and grad‐div stabilized Taylor Hood FE approximations of the Navier‐Stokes equations. SIAM Journal on Numerical Analysis, 49(4):1461-1481_{\mathfrak{i}} 2011..

(8) 43 [7] J. de Frutos. B. García‐Archilla, V. John, and J. Novo. Grad‐div stabilization for the evolu‐ tionary Oseen problem with inf‐sup stable finite elements. Journal of Scientific Computing,. 66(3):991-1024 , 2016. [8] F. Hecht. New development in 4):251-265_{0}. 2012.. FreeFem++ .. Journal of Numerical Mathematics, 20(3 ‐. [9] T.J.R. Hughes, L.P. Franca. and M. Balestra. A new finite element formulation for compu‐ tational fluid dynamics: V. circumventing the Babuška‐Brezzi condition: a stable Petrov‐ Galerkin formulation of the Stokes problem accommodating equal‐order interpolations. Computer Methods in Applied Mechanics and Engineering, 59(1):85-99_{\dot{0}} 1986.. [10] V. John and J. Novo. Analysis of the pressure stabilized Petrov‐Galerkin method for the evolutionary Stokes equations avoiding time step restrictions. SIA Analysis. i53(2):1005-1031_{i} 2015.. M. Journal on Numerical. [11] H. Notsu and M. Tabata. Error estimates of a pressure‐stabilized characteristics finite element scheme for the Oseen equations. Journal of Scientific Computing.. 65(3):940-955_{:}. 2015.. [12] H. Notsu and M. Tabata. Error estimates of a stabilized Lagrange Galerkin scheme for the Navier‐Stokes equations.. 50(2):361-380_{\backslash }.. ESAIM: Mathematical Modelling and Numerical Analysis,. 2016.. [13] M. Tabata and S. Uchiumi. A genuinely stable Lagrange‐Galerkin scheme for convection‐ diffusion problems. Japan Journal of Industrial and Applied Mathematics. 33(1):121-143, 2016.. [14] M. Tabata and S. Uchiumi. An exactly computable Lagrange Galerkin scheme for the Navier Stokes equations and its error estimates. Mathematics of Computation.. 87:39-67_{\dot{}}. 2018.. [15] S. Uchiumi. A viscosity‐independent error estimate of a pressure‐stabilized Lagrange‐ Galerkin scheme for the Oseen problem. arXiv:1712.04150 [math.NA]..

(9)