Python×3次元弾性骨組モデルの静解析プログラム

投稿者: | 2022年1月21日

(本記事は、プログラムの改善により随時アップデートしていきます。)

本記事では、「Excelで解く3次元建築構造解析(藤井大地 著)」の、3次元弾性骨組モデルの静解析プログラムをPythonで作成してみました。この記事で扱っているインプットファイルとソースコードは、以下よりダウンロード可能です。

このプログラムを実行することで、以下のような3次元弾性骨組モデルの変位図が出力されます。

解析モデル

本記事では、同著書3章の骨組解析の例題(下図)を用いて、VBAによる解析結果との比較により検証を行います。

まずは、JSONフォーマット(以下記事参照)を用いたモデルファイルを作成します。モデルの節点座標および物性値をモデルに落とし込むと、以下のようになります。今回は、固定端の節点をNODE1~NODE4、柱―梁の交差部の節点をNODE5~NODE8としています。また、柱部分の要素をELEM1~ELEM4、梁部分の要素をELEM5~ELEM8としています。

{
    "Node": [
        {
            "Name": "NODE1",
            "X": "0", "Y": "0", "Z": "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": "600", "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": "NODE3",
            "X": "600", "Y": "600", "Z": "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": "NODE4",
            "X": "0", "Y": "600", "Z": "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": "NODE5",
            "X": "0", "Y": "0", "Z": "400",
            "BCX": "0", "BCY": "0", "BCZ": "0",
            "BCRX": "0", "BCRY": "0", "BCRZ": "0",
            "FX": "10", "FY": "0", "FZ": "0",
            "MX": "0", "MY": "0", "MZ": "0"
        },
        {
            "Name": "NODE6",
            "X": "600", "Y": "0", "Z": "400",
            "BCX": "0", "BCY": "0", "BCZ": "0",
            "BCRX": "0", "BCRY": "0", "BCRZ": "0",
            "FX": "0", "FY": "0", "FZ": "0",
            "MX": "0", "MY": "0", "MZ": "0"
        },
        {
            "Name": "NODE7",
            "X": "600", "Y": "600", "Z": "400",
            "BCX": "0", "BCY": "0", "BCZ": "0",
            "BCRX": "0", "BCRY": "0", "BCRZ": "0",
            "FX": "0", "FY": "0", "FZ": "0",
            "MX": "0", "MY": "0", "MZ": "0"
        },
        {
            "Name": "NODE8",
            "X": "0", "Y": "600", "Z": "400",
            "BCX": "0", "BCY": "0", "BCZ": "0",
            "BCRX": "0", "BCRY": "0", "BCRZ": "0",
            "FX": "0", "FY": "0", "FZ": "0",
            "MX": "0", "MY": "0", "MZ": "0"
        }
    ],
    "Elem": [
        {
            "Name": "ELEM1",
            "INode": "NODE1", "JNode": "NODE5",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        },
        {
            "Name": "ELEM2",
            "INode": "NODE2", "JNode": "NODE6",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        },
        {
            "Name": "ELEM3",
            "INode": "NODE3", "JNode": "NODE7",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        },
        {
            "Name": "ELEM4",
            "INode": "NODE4", "JNode": "NODE8",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        },
        {
            "Name": "ELEM5",
            "INode": "NODE5", "JNode": "NODE6",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        },
        {
            "Name": "BEAM6",
            "INode": "NODE6", "JNode": "NODE7",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        },
        {
            "Name": "BEAM7",
            "INode": "NODE7", "JNode": "NODE8",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        },
        {
            "Name": "BEAM8",
            "INode": "NODE8", "JNode": "NODE5",
            "E": "20600", "G": "7923",
            "A": "100", "Iy": "833", "Iz": "833", "K": "1666"
        }
    ]
}

解析プログラム

解析プログラムは以下の通りとなります。なお、このプログラムは上記のJSONフォーマットで作成された3次元骨組モデルであれば、例題以外の形状のモデルをインプットとしても解析可能なものとなっています。

import json
import numpy as np
import math
import matplotlib.pyplot as plt
import os

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'])
        self.DX, self.DY, self.DZ = 0., 0., 0.
        self.THX, self.THY, self.THZ = 0., 0., 0.

    # 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 = float(item['K'])
        self.NodeI = dic_nodes[item['INode']]
        self.NodeJ = dic_nodes[item['JNode']]
        self.Lx = abs(self.NodeI.X - self.NodeJ.X)
        self.Ly = abs(self.NodeI.Y - self.NodeJ.Y)
        self.Lz = abs(self.NodeI.Z - self.NodeJ.Z)
        self.L = math.sqrt(self.Lx**2 + self.Ly**2 + self.Lz**2)
        self.Ke = self.get_elem_stiff()
        self.T = self.get_rotation_matrix()
        self.Kge = self.trans_elem_stiff_to_global()
    
    def get_elem_stiff(self):
        eA_L = self.E * self.A / self.L
        gK_L = self.G * self.K / self.L
        eIy_L1 = self.E * self.Iy / self.L
        eIz_L1 = self.E * self.Iz / self.L
        eIy_L2 = eIy_L1 / self.L
        eIz_L2 = eIz_L1 / self.L
        eIy_L3 = eIy_L2 / self.L
        eIz_L3 = eIz_L2 / self.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
    
    def trans_elem_stiff_to_global(self):
        kge = np.dot(np.dot(np.transpose(self.T), self.Ke), self.T)
        return kge
    
    def get_rotation_matrix(self):
        zero = 1.0E-10
        horz_L = math.sqrt(self.Lx**2 + self.Ly**2)
        lx = (self.NodeJ.X - self.NodeI.X) / self.L
        mx = (self.NodeJ.Y - self.NodeI.Y) / self.L
        nx = (self.NodeJ.Z - self.NodeI.Z) / self.L
        ram = math.sqrt(lx**2 + mx**2)
        if horz_L < zero:
            ly, my, ny = nx, 0.0, 0.0
            lz, mz, nz = 0.0, 1.0, 0.0
        else:
            ly, my, ny = -mx/ram, lx/ram, 0.0
            lz, mz, nz = -lx*nx/ram, -mx*nx/ram, ram
        t = np.zeros(shape=(12,12))
        t[0]  = [ lx,  mx,  nx, 0.0, 0.0, 0.0,\
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
        t[1]  = [ ly,  my,  ny, 0.0, 0.0, 0.0,\
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
        t[2]  = [ lz,  mz,  nz, 0.0, 0.0, 0.0,\
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
        t[3]  = [0.0, 0.0, 0.0,  lx,  mx,  nx,\
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,]
        t[4]  = [0.0, 0.0, 0.0,  ly,  my,  ny,\
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
        t[5]  = [0.0, 0.0, 0.0,  lz,  mz,  nz,\
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
        t[6]  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
                  lx,  mx,  nx, 0.0, 0.0, 0.0]
        t[7]  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
                  ly,  my,  ny, 0.0, 0.0, 0.0]
        t[8]  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
                  lz,  mz,  nz, 0.0, 0.0, 0.0]
        t[9]  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
                 0.0, 0.0, 0.0,  lx,  mx,  nx]
        t[10] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
                 0.0, 0.0, 0.0,  ly,  my,  ny]
        t[11] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\
                 0.0, 0.0, 0.0,  lz,  mz,  nz]
        return t

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])

        size_K = len(self.Nodes) * 6
        self.K = np.array([[0] * size_K for i in range(size_K)])
        
        for item in df['Elem']:
            elem = Elem(item, self.dic_nodes)
            self.K = self.map_on_global_stiff(self.K, elem)
            self.Elems.append(elem)
        
    def map_on_global_stiff(self, k, elem):
        idx = []
        kge_map = [[0] * len(k) for i in range(len(k))]
        i_idx = int(elem.NodeI.Name.replace('NODE',''))
        j_idx = int(elem.NodeJ.Name.replace('NODE',''))
        idx = [6*i_idx-6, 6*i_idx-5, 6*i_idx-4,\
               6*i_idx-3, 6*i_idx-2, 6*i_idx-1,\
               6*j_idx-6, 6*j_idx-5, 6*j_idx-4,\
               6*j_idx-3, 6*j_idx-2, 6*j_idx-1]
        i = 0
        for ikeg in idx:
            j = 0
            for jkeg in idx:
                kge_map[ikeg][jkeg] = elem.Kge[i][j]
                j+=1
            i+=1
        kge_map = np.array(kge_map)
        k = k + kge_map
        return k
    
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
        u = np.linalg.solve(kr, fr)
        i = 0
        for bc in model.BC:
            if bc == 1:
                u = np.insert(u, i, 0.)
                i += 1
        return u

class Output:
    def __init__(self, model):
        self.U = []
        self.U = Solver.solve(model)
        self.store_output(self.U, model)
    
    def store_output(self, u, model):
        i = 0
        for node in model.Nodes:
            node.DX = u[i]
            node.DY = u[i+1]
            node.DZ = u[i+2]
            node.THX = u[i+3]
            node.THY = u[i+4]
            node.THZ = u[i+5]
            i += 6

class StaticViewer:
    def __init__(self, model):
        self.beamsToLine(model)
        
    def beamToLine_origin(self, ax, elem):
        x0 = elem.NodeI.X
        y0 = elem.NodeI.Y
        z0 = elem.NodeI.Z
        x1 = elem.NodeJ.X
        y1 = elem.NodeJ.Y
        z1 = elem.NodeJ.Z
        return ax.plot([x0, x1], [y0, y1], [z0, z1],\
                   color='blue', linestyle='solid')

    def beamToLine_result(self, ax, elem):
        disp_mag = 50.
        x0 = elem.NodeI.X + disp_mag * elem.NodeI.DX
        y0 = elem.NodeI.Y + disp_mag * elem.NodeI.DY
        z0 = elem.NodeI.Z + disp_mag * elem.NodeI.DZ
        x1 = elem.NodeJ.X + disp_mag * elem.NodeJ.DX
        y1 = elem.NodeJ.Y + disp_mag * elem.NodeJ.DY
        z1 = elem.NodeJ.Z + disp_mag * elem.NodeJ.DZ
        ax.plot([x0, x1], [y0, y1], [z0, z1],\
                color='red', linestyle='solid')

    def beamsToLine(self, model):
        fig = plt.figure(figsize = (6, 6))
        ax = fig.add_subplot(111, projection='3d')
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.set_zlabel("z")
        u_lim, l_lim = self.model_view_range(model)
        ax.set_xlim(l_lim, u_lim)
        ax.set_ylim(l_lim, u_lim)
        ax.set_zlim(l_lim, u_lim)
        
        for elem in model.Elems:
            self.beamToLine_origin(ax, elem)
            self.beamToLine_result(ax, elem)
        filename = '3d_frame_elasto_static_fig.png'
        plt.savefig(filename)
    
    def model_view_range(self, model):
        u_lim, l_lim = 0., 0.
        for node in model.Nodes:
            if u_lim < node.X:
                u_lim = node.X
            elif node.X < l_lim:
                l_lim = node.X
            if u_lim < node.Y:
                u_lim = node.Y
            elif node.Y < l_lim:
                l_lim = node.Y
        u_lim = u_lim + 0.2 * (u_lim - l_lim)
        l_lim = l_lim - 0.2 * (u_lim - l_lim)
        return u_lim, l_lim

class Controller:
    def __init__(self):
        self.Controller = self
    
    def run(self, path):
        model = Model(path)
        output = Output(model)
        StaticViewer(model)
        print(output.U)

if __name__ == "__main__":
    path = '3d_frame_example.json'
    ctl = Controller()
    ctl.run(path)
    os.system('PAUSE')

プログラム実行結果

各節点の変位および変位角を示す、変位マトリクスが以下の通り出力されました(NODE1~NODE8それぞれの変位および回転角$(u, v, w, \theta_x, \theta_y, \theta_z)$)。「Excelで解く3次元建築構造解析」の記載と同じ数値であることが確認できました。

[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  2.09228944e+00 -2.31580600e-01  3.72653975e-04  3.84703517e-04
  3.74343509e-03  1.64257994e-03  2.09083353e+00  2.31580600e-01
 -3.72653975e-04 -3.84703517e-04  3.73967053e-03  1.64098956e-03
  3.95146071e-01  2.31580600e-01 -1.45107040e-04 -3.84703517e-04
  9.21721644e-04  1.64098956e-03  3.95146171e-01 -2.31580600e-01
  1.45107040e-04  3.84703517e-04  9.22052123e-04  1.64257994e-03]

関連TIPS記事

本プログラムの作成に当たっては、以下のPythonプログラミングに関するTIPSを使用しています。

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

Python×基本的な連立方程式の解法
(https://archit-computation.com/2021/06/12/tips-0002/)

Python×3次元梁の要素剛性マトリクスの作成
(https://archit-computation.com/2021/12/24/tips-0004/)

コメントを残す

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

CAPTCHA