このページの最終更新日 2020/10/25

アイソパラメトリック要素

1、四角形アイソパラメトリック要素とベクトル

 有限要素法に基づいた解析を行うとき、モデルの形状を2次元であれば三角形や四角形、 3次元であれば四面体や六面体に形状を分割して、数値計算を行うことになります。

 このとき、四角形や六面体の分割要素について「アイソパラメトリック要素」などという言葉が登場します。アイソパラメトリック要素の説明として、教科書では形状関数やら、自然座標系やらと出てきますが、理解するまでにつまづいてしまう人も多いのではないでしょうか。そこで、ここでは高校の数学のレベルでのベクトルを考えて、アイソパラメトリック要素というものは何なのかを解説してみたいと思います。また、アイソパラメトリック要素を用いた有限要素法の解析の中身についても触れていこうと思います。


まず次の四角形ABCDがあり、その内部に点Pがあったとします。この点Pを位置ベクトルを用いて表すと次のようになります。

pic.1

$\overrightarrow{OE}=(1-u)\overrightarrow{OA}+u\overrightarrow{OB}$

$\overrightarrow{OF}=(1-u)\overrightarrow{OD}+u\overrightarrow{OC}$

$\overrightarrow{OP}=(1-v)\overrightarrow{OE}+v\overrightarrow{OF}$

$\quad=(1-u)(1-v)\overrightarrow{OA}+u(1-v)\overrightarrow{OB}+uv\overrightarrow{OC}+(1-u)v\overrightarrow{OD}$


 ここまでは高校数学で問題ありませんよね?これを念頭に置いて、アイソパラメトリック要素とは何なのかを考えてみます。 さて、アイソパラメトリック要素とは、自然座標系 $(\xi, \eta)$を導入し、この自然座標系を実際の物理座標系に変換して定式化する要素であると説明されます。

 四角形要素の自然座標系とは、pic.2左のように$(-1\leqq \xi \leqq 1, \ -1\leqq \eta \leqq 1)$で表される正方形として定義されます。 これに対して実際のモデルをメッシュ分割するとき、完全に正方形として分割することは出来ません。実際の解析ではpic.2右のようにちょっと歪んだ物理座標系の四角形に分割して解析を行うことになります。

pic.2


さて、このことについて有限要素の教科書を読んでいると、次のような形状関数たるものが登場し、

$\begin{cases} N_1 (\xi, \ \eta) = \frac{1}{4}(1-\xi)(1-\eta) \\ N_2 (\xi, \ \eta) = \frac{1}{4}(1+\xi)(1-\eta) \\ N_3 (\xi, \ \eta) = \frac{1}{4}(1+\xi)(1+\eta) \\ N_4 (\xi, \ \eta) = \frac{1}{4}(1-\xi)(1+\eta) \\ \end{cases}$


自然座標 $(\xi, \eta)$から物理座標(x,y)の変換は次の式で表せます。などと書かれていたりします。

$\displaystyle x = \sum_{\alpha=1}^{4} N_{\alpha} x_\alpha = \frac{1}{4}(1-\xi)(1-\eta) x_1 + \frac{1}{4}(1+\xi)(1-\eta) x_2 + \frac{1}{4}(1+\xi)(1+\eta) x_3 + \frac{1}{4}(1-\xi)(1+\eta) x_4$  (1a)

$\displaystyle y = \sum_{\alpha=1}^{4} N_{\alpha} y_\alpha = \frac{1}{4}(1-\xi)(1-\eta) y_1 + \frac{1}{4}(1+\xi)(1-\eta) y_2 + \frac{1}{4}(1+\xi)(1+\eta) y_3 + \frac{1}{4}(1-\xi)(1+\eta) y_4$  (1b)

こんなことを言われてもサッパリ解らんという人も出てくるので、 まずこの自然座標系から物理座標系の変換を最初に挙げた高校レベルのベクトルと関連づけて説明しようと思う訳です。 ここで、物理座標系における四角形の辺を$(\xi, \eta)$を用いてpic.3のような比で分割し、座標Pの位置ベクトルを$P_1, P_2, P_3, P_4$を用いて表すことを考えると、

pic.3

$\displaystyle \overrightarrow{OQ_1}=\frac{1-\xi}{2}\overrightarrow{OP_1}+\frac{1+\xi}{2}\overrightarrow{OP_2}$

$\displaystyle \overrightarrow{OQ_2}=\frac{1-\xi}{2}\overrightarrow{OP_4}+\frac{1+\xi}{2}\overrightarrow{OP_3}$

$\displaystyle \overrightarrow{OP}=\frac{1-\eta}{2}\overrightarrow{OQ_1}+\frac{1+\eta}{2}\overrightarrow{OQ_2}$

$\displaystyle=\frac{1}{4}(1-\xi)(1-\eta)\overrightarrow{OP_1} + \frac{1}{4}(1+\xi)(1-\eta)\overrightarrow{OP_2} + \frac{1}{4}(1+\xi)(1+\eta)\overrightarrow{OP_3} + \frac{1}{4}(1-\xi)(1+\eta)\overrightarrow{OP_4}$  (2)


(2)式は(1a)(1b)式と全く同じ式です。すなわち自然座標系である$(\xi, \eta)$は$(-1\leqq \xi \leqq 1, \ -1\leqq \eta \leqq 1)$の範囲を取る訳ですから、$(\xi, \eta)$が動けば点Pも四角形の内部を動くことになります。$(\xi, \eta)$の値と物理座標系の関係は図のように辺の比を取って考えれば良い訳です。具体的に例えば $\xi = 0.5, \ \ \eta = -0.5$ のとき点Pがどの場所に居るのかもイメージ出来ますね?


2、四点が同一平面状にない場合の面積

「四角形」といえば数学的には普通は四角形の4つの頂点が同一平面上にあるものを指しています。 しかし実際の解析でモデルを四角形要素で分割するなどと言っている場合は、この四角形というのは必ずしも同一平面上にはない場合があります。 ここで考えてみたいのが、同一平面上にない四角形の面積はどのように求めれば良いか、という問題です。ぱっと思いつくのは四角形を2つの三角形に分ける方法でしょうか。しかしpic.3で言えば $\triangle P_1 P_2 P_3 + \triangle P_3 P_4 P_1$と考えた場合と、$\triangle P_2 P_3 P_4 + \triangle P_4 P_1 P_2$ と考えた場合では答えが違ってきます。
ここでは(2)式に基づいて面積を考えていきたいと思います。$(\xi, \eta)$が$(-1\leqq \xi \leqq 1, \ -1\leqq \eta \leqq 1)$の範囲を動くとき、Pの軌跡は四角形の4辺を含んだ曲面となっているはずです。ここで、この曲面の面積をベクトル解析に基づいて求めていくことにします。


OPベクトルは$\xi$ と $\eta$ の関数であるから、

$\overrightarrow{OP}=\boldsymbol{P}(\xi,\ \eta)$

と書いて、曲面の第1基本微分形式(この辺はWeb上で調べれば色々出てくるので省略)を考えれば、次式から面積が求められます。

$\displaystyle S=\int_{-1}^1 \int_{-1}^1\ \left|\frac{\partial \boldsymbol{P}}{\partial \xi}\times \frac{\partial \boldsymbol{P}}{\partial \eta}\right| d\xi\ d\eta$