morikomorou’s blog

自分が学んだことなどの備忘録的なやつ

トップ > Python > 【python】有限要素法(FEM)をゼロから実装する4_全体剛性マトリクスの導出

【python】有限要素法(FEM)をゼロから実装する4_全体剛性マトリクスの導出

Python 数値解析 有限要素法


はじめに

前回から本格的に4角形1次要素の平面応力状態における静弾性問題の有限要素解析の実装を行い始めました。今回はその続きで有限要素法の肝となる全体剛性マトリクスの導出まで行いたいと思います。




実装の流れ

前回のおさらいで実装の流れは下記のとおりです。

  1. Dマトリクス \boldsymbol{D}の算出
  2. Bマトリクス \boldsymbol{B}の算出
  3. 要素剛性マトリクス \boldsymbol{K}_{e}の算出
  4. 全体剛性マトリクス \boldsymbol{K}の算出
  5. 連立方程式を解き、全節点の変位 \boldsymbol{d}を求める

このうち、Bマトリクスの導出まで前回までで完了しました。

今回はその次の要素剛性マトリクスの導出からやっていきましょう

実装

以前の記事で作成したコードと、クラス要件をベースに実装を進めてます。

基本的な計算はすべてSolverクラスが実行するようなことが書かれてたので今回もそっちに実行するようにします。

まずは過去記事のおさらいからです。全コード載せときます。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import patches


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:
        id (int): メッシュのID(重複不可)
        nodes ({Node.id: Node}): Nodeの辞書
        elements (list): Elementのリスト
    """

    def __init__(self, idx, nodes_dict, elements):
        self.id = idx
        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)


class Material:
    """
    材料特性をまとめておくクラス
    
    Attributes:
        E (float): ヤング率
        nu (float): ポアソン比
    """

    def __init__(self, E, nu):
        self.E = E
        self.nu = nu


class Solver:
    """
    有限要素法の計算を行うクラス
    
    Attributes:
        mesh (Mesh): Meshクラス
        material (Material): Materialクラス
        D (array 3x3): Dマトリクス
    """

    def __init__(self, mesh, material):
        self.mesh = mesh
        self.material = material
        self.D = self.calc_D()

    def calc_D(self):
        """
        Dマトリクスを算出。

        Returns:
            D (array 3x3): Dマトリクス
        """
        D = np.array([[1, self.material.nu, 0],
                           [self.material.nu, 1, 0],
                           [0, 0, (1 - self.material.nu) / 2]]) \
                 * self.material.E / (1 - self.material.nu ** 2)
        return D
    
    def calc_B(self, elem, xi, eta):
        """
        Bマトリクス及びヤコビアンを算出。

        Args:
            elem (Element): 対象のElementクラス
            xi, eta (float): Bマトリクス算出用の座標

        Returns:
            B (array 3x8): Bマトリクス
            J (array 2x2): ヤコビアン
        """
        dndxi = np.array([- 1 + eta, 1 - eta, 1 + eta, - 1 - eta]) / 4
        dndeta = np.array([- 1 + xi, - 1 - xi, 1 + xi, 1 - xi]) / 4
        x = elem.xy[:, 0]
        y = elem.xy[:, 1]
        dxdxi = np.dot(dndxi, x)
        dydxi = np.dot(dndxi, y)
        dxdeta = np.dot(dndeta, x)
        dydeta = np.dot(dndeta, y)
        J = np.array([[dxdxi, dydxi], [dxdeta, dydeta]])
        B = np.zeros((3, 8))
        for i in range(4):
            Bi = np.dot(np.linalg.inv(J), np.array([dndxi[i], dndeta[i]]).T)
            B[0, 2 * i] = Bi[0]
            B[2, 2 * i + 1] = Bi[0]
            B[2, 2 * i] = Bi[1]
            B[1, 2 * i + 1] = Bi[1]
        return B, J


MESH_SIZE = 5.
LENGTH = 20.
HEIGHT = 10.
THICKNESS = 1.
E = 208000 # ヤング率
NU = 0.27 # ポアソン比

nodes_dict, elems = create_mesh(LENGTH, HEIGHT, MESH_SIZE, THICKNESS)
material = Material(E, NU)
mesh = Mesh(0, nodes_dict, elems, material)
solver = Solver(mesh, material)
B, J = solver.calc_B(solver.mesh.elements[0], -1, -1) # xi, etaはとりあえず-1で算出
print(B)
print(J)




要素剛性マトリクス \boldsymbol{K}_{e}の算出

要素剛性マトリクスの導出はそこまで難しくなくて、以下の式を要素ごとに計算するだけです。

{\displaystyle \boldsymbol{K}_{e} = h \sum_{i=1}^{2} \sum_{j=1}^{2} w_{i}w_{j}\boldsymbol{B}^{T}(\xi_{i}, \eta_{i})\boldsymbol{DB}(\xi_{i}, \eta_{i})\mathrm{det}\boldsymbol{J}d\xi d\eta \ }

そのうえでポイントとなるのは \boldsymbol{B}や \boldsymbol{J}算出における (\xi_{i}, \eta_{i})や w_{i}w_{j}をどう決めるか？というところになると思います。

上の式はガウスルジャンドル積分と呼ばれるもので積分を近似的に解く手法になります。
簡単に言うと、要素全域で積分をしたいけど計算が大変なので、ガウス点 (\xi_{i}, \eta_{i})と呼ばれる要素内の数点における値を計算して、それを重みづけ w_{i}w_{j}して足し合わせることで要素全域の積分近似ができるという計算方法です。

基本的にどの点をガウス点とするか？やその重みは決められています。
Wikiによると以下の通りです

ガウス求積 - Wikipedia

4角形1次のアイソパラメトリック要素の場合は以下のような点になります。

今回はよく使われるn=2のガウス点で実装しましょう。この時の重み w_{i}w_{j}はすべて1です。
ガウス点は以下の4点であり、その時の (\xi_{i}, \eta_{i})を \boldsymbol{K}_{e}の式に代入すればOKです。

(\xi_{i}, \eta_{i}) = (-1/\sqrt{3}, -1/\sqrt{3}), (1/\sqrt{3}, -1/\sqrt{3}), (1/\sqrt{3}, 1/\sqrt{3}), (-1/\sqrt{3}, 1/\sqrt{3})

class Solver:
    """
    有限要素法の計算を行うクラス
    
    Attributes:
        mesh (Mesh): Meshクラス
        material (Material): Materialクラス
        D (array 3x3): Dマトリクス
    """

    def __init__(self, mesh, material):
        self.mesh = mesh
        self.material = material
        self.D = self.calc_D()

    def calc_D(self):
        """
        Dマトリクスを算出。

        Returns:
            D (array 3x3): Dマトリクス
        """
        D = np.array([[1, self.material.nu, 0],
                           [self.material.nu, 1, 0],
                           [0, 0, (1 - self.material.nu) / 2]]) \
                 * self.material.E / (1 - self.material.nu ** 2)
        return D
    
    def calc_B(self, elem, xi, eta):
        """
        Bマトリクス及びヤコビアンを算出。

        Args:
            elem (Element): 対象のElementクラス
            xi, eta (float): Bマトリクス算出用の座標

        Returns:
            B (array 3x8): Bマトリクス
            J (array 2x2): ヤコビアン
        """
        dndxi = np.array([- 1 + eta, 1 - eta, 1 + eta, - 1 - eta]) / 4
        dndeta = np.array([- 1 + xi, - 1 - xi, 1 + xi, 1 - xi]) / 4
        x = elem.xy[:, 0]
        y = elem.xy[:, 1]
        dxdxi = np.dot(dndxi, x)
        dydxi = np.dot(dndxi, y)
        dxdeta = np.dot(dndeta, x)
        dydeta = np.dot(dndeta, y)
        J = np.array([[dxdxi, dydxi], [dxdeta, dydeta]])
        B = np.zeros((3, 8))
        for i in range(4):
            Bi = np.dot(np.linalg.inv(J), np.array([dndxi[i], dndeta[i]]).T)
            B[0, 2 * i] = Bi[0]
            B[2, 2 * i + 1] = Bi[0]
            B[2, 2 * i] = Bi[1]
            B[1, 2 * i + 1] = Bi[1]
        return B, J
    
    def calc_Ke(self, elem):
        """
        Keマトリクスを算出。

        Args:
            elem (Element): 対象のElementクラス

        Returns:
            Ke (array 8x8): 要素剛性マトリクス
        """
        gps = ((-1 / np.sqrt(3), -1 / np.sqrt(3), 1, 1), 
              (1 / np.sqrt(3), -1 / np.sqrt(3), 1, 1), 
              (1 / np.sqrt(3), 1 / np.sqrt(3), 1, 1), 
              (-1 / np.sqrt(3), 1 / np.sqrt(3), 1, 1), 
              )
        
        Ke = np.zeros((8, 8))
        for xi, eta, wi, wj in gps:
            B, J = self.calc_B(elem, xi, eta)
            Ke += wi * wj * np.dot(B.T, np.dot(self.D, B)) * np.linalg.det(J) * elem.thickness

        return Ke

結果は8x8の行列になります。




全体剛性マトリクス \boldsymbol{K}の算出

要素剛性マトリクスを求められたので次は全体剛性マトリクスを算出します。最終的に求めたい式は下記の形です。

\boldsymbol{f} = \boldsymbol{K} \, \boldsymbol{d}

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

以下に示す要素ごとのつり合いの式を全要素分まとめて連立させた形になっております。

\boldsymbol{f_{e}} = \boldsymbol{K}_{e} \, \boldsymbol{d_{e}}

簡単のために以下のようなメッシュを考えます。

ノード iの変位を (u_{i}, v_{i})とすると、要素1のつり合いの式は下記のとおりです。

要素2のつり合いの式は下記のとおり。

この2つの式を連立すると以下の形になり、全体剛性マトリクスは次のように計算できます。

実際にソルバークラスの中にcalc_Kという関数を作って全体剛性マトリクスを求める機能を実装してみます。

class Solver:
    """
    有限要素法の計算を行うクラス
    
    Attributes:
        mesh (Mesh): Meshクラス
        material (Material): Materialクラス
        D (array 3x3): Dマトリクス
    """

    def __init__(self, mesh, material):
        self.mesh = mesh
        self.material = material
        self.D = self.calc_D()

    def calc_D(self):
        """
        Dマトリクスを算出。

        Returns:
            D (array 3x3): Dマトリクス
        """
        D = np.array([[1, self.material.nu, 0],
                           [self.material.nu, 1, 0],
                           [0, 0, (1 - self.material.nu) / 2]]) \
                 * self.material.E / (1 - self.material.nu ** 2)
        return D
    
    def calc_B(self, elem, xi, eta):
        """
        Bマトリクス及びヤコビアンを算出。

        Args:
            elem (Element): 対象のElementクラス
            xi, eta (float): Bマトリクス算出用の座標

        Returns:
            B (array 3x8): Bマトリクス
            J (array 2x2): ヤコビアン
        """
        dndxi = np.array([- 1 + eta, 1 - eta, 1 + eta, - 1 - eta]) / 4
        dndeta = np.array([- 1 + xi, - 1 - xi, 1 + xi, 1 - xi]) / 4
        x = elem.xy[:, 0]
        y = elem.xy[:, 1]
        dxdxi = np.dot(dndxi, x)
        dydxi = np.dot(dndxi, y)
        dxdeta = np.dot(dndeta, x)
        dydeta = np.dot(dndeta, y)
        J = np.array([[dxdxi, dydxi], [dxdeta, dydeta]])
        B = np.zeros((3, 8))
        for i in range(4):
            Bi = np.dot(np.linalg.inv(J), np.array([dndxi[i], dndeta[i]]).T)
            B[0, 2 * i] = Bi[0]
            B[2, 2 * i + 1] = Bi[0]
            B[2, 2 * i] = Bi[1]
            B[1, 2 * i + 1] = Bi[1]
        return B, J
    
    def calc_Ke(self, elem):
        """
        Keマトリクスを算出。

        Args:
            elem (Element): 対象のElementクラス

        Returns:
            Ke (array 8x8): 要素剛性マトリクス
        """
        gps = ((-1 / np.sqrt(3), -1 / np.sqrt(3), 1, 1), 
              (1 / np.sqrt(3), -1 / np.sqrt(3), 1, 1), 
              (1 / np.sqrt(3), 1 / np.sqrt(3), 1, 1), 
              (-1 / np.sqrt(3), 1 / np.sqrt(3), 1, 1), 
              )
        
        Ke = np.zeros((8, 8))
        for xi, eta, wi, wj in gps:
            B, J = self.calc_B(elem, xi, eta)
            Ke += wi * wj * np.dot(B.T, np.dot(self.D, B)) * np.linalg.det(J) * elem.thickness

        return Ke

    def calc_K(self):
        """
        Kマトリクスを算出。

        Returns:
            K (array 2node_numx2node_num): 全体剛性マトリクス
        """
        K = np.zeros((len(self.mesh.nodes) * 2, len(self.mesh.nodes) * 2))
        for elem in self.mesh.elements:
            Ke = self.calc_Ke(elem)
            for i in range(4):
                for j in range(4):
                    K[2 * elem.node[i]:2 * elem.node[i]+2,
                      2 * elem.node[j]:2 * elem.node[j]+2] += Ke[2 * i:2 * i+2, 2 * j:2 * j+2]
        return K

とりあえずエラーは吐いてないので大丈夫そうです。




終わりに

次回はついに最後のステップである連立方程式を解いて変形を求めるところまでやっていきましょう。

ランキング参加中
プログラミング

参考文献