Salome-Meca2019による練習3：書籍『例題で極める非線形有限要素法』より

Salome-MecaはDEXCS2014版を長らく使っています．最近，Salome-Mecaの2019 を使い始めました．そのための練習用に，以下の本に載っているMarcでの事例をSalome-Mecaで計算しようとしています．

Markdown記法に慣れていないため読み難いかもしれません．頑張ってMarkdown記法を覚えていきます．

例題で極める非線形有限要素法
CAEで正しい結果を導くための理論トレーニング
発行元：丸善出版　　発行年月日：2020年10月
https://www.maruzen-publishing.co.jp/item/?book_no=304000

【CS08-01】 矩形断面長柱の座屈

形状：矩形断面1mm×1mmで長さ20mmの部材．
境界条件：長柱の一端を固定し，他端に荷重[N]を掛けます．
縦弾性係数は10,000[MPa]，ポアソン比は0です．
$a=1$，$L=20$，$E=10,000$，断面二次モーメントは，$I=a^4/12$　です．

オイラーの座屈の公式より座屈過重は $Pcr=4π^2EI/L^2=82.16$ [MPa] です．
この公式はポアソン比が0という仮定のもとで導出されているため，CAE計算でもポアソン比を0にしました．

Salome-Mecaのメッシュ作成の設定

 ６面体要素
   3D Hexahedron(I.J,K)
   2D Quadrangle:Mapping
   1D Wire Discreation         Local Length 0.2 

Zの小さいほうを固定面，Zの大きいほうを荷重面としました．形状の初期不整を入れるためスパン中央でパーティションを入れてエッジに名前を付けました．このエッジについてY方向に荷重を加えて初期不整とします．参考にした書籍では，スパン中央が撓んだ正弦波の形に予めしておくことで初期不整を作り出しています．しかし,同じ形状をSalome-Mecaで作成する方法が分からなかったため上記の設定としました．

参考にした書籍では,スパン中央での正弦波形状の初期不整の量$\delta$について0.1mmと0.01mmの２通りの計算を行って，0.01mmでは座屈が発生せず0.1mmで座屈が発生したと記載されています．

両端固定でスパン中央に集中荷重を受けるオイラー梁の中央部の撓みは
$\delta=FL^3/192EI$
となります．したがって，$\delta=0.1$[mm] と $0.01$[mm] を生じさせる荷重はそれぞれ $F=2$ [N] と $0.2$ [N] です．

最大荷重の 100 [N] まで 100 [s] でランプ状に 100ステップで変化する条件で計算しました．荷重面の移動はZ軸方向のみ，かつ荷重面が平面状態を保つようにしました．
DEFI_LIST_REELコマンドについては， 0 [s] から 100 [s] まで100ステップという設定でも良さそうに思えるかもしれません．しかし，実際に行ってみると 0 [s] から 1 [s] までが 0.01 [s] 間隔となり 1 [s] から 100 [s] までが 1 [s] 間隔となってしまいました．

funcload = DEFI_FONCTION(
  NOM_PARA='INST', 
  VALE=(0.0, 0.0, 100.0, 1.0)
) 
listr = DEFI_LIST_REEL(
  DEBUT=0.0, 
  INTERVALLE=(_F(
      JUSQU_A=1.0, 
      NOMBRE=1
    ), _F(
      JUSQU_A=100.0, 
      NOMBRE=99
    ))
)

mecaFix = AFFE_CHAR_MECA(
  DDL_IMPO=(_F(
      DX=0.0, 
      DY=0.0, 
      DZ=0.0, 
      GROUP_MA=('fix', )
    ), _F(
      DX=0.0, 
      DY=0.0, 
      GROUP_MA=('load', )
    )), 
  LIAISON_UNIF=_F(
    DDL=('DZ', ), 
    GROUP_MA=('load', )
  ), 
  MODELE=model
) 
mecaLoad = AFFE_CHAR_MECA(
  FORCE_FACE=_F(
    FZ=-100.0, 
    GROUP_MA=('load', )
  ), 
  MODELE=model
) 
mecaShif = AFFE_CHAR_MECA(
  FORCE_ARETE=_F(
    FY=2.0, 
    GROUP_MA=('shift', )
  ), 
  MODELE=model
) 

さきほど記載したように，$\delta=0.1$[mm] と $\delta=0.01$[mm] を生じさせる荷重はそれぞれ $Fy=2$ [N] と $Fy=0.2$ [N] です．
参考書籍では，$\delta=0.1$[mm] のときに座屈が発生し，$\delta=0.01$[mm] のときは生じません．同様の結果を期待しましたが， $Fy=0.2$ [N] のときも座屈が発生しました． そこで，$Fy$ を１桁づつ小さくしていきましたが $Fy=0.0002$ [N] でも座屈が発生したため，ここで止めました．

横軸にZ方向の変位量，縦軸に荷重$P$としたグラフを下記に記載します．初期不整の荷重$Fy$が小さくなるつれて理論式に近づいていることが判ります．$Fy=0.0002$ [N] のときの折れ線の頂点の荷重の値は82 [N] です．計算での荷重増分ステップが1[N]であるので，理論式通りに再現していると言えるでしょう．

初期不整荷重の大きさによる座屈経路の違い




それぞれの変形状態を動画にしましたが，動画から座屈荷重を読み取るのは無理そうです．
それから，実際は座屈のあとに塑性変形状態になりますが，この計算では塑性を考慮していません．

初期不整荷重：$Fy=20$ [N]

初期不整荷重：$Fy=2$ [N]

初期不整荷重：$Fy=0.0002$ [N]