関連記事で記載したとおり、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