数理解析研究所講究録 1441 巻 2005 年 159-164

159

特性曲線法に基づく有限要素スキームと差分スキーム
九州大学・大学院数理学研究院
野津 裕史 (Hirofumi Notsu) 2
田端 正久 (Masahisa Tabata) 1
Faculty of Mathematics, Kyushu University

1

はじめに時刻を , 空間に関する微分を で表し,は物質微分項と呼ばれる. 流れ問題では,$\nabla$$t$$u$を与えられた流速ベクトルとする.

$\phi$に,$( \frac{\partial}{\partial t}+u\cdot\nabla)\phi$

密度, 運動量, エネルギーが入る.

この項で記述される移流効果が拡散効果より支配的になれば, 有限要素法におけるガレルキン近似や差分法における中心差分近似から導かれるスキームは不安定になる. その処方とし特性曲線に基づく近似が開発されてきた. 本稿では, 特性曲線に基づく近似スキームの性質について考察する.
特性曲線に基づく有限要素法は, 時間 1 次精度スキームは良く知られていた [1]. 時間2 次精度スキームを物質微分項に対して実現することは難しくない. 常微分方程式の数値解法に現れる 2 次ルンゲ・クッタ法を使えば実現できる. 物質微分項以外の項は, クラン
ク. ニコルソン型の近似を用いれば各項時間 2 次精度の近似ができる. しかし、 このようにして得られたスキームは全体として, 時間 2 次精度になっていない [2]. スキーム全体として 2 次精度を維持するには, どのようにすればよいかを, 有限要素法 差分法の場合について考察し, それらのスキームの安定性と収束性を示す.に依存して決まる正定数であること
が $A,$ $B,$
本稿を通して, 記号 $c(A, B, \cdots)$ を,て、 風上近似,$c$$\cdots$を示すために用いる.

2

特性曲線に基づく近似$u=u(x, t)$ を既知の流速とする. $X=X(t)$ を時刻$t$における粒子の位置とする.

$X$は常微分方程式系$\frac{dX}{dt}(t)=u(X(t), t)$を満たしている.

このとき, 物質微分項は $X$ を使って$( \frac{\partial\phi}{\partial t}+u\cdot\nabla\phi)(X(t), t)=\frac{d}{dt}\phi(X(t), t)$

lE-mail : tabata@math.kyushu-u.ac.jp$2\mathrm{E}$(1)-mail : notsu@math.kyushu-u. ac.jp(2)

$1\mathrm{B}0$と書ける.

$\triangle t$

を時間刻みとして, (2) を$t$に関して差分近似した(3)$\frac{\phi(X(t),t)-\phi(X(t-\triangle t),t-\triangle t)}{\triangle t}$は (2) 式左辺の物質微分項の近似になる.

これが特性曲線に基づく近似の基本的な考え方である. (3) の $X(t-\triangle t)$ を$X(t-\triangle t)\simeq X(t)-u(X(t), t)\triangle t$としたものは, 常微分方程式 (1) を後退オイラー近似したものであり, 時間刻み 1 次精度の特性曲線に基づくスキームが得られる.
常微分方程式 (1 戸こ 2 次ルンゲ・クッタ法$X(t- \triangle t)\simeq X(t)-u(X(t)-u(X(t), t)\frac{\triangle t}{2},$ $t- \frac{\triangle t}{2})\triangle t$を適用して (3) を用いると, 物質微分項は時間刻み 2 次精度で近似できる. 他の項に関しに関して 2 次精度近似が実現できるては, クランク. ニコルソン型の近似を用いると$\triangle t$

それだけでは全体として時間刻み 2 次精度近似スキームは得られない [2]. 第 4 節で真に 2 次精度近似を得るために必要な付加項が示される.
が,3

移流拡散方程式$\Omega$

を $\mathrm{R}^{d}(d=2,3)$ の有界領域,$\Sigma_{T}=\Gamma \mathrm{x}(0, T)$と置く.

$\Gamma$をその周とする. $T(>0)$ を時刻とする, $Q_{T}=\Omega \mathrm{x}(0, T)$ ,を未知関数とする移流拡散問題$\phi:Q_{T}arrow \mathrm{R}$,$(x, t)\in Q_{T}$(4a),$(x, t)\in\Sigma_{T}$(4b)$\frac{\partial\phi}{\partial t}+u\cdot\nabla\phi-\nu\triangle\phi=f$$\phi=0$$\phi=\phi^{0}$を考える.
$\phi^{0}$

:

ここに,

$\Omega\prec \mathrm{R}$$u:Q\tau\prec \mathrm{R}^{d}$$x\in\Omega,$は与えられた流速,は初期値である. 流速$\nabla\cdot u=0$,$u$$\iota/$

(4c)

$t=0$

は拡散係数,$f$:$Q_{T}arrow \mathrm{R}$は外力,lま$((x, t)\in Q_{T})$

,

$u=0$

$((x, t)\in\Sigma_{T})$を満たしているとする. この仮定は本質的な仮定ではなく, そうでない場合にも適切な変更 [2] の下で本稿の結果は成立する. 境界条件に関しても同様である.

4

時間刻み 2 次精度近似
時刻での流速を,$n\triangle t$$u^{n}(x)\equiv u(x, n\triangle t)$で表す. 常微分方程式系 (1) の後退オイラー

近似と 2 次ルンゲ・クッタ近似$X_{1}^{n}(x)\equiv$を定義する,$x-u^{n}(x)\triangle t$,$X_{2}^{n}(x)\equiv$$x-u^{n-\frac{1}{2}}(x-u^{n}(x) \frac{\triangle t}{2})\triangle t$

161
4.1

有限要素法での近似

領域 を三角形 (四面体) 分割して得られる近似領域を, 斉次境界条件 (4b) に対応する有限要素空間を限要素空間を, を次式で求める.
$n=1,$
朔$\Omega_{h},$$\Omega$$X_{h}$$\in V_{h},$$\cdots,$有を近似するとする. 有限要素解$H^{1}(\Omega_{h})$$V_{h}$$P_{1}$$N_{T}$$( \frac{\phi_{h}^{n}-\phi_{h}^{n-1}\mathrm{o}X_{2}^{n}}{\triangle t},$$\psi_{h})+\frac{t/}{2}(\nabla\phi_{h}^{n}+\nabla\phi_{h}^{n-1}\mathrm{o}X_{1}^{n}, \nabla\psi_{h})+\frac{\nu\triangle}{2}(J^{n}\nabla\phi_{h}^{n-1}\mathrm{o}X_{1}^{n}, \nabla\psi_{h})$,$= \frac{1}{2}(f^{n}+f^{n-1}\mathrm{o}X_{1}^{n}, \psi_{h})$(5a)$\psi_{h}\in V_{h}$(5b)$\phi_{h}^{0}=\Pi_{h}\phi^{0}$

ここに,$(\cdot, \cdot)$

補間作用素,$0$での内積,はヤコビ行列は関数の合成を示している. (5a) をは$L^{2}(\Omega_{h})$$J$$(J_{ij}^{n}=\partial u_{i}^{n}/\partial xj),$$\mathrm{n}_{h}$:$C(\overline{\Omega})arrow X_{h}$は$\langle A_{h}^{n-1/2}\phi_{h}, \psi_{h}\rangle=\langle \mathcal{F}_{h}^{n-1/2}, \psi_{h}\rangle$と書くことにする.

4.2

差分法での近似を 2 次元長方形領域とする. 3 次元直方体領域のときにも以下の議論は自然に拡張ででの格子点関数とする. 差分解 ,を時刻きる. 空間格子間隔を $h(>0)$ とする.
$\Omega$$n\triangle t$$\phi_{h}^{n}$$n=1,$ $\cdots,$$N_{T}$$\phi_{h}^{n}$, を次式で求める.

$\frac{\phi_{h}^{n}-(\Pi_{h1}\phi_{h}^{n-1})\mathrm{o}X_{2}^{n}}{\triangle t}-\frac{\nu}{2}(\triangle_{h}\phi_{h}^{n}+\triangle_{h}\phi_{h}^{n-1})\sim-\frac{\iota/\triangle t}{2}(u_{1,2}^{n}+u_{2_{1}1}^{n})\nabla_{h,12}\phi_{h}^{n-1}$$= \frac{1}{2}$(は格子点集合,$f^{n-1}\mathrm{o}X_{1}^{n}$)

(6b)
(6c)

翫はその境界となる格子点集合である.

らなる離散ラプラス作用素, 入 は$h$$\triangle_{h}\phi_{h}=\sim(1+u_{1,1}^{n}\triangle t)\nabla_{h,1}$((垣$h \{\frac{1}{2,1},0$)$\nabla_{h,1}\phi_{h})\mathrm{o}X_{1}^{n})$$+(1+u_{2,2}^{n}\triangle t)\nabla_{h,2}((\Pi_{h1}^{\langle 0,\frac{1}{2})}\nabla_{h,2}\phi_{h})\mathrm{o}X_{1}^{n})$で. される変形ラプラス作$\text{定}\ovalbox{\tt\small REJECT}$$\mathrm{f}\mathrm{f}\mathrm{l}*_{\backslash }\not\equiv$(6a)$(x\in\Omega_{h})$$(x\in\overline{\Omega}_{h})$$\phi_{h}^{0}=\phi^{0}$$\Omega_{h}$十$(x\in\Gamma_{h})$$\phi_{h}^{n}=0$

ここに,$f^{n}$である.

$\Pi_{h}(\frac{1}{2,1},0)$は$\Omega_{h}^{(\frac{1}{2},0)}=\{((i+\frac{1}{2}h_{7}jh)\in\Omega;\mathrm{i},j\in \mathrm{Z}\}$$\triangle_{h}$は 5 点差分か$1\mathrm{B}2$

は $((\mathrm{i}\pm 1)h, (j\pm 1)h)$ と $((\mathrm{i}\pm 1)h, (j\mp 1)h)$の値を使う双一次補間作用素である.
は
方向 (中心) 差分作用素の 4 点を使うを近似する差分作用素である.
である. (6a) 式を$\nabla_{h,12}$$\nabla_{h,i}$$\partial^{\prime 2}/\partial x_{1}.\partial x_{2}$$x_{i}$$A_{h}^{n-1/2}\phi_{h}=\mathcal{F}_{h}^{n-1/2}$と書くことにする.

注意 1は$\nabla_{h,1}\phi_{h}$$\Omega_{h}^{(\frac{1}{2},0\rangle}$

で$\text{定}\ovalbox{\tt\small REJECT}$の点にならないので, 双一次される.

$\ovalbox{\tt\small REJECT}\ovalbox{\tt\small REJECT}_{\mathrm{B}}7$

作用$\Omega_{h}^{(\frac{1}{2},0)}$$\not\equiv*\backslash$の点を$X_{1}^{n}$で移動した点は, 一般に\Pi h(2-11 切が必要となる. 双一次$\ovalbox{\tt\small REJECT} 7\mathrm{a}5$

作用素$\Omega_{h}^{(\frac{1}{2},0)}$$\Pi_{h1}^{(0_{\mathrm{t}}\frac{1}{2})}$も同様の理由で必要となる.

5

近似精度

有限要素スキーム (5a), 差分スキーム (6a) の第 3 項は, スキーム全体が精度になるために, 必要な項である.
$A^{n-1/2}\phi$とお$\langle$.

十$u\cdot\nabla\phi-U\triangle\phi)^{n-1/\mathrm{z}}\circ Y_{1}^{h}’$,$O(\triangle t^{2})$の近似$\mathcal{F}^{n-1/2}\equiv f^{n-1/2}\circ Y_{1}^{n}$

ここに,$Y_{1}^{n}(x) \equiv x-\frac{1}{2}u^{n}(x)\triangle t$である.

補題 1(有限要素法, Lemmas $3,4_{7}[2]$ ) 滑らかな関数$|\langle$($A^{n-1/2}$

一護$n-\mathit{1}/\mathit{2}$$h$)$\phi,$$\phi,$$u,$ $f$と$\psi_{h}\in V_{h}$$\psi_{h}\rangle|\leq c(u, \phi)\triangle t^{2}.||\psi_{h}||_{I^{2}(\Omega_{h})}$」$|\langle \mathcal{F}^{n-1/2}-\mathcal{F}_{h}^{n-1/2},$ $\psi_{h\rangle}|\leq c(f)\triangle t^{2}||\psi_{h}||_{L^{2}(\Omega_{h})}$

に対して

(7a)
(7b)が成立する.
補題 2(差分法) 滑らかな関数$\phi,$$u,$ $f$

に対して$\Omega_{h}$で$|(A^{n-^{\tau}}[perp]/2-A_{h}^{n-1/2})\phi|\leq c(u, \phi)(\triangle t^{2}+h)$$|\mathcal{F}^{n-1/2}-F_{h}^{-1/2}|\leq c(f)\triangle t^{2}$(8a)(8b)が成立する.

6

安定性と収束性$\phi$を (4) の解する.

$\phi_{h}$を (5) の有限要素解 あるいは, (6) の差分解とする. 次の結果が成立

163
定理 1(有限要素法, Theorem 2, [2]) スキーム (5) は無条件安定であり,$||\phi-\phi_{h}||_{l^{\infty}(L^{2})}$が成立する.

,$\sqrt{l/}|\phi-\phi_{h}|_{l^{2}(H^{1})}’\leq c(u, \phi^{0}, f, \phi)(\triangle t^{2}+h)$

ここに,$||\psi_{h}||_{l^{\infty}(L^{2})}\equiv \mathrm{n}1\mathrm{a}\mathrm{x}\{||\psi_{h}^{n}||_{L^{2}(\Omega_{h});}n=0, \cdots, N_{T}\}$$| \phi_{h}|_{l^{2}(H^{1})}’\equiv\{\triangle t\sum_{n=1}^{N_{T}}||\frac{\nabla\phi_{h}^{n}+\nabla\phi_{h}^{n-1}\circ X_{1^{\hslash}}’}{2}||_{L^{2}(\Omega_{h})}^{2}\}^{1/2}$である.
定理 2(差分法) スキーム (6) は無条件安定であり,$||\phi-\phi_{h}||_{f(L^{2})}\infty$,$\sqrt{\nu}|\phi-\phi_{h}|_{\acute{\ell}^{2}(H^{1})}\leq c(u, \phi^{0}, f, \phi)(\triangle t^{2}+h)$が成立する. ここに,$|| \psi_{h}||_{\ell^{\infty}(L^{2})}\equiv\max\{||\psi_{h}^{a}r||_{\Omega_{h};}n=0, \cdots, N_{T}\}$,I$\psi_{h}||_{\Omega_{h}}\equiv$$| \phi_{h}|_{\acute{l}^{2}(H^{1})}\equiv.\{\triangle t\sum_{\prime n=1}^{N_{T}}||\frac{1}{2}\tilde{\nabla}_{h}\phi_{h}^{n}||^{2}L^{2}(\Omega_{h}^{\langle 1/2,0\rangle})\cross L^{2}(\Omega_{h}^{(0.1/2)}.)\}^{1/2}$,$\tilde{\nabla}_{h}\phi_{h}^{n}\equiv(\nabla_{h,1}\phi_{h}^{n}+(\Pi_{h1}^{(1/2,0)}\nabla_{h,1}\phi_{h}^{n-1})\circ X_{\mathrm{J}}^{n}$\nabla h,2\phi nh+(垣 (c1’$1/2$)$\nabla_{h,2}\phi_{h}n-1)\circ X_{1}^{\prime n})$である.
定理 2 は, 補題 2 の結果を使って, [2] の方針に従って離散 Gronwa 垣の不等式を用いて証明される.

注意 2 有限要素スキーム (5a) では合成関数を積分するために数値積分を用いることが多い. その際 精度の悪い数値積分法を用いると, 搬入する誤差の影響で計算が不安定になることがある [5] ので注意が必要である. 時間 2 次精度スキームは時間 1 次精度スキー
ムより, 搬入誤差に関して強靭であることが調べられている [3]. 差分スキーム (6a) では数値積分は必要ない.になる [2].
要素 $(k>1)$ に拡張でき, そのとき右辺は,注意 3 定理 1 は要素のときは, 集中質量近似を用いると差分的なスキームを作ることができる.
$P_{k}$$\ovalbox{\tt\small REJECT}$$c(\triangle t^{2}+h^{k})$謝辞
この研究は, 科学研究費 (S), 課題番号 16104001 により日本学術振興会からの支援と九州大学 21 世紀 COE プログラム 「機能数理学の構築と展開」により文部科学省からの支援とを受けて遂行された.

164
参考文献
1. O. Pironneau, Finite Element Methods for Fluids, John Wiley&Sons,

1989.

2. H. Rui and M. Tabata, A second order characteristic finite element
convection-diffusion problems, Numer. Math., 92(2002):161-177.

scheme forscheme of
3. M. Tabata and S. Fujima, Robustness of a characteristic finite element
second order ni time increment, to appear in Proceedings of the Third International
Conference on Computational Fluid Mechanics.

methods for convection-dominated dif4. J. Douglas Jr. and T. F.
fusion problems based on combining the method of characteristic with finite elementof finite difference procedures, SIAM J. Numer. Anal., 19(1982):871-885.
$\mathrm{R}\mathrm{u}\mathrm{s}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{l}:\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}$5. M. Tabata, Discrepancy between Theory and Real Computation on the Stability ofSome Finite Element Schemes, Kyushu Univ. Preprint Series, MHF2005-2, 2005.