有限要素法による二次元ポアソン方程式の解法(with 三角形二次要素)（その２）

前回のおさらい

前回の記事→（その１）要素行列の作成と局所座標系の導入まで
次の記事→　(その３) 境界条件全体行列作成
要素ごとに離散化した式

{\begin{align}
\sum_{e=1} \left[
w_j^e u_i  \mathbf{K^e} 
+w_j^e   \mathbf{f^e}
\right] 
=
\sum_{l=1} w_j^e q_i^e \int_{\Gamma_{N_l}} \psi_i^e \psi_j^e \mathrm{d\Gamma}\\
\mathbf{K^e} = \int_{\Omega_e}  \nabla \phi_i^e \nabla \phi_j^e \mathrm{d}\Omega, ~~~
\mathbf{f^e} = f_i^e \int_{\Omega_e}  \phi_i^e \phi_j^e \mathrm{d}\Omega 
\end{align}
}

局所座標系へ変換

\begin{align}
&J_e^{-1} (J_e^{-1})^T |J_e|= 
\begin{pmatrix}
a&b\\
c&d
\end{pmatrix}\\

&\mathbf{K^e} = a \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
+ b \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega
+ c \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
+ d \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega

\end{align}

ところで、

J_e^{-1} (J_e^{-1})^T |J_e|

は対称行列なんで$b=c$になりますね。

内挿関数(形状関数)を求めていく

この章から計算がえげつないものになっていきます。皆さん関数電卓なり、mathmaticaなりを使って楽をしましょう。私はpythonのsympyというライブラリーをよく使っています。
また、大きく計算がずれているかもしれないので間違っていたら指摘してください。

図のように局所座標系に変換した直角二等辺三角形内の物理量$u(\xi,\eta)$を６つの節点（上の図の赤い丸）での値で表すことを考える。
$i$番目の節点での値を$u_i(i=1\ldots6)$とすると、要素内の物理量の分布は節点の数と同じ数の内挿関数$\phi_i(\xi,\eta)$を用いて以下のように表せる。

u(\xi,\eta) = u_i \phi_i

ここで、

u(\xi,\eta) = a\xi^2 + b\xi\eta + c\eta^2 + d\xi + e\eta + f

と表されるものとする。
更に、$u(0,0) = u_1, u(1,0) = u_2, \ldots , u(0,0.5)=u_6$
である。
そのため、これらを適用すると、以下のような方程式が成り立つ。

\begin{bmatrix}
u_1\\u_2\\u_3\\u_4\\u_5\\u_6
\end{bmatrix}
=
 \left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 1\\1 & 0 & 0 & 1 & 0 & 1\\0 & 0 & 1 & 0 & 1 & 1\\\frac{1}{4} & 0 & 0 & \frac{1}{2} & 0 & 1\\\frac{1}{4} & \frac{1}{4} & \frac{1}{4} & \frac{1}{2} & \frac{1}{2} & 1\\0 & 0 & \frac{1}{4} & 0 & \frac{1}{2} & 1\end{matrix}\right]
\begin{bmatrix}
a\\b\\c\\d\\e\\f
\end{bmatrix}

$a,b,\ldots,f$を求めるために逆行列を求める。すると、以下のようになる。

\begin{align}
\begin{bmatrix}
a\\b\\c\\d\\e\\f
\end{bmatrix}
&=
\left[\begin{matrix}2 & 2 & 0 & -4 & 0 & 0\\4 & 0 & 0 & -4 & 4 & -4\\2 & 0 & 2 & 0 & 0 & -4\\-3 & -1 & 0 & 4 & 0 & 0\\-3 & 0 & -1 & 0 & 0 & 4\\1 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]
\begin{bmatrix}
u_1\\u_2\\u_3\\u_4\\u_5\\u_6
\end{bmatrix}\\
\end{align}

これを、

u(\xi,\eta) = a\xi^2 + b\xi\eta + c\eta^2 + d\xi + e\eta + f

に代入すると、

\begin{align}
u(\xi,\eta) &=
\begin{bmatrix}
a&b&c&d&e&f
\end{bmatrix}
 \begin{bmatrix}
\xi^2 \\ \xi\eta \\ \eta^2 \\ \xi \\ \eta \\ 1
\end{bmatrix}\\
&=
\begin{bmatrix}
u_1&u_2&u_3&u_4&u_5&u_6
\end{bmatrix}
\left[\begin{matrix}2 & 2 & 0 & -4 & 0 & 0\\4 & 0 & 0 & -4 & 4 & -4\\2 & 0 & 2 & 0 & 0 & -4\\-3 & -1 & 0 & 4 & 0 & 0\\-3 & 0 & -1 & 0 & 0 & 4\\1 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]^T
 \begin{bmatrix}
\xi^2 \\ \xi\eta \\ \eta^2 \\ \xi \\ \eta \\ 1
\end{bmatrix}
\end{align}

そして、内挿関数は以下の式を満たす。

u(\xi,\eta) = u_i \phi_i = \begin{bmatrix}
u_1&u_2&u_3&u_4&u_5&u_6
\end{bmatrix}
\begin{bmatrix}
\phi_1\\\phi_2\\\phi_3\\\phi_4\\\phi_5\\\phi_6
\end{bmatrix}\\

以上から、6つの内挿関数が求まる。

\begin{align}
\phi_1 &=  2 \xi^{2} + 4 \xi \eta + 2 \eta^{2} - 3 \xi  - 3 \eta + 1 \\
\phi_2 &= 2\xi^2-\xi\\
\phi_3 &= 2\eta^2 -\eta\\
\phi_4 &= 4 \xi^{2} - 4 \xi \eta + 4 \xi\\
\phi_5 &= 4\xi \eta\\
\phi_6 &= 4 \xi \eta - 4 \eta^{2} + 4 \eta\\
\end{align}

要素内積分

{\begin{align}
\mathbf{K^e} &= a \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
+ b \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega
+ c \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
+ d \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega
\end{align}
}

$i,j = 1...6$であるため、

\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega

等は$6\times6$の行列になります。

ここで、前の章で導出した内挿関数を行列を使って表現すると、

\begin{bmatrix}
\phi_1\\\phi_2\\\phi_3\\\phi_4\\\phi_5\\\phi_6
\end{bmatrix}=
\left[\begin{matrix}2 & 4 & 2 & -3 & -3 & 1\\2 & 0 & 0 & -1 & 0 & 0\\0 & 0 & 2 & 0 & -1 & 0\\-4 & -4 & 0 & 4 & 0 & 0\\0 & 4 & 0 & 0 & 0 & 0\\0 & -4 & -4 & 0 & 4 & 0\end{matrix}\right]
 \begin{bmatrix}
\xi^2 \\ \xi\eta \\ \eta^2 \\ \xi \\ \eta \\ 1
\end{bmatrix}
=
\mathbf{\Phi} \begin{bmatrix}
\xi^2 \\ \xi\eta \\ \eta^2 \\ \xi \\ \eta \\ 1
\end{bmatrix}

ここでは、行列部分を$\mathbf{\Phi}$と置きました。

これを使って、$\frac{\partial \phi_i^e}{\partial \xi}$を表すと、

\frac{\partial }{\partial \xi}
\begin{bmatrix}
\phi_1\\\phi_2\\\phi_3\\\phi_4\\\phi_5\\\phi_6
\end{bmatrix}=
\mathbf{\Phi}
 \begin{bmatrix}
2\xi \\ \eta \\ 0 \\ 1 \\ 0 \\ 0
\end{bmatrix}

そうすると、求めたい$\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega$は、以下のように表せる。

\begin{align}
\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
&=
\mathbf{\Phi}
\int_{\Omega_e}
 \begin{bmatrix}
2\xi \\ \eta \\ 0 \\ 1 \\ 0 \\ 0
\end{bmatrix}
 \begin{bmatrix}
2\xi & \eta & 0 & 1 & 0 & 0
\end{bmatrix}
\mathrm{d}\Omega
\mathbf{\Phi}^T
=
\mathbf{\Phi}
\int_{\Omega_e}
\begin{bmatrix}
4\xi^2  &2\xi\eta & 0 & 2\xi & 0 & 0 \\
2\xi\eta & \eta^2 & 0 & \eta & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
2\xi & \eta & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
\end{bmatrix}
\mathrm{d}\Omega
\mathbf{\Phi}^T
\end{align}

ここで、積分領域は直角三角形であったことを思い出し、$\int_{\Omega_e} \mathrm{d}\Omega$を具体的に考えると、

\int_{\Omega_e} \mathrm{d}\Omega
= \int_0^1 \left[\int_0^{1-\xi}  \mathrm{d}\eta\right] \mathrm{d}\xi

となります。

そして、積分を実際に計算すると、非積分関数は$\xi,\eta$の関数であるため、こんな感じの式が良く出てきます。

\int_0^1 \left[\int_0^{1-\xi} \xi^m \eta^n \mathrm{d}\eta\right] \mathrm{d}\xi = \frac{1}{n+1} \int_0^1 \xi^m (1-\xi)^{n+1} \mathrm{d}\xi = \cdots

右側の式は、大学一年生の頃に見た人が多いのではないでしょうか。

ベータ関数です。
これを使うと以下のように書けます。($\xi,\eta$が$x,y$にかわっていますが同じことです)

\begin{align}
\int_0^1 \left[\int_0^{1-x} x^m y^n \mathrm{d}y\right] \mathrm{d}x &= \frac{1}{n+1} B(m+1,n+2)\\
&= \frac{1}{n+1} \frac{m! (n+1)!}{(m+n+2)!}\\
&= \frac{m! n!}{(m+n+2)!}
\end{align}

これを使えば多少は楽に計算できるんじゃないでしょうか。

\begin{align}
\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
&=
\mathbf{\Phi}
\left[\begin{matrix}\frac{1}{3} & \frac{1}{12} & 0 & \frac{1}{3} & 0 & 0\\\frac{1}{12} & \frac{1}{12} & 0 & \frac{1}{6} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\\frac{1}{3} & \frac{1}{6} & 0 & \frac{1}{2} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]
\mathbf{\Phi}^T\\
&=
\left[\begin{matrix}\frac{1}{2} & \frac{1}{6} & 0 & - \frac{2}{3} & 0 & 0\\\frac{1}{6} & \frac{1}{2} & 0 & - \frac{2}{3} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\- \frac{2}{3} & - \frac{2}{3} & 0 & \frac{4}{3} & 0 & 0\\0 & 0 & 0 & 0 & \frac{4}{3} & - \frac{4}{3}\\0 & 0 & 0 & 0 & - \frac{4}{3} & \frac{4}{3}\end{matrix}\right]
\end{align}

これまでやったことを繰り返して、

\mathbf{K^e} = a \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
+ b \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega
+ c \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
+ d \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega

を求めるんです...

\begin{align}
 \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega
&= 
\mathbf{\Phi}
\int_{\Omega_e}
\begin{bmatrix}
2\xi \\ \eta \\ 0 \\ 1 \\ 0 \\ 0
\end{bmatrix}
 \begin{bmatrix}
0 & \xi & 2\eta & 0 & 1 & 0
\end{bmatrix}
\mathrm{d}\Omega
\mathbf{\Phi}^T\\
&=
\mathbf{\Phi}
\int_{\Omega_e}
 \begin{bmatrix}
0 & 2\xi^2 & 4\xi\eta & 0 & 2\xi & 0 \\
0 & \xi\eta & 2\eta^2 & 0 & \eta & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & \xi & 2\eta & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
\end{bmatrix}
\mathrm{d}\Omega
\mathbf{\Phi}^T\\
&=
\mathbf{\Phi}
\left[\begin{matrix}0 & \frac{1}{6} & \frac{1}{6} & 0 & \frac{1}{3} & 0\\0 & \frac{1}{24} & \frac{1}{6} & 0 & \frac{1}{6} & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{1}{6} & \frac{1}{3} & 0 & \frac{1}{2} & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]
\mathbf{\Phi}^T\\
&=
 \left[\begin{matrix}
\frac{1}{2} & 0 & \frac{1}{6} & 0 & 0 & - \frac{2}{3}\\\frac{1}{6} & 0 & - \frac{1}{6} & - \frac{2}{3} & \frac{2}{3} & 0\\0 & 0 & 0 & 0 & 0 & 0\\- \frac{2}{3} & 0 & 0 & \frac{2}{3} & - \frac{2}{3} & \frac{2}{3}\\0 & 0 & \frac{2}{3} & - \frac{2}{3} & \frac{2}{3} & - \frac{2}{3}\\0 & 0 & - \frac{2}{3} & \frac{2}{3} & - \frac{2}{3} & \frac{2}{3}
\end{matrix}\right]
\end{align}

 \int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
=  \left(\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega\right)^T

$\displaystyle
\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega
$に関しては、$\displaystyle
\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \xi} \frac{\partial \phi_j^e}{\partial \xi} \mathrm{d}\Omega
$と$\xi$と$\eta$が入れ替わっただけなので、同様に計算する。

\int_{\Omega_e} \frac{\partial \phi_i^e}{\partial \eta} \frac{\partial \phi_j^e}{\partial \eta} \mathrm{d}\Omega=
\left[\begin{matrix}\frac{1}{2} & 0 & \frac{1}{6} & 0 & 0 & - \frac{2}{3}\\0 & 0 & 0 & 0 & 0 & 0\\\frac{1}{6} & 0 & \frac{1}{2} & 0 & 0 & - \frac{2}{3}\\0 & 0 & 0 & \frac{4}{3} & - \frac{4}{3} & 0\\0 & 0 & 0 & - \frac{4}{3} & \frac{4}{3} & 0\\- \frac{2}{3} & 0 & - \frac{2}{3} & 0 & 0 & \frac{4}{3}\end{matrix}\right]

さて、$\mathbf{f^e}$についても同様です。

\mathbf{f^e} = 
f_i^e \int_{\Omega_e}  \phi_i^e(x,y) \phi_j^e(x,y) \mathrm{d}\Omega 
=
f_i^e \mathbf{\Phi} \int_{\Omega_e}
\begin{bmatrix}
\xi^2 \\ \xi\eta \\ \eta^2 \\ \xi \\ \eta \\ 1
\end{bmatrix}
\begin{bmatrix}
\xi^2 & \xi\eta & \eta^2 & \xi & \eta & 1
\end{bmatrix}
\mathrm{d}\Omega \mathbf{\Phi}^T |J_e|

(2020/2/24 ヤコビアンを入れ忘れていたことに気が付き加筆)
さて、積分の中身の行列は、ベータ関数を使って頑張って求められますね...

全部計算するとこれになります。

\left[\begin{matrix}\frac{1}{60} & - \frac{1}{360} & - \frac{1}{360} & 0 & - \frac{1}{90} & 0\\- \frac{1}{360} & \frac{1}{60} & - \frac{1}{360} & 0 & 0 & - \frac{1}{90}\\- \frac{1}{360} & - \frac{1}{360} & \frac{1}{60} & - \frac{1}{90} & 0 & 0\\0 & 0 & - \frac{1}{90} & \frac{4}{45} & \frac{2}{45} & \frac{2}{45}\\- \frac{1}{90} & 0 & 0 & \frac{2}{45} & \frac{4}{45} & \frac{2}{45}\\0 & - \frac{1}{90} & 0 & \frac{2}{45} & \frac{2}{45} & \frac{4}{45}\end{matrix}\right]

骨が折れますね…