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

投稿者: | 2021年12月24日

関連記事で記載したとおり、3次元梁の要素剛性マトリクスは12×12のサイズとなります。今回も自分用の保管を兼ねて、 Pythonで3次元梁モデルの要素剛性マトリクスを作成するプログラムを、インターネットの海に投下しておきます。

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

2次元梁の要素剛性マトリクスは、梁要素の縦弾性係数:$E$、せん断弾性係数:$G$、断面積:$A$、$y$軸まわりの断面二次モーメント:$I_y$、$z$軸まわりの断面二次モーメント:$I_z$、サンブナンのねじり定数:$K$、および長さ:$L$を用いた以下の形となります。

$$\begin{bmatrix}
\frac{EA}{L} & & & & & & & & & & & \\
0 & \frac{12EI_z}{L^3} & & & & & & & & & & \\
0 & 0 & \frac{12EI_y}{L^3} & & & & & & & & & \\
0 & 0 & 0 & \frac{GK}{L} & & & & sym. & & & & \\
0 & 0 & \frac{-6EI_y}{L^2} & 0 & \frac{4EI_y}{L} & & & & & & & \\
0 & \frac{6EI_z}{L^2} & 0 & 0 & 0 & 0 & & & & & & \\
\frac{-EA}{L} & 0 & 0 & 0 & 0 & 0 & \frac{EA}{L} & & & & & \\
0 & \frac{-12EI_z}{L^3} & 0 & 0 & 0 & \frac{-6EI_z}{L^2} & 0 & \frac{12EI_z}{L^3} & & & & \\
0 & 0 & \frac{-12EI_y}{L^3} & 0 & \frac{6EI_y}{L^2} & 0 & 0 & 0 & \frac{12EI_y}{L^3} & & & \\
0 & 0 & 0 & \frac{-GK}{L} & 0 & 0 & 0 & 0 & 0 & \frac{GK}{L} & & \\
0 & 0 & \frac{-6EI_y}{L^2} & 0 & \frac{2EI_y}{L} & 0 & 0 & 0 & \frac{6EI_y}{L^2} & 0 & \frac{4EI_y}{L} & \\
0 & \frac{6EI_z}{L^2} & 0 & 0 & 0 & \frac{EI_z}{L} & 0 & \frac{-6EI_z}{L^2} & 0 & 0 & 0 & \frac{4EI_z}{L} \\
\end{bmatrix}$$

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

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

import numpy as np

def get_elem_stiff(elem):
    eA_L = elem.E * elem.A / elem.L
    gK_L = elem.G * elem.K / elem.L
    eIy_L1 = elem.E * elem.Iy / elem.L
    eIz_L1 = elem.E * elem.Iz / elem.L
    eIy_L2 = eIy_L1 / elem.L
    eIz_L2 = eIz_L1 / elem.L
    eIy_L3 = eIy_L2 / elem.L
    eIz_L3 = eIz_L2 / elem.L
    ke = np.zeros((12, 12))
    ke[0]  = [eA_L, 0.0, 0.0, 0.0, 0.0, 0.0,\
              -eA_L, 0.0, 0.0, 0.0, 0.0, 0.0]
    ke[1]  = [0.0, 12*eIz_L3, 0.0, 0.0, 0.0, 6*eIz_L2,\
              0.0, -12*eIz_L3, 0.0, 0.0, 0.0, 6*eIz_L2]
    ke[2]  = [0.0, 0.0, 12*eIy_L3, 0.0, -6*eIy_L2, 0.0,\
              0.0, 0.0, -12*eIy_L3, 0.0, -6*eIy_L2, 0.0]
    ke[3]  = [0.0, 0.0, 0.0, gK_L, 0.0, 0.0,\
              0.0, 0.0, 0.0, -gK_L, 0.0, 0.0]
    ke[4]  = [0.0, 0.0, -6*eIy_L2, 0.0, 4*eIy_L1, 0.0,\
              0.0, 0.0, 6*eIy_L2, 0.0, 2*eIy_L1, 0.0]
    ke[5]  = [0.0, 6*eIz_L2, 0.0, 0.0, 0.0, 4*eIz_L1,\
              0.0, -6*eIz_L2, 0.0, 0.0, 0.0, 2*eIz_L1]
    ke[6]  = [-eA_L, 0.0, 0.0, 0.0, 0.0, 0.0,\
              eA_L, 0.0, 0.0, 0.0, 0.0, 0.0]
    ke[7]  = [0.0, -12*eIz_L3, 0.0, 0.0, 0.0, -6*eIz_L2,\
              0.0, 12*eIz_L3, 0.0, 0.0, 0.0, -6*eIz_L2]
    ke[8]  = [0.0, 0.0, -12*eIy_L3, 0.0, 6*eIy_L2, 0.0,\
              0.0, 0.0, 12*eIy_L3, 0.0, 6*eIy_L2, 0.0]
    ke[9]  = [0.0, 0.0, 0.0, -gK_L, 0.0, 0.0,\
              0.0, 0.0, 0.0, gK_L, 0.0, 0.0]
    ke[10] = [0.0, 0.0, -6*eIy_L2, 0.0, 2*eIy_L1, 0.0,\
              0.0, 0.0, 6*eIy_L2, 0.0, 4*eIy_L1, 0.0]
    ke[11] = [0.0, 6*eIz_L2, 0.0, 0.0, 0.0, 2*eIz_L1,\
              0.0, -6*eIz_L2, 0.0, 0.0, 0.0, 4*eIz_L1]
    return ke

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

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

$j$節点における変位および変位角の理論解は、下式を用いて、$\{U_j\}=\{u_j,v_j,w_j,\theta_{jx},\theta_{jy},\theta_{jz}\}=\{4.10mm,21.7mm,46.7mm,0.063rad,-0.066rad,0.0392rad\}$となります。

外力$u$$v$$w$$\theta_x$$\theta_y$$\theta_z$
$F_x$$\dfrac{F_xL}{EA}$
$F_y$$\dfrac{F_yL^3}{3EI_z}$$\dfrac{F_yL^2}{2EI_z}$
$F_z$$\dfrac{F_zL^3}{3EI_y}$$\dfrac{F_zL^2}{2EI_y}$
$M_x$$\dfrac{M_xL}{GK}$
$M_y$$\dfrac{M_yL^2}{2EI_y}$$\dfrac{M_yL}{EI_y}$
$M_z$$\dfrac{M_zL^2}{2EI_z}$$\dfrac{M_zL}{EI_z}$

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

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

{
    "Node": [
        {
            "Name": "NODE1",
            "X": "0.0", "Y": "0.0", "Z": "0.0",
            "BCX": "1", "BCY": "1", "BCZ": "1",
            "BCRX": "1", "BCRY": "1", "BCRZ": "1",
            "FX": "0", "FY": "0", "FZ": "0",
            "MX": "0", "MY": "0", "MZ": "0"
        },
        {
            "Name": "NODE2",
            "X": "1000", "Y": "0.0", "Z": "0.0",
            "BCX": "0", "BCY": "0", "BCZ": "0",
            "BCRX": "0", "BCRY": "0", "BCRZ": "0",
            "FX": "1000000", "FY": "10000", "FZ": "5000",
            "MX": "10000000", "MY": "500000", "MZ": "10000000"
        }
    ],
    "Elem": [
        {
            "Name": "ELEM1",
            "INode": "NODE1", "JNode": "NODE2",
            "E": "205000", "G": "78800", "A": "1190",
            "IY": "148000", "IZ": "1870000"
        }
    ]
}

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

Python×2次元梁の要素剛性マトリクスの作成
(https://archit-computation.com/2021/06/15/tips-0003/)

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.Z = float(item['Z'])
        self.BCX = int(item['BCX'])
        self.BCY = int(item['BCY'])
        self.BCZ = int(item['BCZ'])
        self.BCRX = int(item['BCRX'])
        self.BCRY = int(item['BCRY'])
        self.BCRZ = int(item['BCRZ'])
        self.FX = float(item['FX'])
        self.FY = float(item['FY'])
        self.FZ = float(item['FZ'])
        self.MX = float(item['MX'])
        self.MY = float(item['MY'])
        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.G = float(item['G'])
        self.A = float(item['A'])
        self.IY = float(item['IY'])
        self.IZ = float(item['IZ'])
        self.K = self.IY + self.IZ
        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)
        lz = abs(self.NodeI.Z - self.NodeJ.Z)
        l = math.sqrt(lx**2 + ly**2 + lz**2)
        return l
    
    def get_elem_stiff(elem):
        eA_L = elem.E * elem.A / elem.L
        gK_L = elem.G * elem.K / elem.L
        eIy_L1 = elem.E * elem.IY / elem.L
        eIz_L1 = elem.E * elem.IZ / elem.L
        eIy_L2 = eIy_L1 / elem.L
        eIz_L2 = eIz_L1 / elem.L
        eIy_L3 = eIy_L2 / elem.L
        eIz_L3 = eIz_L2 / elem.L
        ke = np.zeros((12, 12))
        ke[0]  = [eA_L, 0.0, 0.0, 0.0, 0.0, 0.0,\
                  -eA_L, 0.0, 0.0, 0.0, 0.0, 0.0]
        ke[1]  = [0.0, 12*eIz_L3, 0.0, 0.0, 0.0, 6*eIz_L2,\
                  0.0, -12*eIz_L3, 0.0, 0.0, 0.0, 6*eIz_L2]
        ke[2]  = [0.0, 0.0, 12*eIy_L3, 0.0, -6*eIy_L2, 0.0,\
                  0.0, 0.0, -12*eIy_L3, 0.0, -6*eIy_L2, 0.0]
        ke[3]  = [0.0, 0.0, 0.0, gK_L, 0.0, 0.0,\
                  0.0, 0.0, 0.0, -gK_L, 0.0, 0.0]
        ke[4]  = [0.0, 0.0, -6*eIy_L2, 0.0, 4*eIy_L1, 0.0,\
                  0.0, 0.0, 6*eIy_L2, 0.0, 2*eIy_L1, 0.0]
        ke[5]  = [0.0, 6*eIz_L2, 0.0, 0.0, 0.0, 4*eIz_L1,\
                  0.0, -6*eIz_L2, 0.0, 0.0, 0.0, 2*eIz_L1]
        ke[6]  = [-eA_L, 0.0, 0.0, 0.0, 0.0, 0.0,\
                  eA_L, 0.0, 0.0, 0.0, 0.0, 0.0]
        ke[7]  = [0.0, -12*eIz_L3, 0.0, 0.0, 0.0, -6*eIz_L2,\
                  0.0, 12*eIz_L3, 0.0, 0.0, 0.0, -6*eIz_L2]
        ke[8]  = [0.0, 0.0, -12*eIy_L3, 0.0, 6*eIy_L2, 0.0,\
                  0.0, 0.0, 12*eIy_L3, 0.0, 6*eIy_L2, 0.0]
        ke[9]  = [0.0, 0.0, 0.0, -gK_L, 0.0, 0.0,\
                  0.0, 0.0, 0.0, gK_L, 0.0, 0.0]
        ke[10] = [0.0, 0.0, -6*eIy_L2, 0.0, 2*eIy_L1, 0.0,\
                  0.0, 0.0, 6*eIy_L2, 0.0, 4*eIy_L1, 0.0]
        ke[11] = [0.0, 6*eIz_L2, 0.0, 0.0, 0.0, 2*eIz_L1,\
                  0.0, -6*eIz_L2, 0.0, 0.0, 0.0, 4*eIz_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.FZ, \
                           node.MX, node.MY, node.MZ])
            self.BC.extend([node.BCX, node.BCY, node.BCZ, \
                            node.BCRX, node.BCRY, 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], '\nw = ', u[2], \
              '\ntheta_x = ', u[3], \
              '\ntheta_y = ', u[4], \
              '\ntheta_z = ', u[5])

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_3d.json'
    ctl = Controller()
    ctl.run(path)

検証結果

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

u =  4.099200655872105
v =  21.73818529629144
w =  46.69303449791255
theta_x =  0.06288580441005569
theta_y =  -0.06591957811470006
theta_z =  0.0391287335333246

カテゴリー: Tips

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

CAPTCHA