（本記事は、プログラムの改善により随時アップデートしていきます。）
本記事では、「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/)