線形応力解析 ( Stress Analysis )

線形応力解析の例を以下に示す．

• 線形問題なので、外力がかかった際の変位を、構造物の変形を考えずにどの程度曲がるか( 十分小さくないと線形近似が破綻する )を計算している．

• 非線形問題では、変形により負荷のかかり方も異なってくることが考慮に入る．

問題設定

問題としては、以下に示すような片持ち梁にかかる応力変形を計算する．

• 手前側の側面を固定とし、奥側の面を -z方向に 1 (tF) で引っ張ることを考える (境界条件)．

• 構造物は、ヤング率 200 (GPa)、ポアソン比 0.3 の材質( 構造鉄を意識 )からなるとする(材料条件)．

片持ち梁を生成する gmsh-API python プログラム

gmsh-API pythonを利用した片持ち梁モデルの生成プログラム． 自分で作成した箱生成用関数を内部で利用している． 中身については、 gmsh もしくは github で解説されているはず．たぶん．

モデル生成

$ python beam.py
$ ElmerGrid 14 2 model.msh

変換中の数字は、( 14 : gmshの.mshファイル、 2 : ElmerMeshファイル4つを含んだディレクトリ )である．同ディレクトリに model というディレクトリができて、その中に model.header / model.element / model.node / model.boundary が生成されているはず．

片持ち梁生成用プログラム( beam.py )

import numpy as np
import sys
import gmsh

gmshlib = "/Users/kent/gmsh/pygmshLibrary/"
sys.path.append( gmshlib )
import generateRectangularBox as box

# ------------------------------------------------- #
# --- [1] initialization of the gmsh            --- #
# ------------------------------------------------- #
gmsh.initialize()
gmsh.option.setNumber( "General.Terminal", 1 )
gmsh.model.add( "model" )

# ------------------------------------------------- #
# --- [2] initialize settings                   --- #
# ------------------------------------------------- #
ptsDim , lineDim , surfDim , voluDim  =  0,  1,  2,  3
pts    , line    , surf    , volu     = {}, {}, {}, {}
ptsPhys, linePhys, surfPhys, voluPhys = {}, {}, {}, {}
lc                                    = 0.01

# ------------------------------------------------- #
# --- [3] Modeling                              --- #
# ------------------------------------------------- #
v0  = [    0.0,   0.0,   0.0   ]
v1  = [    1.0,   0.0,   0.0   ]
v2  = [    0.0,   0.1,   0.0   ]
v3  = [    0.0,   0.0,   0.030 ]
ret = box.generateRectangularBox( origin=v0, v1=v1, v2=v2, v3=v3, lc=lc )
lft = ret["surfs"]["side2"]
rgt = ret["surfs"]["side4"]
bot = ret["surfs"]["bottom"]
box = ret["volus"]["box"]

# ------------------------------------------------- #
# --- [4] Physical Grouping                     --- #
# ------------------------------------------------- #
gmsh.model.occ.synchronize()
surfPhys["fixed"] = gmsh.model.addPhysicalGroup( surfDim, [ lft ], tag=201 )
surfPhys["free" ] = gmsh.model.addPhysicalGroup( surfDim, [ rgt ], tag=202 )
surfPhys["bot"  ] = gmsh.model.addPhysicalGroup( surfDim, [ bot ], tag=203 )
voluPhys["box"  ] = gmsh.model.addPhysicalGroup( voluDim, [ box ], tag=301 )

# ------------------------------------------------- #
# --- [5] post process                          --- #
# ------------------------------------------------- #
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)
gmsh.write( "model.geo_unrolled" )
gmsh.write( "model.msh" )
gmsh.finalize()

• 201 / 202 : 手前、及び、奥側のSurface

• 301 : 梁のVolume

線形応力解析用のElmer入力ファイル ( .sif ファイル)

以下にElmer入力ファイルのサンプルを示す．

Elmer 線形応力解析 ( elastic_linear.sif )

CHECK KEYWORDS "Warn"

Header
  Mesh DB "." "model"
  Include Path ""
  Results Directory ""
End

Simulation
  coordinate system = "Cartesian 3D"
  Coordinate Mapping(3) = 1 2 3

  Simulation Type = "Steady State"
  Steady State Max Iterations = 20
  
  Solver Input File = "elastic_linear.sif"
  Output File = "elastic_linear.dat"
  Post File = "elastic_linear.vtu"
End


Constants
  Gravity(4) = 0 0 -1 9.82
End


Body 1
  Target Bodies(1) = 301
  Name = "beam"

  Equation = 1
  Material = 1
  Body Force = 1
End


Solver 1
  Exec Solver   = "Always"
  Equation      = "Stress Analysis"
  Variable      = "Displacement"
  Variable Dofs = 3
  Linear System Solver = "Iterative"
  Linear System Iterative Method = "BiCGStab"
  Linear System Max Iterations = 10000
  Linear System Convergence Tolerance = 1.0e-5
  Linear System Preconditioning = "ILU0"
  Steady State Convergence Tolerance = 1.0e-5
End


Equation 1
  Name = "StressAnalysis"
  Stress Analysis   = True
  Calculate Stresses = Logical False
End


Material 1
  Name           = "Iron_SS400"
  Youngs Modulus = 200.0e9
  Poisson Ratio  = 0.3
End


Body Force 1
  Pressure = 5.0e4
End


Boundary Condition 1
  Name = "Fixed_Surface"
  Target Boundaries(1) = 202

  Displacement 1 = 0
  Displacement 2 = 0
  Displacement 3 = 0
End


Boundary Condition 2
  Name = "Load_Surface"
  Target Boundaries(1) = 201
  Force 3 = -1.00e4

End

• 最も重要な箇所は Equation Section の "Stress Analysis = True" の指定

線形応力解析の結果

解析実行と結果の表示は以下の通り．

解析実行と結果の表示

$ ElmerSolver elastic_linear.sif
$ paraview model/elastic_linear_t0001.vtu

結果は次のようになった．

1 (m) の梁に対して、 0.17 (mm) 程度曲がっているらしい．x正方向に進むほどより大きな変位がある．カラーバー上では正値となっているのは、絶対値をとっているからであって、本来はz方向負側に曲がるはずである．曲がっているのか、よく見えないので、加える外力を増やしてみる．試しに 10 倍の 10 (tF) をかけてみる．結果は次のようになった．

カラーバーが示す変位量が10倍の 1.7 (mm) になった．ちょっとぐらい曲がっているのだろうか．変位しているのかが、まだはっきりとはわからないので、さらに10倍の 100 (tF) をかけてみる．結果は次のようになった．

また10倍の変位量となり、17 (mm) 曲がることがわかった． 根性がひねくれている人でも曲がっていることがはっきりわかるように軸をつけている． ここでは、 1 (m) の片持ち梁に 100 (tF) かけた構造鉄の変位量が 17 (mm) であるという答えを得たが、 17 (mm) は線形解析にとって無視できない変位量かもしれない．線形解析の適用範囲から逸脱し、この値が正しくない可能性があることを付記しておく．