要素剛性マトリックス［K］ | 有限要素法の数値計算第6回

今回は、「素剛性マトリックス［K］」について説明します。

本コラムでは、第1回の「有限要素法の歴史」から始まり、CAEについて、有限要素法（FEM）を中心に解説してきました。第2回の「CAEの計算の流れ」では、CAEの基本的な計算ステップを紹介し、第3回から第5回にかけて、「形状関数」、「変位-ひずみマトリックス［B］」、「応力-ひずみマトリックス［D］」について詳しく説明しました。今回はそれらを統合し、「要素剛性マトリックス［K］」について説明します。そして、この回を最終回として、コラムの総まとめを行います。

これまでのおさらい

第1回：有限要素法の歴史
第2回：CAEの計算の流れ
第3回：形状関数
第4回：変位－ひずみマトリックス［B］
第5回：応力－ひずみマトリックス［D］

今回はシリーズの最終回として、これらの項目がどのように連携し、有限要素解析の「心臓部」である剛性マトリックス［K］の構築につながるのかをまとめます。

形状関数による変位の補間

形状関数は、要素内部の任意点の変位を節点の変位から補間するための関数です。四角形一次要素では、自然座標系を用いて変位を \( \xi,\ \eta \) の関数として表現し、各節点の変位から式を構築します。

変位-ひずみマトリックス［B］の導出

ひずみは変位の空間微分で表され、\( \varepsilon_x = \dfrac{\partial u}{\partial x} \)、\( \varepsilon_y = \dfrac{\partial v}{\partial y} \)、\( \gamma_{xy} = \dfrac{\partial u}{\partial y} + \dfrac{\partial v}{\partial x} \) のように定義されます。形状関数を微分することで、変位-ひずみマトリックス［B］が構築され、節点変位から要素内のひずみを求めることができます。

応力-ひずみマトリックス［D］の導出

材料の力学特性（ヤング率E、ポアソン比ν）を用いて、ひずみから応力を求める関係式を構築します。一般化されたフックの法則に基づき、応力-ひずみマトリックス［D］が導出され、\( \sigma = D\varepsilon \) の形で表現されます。

これらが組み合わさることで、要素剛性マトリックス［K］が構築されます。

要素剛性マトリックス［K］の構築

形状関数から［B］を、材料特性から［D］を導出し、要素剛性マトリックス［K］は以下の積分式で構築されます。

\[ [K] = \int_V [B]^{T}\,[D]\,[B]\,dV \qquad \cdots ① \]

［K］：要素剛性マトリックス
［B］：変位－ひずみマトリックス
［D］：応力－ひずみマトリックス
V ：要素体積

この積分は通常、数値積分（ガウス積分など）によって計算されます。要素ごとにつり合い方程式を作成します。

構造全体の挙動を解析するためには、全要素の剛性マトリックスを組み合わせて、構造全体のグローバル剛性マトリックスを構築する必要があります。

\[ [K]\{U\} = \{F\} \qquad \cdots ② \]

[K]：全体剛性マトリックス
{U}：全体節点変位ベクトル
{F}：外力ベクトル

境界条件を割り当てた後に、この全体剛性マトリックスの全体を表す連立方程式を解きます。

\[ \{U\} = [K]^{-1}\{F\} \qquad \cdots ③ \]

変位ベクトル{U}が得られれば、［B］マトリックスにより要素のひずみが求まります。（※1）さらに［D］マトリックスによって要素の応力が求められます。（※2）このように有限要素法において、変位、応力、ひずみが計算されます。

※1：第4回：変位－ひずみマトリックス［B］　式⑥参照

\[ \begin{Bmatrix} \varepsilon_x \\ \varepsilon_y \\ \gamma_{xy} \end{Bmatrix} = \begin{bmatrix} \dfrac{\partial N_1}{\partial x} & 0 & \dfrac{\partial N_2}{\partial x} & 0 & \dfrac{\partial N_3}{\partial x} & 0 & \dfrac{\partial N_4}{\partial x} & 0 \\[6pt] 0 & \dfrac{\partial N_1}{\partial y} & 0 & \dfrac{\partial N_2}{\partial y} & 0 & \dfrac{\partial N_3}{\partial y} & 0 & \dfrac{\partial N_4}{\partial y} \\[6pt] \dfrac{\partial N_1}{\partial y} & \dfrac{\partial N_1}{\partial x} & \dfrac{\partial N_2}{\partial y} & \dfrac{\partial N_2}{\partial x} & \dfrac{\partial N_3}{\partial y} & \dfrac{\partial N_3}{\partial x} & \dfrac{\partial N_4}{\partial y} & \dfrac{\partial N_4}{\partial x} \end{bmatrix} \begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\ v_4 \end{Bmatrix} \qquad \cdots ⑥ \]

※2：第5回：応力－ひずみマトリックス［D］　式⑨参照

\[ \{\sigma\} = [D]\{\varepsilon\} \qquad \cdots ⑨ \]

上記のように第2回「CAEの計算の流れ」で紹介したCAEの計算ステップ（弱形式化、離散化、補間、方程式の構築、解法）は、今回の説明で理論的に裏付けられました。

形状関数、［B］、［D］を通じて、CAEソフトウェアがどのようにして物体の挙動を数値的に予測しているかが明らかになりました。

まとめ

本シリーズを通じて、有限要素法の数値計算の仕組みを理解する一助となれば幸いです。今回が「有限要素法解析の数値計算」最終回となります。長い間ご精読ありがとうございました。