2次元Poisson方程式に対する有限要素法

4年16組48番 南木集 2005年2月21日

– Typeset by FoilTEX –

問題： 2次元Poisson方程式 (Dirichlet境界条件＆Neumann境界条件)

次式を満たすuを求めよ。

−Δu = f ( (x, y) ∈ Ω)  (1)
u(x, y) = g_{1} ( (x, y) ∈ Γ_{1})  (2)
\frac{∂u}{∂n} = g_{2} ( (x, y) ∈ Γ_{2})  (3)

ただし、 Ωは２次元領域、 ΓをΩの境界、 Γ_{1}, Γ_{2} はΓの一部でΓ_{1} ∪ Γ_{2} = Γかつ Γ_{1} ∩ Γ_{2} = ∅とし、 f, g_{1}, g_{2} はそれぞれΩ, Γ_{1}, Γ_{2} で与えられた関数とする。また nはΓ_{2} 上での外向き単位法線ベクトルである

弱形式の導出

v は v = 0 (on Γ_{1})を満たす任意関数とする Poisson方程式(1)の両辺にv をかけ両辺をΩで積分する

− ∫∫_{Ω}(Δu)vdx = ∫∫_{Ω} fvdx

左辺をGreenの公式を用いて計算すると

∫∫_{Ω}(\frac{∂u}{∂x}\frac{∂v}{∂x} + \frac{∂u}{∂y}\frac{∂v}{∂y})dxdy − ∫_{Γ} v\frac{∂u}{∂n}dγ = ∫∫_{Ω} fvdx

– Typeset by FoilTEX – 2

Neumann境界条件とv = 0 (on Γ_{1}) より次式の弱形式を得る。

∫∫_{Ω}(\frac{∂u}{∂x}\frac{∂v}{∂x} + \frac{∂u}{∂y}\frac{∂v}{∂y})dxdy = ∫∫_{Ω} fvdx + ∫_{Γ_{2}} g_{2}vdγ (dγは線要素)

< u, v > = ∫∫_{Ω}(\frac{∂u}{∂x}\frac{∂v}{∂x} + \frac{∂u}{∂y}\frac{∂v}{∂y})dxdy
(u, v) = ∫∫_{Ω} uvdxdy
[g_{2}, v] = ∫_{Γ_{2}} g_{2}vdγ

という記号を導入すれば弱形式は

< u, v >= (f, v) + [g_{2}, v] となる

有限要素分割

領域を有限個の小領域に分割する各小領域を要素と呼ぶ要素の頂点を節点と呼ぶ

右図の場合各要素の節点の個数は3個

近似関数

各要素内における近似関数の形をα_{1}, α_{2}, α_{3} をパラメータとして

û = α_{1} + α_{2}x + α_{3}y

とする。すなわちûはx, y の1次式である。

要素間の連続性を保障するため要素[ e ]の節点における近似関数ûの値u_{i} を節点パラメータとして用いる。このとき次式が成立する

節点 i で û = u_{j} = α_{1} + α_{2}x_{j} + α_{3}y_{j} (1 ≤ j ≤ 3)

すなわち

(\begin{matrix}1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix}) (\begin{matrix}α_{1} \\ α_{2} \\ α_{3} \end{matrix}) = (\begin{matrix}u_{1} \\ u_{2} \\ u_{3} \end{matrix})

(\begin{matrix}α_{1} \\ α_{2} \\ α_{3} \end{matrix}) = (\begin{matrix}1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix})^{−1} (\begin{matrix}u_{1} \\ u_{2} \\ u_{3} \end{matrix})

= \frac{1}{|\begin{matrix}1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix}|} (\begin{matrix}x_{2}y_{3} − x_{3}y_{2} & x_{3}y_{1} − x_{1}y_{3} & x_{1}y_{2} − x_{2}y_{1} \\ y_{2} − y_{3} & y_{3} − y_{1} & y_{1} − y_{2} \\ x_{3} − x_{2} & x_{1} − x_{3} & x_{2} − x_{1} \end{matrix}) (\begin{matrix}u_{1} \\ u_{2} \\ u_{3} \end{matrix})

(i, j, k)を(1, 2, 3)の偶置換としてa_{i}, b_{i}, c_{i} を次のように置く

a_{i} = \frac{x_{j}y_{k} − x_{k}y_{j}}{|\begin{matrix}1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix}|} , b_{i} = \frac{y_{j} − y_{k}}{|\begin{matrix}1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix}|} , c_{i} = \frac{x_{k} − x_{j}}{|\begin{matrix}1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix}|}

– Typeset by FoilTEX – 6

a_{i}, b_{i}, c_{i} を用いるとûは次式で書ける

要素[ e ]で û = ∑_{i=1}^{3}(a_{i} + b_{i}x + c_{i}y)u_{i}

= ∑_{i=1}^{3} L_{i}u_{i}

ただし、 L_{i} = a_{i} + b_{i}x + c_{i}y 各要素内ではx, y の1次式で、領域Ωでは連続な区分多項式の近似関数ûが得られた

Dirichlet 境界条件の近似については対応する境界上の節点での節点パラメータを関数g_{1} (on Γ_{1})の値に等しくする。

v̂も同様にv̂ = 0 (on Γ_{1})を満たすように構成する

要素係数行列の計算

< u, v >_{e_{k}} =^{def} ∫∫_{e_{k}}(\frac{∂u}{∂x}\frac{∂v}{∂x} + \frac{∂u}{∂y}\frac{∂v}{∂y})dxdy , (f, v)_{e_{k}} =^{def} ∫∫_{e_{k}} fvdxdy

という記号を用いれば要素数N、要素をe_{k} (1 ≤ k ≤ N)として弱形式 < û, v̂ >= (f, v̂) + [g_{2}, v] は次式になる

∑_{k=1}^{N} < û, v̂ >_{e_{k}}= ∑_{k=1}^{N}(f, v̂)_{e_{k}} + [g_{2}, v̂]

< û, v̂ >_{e_{k}}, (f, v̂)_{e_{k}} を計算する

< û, v̂ >_{e_{k}} = < ∑_{j=1}^{3} u_{j}L_{j}, ∑_{i=1}^{3} v_{i}L_{i} >_{e_{k}}= ∑_{i=1}^{3} ∑_{j=1}^{3} v_{i} < L_{j}, L_{i} >_{e_{k}} u_{j}

= ∑_{i=1}^{3} ∑_{j=1}^{3} v_{i}A_{ij}^{(e_{k})}u_{j}

(f, v̂)_{e_{k}} = (f, ∑_{i=1}^{3} v_{i}L_{i})_{e_{k}} = ∑_{i=1}^{3} v_{i}(f, L_{i})_{e_{k}} = ∑_{i=1}^{3} v_{i}f_{i}^{(e_{k})}

ただし、 A_{ij}^{(e_{k})} =< L_{j}, L_{i} >_{e_{k}} , f_{i}^{(e_{k})} = (f, L_{i})_{e_{k}} (1 ≤ i, j ≤ 3)

次のベクトルと行列を定義する

u_{e_{k}} = (\begin{matrix}u_{1} \\ u_{2} \\ u_{3} \end{matrix}) , v_{e_{k}} = (\begin{matrix}v_{1} \\ v_{2} \\ v_{3} \end{matrix}) , f_{e_{k}} = (\begin{matrix}f_{1}^{(e_{k})} \\ f_{2}^{(e_{k})} \\ f_{3}^{(e_{k})} \end{matrix})

A_{e_{k}} = (\begin{matrix}A_{11}^{(e_{k})} & A_{12}^{(e_{k})} & A_{13}^{(e_{k})} \\ A_{21}^{(e_{k})} & A_{22}^{(e_{k})} & A_{23}^{(e_{k})} \\ A_{31}^{(e_{k})} & A_{32}^{(e_{k})} & A_{33}^{(e_{k})} \end{matrix})

するとA_{e_{k}} は対称行列で

< û, v̂ >_{e_{k}}= v_{e_{k}}^{T}A_{e_{k}}u_{e_{k}} , (f, v̂)_{e_{k}} = v_{e_{k}}f_{e_{k}} となる

直接剛性法

Ωに節点番号を与え節点数をm個とする。前ページで与えたベクトル、行列を m次元に拡大する

u = ( u_{1} ; u_{2} ; ⋮ ; u_{m}) , v = ( v_{1} ; v_{2} ; ⋮ ; v_{m}) とおく

有限要素e_{k} (1 ≤ k ≤ N)を１つ取って考える。 e_{k} の局所的な節点番号1^{∗}, 2^{∗}, 3^{∗} に対してΩの全体節点番号をそれぞれ i, j, k が対応しているとする。次のベクトル、行列を定義する。

f_{k}^{∗} = (0 … 0f_{1}0 … 0f_{2}0 … 0f_{3}0 … 0)^{T} f_{1}は第i成分、 f_{2}は第j 成分、 f_{3}は第k 成分である

A_{k}^{∗} = ( \begin{matrix}… & & & & & & … \\ & A_{11}^{(e_{k})} & … & A_{12}^{(e_{k})} & … & A_{13}^{(e_{k})} \\ … & … & … & … & … & … & … \\ … & A_{21}^{(e_{k})} & … & A_{22}^{(e_{k})} & … & A_{23}^{(e_{k})} & … \\ … & … & … & … & … & … & … \\ & A_{31}^{(e_{k})} & … & A_{31}^{(e_{k})} & … & A_{33}^{(e_{k})} \\ … & & & & & & … \end{matrix} )

示していない成分はすべて0 この記号を用いると < û, v̂ >= v^{T} A_{k}^{∗}u (1 ≤ k ≤ N) よって弱形式は

∑_{k=1}^{N} v^{T} A_{k}^{∗}u = ∑_{k=1}^{N} v^{T} f_{k}^{∗}

v がv = 0 (on Γ_{1})を満たす任意関数であることを考えると上式においてΓ_{1} 上にくる節点番号の行を除いた

A^{∗∗}u = f^{∗∗} が成立する

ここでN_{k} (1 ≤ k ≤ N)を節点番号とすると

A^{∗∗} = A^{∗}の第k 行(N_{k} ∈ Γ_{1}なるk)を除いた行列
f^{∗∗} = ∑_{k=1}^{N} f_{k} の第k 成分(N_{k} ∈ Γ_{1}なるk)を除いた縦ベクトル

またΓ_{1} 上ではû = g_{1} （既知）であるから、これを代入して消去できる。 A^{∗∗} を列ベクトルでA^{∗∗} = (a_{1}a_{2} … a_{N})と表示すると

∑_{k=1}^{N} a_{k}u_{k} = f^{∗∗}

左辺をΓ_{1} に属す部分と属さない部分に分解して移行すると

∑_{N_{k}∉Γ_{1}なるi} a_{i}u_{i} = f^{∗∗} − ∑_{N_{k}∈Γ_{1}なるi} a_{i}u_{i}

これより

Au^{∗} = f

ここで

A = A^{∗∗}から第i列(N_{k} ∈ Γ_{1}なるk)を除いた正方行列
u^{∗} = uから第k 成分(N_{k} ∈ Γ_{1}なるk)を除いた縦ベクトル
f = f^{∗∗} − ∑_{N_{k}∈Γ_{1}なるk} a_{i}u_{i}

この方程式が解くべき方程式である

– Typeset by FoilTEX – 15

数値実験

実際に数値実験をした問題は

−Δu = 1 ( (x, y) ∈ Ω = (0, 1) × (0, 1))
u(x, y) = 0 ( (x, y) ∈ Γ_{1})
\frac{∂u}{∂n} = 0 ( (x, y) ∈ Γ_{2})

ただし、 Γ_{1} はx = 0またはy = 0で与えられる2辺、 Γ_{2} はx = 1またはy = 1で与えられる２辺とする

Ω Γ_{1}  Γ_{2} (0,0) (0,1)  (1,0) (1,1)

– Typeset by FoilTEX – 16

実験結果:等高線

– Typeset by FoilTEX – 17

有限要素法とは

• 問題の微分方程式に対する弱形式を求める

• 領域Ωを任意の有限個(ｎ個)の小領域(要素)Ω̂_{i} に分割

Ω ≅ ⋃_{i=1}^{n} Ω̂_{i}

• 各要素において近似関数を構成し弱形式への寄与分を求め、直接剛性法により全体的な近似方程式を構成する

• 得られた連立1次方程式を数値的に解く。得られた値が有限要素解である