(本記事は、プログラムの改善により随時アップデートしていきます。)
本記事では、「Excelで解く構造力学(藤井大地 著)」の、2次元弾性骨組モデルの静解析プログラムをPythonで作成してみました。この記事で扱っているインプットファイルとソースコードは、以下よりダウンロード可能です。
先に結果からお示ししますが、このプログラムを実行することで、以下のような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]
また冒頭で述べた通り、変位図も出力されました。

関連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×2次元梁の要素剛性マトリクスの作成
(https://archit-computation.com/2021/06/15/tips-0003/)