Python×2次元梁の要素剛性マトリクスの作成

投稿者: | 2021年6月15日

有限要素法構造解析プログラムを作成するにあたって、マトリクス(行列)の記述は骨の折れる作業です。たとえば要素剛性マトリクスであれば、2次元梁なら6×63次元梁なら12×12もの行列になりますが、もちろん1か所のタイプミスでも答えが合わなくなるので、エラーの発見が大変です。

そこで自分用の保管を兼ねて、Pythonで2次元梁モデルの要素剛性マトリクスを作成するプログラムを、インターネットの海に投下しておきます。

2次元梁の要素剛性マトリクス

要素剛性マトリクスについて大雑把に解説すると、ある1つの要素内で節点変位ベクトル$\{U_e\}$と節点力ベクトル$\{F_e\}$の間に剛性方程式$[K_e]\{U_e\}=\{F_e\}$が成立するときの、マトリクス$[K_e]$のことです。2次元梁の要素剛性マトリクスは、梁要素の縦弾性係数:$E$、断面積:$A$、断面2次モーメント:$I$および長さ:$L$を用いた以下の形となります。

要素剛性マトリクス作成プログラム

以下が、Python上で2次元梁の要素剛性マトリクスを作成する関数です。外部から1つの要素のインスタンスであるelemを読み込んで、そのパラメータである縦弾性係数、断面積、断面2次モーメントおよび長さを用いて$K_e$マトリクスを作成します。

import numpy as np

def get_elem_stiff(elem):
    eA_L = elem.E * elem.A / elem.L
    eI_L1 = elem.E * elem.I / elem.L
    eI_L2 = eI_L1 / elem.L
    eI_L3 = eI_L2 / elem.L
    ke = np.zeros((6, 6))
    ke[0] = [eA_L, 0.0, 0.0, -eA_L, 0.0, 0.0]
    ke[1] = [0.0, 12.0*eI_L3, 6.0*eI_L2, 0.0, -12.0*eI_L3, 6.0*eI_L2]
    ke[2] = [0.0, 6.0*eI_L2, 4.0*eI_L1, 0.0, -6.0*eI_L2, 2.0*eI_L1]
    ke[3] = [-eA_L, 0.0, 0.0, eA_L, 0.0, 0.0]
    ke[4] = [0.0, -12.0*eI_L3, -6.0*eI_L2, 0.0, 12.0*eI_L3, -6.0*eI_L2]
    ke[5] = [0.0, 6.0*eI_L2, 2.0*eI_L1, 0.0, -6.0*eI_L2, 4.0*eI_L1]
    return ke

要素剛性マトリクスの検証計算

1つの要素からなる片持ち梁モデルで検証してみます。下図のとおり、図左側の節点$i$を固定端、右側の$j$節点を自由端とし、$j$節点に集中荷重$F_x$および$F_y$、集中曲げモーメント$M_z$を載荷します。

変位および変位角の理論解は、下式を用いて、$\{U_j\}=\{u_j,v_j,\theta_j\}=\{4.10mm,-21.7mm,-0.0392rad\}$となります。

変位$u$変位$v$変位角$\theta$
軸方向荷重$\dfrac{PL}{EA}$
軸直方向荷重$\dfrac{PL^3}{3EI}$$\dfrac{PL^2}{2EI}$
集中曲げモーメント$\dfrac{ML^2}{2EI}$$\dfrac{ML}{EI}$

モデルは以下の別記事でご紹介した、JSONフォーマットで作成します。なお、数値計算上の単位については、力をN、長さをmmとします。

Python×JSONファイルの読込
(https://archit-computation.com/2021/06/06/tips-0001/)

{
    "Node": [
        {
            "Name": "NODE1",
            "X": "0.0", "Y": "0.0",
            "BCX": "1", "BCY": "1", "BCZ": "1",
            "BCRX": "1", "BCRY": "1", "BCRZ": "1"
            "FX": "0", "FY": "0", "M": "0"
        },
        {
            "Name": "NODE2",
            "X": "1000", "Y": "0.0",
            "BCX": "0", "BCY": "0", "BCZ": "0",
            "BCRX": "0", "BCRY": "0", "BCRZ": "0"
            "FX": "1000000", "FY": "10000", "M": "10000000"
        }
    ],
    "Elem": [
        {
            "Name": "ELEM1",
            "INode": "NODE1", "JNode": "NODE2",
            "E": "205000", "A": "1190", "I": "1870000"
        }
    ]
}

この例題を解くためのプログラムは以下の通りです。$K_e$マトリクスは各要素固有のものなので、要素(Elem)クラスのオブジェクト生成時に$K_e$マトリクスが作成されるようにしました。今回は1要素のみのモデルであるため、全体剛性マトリクスは$K_e$マトリクスと同一です。また、要素は$x$軸と平行(傾きはゼロ)であるため、マトリクスの回転は考慮不要です。

import json
import numpy as np
import math

class Node:
    def __init__(self, item):
        self.Name = item['Name']
        self.X = float(item['X'])
        self.Y = float(item['Y'])
        self.BCX = int(item['BCX'])
        self.BCY = int(item['BCY'])
        self.BCRZ = int(item['BCRZ'])
        self.FX = float(item['FX'])
        self.FY = float(item['FY'])
        self.MZ = float(item['MZ'])

    # Define dictionary type to be referred in Elem class    
    def get_dic_node(self):
        dic_node = {}
        dic_node[self.Name] = self
        return dic_node

class Elem:
    def __init__(self, item, dic_nodes):
        self.Name = item['Name']
        self.E = float(item['E'])
        self.A = float(item['A'])
        self.I = float(item['I'])
        self.NodeI = dic_nodes[item['INode']]
        self.NodeJ = dic_nodes[item['JNode']]
        self.L = self.get_elem_length()
        self.Ke = self.get_elem_stiff()
    
    def get_elem_length(self):
        lx = abs(self.NodeI.X - self.NodeJ.X)
        ly = abs(self.NodeI.Y - self.NodeJ.Y)
        l = math.sqrt(lx**2 + ly**2)
        return l
    
    def get_elem_stiff(self):
        eA_L = self.E * self.A / self.L
        eI_L1 = self.E * self.I / self.L
        eI_L2 = eI_L1 / self.L
        eI_L3 = eI_L2 / self.L
        ke = np.zeros((6, 6))
        ke[0] = [eA_L, 0.0, 0.0, -eA_L, 0.0, 0.0]
        ke[1] = [0.0, 12.0*eI_L3, 6.0*eI_L2, 0.0, -12.0*eI_L3, 6.0*eI_L2]
        ke[2] = [0.0, 6.0*eI_L2, 4.0*eI_L1, 0.0, -6.0*eI_L2, 2.0*eI_L1]
        ke[3] = [-eA_L, 0.0, 0.0, eA_L, 0.0, 0.0]
        ke[4] = [0.0, -12.0*eI_L3, -6.0*eI_L2, 0.0, 12.0*eI_L3, -6.0*eI_L2]
        ke[5] = [0.0, 6.0*eI_L2, 2.0*eI_L1, 0.0, -6.0*eI_L2, 4.0*eI_L1]
        return ke

class Model:
    def __init__(self, path):
        # Store all nodes, elements, force
        # and boundary conditions to each list
        self.Nodes = []
        self.Elems = []
        self.F = []
        self.BC = []

        with open(path) as f:
            df = json.load(f)

        self.dic_nodes = {}        
        for item in df['Node']:
            node = Node(item)
            self.Nodes.append(node)
            self.dic_nodes.update(node.get_dic_node())
            self.F.extend([node.FX, node.FY, node.MZ])
            self.BC.extend([node.BCX, node.BCY, node.BCRZ])

        for item in df['Elem']:
            elem = Elem(item, self.dic_nodes)
            self.Elems.append(elem)
        
        # In case there is only one element in model
        self.K = self.Elems[0].Ke

class Solver:
    @staticmethod
    def solve(model):
        i = 0
        kr = model.K
        fr = model.F
        
        # Contraction of K matrix
        for bc in model.BC:
            if bc == 1:
                kr = np.delete(kr, i, 0)
                kr = np.delete(kr, i, 1)
                fr = np.delete(fr, i, 0)
            else:
                i += 1
        
        kr_inv = np.linalg.inv(kr)
        u = np.dot(kr_inv, fr)
        print('u = ', u[0], '\nv = ', u[1], '\ntheta = ', u[2])

class Controller:
    def __init__(self):
        self.Controller = self
    
    def run(self, path):
        model = Model(path)
        Solver.solve(model)

if __name__ == "__main__":
    path = 'cantilever_w_1elem.json'
    ctl = Controller()
    ctl.run(path)

検証結果

プログラム実行の結果、理論解と一致する以下の結果が出力されました。このほか、同様に$i$端を自由端、$j$端を固定端とするケースなどの検証解析を実施することで、誤りなく$K_e$マトリクスが作成できていることを確認できます。

u =  4.099200655872105 
v =  21.73818529629144 
theta =  0.0391287335333246

カテゴリー: Tips