【python】有限要素法(FEM)をゼロから実装する1_概要とメッシュ生成

はじめに

有限要素法(FEM)による解析は各種設計を行う際には必須です。

市販のFEMの解析ソフトで変形とか、振動とか特に知識がなくても簡単に解析できる時代ですが、やはりFEM解析を行う上では、どういう原理なのかを知っておきたいですよね。

なので、有限要素法を用いた2次元弾性体の静的変形解析のプログラムをゼロから作っていきたいと思います。長いシリーズになりそうなので今回はメッシュの生成まで行います。

はじめに
有限要素法(Finite Element Method)とは
- プログラムの流れ
- FEMプログラムにおけるクラス要件
メッシュの生成
終わりに
参考文献

有限要素法(Finite Element Method)とは

偏微分方程式の近似解を数値的に解く方法の一つです。解析対象モデルを小さなの要素(element)に分割してそれらの要素ごとで計算した結果を重ね合わせて全体の解を求めます。

要素は節点(node)に囲まれた領域で、弾性体の場合、それぞれの要素は以下の式に従います。
$\boldsymbol{K}_{e}\boldsymbol{u}_{e} = \boldsymbol{F}_{e}$

ここで $\boldsymbol{K}_{e}$ は要素剛性マトリクス、 $\boldsymbol{u}_{e}$ は要素の節点変位ベクトル、 $\boldsymbol{F}_{e}$ は等価節点力と呼ばれます。

これを全要素分足し合わせてまとめると以下の連立方程式の形になります。
$\boldsymbol{K}\boldsymbol{u} = \boldsymbol{F}$

ここで $\boldsymbol{K}$ は全体剛性マトリクス、 $\boldsymbol{u}$ は全節点の変位をまとめた全体変位ベクトル、 $\boldsymbol{F}$ は全節点の等価節点力をまとめた全体荷重ベクトルです。

この式に各節点変位の境界条件や、荷重の境界条件を代入して $\boldsymbol{u}$ に対して連立方程式を解くことで各節点の変位を求めることができます。

変位が求まればひずみや応力も簡単に求まります。

プログラムの流れ

FEM解析プログラムの流れとしては以下のような流れになります。

モデル作成
モデルを細分化し要素作成
要素に物性値や境界条件等を指定
要素ごとに要素剛性マトリクス $\boldsymbol{K}_{e}$ を求める
全要素の要素剛性マトリクス $\boldsymbol{K}_{e}$ をまとめて全体剛性マトリクス $\boldsymbol{K}$ を作成
境界条件から全体変位ベクトル $\boldsymbol{u}$ 、全体荷重ベクトル $\boldsymbol{F}$ を作成
連立方程式を解き、全体変位ベクトル $\boldsymbol{u}$ を求める
各要素の応力、ひずみを計算する
結果を出力

FEMプログラムにおけるクラス要件

せっかくなのでオブジェクト指向でコード作成したいと思います。クラスの要件をChatGPTさんに作ってもらいました。

これベースで開発していきたいと思います。

メッシュの生成

今回はモデルを作成し、要素に細分化してみるところから始めます。

要素の形状

計算を行う上で要素の形状は同じにしておく必要があります。2次元解析においてよく使われるのは3角形要素と4角形要素です。
今回は4角形要素を使用します。

Nodeクラス、Elementクラス、Meshクラスの作成

まずはクラスを定義しておきます。

各クラスには一意のidを付与して置き、識別できるようにしておきます。

Nodeクラスには節点のx, y座標の情報とそのノードのIDを定義します。
Elementクラスには要素を構成する4つのノードIDを入れられるリストと、その要素のIDを定義します。2次元要素なので厚みも定義します。
Elementの各節点の座標も必要になる気がするのでノードIDのリストから座標を取得できるメソッドも定義しておきます。
MeshクラスにはそのメッシュのID及び、一つのモデルから作成した全Nodeクラスの辞書と、全Elementクラスのリストを定義します。
Meshクラス生成時にElementの各節点座標を取得できるようにしておきます。

import numpy as np


class Node:
    """
    Nodeのインデックス、座標および境界条件をまとめておくクラス(2次元用)
    
    Attributes:
        id (int): Nodeのインデックス番号(重複不可)
        x, y (float): Nodeのx, y座標
    """
    def __init__(self, idx, x, y):
        self.id = idx
        self.x = x
        self.y = y
    
    def __str__(self):
        return 'Node {}: x {}, y {}'.format(self.id, self.x, self.y)


class Element:
    """
    4角形1次要素Elementのインデックス、要素内のNode情報および要素の厚みをまとめておくクラス
    
    Attributes:
        id (int): Elementのインデックス番号(重複不可)
        node ([int, int, int, int]): 要素内のNodeのインデックスのリスト(左下から反時計回りで指定)
        thickness (float): Elementの厚み
        xy (ndarray(size: 2, 4)): Element内のNodeのxy座標まとめ
    """
    def __init__(self, idx, node_idxs, thickness):
        self.id = idx
        self.node = node_idxs
        self.thickness = thickness
        self.xy = None

    def __str__(self):
        return 'Element {}: Nodes {}, Thickness: {}'.format(self.id, self.node, self.thickness)
    
    def get_coordination(self, nodes_dict):
        """
        要素の各節点のx,y座標を求めてself.xyに代入。

        Args:
            nodes_dict ({Node.id: Node}): Nodeの辞書
        """

        res = []
        for node_idx in self.node:
            res.append([nodes_dict[node_idx].x, nodes_dict[node_idx].y])
        self.xy = np.array(res)


class Mesh:
    """
    NodeオブジェクトとElementオブジェクトをまとめておくクラス
    
    Attributes:
        nodes ({Node.id: Node}): Nodeの辞書
        elements (list): Elementのリスト
    """

    def __init__(self, nodes_dict, elements):
        self.nodes = nodes_dict
        self.elements = elements
        self.get_element_coord()
    
    def get_element_coord(self):
        for elm in self.elements:
            elm.get_coordination(self.nodes)

モデル作成

モデルの形状と細分化は下記の図のようにします。

ノードの作成

ノードを作成してみましょう。モデルやメッシュのサイズは可変できるようにしておきます。

mesh_size = 5.
length = 10.
height = 10.

# Nodeの作成
num_mesh_len = int(length / mesh_size)
num_mesh_hei = int(height / mesh_size)
if num_mesh_len == 0:
    num_mesh_len = 1
if num_mesh_hei == 0:
    num_mesh_hei = 1
x = np.linspace(0, length, num_mesh_len + 1)
y = np.linspace(0, height, num_mesh_hei + 1)
X, Y = np.meshgrid(x, y)
X = X.ravel()
Y = Y.ravel()

nodes_dict = {}
for i, coord in enumerate(zip(X, Y)):
    nodes_dict[i] = Node(i, coord[0], coord[1])

for node in nodes_dict.values():
    print(node)

結果は下記のとおりです。ちゃんと作成できてます。

Node 0: x 0.0, y 0.0
Node 1: x 5.0, y 0.0
Node 2: x 10.0, y 0.0
Node 3: x 0.0, y 5.0
Node 4: x 5.0, y 5.0
Node 5: x 10.0, y 5.0
Node 6: x 0.0, y 10.0
Node 7: x 5.0, y 10.0
Node 8: x 10.0, y 10.0

4角形要素の作成

先ほどノードを作成しましたが、それらの隣り合う4つからそれぞれ1つの要素を作成します。
Elementオブジェクトを作る際に、構成するノードIDのリストの指定が必要ですが、このリストを作る際、ノードの順番も今後の計算に重要になるので注意です。

以下の図を例に示しますが、左下から反時計回りにノードを指定するようにします。この場合だと構成するノードIDのリストは[1,2,4,3]と指定します。

thickness = 1.

# Elementの作成
node_idx = 0
elem_idx = 0
elems = []
for i in range(num_mesh_hei):
    for j in range(num_mesh_len + 1):
        if (node_idx + 1) % (num_mesh_len + 1) == 0:
            node_idx += 1
            continue
        else:
            node_idxs = [node_idx, node_idx + 1,
                        node_idx + num_mesh_len + 2, node_idx + num_mesh_len + 1]
            elems.append(Element(elem_idx, node_idxs, thickness))
            node_idx += 1
            elem_idx += 1

for elem in elements:
    print(elem)

結果は以下。

Element 0: Nodes [0, 1, 4, 3], Thickness: 1.0
Element 1: Nodes [1, 2, 5, 4], Thickness: 1.0
Element 2: Nodes [3, 4, 7, 6], Thickness: 1.0
Element 3: Nodes [4, 5, 8, 7], Thickness: 1.0

可視化

コマンドラインの出力だけだとわかりづらいので可視化しておきます。ついでに先ほどのメッシュ生成も関数にまとめておきます。

import matplotlib.pyplot as plt
from matplotlib import patches


def create_mesh(length, height, mesh_size, thickness):
    """
    モデルを格子点に分割し、全ての点のNodeクラスおよびそれらを使用した4角形1次要素を作成する.
    ただし、モデル形状は四角形とする。
    
    Args:
        length, height (float): モデルの長さ、高さ
        mesh_size (float): 分割したいメッシュのサイズ
        thickness (float): 要素の厚み(今回は全て一定とする)
        
    Returns:
        nodes_dict ({Node.id: Node}): 作成したNodeの辞書
        elsms (list): 作成したElementクラスをまとめたリスト
    """
    # Nodeの作成
    num_mesh_len = int(length / mesh_size)
    num_mesh_hei = int(height / mesh_size)
    if num_mesh_len == 0:
        num_mesh_len = 1
    if num_mesh_hei == 0:
        num_mesh_hei = 1
    x = np.linspace(0, length, num_mesh_len + 1)
    y = np.linspace(0, height, num_mesh_hei + 1)
    X, Y = np.meshgrid(x, y)
    X = X.ravel()
    Y = Y.ravel()

    nodes_dict = {}
    for i, coord in enumerate(zip(X, Y)):
        nodes_dict[i] = Node(i, coord[0], coord[1])
    
    # Elementの作成
    node_idx = 0
    elem_idx = 0
    elems = []
    for i in range(num_mesh_hei):
        for j in range(num_mesh_len + 1):
            if (node_idx + 1) % (num_mesh_len + 1) == 0:
                node_idx += 1
                continue
            else:
                node_idxs = [node_idx, node_idx + 1,
                         node_idx + num_mesh_len + 2, node_idx + num_mesh_len + 1]
                elems.append(Element(elem_idx, node_idxs, thickness))
                node_idx += 1
                elem_idx += 1
    return nodes_dict, elems


def plot_mesh(mesh):
    fig, ax = plt.subplots()
    for elem in mesh.elements:
        patch = patches.Polygon(xy=elem.xy, ec='black')
        ax.add_patch(patch)
        text_xy = np.mean(elem.xy, axis=0)
        ax.text(text_xy[0], text_xy[1], elem.id, fontsize=12, va='center', ha='center')
    for node in mesh.nodes.values():
        ax.scatter(node.x, node.y, fc='black', s=100)
        ax.text(node.x, node.y, node.id, fontsize=8, color='white', va='center', ha='center')
    ax.autoscale()
    ax.set_aspect('equal', 'box')
    plt.show()


# モデル、要素サイズ定義
MESH_SIZE = 5.
LENGTH = 10.
HEIGHT = 10.
THICKNESS = 1.

nodes_dict, elems = create_mesh(LENGTH, HEIGHT, MESH_SIZE, THICKNESS)
mesh = Mesh(0, nodes_dict, elems)
for node in mesh.nodes.values():
    print(node)
for elem in mesh.elements:
    print(elem)

plot_mesh(mesh)

結果は以下。

Node 0: x 0.0, y 0.0
Node 1: x 5.0, y 0.0
Node 2: x 10.0, y 0.0
Node 3: x 0.0, y 5.0
Node 4: x 5.0, y 5.0
Node 5: x 10.0, y 5.0
Node 6: x 0.0, y 10.0
Node 7: x 5.0, y 10.0
Node 8: x 10.0, y 10.0
Element 0: Nodes [0, 1, 4, 3], Thickness: 1.0
Element 1: Nodes [1, 2, 5, 4], Thickness: 1.0
Element 2: Nodes [3, 4, 7, 6], Thickness: 1.0
Element 3: Nodes [4, 5, 8, 7], Thickness: 1.0