有限要素法構造解析プログラムを作成するにあたって、マトリクス（行列）の記述は骨の折れる作業です。たとえば要素剛性マトリクスであれば、2次元梁なら6×6、3次元梁なら12×12もの行列になりますが、もちろん1か所のタイプミスでも答えが合わなくなるので、エラーの発見が大変です。
そこで自分用の保管を兼ねて、Pythonで2次元梁モデルの要素剛性マトリクスを作成するプログラムを、インターネットの海に投下しておきます。
2次元梁の要素剛性マトリクス
要素剛性マトリクスについて大雑把に解説すると、ある1つの要素内で節点変位ベクトル$\{U_e\}$と節点力ベクトル$\{F_e\}$の間に剛性方程式$[K_e]\{U_e\}=\{F_e\}$が成立するときの、マトリクス$[K_e]$のことです。2次元梁の要素剛性マトリクスは、梁要素の縦弾性係数：$E$、断面積：$A$、断面2次モーメント：$I$および長さ：$L$を用いた以下の形となります。
要素剛性マトリクス作成プログラム
以下が、Python上で2次元梁の要素剛性マトリクスを作成する関数です。外部から1つの要素のインスタンスであるelemを読み込んで、そのパラメータである縦弾性係数、断面積、断面2次モーメントおよび長さを用いて$K_e$マトリクスを作成します。
import numpy as np
def get_elem_stiff(elem):
eA_L = elem.E * elem.A / elem.L
eI_L1 = elem.E * elem.I / elem.L
eI_L2 = eI_L1 / elem.L
eI_L3 = eI_L2 / elem.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
要素剛性マトリクスの検証計算
1つの要素からなる片持ち梁モデルで検証してみます。下図のとおり、図左側の節点$i$を固定端、右側の$j$節点を自由端とし、$j$節点に集中荷重$F_x$および$F_y$、集中曲げモーメント$M_z$を載荷します。
変位および変位角の理論解は、下式を用いて、$\{U_j\}=\{u_j,v_j,\theta_j\}=\{4.10mm,-21.7mm,-0.0392rad\}$となります。
| 変位$u$ | 変位$v$ | 変位角$\theta$ | |
| 軸方向荷重 | $\dfrac{PL}{EA}$ | ― | ー |
| 軸直方向荷重 | ― | $\dfrac{PL^3}{3EI}$ | $\dfrac{PL^2}{2EI}$ |
| 集中曲げモーメント | ― | $\dfrac{ML^2}{2EI}$ | $\dfrac{ML}{EI}$ |
モデルは以下の別記事でご紹介した、JSONフォーマットで作成します。なお、数値計算上の単位については、力をN、長さをmmとします。
Python×JSONファイルの読込
(https://archit-computation.com/2021/06/06/tips-0001/)
{
"Node": [
{
"Name": "NODE1",
"X": "0.0", "Y": "0.0",
"BCX": "1", "BCY": "1", "BCZ": "1",
"BCRX": "1", "BCRY": "1", "BCRZ": "1"
"FX": "0", "FY": "0", "M": "0"
},
{
"Name": "NODE2",
"X": "1000", "Y": "0.0",
"BCX": "0", "BCY": "0", "BCZ": "0",
"BCRX": "0", "BCRY": "0", "BCRZ": "0"
"FX": "1000000", "FY": "10000", "M": "10000000"
}
],
"Elem": [
{
"Name": "ELEM1",
"INode": "NODE1", "JNode": "NODE2",
"E": "205000", "A": "1190", "I": "1870000"
}
]
}この例題を解くためのプログラムは以下の通りです。$K_e$マトリクスは各要素固有のものなので、要素(Elem)クラスのオブジェクト生成時に$K_e$マトリクスが作成されるようにしました。今回は1要素のみのモデルであるため、全体剛性マトリクスは$K_e$マトリクスと同一です。また、要素は$x$軸と平行（傾きはゼロ）であるため、マトリクスの回転は考慮不要です。
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.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'])
# 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.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)
l = math.sqrt(lx**2 + ly**2)
return l
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
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])
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], '\ntheta = ', u[2])
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.json'
ctl = Controller()
ctl.run(path)
検証結果
プログラム実行の結果、理論解と一致する以下の結果が出力されました。このほか、同様に$i$端を自由端、$j$端を固定端とするケースなどの検証解析を実施することで、誤りなく$K_e$マトリクスが作成できていることを確認できます。
u = 4.099200655872105
v = 21.73818529629144
theta = 0.0391287335333246