数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。noteにも書いています。https://note.mu/kumag_min

【有限要素法】2次元有限要素法における面積座標(局所座標)の積分公式

今回は、2次元の有限要素法における面積座標(局所座標)の積分公式


 \displaystyle \iint_T L_1^l L_2^m L_3^n dxdy = 2S\frac{l! m! n!}{(l+m+n+2)!}

を示します。ここで  T は2次元の三角形要素、 S はその面積、 L_1 L_2 L_3 は形状関数で、それぞれ

 
\displaystyle 
\begin{eqnarray} 

L_1 &=& A_1+B_1x+C_1y \\
L_2 &=& A_2+B_2x+C_2y \\
L_3 &=& A_3+B_3x+C_3y \\

\end{eqnarray}

とあらわされます。係数は

 
\displaystyle 
\begin{eqnarray} 

&A_1& = \frac{x_2y_3-x_3y_2}{2S}, &B_1& = \frac{y_2-y_3}{2S}, &C_1& = \frac{x_3-x_2}{2S} \\
&A_2& = \frac{x_3y_1-x_1y_3}{2S}, &B_2& = \frac{y_3-y_1}{2S}, &C_2& = \frac{x_1-x_3}{2S} \\
&A_3& = \frac{x_1y_2-x_2y_1}{2S}, &B_3& = \frac{y_1-y_2}{2S}, &C_3& = \frac{x_2-x_1}{2S} \\

\end{eqnarray}

で与えられます。 x_i y_j はそれぞれ三角形要素の頂点の座標です。 これを知っていると積分がサクサクできます。というよりも2次元の場合この公式無しに積分を求めるのは無理です。


まず、 X=L_1 Y=L_2 とおきます。すると、 L_1+L_2+L_3=1 より  L_3=1-L_1-L_2=1-X-Y とかけるので上式に代入して


 \displaystyle \iint_T L_1^l L_2^m L_3^n dxdy = \iint_T X^l Y^m (1-X-Y)^n dxdy

を得ます。次に  x y に関する積分 X Y に関する積分に変更します。 (x,y) から  (X,Y) への変数変換なのでまずヤコビアンを計算します。すると

 \displaystyle 
\left| \begin{array}{cc}
      \frac{\partial x}{\partial X} & \frac{\partial x}{\partial Y} \\
      \frac{\partial y}{\partial X} & \frac{\partial y}{\partial Y} \\
    \end{array}
  \right|
=
\left| \begin{array}{cc}
      x_1-x_3 & x_2-x_3 \\
      y_1-y_3 & y_2-y_3 \\
    \end{array}
  \right|
=
2S

となることがわかります。 L_1=X は頂点1でのみ1をとりその他の頂点で0をとります。他の形状関数も同様です。なので三角形領域  T は原点に直角部分がくる  1:1:\sqrt{2} の直角三角形の領域になります(頂点1は  X 軸上長さ1の点に、頂点2は  Y 軸上長さ1の点に、頂点3は原点へと変換されます)。よって

 \displaystyle 
\begin{eqnarray} 

\iint_T X^l Y^m (1-X-Y)^n dxdy &=& \iint_T X^l Y^m (1-X-Y)^n 2S dXdY \\
                                                    &=& 2S \int_0^1 X^l \int_0^{1-X} Y^m (1-X-Y)^n dY dX

\end{eqnarray}

となります。さらに変数変換

 \displaystyle \tilde{Y} = \frac{Y}{1-X}

を用いて

 \displaystyle \int_0^{1-X} Y^m (1-X-Y)^n dY = (1-X)^{m+n+1} \int_0^1 \tilde{Y}^m (1-\tilde{Y})^n d\tilde{Y}

を得ます。ここで、こちらの記事

で紹介している公式

 \displaystyle \int_0^a x^m (a-x)^n dx = \frac{m! n!}{(m+n+1)!} a^{m+n+1}

を使うと

 \displaystyle \int_0^1 \tilde{Y}^m (1-\tilde{Y})^n d\tilde{Y} = \frac{m! n!}{(m+n+1)!}

を得ます。これらの式を用いると

 \displaystyle 
\begin{eqnarray}

&2S& \int_0^1 X^l \int_0^{1-X} Y^m (1-X-Y)^n dY dX \\
&=& 2S \frac{m! n!}{(m+n+1)!} \int_0^1 X^l (1-X)^{m+n+1}dX

\end{eqnarray}

を得ます。もう一度上の公式を  X積分に対して適用して

 \displaystyle  
\begin{eqnarray}

&2S& \frac{m! n!}{(m+n+1)!} \int_0^1 X^l (1-X)^{m+n+1}dX \\
&=& 2S \frac{m! n!}{(m+n+1)!} \frac{l! (m+n+1)!}{(l+m+n+2)!} \\
&=& 2S \frac{m! n! l!}{(l+m+n+2)!}

\end{eqnarray}

となり目的の式を得ます。


1次元有限要素法における面積座標(局所座標)の積分公式はこちらです。


参考
菊地文雄『有限要素法概説』のp.56とpp.148-150

有限要素法概説―理工学における基礎と応用 (FEM+BEM (3))

有限要素法概説―理工学における基礎と応用 (FEM+BEM (3))

中山司『流れ解析のための有限要素法入門』のpp.47-49とpp.73-85

流れ解析のための有限要素法入門

流れ解析のための有限要素法入門