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

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

2d_frame_elasto_static

先に結果からお示ししますが、このプログラムを実行することで、以下のような2次元弾性骨組モデルの変位図が出力されます。

解析モデル

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

まずは、JSONフォーマット（以下記事参照）を用いたモデルファイルを作成します。モデルの節点座標および物性値をモデルに落とし込むと、以下のようになります。今回は、固定端をNODE1、折れ点をNODE2、自由端をNODE3とし、NODE1―NODE2間の要素をELEM1、NODE2―NODE3間の要素をELEM1としています。

{
    "Node": [
        {
            "Name": "NODE1",
            "X": "0.0", "Y": "0.0",
            "BCX": "1", "BCY": "1", "BCRZ": "1",
            "FX": "0", "FY": "0", "MZ": "0"
        },
        {
            "Name": "NODE2",
            "X": "0.0", "Y": "200.0",
            "BCX": "0", "BCY": "0", "BCRZ": "0",
            "FX": "0", "FY": "0", "MZ": "0"
        },
        {
            "Name": "NODE3",
            "X": "100.0", "Y": "200.0",
            "BCX": "0", "BCY": "0", "BCRZ": "0",
            "FX": "0", "FY": "-10", "MZ": "0"
        }
    ],
    "Elem": [
        {
            "Name": "ELEM1",
            "INode": "NODE1", "JNode": "NODE2",
            "E": "2100", "A": "100", "I": "833"
        },
        {
            "Name": "ELEM2",
            "INode": "NODE2", "JNode": "NODE3",
            "E": "2100", "A": "100", "I": "833"
        }
    ]
}

解析プログラム

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

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.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'])
        self.DX, self.DY, 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.A = float(item['A'])
        self.I = float(item['I'])
        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.L = math.sqrt(self.Lx**2 + self.Ly**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
        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
    
    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):
        cosa = self.Lx / self.L
        sina = self.Ly / self.L
        t = ([ [ cosa, sina, 0.0,   0.0,  0.0, 0.0],
               [-sina, cosa, 0.0,   0.0,  0.0, 0.0],
               [  0.0,  0.0, 1.0,   0.0,  0.0, 0.0],
               [  0.0,  0.0, 0.0,  cosa, sina, 0.0],
               [  0.0,  0.0, 0.0, -sina, cosa, 0.0],
               [  0.0,  0.0, 0.0,   0.0,  0.0, 1.0]
               ])
        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.MZ])
            self.BC.extend([node.BCX, node.BCY, node.BCRZ])

        size_K = len(self.Nodes) * 3
        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 = int(elem.NodeI.Name.replace('NODE',''))
        j = int(elem.NodeJ.Name.replace('NODE',''))
        idx = [3*i-3, 3*i-2, 3*i-1, 3*j-3, 3*j-2, 3*j-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.THZ = u[i+2]
            i += 3

class Viewer:
    def __init__(self, model):
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        u_lim, l_lim = Utils.model_view_range(model)
        ax.set_xlim(l_lim, u_lim)
        ax.set_ylim(l_lim, u_lim)
        ax.set_aspect('equal')
        for elem in model.Elems:
            Viewer.lines_origin(ax, elem)
            Viewer.lines_result(ax, elem)
        plt.show()

    def lines_origin(ax, elem):
        x0 = elem.NodeI.X
        y0 = elem.NodeI.Y
        x1 = elem.NodeJ.X
        y1 = elem.NodeJ.Y
        ax.plot([x0, x1], [y0, y1], color='blue', linestyle='dashed')
    
    def lines_result(ax, elem):
        x0 = elem.NodeI.X + elem.NodeI.DX
        y0 = elem.NodeI.Y + elem.NodeI.DY
        x1 = elem.NodeJ.X + elem.NodeJ.DX
        y1 = elem.NodeJ.Y + elem.NodeJ.DY
        ax.plot([x0, x1], [y0, y1], color='red', linestyle='solid')

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

class Utils:
    @staticmethod    
    def model_view_range(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

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

プログラム実行結果

各節点の変位および変位角が以下の通り出力されました（NODE1、NODE2それぞれの変位および回転角$(u_1, v_1, \theta_1, u_2, v_2, \theta_2)$）。「Excelで解く構造力学」の記載と同じ数値であることが確認できました。

[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  1.14331447e+01
 -9.52380952e-03 -1.14331447e-01  1.14331447e+01 -1.33481926e+01
 -1.42914309e-01]

また冒頭で述べた通り、変位図も出力されました。

月	火	水	木	金	土	日
		1	2	3	4	5
6	7	8	9	10	11	12
13	14	15	16	17	18	19
20	21	22	23	24	25	26
27	28	29	30

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

解析モデル

解析プログラム

プログラム実行結果

関連TIPS記事