14
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

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

Last updated at Posted at 2020-01-31

前回のおさらい

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

{\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というライブラリーをよく使っています。
また、大きく計算がずれているかもしれないので間違っていたら指摘してください。

quad_element.png

図のように局所座標系に変換した直角二等辺三角形内の物理量$u(\xi,\eta)$を6つの節点(上の図の赤い丸)での値で表すことを考える。
$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]

骨が折れますね…

14
6
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
14
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?