はじめに
仕事で部署移動したこともあり、忙しくて更新の間がかなり空いてしまいました…
前回は有限要素法に関する関係式を整理まで行ったので今回は実際に実装してみます。
対象として扱うのは4角形1次要素の平面応力状態における静弾性解析です。
実装の流れ
まずは実装の流れを整理してみます。
先回の記事にて関係式を整理しました。
最終的にやりたいことは以下の方程式を解くことです
ここでは全節点の荷重をまとめたベクトル、
は全体剛性マトリクス、
は全節点の変位をまとめたベクトルです。
このためには要素ごとに要素剛性マトリクスを算出する必要があります。
要素剛性マトリクスは以下の式で表せます。
したがって以下のステップで実装していきます。
- Dマトリクス
の算出
- Bマトリクス
の算出
- 要素剛性マトリクス
の算出
- 全体剛性マトリクス
の算出
- 連立方程式を解き、全節点の変位
を求める
初めから順番にやっていきます。今回はBマトリクス算出まで行います。
実装
以前の記事で作成したコードと、クラス要件をベースに実装を進めます。
以前の記事では基本的な計算はすべてSolverクラスが実行するようなことが書かれてたのでそっちに実行するようにします。
まずは過去記事のおさらいからです。過去記事ではメッシュの生成までやってました。全コード載せときます。
import numpy as np import matplotlib.pyplot as plt from matplotlib import patches class Node: """ Nodeのインデックス、座標および境界条件をまとめておくクラス(2次元用) Attributes: id (int): Nodeのインデックス番号(重複不可) x, y (float): Nodeのx, y座標 """ def __init__(self, idx, x, y): self.id = idx self.x = x self.y = y def __str__(self): return 'Node {}: x {}, y {}'.format(self.id, self.x, self.y) class Element: """ 4角形1次要素Elementのインデックス、要素内のNode情報および要素の厚みをまとめておくクラス Attributes: id (int): Elementのインデックス番号(重複不可) node ([int, int, int, int]): 要素内のNodeのインデックスのリスト(左下から反時計回りで指定) thickness (float): Elementの厚み xy (ndarray(size: 2, 4)): Element内のNodeのxy座標まとめ """ def __init__(self, idx, node_idxs, thickness): self.id = idx self.node = node_idxs self.thickness = thickness self.xy = None def __str__(self): return 'Element {}: Nodes {}, Thickness: {}'.format(self.id, self.node, self.thickness) def get_coordination(self, nodes_dict): """ 要素の各節点のx,y座標を求めてself.xyに代入。 Args: nodes_dict ({Node.id: Node}): Nodeの辞書 """ res = [] for node_idx in self.node: res.append([nodes_dict[node_idx].x, nodes_dict[node_idx].y]) self.xy = np.array(res) class Mesh: """ NodeオブジェクトとElementオブジェクトをまとめておくクラス Attributes: nodes ({Node.id: Node}): Nodeの辞書 elements (list): Elementのリスト """ def __init__(self, nodes_dict, elements): self.nodes = nodes_dict self.elements = elements self.get_element_coord() def get_element_coord(self): for elm in self.elements: elm.get_coordination(self.nodes) def create_mesh(length, height, mesh_size, thickness): """ モデルを格子点に分割し、全ての点のNodeクラスおよびそれらを使用した4角形1次要素を作成する. ただし、モデル形状は四角形とする。 Args: length, height (float): モデルの長さ、高さ mesh_size (float): 分割したいメッシュのサイズ thickness (float): 要素の厚み(今回は全て一定とする) Returns: nodes_dict ({Node.id: Node}): 作成したNodeの辞書 elsms (list): 作成したElementクラスをまとめたリスト """ # Nodeの作成 num_mesh_len = int(length / mesh_size) num_mesh_hei = int(height / mesh_size) if num_mesh_len == 0: num_mesh_len = 1 if num_mesh_hei == 0: num_mesh_hei = 1 x = np.linspace(0, length, num_mesh_len + 1) y = np.linspace(0, height, num_mesh_hei + 1) X, Y = np.meshgrid(x, y) X = X.ravel() Y = Y.ravel() nodes_dict = {} for i, coord in enumerate(zip(X, Y)): nodes_dict[i] = Node(i, coord[0], coord[1]) # Elementの作成 node_idx = 0 elem_idx = 0 elems = [] for i in range(num_mesh_hei): for j in range(num_mesh_len + 1): if (node_idx + 1) % (num_mesh_len + 1) == 0: node_idx += 1 continue else: node_idxs = [node_idx, node_idx + 1, node_idx + num_mesh_len + 2, node_idx + num_mesh_len + 1] elems.append(Element(elem_idx, node_idxs, thickness)) node_idx += 1 elem_idx += 1 return nodes_dict, elems def plot_mesh(mesh): fig, ax = plt.subplots() for elem in mesh.elements: patch = patches.Polygon(xy=elem.xy, ec='black') ax.add_patch(patch) text_xy = np.mean(elem.xy, axis=0) ax.text(text_xy[0], text_xy[1], elem.id, fontsize=12, va='center', ha='center') for node in mesh.nodes.values(): ax.scatter(node.x, node.y, fc='black', s=100) ax.text(node.x, node.y, node.id, fontsize=8, color='white', va='center', ha='center') ax.autoscale() ax.set_aspect('equal', 'box') plt.show() # モデル、要素サイズ定義 MESH_SIZE = 5. LENGTH = 10. HEIGHT = 10. THICKNESS = 1. nodes_dict, elems = create_mesh(LENGTH, HEIGHT, MESH_SIZE, THICKNESS) mesh = Mesh(0, nodes_dict, elems) for node in mesh.nodes.values(): print(node) for elem in mesh.elements: print(elem) plot_mesh(mesh)
Dマトリクスの算出
Dマトリクスは各要素ごとに以下の式で決まってます
はヤング率、
はポアソン比です。
基本的には材料特性で決まるものなので、まずは材料クラスを作成してヤング率、ポアソン比を定義しておきます
class Material: """ 材料特性をまとめておくクラス Attributes: E (float): ヤング率 nu (float): ポアソン比 """ def __init__(self, E, nu): self.E = E self.nu = nu
新たにソルバークラスを作成して、メッシュクラス及び上記の材料クラスを参照できるようにしておき、Dマトリクスもここで算出します
class Solver: """ 有限要素法の計算を行うクラス Attributes: mesh (Mesh): Meshクラス material (Material): Materialクラス D (array 3x3): Dマトリクス """ def __init__(self, mesh, material): self.mesh = mesh self.material = material self.D = self.calc_D() def calc_D(self): """ Dマトリクスを算出。 Returns: D (array 3x3): Dマトリクス """ D = np.array([[1, self.material.nu, 0], [self.material.nu, 1, 0], [0, 0, (1 - self.material.nu) / 2]]) \ * self.material.E / (1 - self.material.nu ** 2) return D
ちゃんとできてるか確認しておきましょう
MESH_SIZE = 5. LENGTH = 20. HEIGHT = 10. THICKNESS = 1. E = 208000 # ヤング率 NU = 0.27 # ポアソン比 nodes_dict, elems = create_mesh(LENGTH, HEIGHT, MESH_SIZE, THICKNESS) material = Material(E, NU) mesh = Mesh(0, nodes_dict, elems) solver = Solver(mesh, material) print(solver.D)
[[224355.51720419 60575.98964513 0. ] [ 60575.98964513 224355.51720419 0. ] [ 0. 0. 81889.76377953]]
まあちゃんとできてそうです。
Bマトリクスの算出
続いてBマトリクスを求めていきましょう。Bマトリクスは以下の2式から求められます。
まずヤコビアン以外の部分を整理すると以下のようになります
ヤコビアンの各成分は以下のように整理できます。
ここでは要素を構成する各ノードの
座標です。これを見るにBマトリクスは
の関数となります。
に何を入れるのかはさておきソルバークラスにBマトリクス計算用の関数を作成しておきます。
class Solver: """ 有限要素法の計算を行うクラス Attributes: mesh (Mesh): Meshクラス material (Material): Materialクラス D (array 3x3): Dマトリクス """ def __init__(self, mesh, material): self.mesh = mesh self.material = material self.D = self.calc_D() def run(self): self.f, self.bc_dist_dict = self.read_boundary_cond() self.K = self.calc_K() self.d = self.solve() def calc_D(self): """ Dマトリクスを算出。 Returns: D (array 3x3): Dマトリクス """ D = np.array([[1, self.material.nu, 0], [self.material.nu, 1, 0], [0, 0, (1 - self.material.nu) / 2]]) \ * self.material.E / (1 - self.material.nu ** 2) return D def calc_B(self, elem, xi, eta): """ Bマトリクス及びヤコビアンを算出。 Args: elem (Element): 対象のElementクラス xi, eta (float): Bマトリクス算出用の座標 Returns: B (array 3x8): Bマトリクス J (array 2x2): ヤコビアン """ dndxi = np.array([- 1 + eta, 1 - eta, 1 + eta, - 1 - eta]) / 4 dndeta = np.array([- 1 + xi, - 1 - xi, 1 + xi, 1 - xi]) / 4 x = elem.xy[:, 0] y = elem.xy[:, 1] dxdxi = np.dot(dndxi, x) dydxi = np.dot(dndxi, y) dxdeta = np.dot(dndeta, x) dydeta = np.dot(dndeta, y) J = np.array([[dxdxi, dydxi], [dxdeta, dydeta]]) B = np.zeros((3, 8)) for i in range(4): Bi = np.dot(np.linalg.inv(J), np.array([dndxi[i], dndeta[i]]).T) B[0, 2 * i] = Bi[0] B[2, 2 * i + 1] = Bi[0] B[2, 2 * i] = Bi[1] B[1, 2 * i + 1] = Bi[1] return B, J
ちゃんと計算できているか確認しておきましょう。
MESH_SIZE = 5. LENGTH = 20. HEIGHT = 10. THICKNESS = 1. E = 208000 # ヤング率 NU = 0.27 # ポアソン比 nodes_dict, elems = create_mesh(LENGTH, HEIGHT, MESH_SIZE, THICKNESS) material = Material(E, NU) mesh = Mesh(0, nodes_dict, elems, material) solver = Solver(mesh, material) B, J = solver.calc_B(solver.mesh.elements[0], -1, -1) # xi, etaはとりあえず-1で算出 print(B) print(J)
[[-0.2 0. 0.2 0. 0. 0. 0. 0. ] [ 0. -0.2 0. 0. 0. 0. 0. 0.2] [-0.2 -0.2 0. 0.2 0. 0. 0.2 0. ]] [[2.5 0. ] [0. 2.5]]
まぁよさそうですね。
参考文献
リンク
リンク
