はじめに

本記事では、FEM解析の概要を説明し、三角形・四角形1次要素の平面歪要素の定式化について説明します。参考にしたページのリンクを貼っておきます。

本記事では、読者の方の工学的判断の役に立つように、FEM解析においてどこにどのような仮定が含まれているのか、コンピューターによる計算でどのような工夫がなされているのかを中心に解説していきます。

FEM解析の原理

　有限要素法では、対象の物体を有限個の要素に分割し、要素の変位を、1次式、2次式などで仮定し、物体に生じる変位、歪み、応力を算出します。本記事では、変位を各節点座標の1次式で仮定(＝要素内で歪み一定で仮定)した、1次要素について記述します。
　また、有限要素法で扱う要素には、ライン要素、平面応力要素、シェル要素、3次元要素など、扱う次元や、節点の自由度に応じて様々な要素がありますが、今回は面外方向の応力を考慮しない並進2成分のみを節点自由度とする平面応力要素について扱います。

四角形1次要素

　前項でも述べたように、四角形1次要素は「要素内の変位が線形に変化する」と仮定して計算を行います。しかし、実際には荷重下での変位は非線形であるため、一定の誤差が生じます。そのため、精度はあまり高くありません。これを解決するために、「2次要素」や「1次適合要素」などがなどがあります。こちらについては後日記事化します。

誤差を小さくするためには、変位（あるいは応力）の勾配が急な領域で、線形近似が成立するくらい要素を細分化する必要があります。「要素内で定ひずみ」とする仮定は実際にはかなり無理があって、多くの商用FEMソフトではこの仮定を採用していません。

定式化

　ここからは、定式化のキーとなる式を紹介します。具体的な展開は、冒頭で紹介した参考リンクをご覧ください。この記事では、コンピューターに計算させるために、どこにどんな仮定が含まれているのかという視点で解説していきます。

FEMの根幹をなす式は、主に以下の3つです：

荷重-変位関係
$$
\mathbf{F} = \mathbf{K} \mathbf{U}
$$

変位-歪み関係
$$
\mathbf{\epsilon} = \mathbf{B} \mathbf{U}
$$
歪み-応力関係
$$
\mathbf{\sigma} = \mathbf{D} \mathbf{\epsilon}
$$

ここで、

$\mathbf{K} :剛性マトリクス$
$\mathbf{B} :変位ひずみマトリクス$
$\mathbf{D}:歪み応力マトリクス$
がわかれば、$\mathbf{U}$がもとまり、$\mathbf{\epsilon}$がまとまり、$\mathbf{\sigma}$がまとまるわけです。

それぞれのマトリクスについてみていきましょう。

$\mathbf{D}$マトリクスについて
$\mathbf{D}$は、要素の材料特性からまとまり、ヤング率$\mathbf{E}$と、ポアソン比$\mathbf{\nu}$を用いて、下記で表されます。
$$
\mathbf{D} =
\frac{E}{(1 - \nu^2)}
\begin{bmatrix}
1 - \nu & \nu & 0
\\ \nu & 1 - \nu & 0
\\ 0 & 0 & \frac{1 - \nu}{2}
\end{bmatrix}
$$

$$
$\mathbf{B}$マトリクスについて
$\mathbf{B}$は、変位歪みマトリクスであり、変位の偏微分として、
$$
\mathbf{B} =
\begin{bmatrix}
\frac{\partial u}{\partial x}
\\ \frac{\partial u}{\partial y}
\\ \frac{\partial u}{\partial y} + \frac{\partial u}{\partial x}
\end{bmatrix}
$$
で表されます。

1次要素の場合要素内の変位を1次式で仮定するので、要素内の任意の点における変位は以下で表されます。
$$
u(x,y) = \alpha_0 + \alpha_1x+\alpha_2y +\alpha_3xy
$$
$$
v(x,y) = \beta_0 + \beta_1x+\beta_2y +\beta_3xy
$$
各節点において、上式が成り立つので、
$$
u_i(x,y)=\alpha_0 + \alpha_1x_i+\alpha_2y_i +\alpha_3x_iy_i
$$
$$
v_i(x,y) = \beta_0 + \beta_1x_i+\beta_2y_i +\beta_3x_iy_i
$$
これを、$\alpha$について解くと、要素内の任意の点での変位が表現できます。

ただ、この式を解くとかなり複雑になってしまうのと、この後出てくる積分の処理が煩雑になってしまうので、別の方法を考えます。

$\xi$,$\eta$ の正規化座標系を導入して、座標変換を行い、変位$u,v$を以下のように表します。

$$
u(\xi,\eta) = \sum_{i=1}^4 N_i(\xi,\eta)u_i
$$
$$
v(\xi,\eta) = \sum_{i=1}^4 N_i(\xi,\eta)v_i
$$
ここで、$N_i(\xi,\eta)$は形状関数と呼ばれ、
$$
N_1 = \frac14(1-\xi)(1-\eta)
$$
$$
N_2 = \frac14(1+\xi)(1-\eta)
$$
$$
N_3 = \frac14(1+\xi)(1+\eta)
$$
$$
N_4 = \frac14(1-\xi)(1+\eta)
$$

で表されます。要素内変位を各節点の重みで表現した形になっています。

$B$マトリクスを求めるのに必要なのは、$u,v$の$x,y$微分です。

$$
\frac{\partial u}{\partial x}
= \sum_{i=1}^4 \frac{\partial N_i}{\partial x}
$$

$$
\frac{\partial v}{\partial x}
= \sum_{i=1}^4 \frac{\partial N_i}{\partial y}
$$

となるため、形状関数の微分を求める必要があり、以下の式で求める。

$$
\begin{Bmatrix}
\frac{\partial N_i}{\partial x}
\\ \frac{\partial N_i}{\partial y}
\end{Bmatrix} =
\begin{bmatrix}
\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x}
\\ \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y}
\end{bmatrix}
\begin{Bmatrix}
\frac{\partial N_i}{\partial \xi}
\\ \frac{\partial N_i}{\partial \eta}
\end{Bmatrix}
=[\mathbf{J}]^{-1}
\begin{Bmatrix}
\frac{\partial N_i}{\partial \xi}
\\ \frac{\partial N_i}{\partial \eta}
\end{Bmatrix}
$$

これらより$\mathbf{B}$がまとまる。
この時、$\mathbf{J}$はヤコビ行列で、下記のように定義される。このヤコビ行列は、座標変換に用いられ、この後もよく出てきます。

$$
J =
\begin{bmatrix}
\frac{\partial x}{\partial\xi} & \frac{\partial y}{\partial\xi}\
\\ \frac{\partial x}{\partial\eta} & \frac{\partial y}{\partial\eta}\
\end{bmatrix}
$$

ここまでの記述から$\mathbf{B}$は下記のように表されます。

$$
\mathbf{B} =
\begin{bmatrix}
\frac{\partial N_1}{\partial x} & 0 &\cdots& \frac{\partial N_4}{\partial x} & 0
\\ 0 & \frac{\partial N_1}{\partial y} &\cdots& 0 & \frac{\partial N_4}{\partial y}
\\ \frac{\partial N_1}{\partial y} & \frac{\partial N_1}{\partial x} &\cdots& \frac{\partial N_4}{\partial y} & \frac{\partial N_4}{\partial x}
\end{bmatrix}
$$

$K$マトリクスについて
剛性マトリクスの詳細な式展開は、参考のページをご参照ください。仮想仕事の原理などを使い、下記の式で求めることができます。

$$
\mathbf{K} = \int_v \mathbf{B}^\mathsf{T}\mathbf{D}\mathbf{B}dV=
t\int_v \mathbf{B}^\mathsf{T}\mathbf{D}\mathbf{B}dxdy
\hspace{10pt}
t:要素の厚さ
$$
ここで、積分計算が出てきますが、xy座標系で積分するのが手間なので、下式のように、正規化した$\xi$,$\eta$座標系に変換することにウマミがある訳です。

$$
\mathbf{K}=
t\int_v \mathbf{B}^\mathsf{T}\mathbf{D}\mathbf{B}dxdy=
t \int_{-1}^{1} \int_{-1}^{1} \mathbf{B}(\xi,\eta)^\mathsf{T}\mathbf{D}\mathbf{B}(\xi,\eta)det|\mathbf{J}|d\xi d\eta
$$

これで、$K$が求められる訳ですが、この積分も簡単ではありません。ここで、コンピューターで積分するために、近似的に積分を解く数値積分（ガウス積分）という先人の知恵を使います。関数の積分を特定の点での関数値の重み付きの総和として近似する方法です。
詳しい内容は初めに挙げた参照リンクの記事をご覧ください。

四角形1次要素の場合積分点と、その重み$w_i$は以下4点で表されます。

$$
\xi_i,\eta_i =
\pm\frac{1}{\sqrt{3}},w_i=1
$$

よって$K$マトリクスは、以下のようになります。
$$
\mathbf{K}=
t \sum_{i=1}^{2}\sum_{j=1}^{2} w_iw_j\mathbf{B}(\xi_i,\eta_i)^\mathsf{T}\mathbf{D}\mathbf{B}(\xi_i,\eta_i)det|\mathbf{J(\xi_i,\eta_i)}|
$$

これで、コンピューターで計算可能な形に書き換えられました。
先人の知恵は偉大ですね。

$K$マトリクスが求まったので、これで各節点での変位$U$が求められます。
変位$Ｕ$と、$Ｂ$,$D$から、$\epsilon$,$\sigma$が求められます。

ここで$Ｂ$マトリクスは各積分点で算出するため、$\epsilon=BU$から求まるひずみは、積分点での応力である点に注意が必要です。設計では、通常節点位置での応力を使用するので、４つの積分点での応力から節点位置での応力を外挿（線形補完）します。
また、各節点の応力は要素ごとに算出されるため、同じ節点でも要素ごとに異なる値となります。通常、応力を滑らかに表示するための処理として、節点での各応力を平均してスムージングする処理が行われているので、この点も注意が必要です。

まとめ

本記事では、平面応力要素の四角形１次要素の定式化を行いました。本記事のまとめとして、商用ソフトを扱う上で、設計者でも知っておくべきだと個人的に思う点について示します。

同じ平面応力要素でも、要素内の変位をどのように仮定するかで何通りか方法がある。
積分の処理は数値積分を使用していて、近似解を求めている。
応力は積分点で求まり、節点位置での応力は積分点の応力を外挿して求めている。
節点位置での応力は要素ごとに異なり、平均値を求めてスムージングされてい場合が多い。

次回は、実装編です。Grasshopperの環境で、平面応力要素の四角形１次要素の実装する過程を記事にします。