有限要素法の実装(1)　ひずみ変位行列(1)

　初めまして食卓塩と名乗っているものです。
　今回はひずみ変位行列の式展開について説明します。そんなのどこでもやってるだとと思われますが、大概の書籍はグダグダと式を無意味に書きつけているので、エレガント(実は物ぐさ)な方法を紹介します。
　ちなみにこの知見を活かしたものはGithubの私のスペースにあるので、暇な人はご参照くださいませ。
　初めに歪ベクトルがあり、これを最終的なひずみ変位行列にするのですが、先にも書いたとおり大概の書籍はエレガントの欠片もないグダグダっぷりで、人によっては何が何だか分からない状況になっているものと推察します。
　数学というのは人によって得手不得手が大きく分かれる学問ですが、それを抜きにしても何の考えもなしにいきなり形状関数を全体座標で偏微分したら理解できない人が大勢出るのではないかと思われます。そのあとにヤコビアンが出てきますが、誰が見ても超展開としか言いようのない唐突さで出てこられたら混乱に拍車がかかるでしょう。
　なお本稿では具体的な形状関数については触れません。アイソパラメトリック要素を使う限り同じ次元の要素を使用する限り、同じ式を使うからです。
　前置きはこれくらいにしてさっそく始めます。初めに二次元平面要素のひずみ変位行列の導出を行いましょう。

```math
\left[\begin{array}{c}
\epsilon_x \\ \epsilon_y \\ \epsilon_{xy}
\end{array}\right]=
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{x}} \normalsize + 
\Large\frac{\partial{u}}{\partial{y}}
\end{array}\right] \cdots (1)
```
　ここまでならだれでもご存知の工学ひずみの式です。$(1)$右辺を分解します。右辺最初の項について名称が無いようですが、ここでは符号行列と称しておきます。

```math
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{x}} \normalsize+
\Large\frac{\partial{u}}{\partial{y}}
\end{array}\right] =
\left[\begin{array}{c}
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 1 & 1 & 0 
\end{array}\right] 
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{x}} \\
\Large\frac{\partial{u}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{y}}
\end{array}\right]
\cdots (2)
```

　次に$(2)$右辺の第2項について変形を行います。アイソパラメトリック要素は自然座標系と言われるもので定義されています。ここでは座標系に$r_1\cdot$$r_2$を用います。
　$(2)$の左辺一行目の$\partial{u}/\partial{x}$は偏微分の連鎖則を用いて次のようになります。

```math
\frac{\partial{u}}{\partial{x}} = 
\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}+
\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}}
 \cdots(3)
```

　以下同様に$(2)$の右辺第2項を変形すると次のようになります。

```math
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{x}} \\
\Large\frac{\partial{u}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{y}}
\end{array}\right] \normalsize=
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}}
\end{array}\right]
\cdots (4)
```

　ここで左辺第二項に注目すると$\partial{r_1}/\partial{x}$と$\partial{r_2}/\partial{x}$といった項が出ています。解かっている人は分かるかと思いますが、これらはヤコビアンの逆行列の構成要素になります。
　というわけで式$(4)$第2項を変形すると次のようになります。

```math
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}}
\end{array}\right] =
\left[\begin{array}{c}
\left[\begin{array}{c}
\Large\frac{\partial{r_1}}{\partial{x}}&
\Large\frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{r_1}}{\partial{y}}&
\Large\frac{\partial{r_2}}{\partial{y}} 
\end{array}\right]
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \\
\Large\frac{\partial{v}}{\partial{r_2}}
\end{array}\right] \\
\left[\begin{array}{c}
\Large\frac{\partial{r_1}}{\partial{x}}&
\Large\frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{r_1}}{\partial{y}}&
\Large\frac{\partial{r_2}}{\partial{y}} 
\end{array}\right]
\left[\begin{array}{c}
\Large\frac{\partial{v}}{\partial{r_1}} \\
\Large\frac{\partial{v}}{\partial{r_2}}
\end{array}\right]
\end{array}\right]\cdots (5)
```

　ここで、右辺第2項左側の行列を次のように置き換えます。

```math
\mathbf J^{-1}=
\left[\begin{array}{c}
\Large\frac{\partial{r_1}}{\partial{x}}&
\Large\frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{r_1}}{\partial{y}}&
\Large\frac{\partial{r_2}}{\partial{y}} 
\end{array}\right] \cdots (6)
```

　そうすると、式$(5)$右辺は次のようになります。

```math
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}}
\end{array}\right] =
\left[\begin{array}{c}
\mathbf J^{-1}
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \\
\Large\frac{\partial{v}}{\partial{r_2}}
\end{array}\right] \\
\mathbf J^{-1}
\left[\begin{array}{c}
\Large\frac{\partial{v}}{\partial{r_1}} \\
\Large\frac{\partial{v}}{\partial{r_2}}
\end{array}\right]
\end{array}\right] \cdots (7)
```

　アイソパラメトリック要素の場合、要素の変位は次のように表されます。なお、右上のインデックスは形状関数の番号を表し、アインシュタインの総和規約に基づいて合計記号は省略しています。

```math
u=N^iu^i \cdots(8)
```

偏微分は当然次のようになります。

```math
\frac{\partial{u}}{\partial{r_j}}=\frac{\partial{N^iu^i}}{\partial{r_j}} \cdots(9)
```

　$u^i$は定数ですから、式の外に出すことができるので$(9)$式は次のようになります

```math
\frac{\partial{u}}{\partial{r_j}}=\frac{\partial{N^i}}{\partial{r_j}} u^i\cdots(10)
```

　式$(8)$を利用して式$(5)第二項の行列を次のように変形します。

```math
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \\
\Large\frac{\partial{v}}{\partial{r_2}}
\end{array}\right] =
\left[\begin{array}{c}
\Large\frac{\partial{N^i}}{\partial{r_1}} \normalsize u^i \\
\Large\frac{\partial{N^i}}{\partial{r_2}} \normalsize u^i
\end{array}\right] =
\left[\begin{array}{c}
\Large\frac{\partial{N^i}}{\partial{r_1}} \\
\Large\frac{\partial{N^i}}{\partial{r_2}}
\end{array}\right]  \normalsize \left \{ u^i \right \} \cdots (11)
```

 更に以下のように置き換えます。

```math
\mathbf H=
\left[\begin{array}{c}
\Large\frac{\partial{N^i}}{\partial{r_1}} \\
\Large\frac{\partial{N^i}}{\partial{r_2}}
\end{array}\right] \cdots (12)
```

 これらを利用すると式$(5)$は最終的に次のように変形できます。

```math
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}}
\end{array}\right] =
\left[\begin{array}{c}
\mathbf J^{-1} \mathbf H \normalsize \left \{ u^i \right \} \\
\mathbf J^{-1} \mathbf H \normalsize \left \{ v^i \right \}
\end{array}\right] \cdots (13)
```

　実際の計算では当然のことですが、変位の項は一つにしたいところです。よって式$(13)$はこのような形になります。

```math
\left[\begin{array}{c}
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{u}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{u}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{x}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{v}}{\partial{r_1}} \frac{\partial{r_1}}{\partial{y}}
\normalsize+
\Large\frac{\partial{v}}{\partial{r_2}} \frac{\partial{r_2}}{\partial{y}}
\end{array}\right] =
\left[\begin{array}{c}
\mathbf J^{-1} \mathbf H \normalsize&  0 \\
0& \mathbf J^{-1} \mathbf H
\end{array}\right]
\left[\begin{array}{c}
u^i\\ v^i
\end{array}\right]
 \cdots (13)
```

　最終的に式$(2)$と組み合わせると、ひずみ変位行列を得ることができます。

```math
\left[\begin{array}{c}
\epsilon_x \\ \epsilon_y \\ \epsilon_{xy}
\end{array}\right]=
\left[\begin{array}{c}
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 1 & 1 & 0 
\end{array}\right] 
\left[\begin{array}{c}
\mathbf J^{-1} \mathbf H \normalsize&  0 \\
0& \mathbf J^{-1} \mathbf H
\end{array}\right]
\left[\begin{array}{c}
u^i\\ v^i
\end{array}\right]
\cdots (14)
```

　式の中でヤコビアンの逆行列が出てきていますが、これを数値的に求めるにはヤコビアンの陽な形が必要です。この行列は次のように表されました。

```math
\mathbf J^{-1}=
\left[\begin{array}{c}
\Large\frac{\partial{r_1}}{\partial{x}}&
\Large\frac{\partial{r_2}}{\partial{x}} \\
\Large\frac{\partial{r_1}}{\partial{y}}&
\Large\frac{\partial{r_2}}{\partial{y}} 
\end{array}\right] \cdots (6)
```

　ヤコビアンは次のようになります。

```math
\mathbf J=
\left[\begin{array}{c}
\Large\frac{\partial{x}}{\partial{r_1}}&
\Large\frac{\partial{y}}{\partial{r_1}} \\
\Large\frac{\partial{x}}{\partial{r_2}}&
\Large\frac{\partial{y}}{\partial{r_2}} 
\end{array}\right] \cdots (15)
```

　要素内の座標も変位同様形状関数を使って補間できますので次のようになります。右上のインデックスは例によって形状関数と節点の番号を指し、
右下のインデックスは座標番号をさします。$1=x,2=y$となります。

```math
x_i=N^kx^k_i \cdots(16)
```

偏微分も同じように表せます。また、座標値そのものも定数ですので、偏微分の式の外に出すことができます。

```math
\begin{eqnarray*}
\frac{\partial{x_i}}{\partial{r_j}}&=&\frac{\partial{N^kx^k_i}}{\partial{r_j}} \\
&=&\frac{\partial{N^k}}{\partial{r_j}} x^k_i
\end{eqnarray*} \cdots(17)
```

　式$(17)$を参考にすると式$(15)$は次のようになります。


```math
\begin{eqnarray*}
\mathbf J &=&
\left[\begin{array}{c}
\Large\frac{\partial{x}}{\partial{r_1}}&
\Large\frac{\partial{y}}{\partial{r_1}} \\
\Large\frac{\partial{x}}{\partial{r_2}}&
\Large\frac{\partial{y}}{\partial{r_2}} 
\end{array}\right] \\ &=&

\left[\begin{array}{c}
\Large\frac{\partial{N^i}}{\partial{r_1}} x^i &
\Large\frac{\partial{N^i}}{\partial{r_1}} y^i \\
\Large\frac{\partial{N^i}}{\partial{r_2}} x^i &
\Large\frac{\partial{N^i}}{\partial{r_2}} y^i 
\end{array}\right] \\ &=& 

\left[\begin{array}{c}
\left[\begin{array}{c}
\Large\frac{\partial{N^i}}{\partial{r_1}} \\
\Large\frac{\partial{N^i}}{\partial{r_2}}
\end{array}\right] x^i
\left[\begin{array}{c}
\Large\frac{\partial{N^i}}{\partial{r_1}} \\
\Large\frac{\partial{N^i}}{\partial{r_2}}
\end{array}\right] y^i
\end{array}\right] \\ &=& 

\left[\begin{array}{c}
\Large\frac{\partial{N^i}}{\partial{r_1}} \\
\Large\frac{\partial{N^i}}{\partial{r_2}}
\end{array}\right]
\left[\begin{array}{c}
x^i & y^i
\end{array}\right] \\ &=&
\mathbf H \mathbf X
\end{eqnarray*}
\cdots (18)
```

 ひずみ変位行列の変位の項はX軸方向変位$u$とY軸方向変位$v$が形状関数の数が縦に並びますが、
式$(18)$の座標マトリックス$\mathbf X$は行ごとに各節点座標が入ることになります。
　これらをまとめると、ひずみ変位行列の計算は次のようになります。

* ヤコビアン$\mathbf J$の計算を行う
* ヤコビアンの逆行列$\mathbf J^{-1}$と行列値$det(\mathbf J)$を求める。行列値は逆行列を求めるときの副産物として出る。
* $\mathbf J^{-1}\mathbf H$を計算する
* 符号行列に従って$\mathbf J^{-1}\mathbf H$の要素をひずみ変位行列に配置する

2019.12.30 文書の修正
参考文献
+ 平居 孝之 パソコン3次元有限要素法 1990年 共立出版