【python】自作の有限要素法(FEM)プログラムの検証

はじめに

これまで全6回にわたって静弾性問題の有限要素解析の実装をゼロから実施してきました。
今回はそこで作ったプログラムの確からしさを検証したいと思います。

これまでの記事は以下です。

はじめに
対象の静弾性問題
自作プログラムでの解析
- 理論値との比較
終わりに
参考文献

対象の静弾性問題

今回対象とするのは以下のような片持ち梁のたわみの問題について考えます。

材料はS45Cとし、ヤング率205GPa, ポアソン比0.27とします。
理論値は以下のWebサービスを使って求めます(めんどくさかったので…)。

各種パラメータを入れると結果は以下の通りになりました。

合計のたわみ量としては0.0156mmという結果でした。

自作プログラムでの解析

上記の問題を自作のプログラムで解いてみます。各種パラメータとメッシュを定義すればすぐ解析可能です。
コードは以下の通りです。(クラスの定義は長くなるので省いています。)

MESH_SIZE = 1
LENGTH = 20.
HEIGHT = 10.
THICKNESS = 1.
E = 205000 # ヤング率
NU = 0.27 # ポアソン比

nodes_dict, elems = create_mesh(LENGTH, HEIGHT, MESH_SIZE, THICKNESS)
mesh = Mesh(0, nodes_dict, elems)
material = Material(E, NU)

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

for k, v in mesh.nodes.items():
    # x変位の境界条件を設定
    if v.x == 0.0:
        v.bc_dist[0] = 0.0
    # y変位の境界条件を設定
    if v.x == 0.0:
        v.bc_dist[1] = 0.0
    # 荷重の設定
    if v.x == LENGTH and v.y == HEIGHT:
        v.bc_force = [0.0, -100.0]

solver = Solver(mesh, material)
solver.run()
solver.plot_deform()
solver.plot_stress(3, 50)

# x = 20, y = 5のところの変位をたわみとする
for k,v in solver.mesh.nodes.items():
    if v.x == LENGTH and v.y == HEIGHT / 2:
        idx = [k * 2, k * 2 + 1]
print(solver.d[idx]) # x = 20, y = 5の変位

結果は以下のようになりました。

変形結果

von Mises応力分布

x=20, y=5における変形量

[-9.17314505e-05 -1.80764243e-02]

たわみ量でいうと約0.0181mmでした

理論値との比較

たわみ量の理論値は0.0156mm, 自作プログラムでの解析値は0.0181mmとなり若干ですが自作プログラムのほうが大きく出ました。
ただ、結果のオーダはばっちりあっているのでまともに解析できてそうな結果でした。

差がある理由ですが、有限要素プログラムの要素の選択にあります。
今回私が自作した有限要素プログラムは4角形1次要素を使っておりました。
1次要素は節点間の変位を1次関数で近似するため変位量が要素内で線形に変化するのですが、今回のような曲げのモードの場合は変位量が非線形に分布するためどうしても誤差が生じるということです。

対策としてはメッシュを細かく切って要素数を増やしたり、2次要素を使うということがあげられます。
2次要素は1次要素の節点の間にもう1つ節点を増やした要素で、4角形2次要素は以下のように定義されます。