数理解析研究所講究録
1129 巻 2000 年 9-22

9A domain-decomposition finite-element scheme for flow problems
Choice of elements-

茨城大学理学部1

藤間昌– (Shoichi Fujima)

はじめに
近年、並列計算機はいろいろな分野で計算手法を変革しつつある。 メモリ分散型並列計算機では、 プロセッサ間の情報の共有はメッセージパッシングと呼ばれる通信手続きによる。その実行はプロセッサ内の計算に比べて低速である。 したがって、 メッセージパッシングの少ない計算手法を開発する必要がある。領域

分割法はこのようなメモリ分散型並列計算機で偏微分方程式を解く有力な手法と考えられる。
領域分割法の–種として、部分領域間の弱い連続性のために Lagrange 乗数を導入する手法がある。 こ

.
.

れは、

部分領域の計算の独立性が大きいこと共有面のない部分領域間では局所的な情報の共有の必要がないことの特長があり、並列計算化の手段として有望である。
Navier-Stokes 方程式に対して、 この型の領域分割法がいくつか提案されている [1, 2, 7, 8, 10, 17]。ここでは 2 次元非圧縮 Navier-Stokes 方程式の領域分割型解法 $[7, 8]$ を発展させて、 Lagrange 乗数の基底関数の選択の拡大や、流速・圧力の基底関数の異なる選択と合わせて考えることにより、 より自由度の高い領域分割への適用や、 2 流体問題への適用を検討する。以下では第 2 節で領域分割型解法を総括し、第 3 節で

3 角形キャヴィテイ内流れの計算例を示す。第 4 節で自由度の高い領域分割への適用を、第 5 節で 2 流体問
題への適用を検討する。

2

Navier-Stokes 方程式の領域分割有限要素解法

2.1

Navier-Stokes 方程式

2 次元の有界領域$\Omega$

内で,非圧縮 Navier-Stokes 方程式,$\frac{\partial u}{\partial t}+u\cdot \mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}u+\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}P$$\mathrm{d}\mathrm{i}\mathrm{v}u$,$=$$\frac{1}{Re}\triangle u+f$

$=$

$0$

(1)

,

(2)

: $\Omega\cross(0, T)arrow R^{2}$ と圧力 $p:\Omega\cross(0,T)arrow R$ は未知関数, $f$ は外力, $Re$ は数である。
Reynolds
境界の 2 つの部分,とでそれぞれ, 流速境界条件 u $=u_{D}$ と応力境界条件\mbox{\boldmath $\sigma$}(u, )を課す。ここに, $\sigma_{ij}(u_{)}p)=-p\delta ij+(2/Re)Dij(u)(i,j=1,2)$ は応力テンソル,を考える。 ここに, 流速$u$$\Gamma_{D}(\neq\emptyset)$$\Gamma_{N}$$p$$\cdot n=\tau N$$D_{\dot{\mathrm{a}}j}(u)=(\partial u_{i}/\partial_{X_{j}}+\partial uj/\partial_{X_{i})}/2$$(i,j=1,2)$ は変形速度テンソル,は単位外向き法線ベクトル,とは境界上の流速と表面力で既知のデータである。 初期条件は時刻 $t=0$ で u=uI。とする。 なお以降のスキームの説明では簡単のため、$n$$u_{D}=0$とする。

$u_{D}$$\tau_{N}$

10

2.2

標準的な有限要素解法

221

混合型有限要素近似

空間の離散化には混合型有限要素近似 [9] を用いる。 の 3 角形要素分割を行い, 要素の最大直径を内に求める。各 3 角形を 4 等分とする。圧力は要素分割に基づく(1 次) 要素による有限次元関数空間要素による有限次元関数空間殊内に求めるした細かい要素分割を考え, 流速はその要素分割に基づ(節点数が P2 要素と等しい Pl 要素なので、 P2 同相 Pl 要素と呼ばれる。 また最大直径 $h/2$ の要素分割に基づく Pl 要素なので、 Pl(h/2) 要素とも呼ばれる。 本論文では後者の表記法を用いる。 すなわち、 流速/$h$$\Omega$$Q_{h}$$\mathrm{P}1$$\text{く}\mathrm{P}1$

である) これらの関数空間に属する関数を下付添え字 で示す。 この関数空間の組み合わせば, 後で定義する双 1 次形式 が下限上限条件を満たす [3]。このため, 定常
Stokes 方程式に関して, 有限要素スキームの–意可解性が保証され, また, 解が $u\in(H^{2}(\Omega))^{2},$ $p\in H^{1}(\Omega)$ならば, 誤差評価,$||u-u_{h}||_{V}+||P-p_{h}||Q\leq Ch(||u||_{(}H^{2}(\Omega))2+||p||_{H^{1}(}\Omega))$(3)
圧力の有限要素の組合せは、 Pl$h$$(\mathrm{h}/2)/\mathrm{P}\mathrm{l}(\mathrm{h})$。

$b$が成り立つ。 ここに$c$

は $h$ に依存しない定数である。

$N$を$\phi_{i,u}i=1,$$V_{h}$$\ldots,$の基底関数,

離散化 [11] を用いる。 すなわち,$\psi_{i},i=1,$$\ldots,$$N_{p}$

を$Q_{h}$の基底関数とする。 時間方向には準陽的なを時間刻みとして, 上付添え字$\Delta t$

$n$

は時刻

$n\triangle t$

におけるデータを表わすものとし, 式 (1) を陽的に離散化し, 式 (2) を陰的に扱う。

以上により次の有限要素スキームを得る。
$(u_{h},p_{h}n+1n)\in V_{h}\mathrm{x}Q_{h}$,,$( \frac{u_{h}^{n+1}-u_{h}n}{\Delta t},\phi_{i})h^{+(u_{h}})a_{1}^{h},uh’\emptyset i)+\Omega_{0}(u_{h}^{n},\phi_{i}+nnb(\phi i,p_{h}^{n}) = \langle\hat{f},\phi_{i})$$b(u_{h}^{n+1},\psi_{i})$

$0$

$=$

,

$1\leq i\leq N_{u}$

(4)

,

(5)

$1\leq i\leq N_{p}$

ここに,

で,

$(, )_{h}$,$(u,v)$$=$$\int_{\Omega}u\cdot vdx$$a_{1}(w, u,v)$$=$$\int_{\Omega}(w\cdot \mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}u)\cdot vd_{X}$$a_{0}(u,v)$$=$$\frac{2}{Re}\int_{\Omega}D(u)$$b(v_{\rangle}q.)$$=$$- \int_{\Omega}q\mathrm{d}\mathrm{i}\mathrm{v}vdx$$\langle\hat{f},v\rangle$$=$$\int_{\Omega}f\cdot vdx+\int_{\Gamma_{N}}\tau N.vd_{S}$は $(, )$ の質量集中化近似を,$a_{1}^{h}$は$a_{1}$,: $D(v)dx$ ,,,の上流下流点選択型有限要素近似 [18] を示す。 なお、$(u_{h}^{0}, \phi i)=(u_{Ic,\phi)}i$,$u_{h}^{0}$は、$1\leq i\leq N_{u}$を解いて求める。
各時刻ステップで解くべき連立 1 次方程式は,となる。 ここに,トル,$F^{n}$

は$U^{n}$$\langle\hat{f}\rangle\phi_{i}\rangle$

は$\overline{M}\frac{U^{n+1}-U^{n}}{\Delta t}+A_{1}(U^{n})U^{n}+A_{0}U^{n}+B^{T}P^{n}$$=$$F^{n}$$BU^{n+1}$

$=$

$0$

$u_{h}^{n}$の節点自由度からなるからなる$N_{u^{-}}$$N_{u^{-}}$

ベクトル,$P^{n}$

は$p_{h}^{n}$ベクトルで,$[\overline{M}]_{\iota j}$$=$$(\phi_{\mathrm{j}},\phi_{i})_{h}$,$i,j=1,$$N_{u}$$\ldots,$

,

,

(6)

,

の節点自由度からなる(7)$N_{p^{-}}$

ベク

11

$[A_{1}(U^{n})]_{ij}$$=$$a_{1(u_{h},\emptyset\phi_{i}}^{hn}j,)$$[A_{0}]_{ij}$$=$$a_{0}(\emptyset j’\phi_{i})$,$i,j=1,$ $\ldots,N_{u}$ ,$i,j=1,$ $\ldots,N_{u}$ ,,はそれぞれ, 集中化質量行列, 移流行列, 粘性行列であり,$[B]_{ij}=b(\phi_{j}, \psi_{i})$,$i=1,$$\ldots,$$N_{p}$,$j=1,$ $\ldots,N_{u}$ ,は発散行列である。
整合離散化圧力 Poisson 方程式

222

式 (6) $,(7)$ を解くには, 以下の整合離散化圧力 Poisson 方程式による算法 [11] が計算量を減らすために有効である。 中間変数$\tilde{U}$を導入し, (6) を 2 つに分ける :$\overline{M}\frac{\tilde{U}-U^{n}}{\triangle t}+A_{1}(U^{n})U^{n}+A_{0}U^{n}$$\overline{M}\frac{U^{n+1}-\tilde{U}}{\triangle t}+B^{\tau_{P}n}$$B\overline{M}\cdot(9)$$=$$F^{n}$

$=$

$0$

,

(8)

.

(9)

(7) を代入すると,

に,

(10)

$(B \overline{M}B^{T})P1n=\frac{1}{\Delta t}B\tilde{U}$が得られ,$U^{n}arrow\tilde{U}arrow P(8\rangle(10)n(9arrow U^{n+1})$の手順で 1 つの時間ステップの計算を進められる。 最終的に解くべき連立 1 次方程式のサイズはり, これは$N_{u}$

領域分割型解法

2.3.1

領域分割
$\Omega$

であに比べて小さい。

2.3

領域$N_{p}$

を重なりのない $K$ 個の部分領域,$\overline{\Omega}=\overline{\Omega_{1}}\cup\cdots\cup\overline{\Omega_{K}}$,$\Omega_{k}\cap\Omega_{l}=\emptyset(1\leq k<l\leq K)$に分割する。 部分領域間の界面を$\overline{\mathrm{r}_{m}}=\overline{\Omega_{\kappa(m)^{\cap\overline{\Omega_{\kappa}}}(}-}+m)’(1\leq m\leq M)$とする。 ここで,$\kappa_{-}(m)$と $\kappa+(m)$ は$\Gamma_{m}$の両側の部分領域の番号である $(\kappa_{-}(m)<\kappa_{+}(m))$ 。各部分領域の 3 角形要素分割を行う。 この際, 以下を仮定する。
仮定 1 界面は直線分である。

仮定 2 界面上で両側の部分領域の要素分割の節点位置は–致する。

2.3.2

Lagrange 乗数法

ごとに, 流速の Pl(h/2) 有限要素空間, 圧力の Pl(h) 有限要素空間
部分領域
全領域 に対する有限要素空間はそれらの積空間により作る :$\Omega_{k}$$V_{k,h}$$\Omega$ズ$V_{h}^{*}$$=$$\prod_{k=1}V_{k,h}$$Q_{h}^{\star}$$=$,$k \tau\prod_{-,-}^{K}Q_{k},h$.

$Q_{k,h}$を構成する。

12

全領域\Omega での自然な数値解を得るために, 界面での部分領域間の流速の連続性を拘束条件とする Lagrange
上で定義される Lagrange 乗数の有限要素空間 $W_{m,h}$ には、 $P1(h)$ 要素を用いる。
乗数法を用いる。界面$\Gamma_{m}$Lagrange 乗数全体の有限要素空間はこれらの積空間により構成する:
レ$W_{h}^{\star}= \prod_{m=1}W_{m},h$スキーム (4),(5) に,.

Lagrange 乗数法を適用すると, 次のスキームを得る o$(u_{h}^{n+1n},p_{h}, \lambda_{h}^{n})\in V_{h}\star \mathrm{X}Q\star h\cross W_{h}^{\star}$$( \frac{u_{h}^{n+1}-u^{n}h}{\triangle t},$,$=$$v_{h})_{h^{+}}b^{\star}(vh_{J}p_{h}^{n})+j\star(vh, \lambda_{h}^{n})$$\langle\hat{f}^{\star},$$v_{h}\rangle-a1(h\star n)-a^{*}\mathrm{o}(u_{h},v_{h})uh’ u_{h}^{n},v_{h}n$,(11)$\forall v_{h}\in V_{h}^{\star}$$b^{\star}(u_{h}^{n+}, q_{h}1)$

$=$

$0$

,

$\forall q_{h}\in Q_{h}^{\star}$(12)$j^{\star}(u_{h^{+}}^{n},\mu_{h}1)$

$=$

$0$

,

$\forall\mu_{h}\in W_{h}^{\star}$

(13)

ここで,

$=$

$\sum_{k=1}^{R’}\int\Omega k(w_{k}\cdot \mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}u_{k})\cdot$$a_{0}^{\star}(u, v)$$=$$\frac{2}{Re}\sum_{k=1}^{K}\int_{\Omega}k)D(u_{k}:D(v_{k})d_{X}$$b^{\star}(v,q)$$=$$- \sum_{k=1}^{\mathrm{A}’}\int_{\Omega}kqk\mathrm{d}\mathrm{i}\mathrm{V}v_{k}dx$$j^{\star}(v,\mu.)$$=$$- \sum_{m=1}^{M}\int\Gamma_{m})(v_{\kappa}(m-v-(m\rangle)+\kappa\mu md\mathit{8}$

$v)$

$=$

$\sum_{k=1}^{K}(\int_{\Omega_{h}}f\cdot vkd_{X}+\int_{\partial\Omega_{k}\cap}\gamma_{N}\Gamma_{N}’\cdot v_{k}d_{S})$$\langle\hat{f}^{\star},$である。

$a_{1}^{h\star}$

,

,
,

,

は礎の上流下流点選択型有限要素近似である。

は表面力ベクトル$\lambda_{m,h}$vkdx,$a_{1}^{\star}(w,u,v)$$\sigma\cdot n_{\kappa}+(m)|\mathrm{r}_{m}$の近似に成り得る (実際に近似になるか否かは, 後で見るように,有限要素の組み合わせに依存する)。

$=$

解くべき連立 1 次方程式は次のようになる。

ここで,

行列,

2.3.3

$\overline{M}$$F^{n}$は$(, )_{h}^{\star}$に対応する集中化質量行列,$B$は$b^{\star}$に対応する発散行列,

(14)

$J$

は

$j^{\star}$に対応するジャンプは既知ベクトルである。

整合離散化圧力 Poisson 方程式の領域分割法主$U^{n+1}$を消去することにより, 整合離散化 Poisson 方程式の領域分割法スキーム,$=\backslash (\mathcal{F}\mathcal{G})$

(15)

13

が得られる。 ここで,
$A$

$=$

$J$

$=$

$\mathcal{M}$$\mathcal{F}$$\mathcal{G}$

$=$

$=$

$B\overline{M}B^{T}1$$\overline{JM}B^{T}1$

$=$

,

,

$J\overline{M}J^{T}1$$B\overline{M}F^{n}1$$\overline{JM}F^{n}1$

,

,

,

である。行列 $A$ は, K- ブロック対角行列であり, 各ブロックは部分領域の整合離散化圧力 Poisson 方程式に対応する。
更に$P^{n}$

を消去すると,(16)$(\mathcal{M}-JA^{-}1JT)\Lambda n=\mathcal{G}-JA-1\mathcal{F}$が得られる。 この式に共役勾配法を適用することにより, 並列計算に適した反復計算式が得られる。

24

.
.

すでに得られた結果$[7, 8]$.

.

の結果を以下にまとめる。

流速/ 圧力/Lagrange 乗数に対する $Pl(h/\mathit{2})/Pl(h)/P\mathit{1}(h)$ 要素の組合せでは、試験問題の計算で、数値解 (流速、圧力、 Lagrange 乗数= 界面上の表面力) の厳密解への収束が、数値的に確かめられた。

Lagrange 乗数に Pl(h) でな Pl
要素を用いる場合 : 流速/ 圧力の収束は数値的に観察されるが、 Lagrange 乗数は表面力に収束しない。
$\text{く}$$(\mathrm{h}/2)$

Lagrange 乗数に Pl(h) でなく改良 Pl$(\mathrm{h}/2)$

要素 (改良はモルタル法に基づき, 界面端で連続

$0$

次) を

用いる場合 : 数値解の厳密解への収束が、数値的に確かめられた。 ただし、 誤差の絶対値は Pl(h) 要素の方が小さい。

(上記の 3 事項は第 4 節の図 5,6 でも観察できる。)
並列計算機 (Intel Paragon$K^{0.7}$ノード) への実装により、プロセッサ数 $K$ のとき、計算速度は倍となり、 この手法が並列計算に有効であることが確かめられた。
$\mathrm{X}\mathrm{P}/\mathrm{S},$$56$

数値例: 3 角形キャヴィティ内流れの計算

3

では矩形領域を矩形の部分領域に分割する場合の数値例を示した。 ここでは有限要素法の特長\sim
領域形状の自由性\sim と領域分割有限要素法の適合性を示すため、矩形でない領域の問題の例として、 3 角形$[7, 8]$

キャヴィティ内流れの計算例を示す。
正方形キャヴィティ内の流れ問題は, 単純な領域形状にも関わらず Reynolds 数に応じて変化する流れ場を生じるため, 要素 (格子) 生成が容易なことも相まって, 計算スキームのよい試験問題となっている。
これに対し, 3 角形キャヴィティ内の流れ問題は, 直交格子生成が困難なためにあまり扱われなかった。 しかしこれは例えば流路の角ごとに現れる重要な流れ現象である。 この観点に着目し最近, 計算例が示されている $[12, 16]$

3.1

3 角形キャヴィティ内流れ問題

領域
する。

。

$\Omega$$\Omega$を, 底辺の 2 頂点が $B(0,0),$ $C(.1,0)$ で, 頂点 $A$ の角度がで Navier-Stokes 方程式をみたし, 境界$\mathrm{B}\mathrm{C}$$2\alpha$

上で, $u=(\mathrm{i}, 0)$ ,れ場を考える。 Reynolds 数は底辺長を代表長さとして定義する。

下向き 2 等辺 3 角形の領域とAB,AC 上で, $u=0$ , をみたす流
の,

14
32

数値結果

の領域形状で, $Re=2\sqrt{3},200^{\sqrt{3},0^{\sqrt{3}}}100$ の場合の計算を行った。 要素分割数は 3550, 領域
分割数は 7 である。流線を図 1 に示す。流れ場は $Re=2\sqrt{3}$ ではほぼ左右対称であったが, $Re=200\sqrt{3}$ で$2\alpha=\pi/3$

は特写が少し下流に移動し,$Re=1000\sqrt{3}$ では再び中心方向に移動する。 この挙動は [16] の結果とよく –致している。

図 3 に示す。 $Re=400$ の計算
次に深さが 2 の領域形状の場合を計算した。 流線と圧力等高線を図
では 6 次渦の存在まで捕らえられている。 また, スキームで圧力には部分領域間の連続性を課していないが, ほぼ連続な圧力場が観察できる。 図 4 は中心軸上の水平方向流速成分の分布を, 縦軸は頂角からの距は角付近
離により, 横軸は流速成分の絶対値により, 両対数軸上にプロットしたものである。
の Stokes 流れでは角度に応じて無限渦列が生じ, 渦の大きさ比, 強さ比の漸近値は角度により決まること$2_{\text{、}}$$\mathrm{M}\mathrm{o}\mathrm{f}\mathrm{f}\mathrm{a}\mathrm{t}\mathrm{t}[13]$を解析的に導いた。従って, 図 4 のプロットの包絡線は直線になるはずである。 およその傾向は観察されるが, 小さい渦になるに従って解像度は粗くなる。 このような問題で主渦を適切に捕らえながら, 無限岡列をいくつまで捕らえられるか ? は異なるスケールの現象を数値スキームにより同時に捕らえるという観点からの適切な試験問題であると思われる。
以上のように、領域分割有限要素法により、 3 角キャヴィティ内流れの現象が適切に捕られた。

4

界面に定数要素を用いる場合

ここまでの領域分割有限要素法では、領域分割、 有限要素分割に、 仮定 1, 仮定 2 を考えた。 Lagrange
乗数法を用いる領域分割法であるモルタル法 [4] の特長は、仮定 2 が不要なことである。 しかし、仮定 1 は必要である。仮定 1 の下では界面上の表面力は連続関数と考えられ、 Pl 要素などで近似することが可能になる。

この節では Lagrange 乗数に定数要素を用いることにより、仮定 1 を外すことを検討する。Lagrange 乗数$\mathrm{P}0$4.1($\mathrm{h}$または$\mathrm{h}/2$) 要素の収束性に関する数値的な考察を行う。

界面の形状と要素の選択

領域分割有限要素法では, 解析対象領域の領域分割と (3 角形) 要素分割を行う。そのどちらを先に行う
かにより 2 つに分けられる。

領域分割を先に行う場合領域の形状や現象の解の空間構造を考慮して, または他の目的で, 先に領域分割を行い,各部分領域の要素分割を後で行う。形状が比較的単純な部分領域ごとに要素分割することになり, 一般的には要素分割が容易になる。部分領域間界面では節点位置を整合させる必要がある。 ただし, 解析手法としてモルタル法を使用すれば, その不整合は許容される。
要素分割を先に行う場合 要素分割が先にあり, 領域分割法を適用するために要素をグルーピングするものである。 負荷分散の観点から各グループの要素数がなるべく等しく, かつ界面がなるべく小さい方が良い。 自動グルーピングの手法が研究されている ([15] など)。
両者にはそれぞれ得失があり, また前者で領域分割後に要素分割をしてから、後者に戻して再度領域分割をすることも考えられる。 したがって, 汎用性のある領域分割型解法には, 要素のグルーピングで得られる領域間界面の形状でも可能なことが望まれる。

4.2

精度の比較に関する数値実験

Lagrange 乗数の有限要素空間を次の 5 種類に変えて, 流速, 圧力, Lagrange 乗数の精度を比較する実験を行う。 比較は,

15$(\mathrm{c})Re=1000\sqrt{3}$

図 1: 正三角形内流れ

16

図 2: 深さ 2 の場合の流線と圧力 $(Re=400)$

図 3: 中心軸上の水平方向流速, Re=400(両対数軸)

17

..
.
.
.

Pl(h) 要素,$\mathrm{P}1(\mathrm{h}/2)$

要素,

改良 Pl(h/2) 要素 (改良はモルタル法に基づき, 界面端で連続

$0$

次),

PO(h) 要素,PO(h/2) 要素,

について行った。

試験問題は, 領域\Omega $=(0,1)\cross(0,1)$ , 厳密解が, $u=(x^{2}y+y, -x-x3sy2),$ $p=- \frac{1}{2}+x+y^{s}s$ , Reynolds
数 $Re=400$ で, これらから定まる境界条件と外力をデータとし, 初期条件 $u=0$ から非定常問題を解き,の等分割であり, 要素分割は得られる定常解の誤差を求めた。領域分割はに対して,$r\iota=6,12,24,48$ , の等 3 角形分割である。 なお, この実験では部分領域間界面は直線である。
相対誤差, 圧力の流速の相対誤差を図 4 に, Lagrange 乗数の相対誤差を図 5 に示す。
流速, 圧力の誤差ともに要素サイズに比例して誤差が/J| さくなっている。 Lagrange 乗数の誤差の比較ではPl(h/2) 要素では誤差が大きいままであるが, 他の 4 つの選択では誤差が収束した (この場合も Pl(h) が最$3\cross 2$$H^{1}(\Omega)$$\Omega$$n\cross n(\cross 2)$$L^{2_{-}}$$L^{2}(\Omega)$

良であった)。
以上のように予備的な実験ではあるが, Lagrange 乗数に PO(h) 要素または PO(h/2) 要素を用いる場合についても解の収束性が観察でき, 使用可能性が確認できた。

52 流体閣題への応用
最後に、 2 流体による移動境界問題への領域分割有限要素法の適用を考える。

界面捕捉法による有限要素スキームの 1 つとして、流速/ 圧力/ 擬密度関数に、3 角形 Pl(h/2)/P1(h)/P1(h/2)
要素を用いる方法が提案されている [14]。また界面捕捉法で表面張力の効果を取り入れる手法として, 表面張力をあたかも体積力のように扱う手法が提案されている [6]。
ここでは、 上記のようなスキームによる大規模計算を意識して、領域分割法による並列計算化を考える。 なお、精度の観点から Pl(h/2) 要素の代わりに$\mathrm{b}$要素も双 1 次形式
が下限上限条件をみたす組合せである。 領域分割に基づき各部分領域で流速/ 圧力/擬密度関数に関して$\mathrm{P}2(\mathrm{h})$要素を用いる。

$\mathrm{P}2(\mathrm{h})/\mathrm{P}\mathrm{l}(\mathrm{h})$の有限要素空間を構成する。全領域での整合をとるために、部分領域間の界面における流速と面密度関数の跳びが の制約に関する Lagrange 乗数法を用いる。 従って、界面上に流速の Lagrange 乗数、
擬密度関数の Lagrange 乗数の有限要素空間をとることになる。
以下では領域分割有限要素法の算法を示し、 P2 要素の未知変数の部分領域間の接続に用いる Lagrange$\mathrm{P}2/\mathrm{P}1/\mathrm{P}2$$0$

乗数の有限要素空間の選択について言及する。

5.1

2 流体の流れ問題と有限要素スキーム

混ざらない (immiscible)2 種類の流体 (流体 1, 流体 2 とよぶ) による流れ問題を考える。 流体により密度と粘性係数は異なる。 初期時刻において各流体が占める領域と初期流速が与えられたときに, 流体領域の変化と速度圧力分布を求める。
の$(0\leq t\leq T,\alpha=1,2)$ を流体\alpha が占めるを 2 種類の流体が占める有界領域 (固定),部分領域とする。各流体の密度を \rho \alpha 、粘性係数をとする。各部分領域の内部では, 非圧縮 Navier-Stoke8
方程式を考える。 2 種類の流体間の界面 $I(t)$ 上では運動学的条件から、流速の連続性と界面の移動に関する条件,$\Omega\subset R^{2}$$\Omega_{\alpha}(t)$$\Omega$$\mu_{\alpha}$$u^{(1)}$$u\cdot n$

$=$

$=$

$u^{(2)}$$v\cdot n$

(17)

(18)

18

$\text{図}4$: Error for the velocity(left) and the pressure(right)

図 5: Error for the Lagrange multiplier19を課す。 ここで、$u^{(\alpha)}$は各流体の流速,は界面の移動速度,

$v$

$n$

は界面上の単位法線ベクトルである。 また、

動力学的条件から,(19)$\tau^{(1)}-\tau=\sigma\kappa n(2)$を課す。 ここで,
面力,$p^{(\alpha)}$$\sigma$は表面張力係数,$\kappa$

は界面の曲率,は界面上で流体$\tau^{(\alpha)(\alpha)}=\tau(u^{\alpha\alpha},p)$$\alpha$にはたらく表は各流体の圧力である。

これらの方程式を有限要素近似する。 この際, (19) の右辺の表面張力項を体積力として Navier-Stokes
方程式の外力項に付け替えること [6] を行えば, 界面の条件 (19) は単純な表面力の釣り合い条件になる。密度と粘性係数の違いに注意しながら, 全領域 で Navier-Stokes 方程式を解けば良い。 そこで流速の 3 角$\Omega$

形 P2 要素による有限要素関数空間を式を満たす流速と圧力を求める。

$V_{h}$, 圧力の Pl 要素による空間を$Q_{h}$

として, 各時刻ステップで次の$(u_{h^{+1}}^{n},p_{h}^{n})\in V_{h}\mathrm{x}Q_{h}$$( \frac{u_{h}^{n+1}-u_{h}n}{\Delta t},v_{h})+(u_{h}^{n}\cdot \mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}uh’ vh)n-\frac{1}{\rho}(p_{h}^{n},\mathrm{d}\mathrm{i}\mathrm{v}v_{h})+2\nu(D(u_{h}^{n}) :D(v_{h}))$$(\mathrm{d}\mathrm{i}\mathrm{v}u_{h^{+}}^{n},q_{h})1$

ここで,

$f$

は重力,

$f_{ST}$$=$$=$$(f, v_{h})+(f_{ST}, v_{h})$

$0$

,

,

$v_{h}\in V_{h}$, (20),(21)$q_{h}\in Q_{h}$は体積力としてモデル化した表面張力である。

方, 界面の形状決定のためには界面捕捉法を用いる。すなわち、流体の種類を表す関数 $\phi=\phi$ ( ,t)(こ$x$

こでは擬密度関数と呼ぶ) :$\phi$
を導入し, 擬密度関数に関する移流方程式を, 3 角形 P2 要素による有限要素関数空間$( \frac{\phi_{h}^{n+1}-\phi^{n}h}{\triangle t},$$\xi_{h})+(u^{n}h^{+}1.

\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\phi_{h}^{n},\xi_{h})=0$,$X_{h}$

により,$\xi_{h}\in X_{h}$(22)を解いて求める。
表面張力項の算出で単位法線ベクトルと曲率は, 擬密度関数を用いて,

$n$

$=$

$\kappa$$=$のように求められる。 ここでは

(23)$\frac{\nabla\phi}{|\nabla\phi|}=(\frac{\phi_{x}}{\sqrt{\phi_{x}^{2}+\phi_{y}^{2}}},$ $\frac{\phi_{y}}{\sqrt{\phi_{x}^{2}+\phi_{y}2}})$7$\cdot n=\frac{\phi_{yy}^{22}\phi_{xx}-2\phi_{x}\emptyset\phi xy+\emptyset x\phi yy}{(\phi_{x}^{2}+\phi 2)^{3}y/2}$(24)に付した添字 $x,$ は偏微分を表す。
擬密度関数に 3 角形 P2 要素を用いることの利点は, 単位法線ベクトルと曲率が要素単位で求められ、$\phi$$y$

表面張力効果を有限要素に取り込むのに適していることである。

5.2

領域分割有限要素法スキーム

上記の 2 流体問題の有限要素スキームで並列計算を行うことを目的として領域分割有限要素法を適用する。 重なりのない領域分割下において、 自然な拡張として次のような算法が考えられよう。

521$\Omega_{k}$流速・圧力の解法
毎に, 流速の P2 有限要素空間

空間は積空間,$V_{h}^{\star}= \prod_{k}V_{k,h},$$V_{k,h}$, 圧力の Pl 有限要素空間$Q_{h}^{\star}= \prod_{k}Qk,h$$Q_{k,h}$を作り, 全領域に関する有限要素, による。界面での流速の連続性を拘束条件とする Lagrange

20上で定義する Lagrange 乗数の有限要素空間 $W_{m,h}$ を作り, 全界面に関する有限乗数法を用いる。界面
要素空間は積空間, $W_{h}^{\star}= \prod_{m}W_{m,h}$ , による。
$\Gamma_{m}$(20),(21) に対応して, 各時刻ステップには次の問題を解く。
$(u_{h}^{n+1},p_{h}^{n}, \lambda_{h}^{n})\in V_{hQ_{h}\mathrm{X}W_{h}}^{\star}\mathrm{X}\star\star$$=$$f(v_{h})$$b^{\star}(u_{h^{+}}^{n}, q_{h}1)$

$=$

$0$

$j^{\star}(u_{h^{+}}^{n},\mu_{h}1)$

$=$

$0$

$( \frac{u_{h}^{n+1}-u^{n}h}{\Delta t},$ $v_{h})_{h}+ \frac{1}{\rho}b^{\star}(v_{h,p_{h}^{n}})+j\star(vh, \lambda^{n}h)$

ここで,

$f$

6.2.2

擬密度関数の癬法$\Omega_{k}$,,,$v_{h}\in V_{h}^{\star}$$q_{h}\in Q_{h}^{\star}$

,

(25)

,

$\mu_{h}\in W_{h}^{\star}$

(26)

,

(27)

は既知項である。

毎に, P2 有限要素空間 $X_{k,h}$ を作り, 全領域に関しては積空間, $x_{h}^{\star}= \prod_{k}xk,h$ ’ により構成する。

界面での連続性を拘束条件とする Lagrange 乗数法を用いる。
, による。
を作り, 全界面に関しては積空間,$\Gamma_{m}$上で Lagrange 乗数の有限要素空間 $Y_{m,h}$$Y_{h}^{\star}= \prod_{m}Y_{m,h}$(22) に対応して, 各時刻ステップには次の問題を解く。
$(\phi_{h}^{n+1}, \eta_{h^{+1}})n\in Xh\cross Y^{\star}h\star$

5.3$( \frac{\phi_{h}^{n+1}-\phi_{h}^{n}}{\triangle t})\xi_{h})+j^{*}(\xi_{h,\eta_{h}^{n}})$$=$$-a_{1}^{\star}(u_{h}^{n+1},\phi h’\xi_{h}n)$$j^{\star}(\phi^{n}h^{+}’\mu h)1$

$=$

$0$

,

$\mu_{h}\in \mathrm{Y}_{h}^{\star}$,$\xi_{h}\in x_{h}^{\star}$

,

,

(28)

(29)

有限要素の還択

以上のスキームのように、部分領域ごとに P2 有限要素空間を構成し Lagrange 乗数で接合する場合の有限要素の選択の方法の 1 つは,.

改良$\mathrm{P}2(\mathrm{h})$要素 (改良はモルタル法 [4] に基づき、界面端で連続 1 次)を用いることである。 この選択では、 Stokes 方程式などにおいて P2 要素が持つ最適な収束性が領域分割ス
キームにおいても保持されることが証明されている。 しかし [8] で得られた経験からの類推では, 計算速度
) 要素や PO(h) 要の観点から自由度のより少ない有限要素を選ぶことも可能性がある。 具体的には、$\mathrm{P}1_{\backslash }^{(}\mathrm{h}$素などである。 今後、 検証を進めたいと考えている。

6

おわりに

2 次元非圧縮 Navier-Stokes 方程式の領域分割有限要素解法 $[7, 8]$ を発展させて、 Lagrange 乗数の基底
関数の選択の拡大や、流速・圧力の基底関数の異なる選択と合わせて考えることにより、 より自由度の高い領域分割への適用や、 2 流体問題への適用を検討した。

後半の節のいくつかの内容は研究の中間報告的な内容であるが、今後、 コードの開発及び並列計算機への実装を進めて、 実際の問題に応用していきたいと考えている。

謝辞
本研究は、次の先生方との議論の中でさまざまな有益な示唆を受けている。 ここにお礼を申し上げたい :九州大学田端正久教授、 九州大学金山寛教授、 富山大学大森克史教授。

21

本研究は次の科学研究費補助金の援助を受けた。基盤研究$(\mathrm{A})(1)$「偏微分方程式の無限精度数値シミ「領$\iota$レーションに関する研究」 (代表者徳島大学今井仁司教授)(課題番号 10354001) 、基盤研究$(\mathrm{B})(1)$

域分割法による計算力学の新解法の開発とその数学的検証」(代表者日本大学中村正彰助教授)(課題番号10555023).
この研究に使用した並列計算機は、 広島大学理学部の Intel Paragon(56 プロセッサ) 九州大学
大型計算機センターの富士通 VPP-700(56 プロセッサ) 、茨城大学総合情報処理センターの日立
プロセッサ) である。
$\mathrm{X}\mathrm{P}/\mathrm{s}$、$\mathrm{S}\mathrm{R}2201(16$参考文献[1] Y. Achdou and.Kuznetsov, Algorithm for a non conforming domain decomposition method, Tech. Rep.
296, Ecole Polytechnique, 1994.
$\mathrm{Y}.\mathrm{A}$[2] Y. Achdou and O. Pironneau, A fast solver for Navier-Stokes equations in the laminar regime using mortar
finite element and boundary element methods, SIAM. J. Numer. Anal., 32, 985-1016, 1995.

[3] M. Bercovier and O. Pironneau, Error estimates for finite element method solution of the Stokes problem inthe primitive variable, Numer. Math., 33, 211-224, 1979.

[4] C. Bernardi, Y. Maday, and A. Patera, A new nonconforming approach to domain decomposition: the mortar element method, in: H. Brezis and.Lions, Nonlinear Partial Differential Equations and their
Applications, Vol.XI, pp.13-51, Longman Sci.&Tech.,, 1994.
$\mathrm{J}.\mathrm{L}$$\mathrm{e}\mathrm{d}\mathrm{s}.$$\mathrm{U}\mathrm{K}$[5] C.Bernardi and Y. Maday, Spectral, spectral element and mortar element methods, preprint.

[6].Brackbill,.Kothe and C.Zemach, A continuum method for modeling surface tension, J. Comp. Phys.,
100, pp.335-354, 1992.

$\mathrm{J}.\mathrm{U}$$\mathrm{D}.\mathrm{B}$[7] S. Fujima, An upwind finite element scheme for the Navier-Stokes equations and its domain decomposihon
algorithm, Ph.D. thesis, Hiroshima University, 1997.
[8] S. Fujima, Iso-P2
J. Mandel et al.

-element method for the Navier-Stokesin, Domain Decomposition Methods 10, $AMS$ Contemporary Mathematics Vol.218, pp.246-$\mathrm{P}1/\mathrm{P}1/\mathrm{p}\mathrm{l}\mathrm{d}_{0}\mathrm{m}\mathrm{a}\mathrm{i}\mathrm{n}-\mathrm{d}\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{p}_{\mathrm{o}\mathrm{S}}\mathrm{i}\mathrm{t}\mathrm{i}_{\mathrm{o}\mathrm{n}}/\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{e}$$\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{S}_{)}$$\mathrm{e}\mathrm{d}\mathrm{s}.$253, 1998.

[9] V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms,
Springer, Berlin, 1986.
[10] R. Glowinski,.Pan and J. Periaux, A one shot domain
domain method for the
Navier-Stokes equations, Domain Decomposition Methocls in Scientific and Engineering Computing, Proc. 7th
Int. Conf. on Domain Decomposition,.Keyes and J.Xu
, Contemporary Mathematics, vol.180,.S.,
1994, 211-220.
$\mathrm{T}.\mathrm{W}$$\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}/\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{o}\mathfrak{U}\mathrm{s}$$\mathrm{D}.\mathrm{E}$$\mathrm{e}\mathrm{d}\mathrm{s}.$$\mathrm{A}.\mathrm{M}$[11] P. M. Gresho, S. T. Chan, R. L. Lee and C. D. Upson, A modified finite element method for solving the timedependent, incompressible Navier-Stokes equations, Part 1: Theory, Int. J. $Num$. Meth. Fluids, 4, 557-598,

1984.

[12] R.Jyotsna and.Vanka, Multigrid calculation of steady, viscous flow in a triangular cavity, J. Comp. Phys.,Vol.122, pp.107-117, 1995.
$\mathrm{S}.\mathrm{P}$[13]$\mathrm{H}.\mathrm{K}$.Moffatt, Viscous and resistive eddies near a sharp corner, J. Fluid Mech., Vol.18, pp.1-18, 1964.

[14] K. Ohmori and H. Kawarada, A sharp interface capturing technique in the finite element approximation for
two-fluid flows, in R. Salvi ed., Nanier-Stokes Equations: Theory and Numerical Methods, Pitman Research
Notes in Math. Vol. 388, pp.310-321, 1998
[15] A. Pothen,.Simon and.Liou, Partitioning sparse matrices with eigenvectors of graphs, SIAM .
Matfix Analysis and Appl., Vol.11, pp.430-452, 1990.
$\mathrm{H}.\mathrm{D}$[16].Ribbens and1994.
$\mathrm{C}.\mathrm{J}$$\mathrm{K}.\mathrm{P}$$\mathrm{C}.\mathrm{Y}$$J$.Wang, Steady viscous flow in a triangular cavity, J. Comp. Phys. Vol.112, pp.173-181,

22[17] A. Suzuki, Implementation of domain decomposition methods on parallel computer ADENART, Parallel, Elsevier, 1995,
Computational Fluid Dynamics: New Algorithms and Applications, N. Satofuka et al.
pp.231-238.
$\mathrm{e}\mathrm{d}\mathrm{s}$[18] M. Tabata and S. Fujima, An upwind finite element scheme forMeth. Fluids, Vo1.12, pp.305-322, 1991.

$\mathrm{h}\mathrm{l}\mathrm{g}\mathrm{h}- \mathrm{R}\mathrm{e}\mathrm{y}\mathrm{n}\mathrm{o}$]-number flows, Int. J. $Num$ .

$\mathrm{d}\mathrm{s}$