Code_Aster構造解析

カテゴリ: オープンソースCAE | 2026-10-02 改稿
CAE visualization for code aster structural theory - technical simulation diagram
Code_Aster構造解析

理論:3種類の解析と照合に使う解

概要

🙋

Code_Asterの使い方は一通り分かりました。構造解析では、静解析のほかに何ができるのですか?

🎓

よく使うのは、静解析・固有値解析・線形座屈の3つだ。静解析(MECA_STATIQUE)は荷重に対する変形と応力を求める。固有値解析(CALC_MODES)は、構造がどの振動数で揺れやすいかを求める。線形座屈も同じCALC_MODESで解くが、質量の代わりに「応力がかかった状態の剛性の変化(幾何剛性)」を使い、荷重を何倍にすると構造が座屈するかを求める。3つとも同じメッシュとモデルの指定(AFFE_MODELE)を使い回せるのがCode_Asterの便利なところだ。ただし、それぞれ材料に必要な値と、間違えやすい点が違う。この記事では、理論解のある3つの例題を通して、その違いを確かめていく。

3種類の解析の違い

解析主なコマンド材料に必要な値得られるもの
静解析MECA_STATIQUEE, NU変位・応力
固有値解析ASSEMBLAGE, CALC_MODESE, NU, RHO固有振動数・モード形状
線形座屈MECA_STATIQUE → RIGI_GEOM → CALC_MODES (MODE_FLAMB)E, NU基準荷重に対する座屈倍率

照合に使う理論解

穴あき板の引張(有限幅の板の応力集中係数、Heywoodの近似式):

$$ \sigma_{net} = \frac{P}{(W-d)\,t},\qquad K_{t,net} = 2 + \left(1 - \frac{d}{W}\right)^3,\qquad \sigma_{max} = K_{t,net}\,\sigma_{net} $$

片持ち梁の曲げ振動(オイラー・ベルヌーイ梁):

$$ f_n = \frac{(\beta_n L)^2}{2\pi L^2}\sqrt{\frac{E I}{\rho A}},\qquad \beta_1 L = 1.875,\ \beta_2 L = 4.694,\ \beta_3 L = 7.855 $$

柱の座屈(オイラーの式、K は有効長さ係数):

$$ P_{cr} = \frac{\pi^2 E I}{(K L)^2},\qquad K = 2\ \text{(fixed-free)},\quad K = 1\ \text{(pinned-pinned)} $$

P は荷重、W・d・t は板の幅・穴の直径・板厚、E はヤング率、I は断面二次モーメント、ρ は密度、A は断面積、L は長さ。

Coffee Break よもやま話

穴のまわりの応力が3倍になると示した論文

無限に広い板に小さな円孔があると、穴の縁の応力が遠くの応力の3倍になることは、ドイツの技術者キルシュが1898年に弾性論から導いた。この「3倍」は、有限要素法が生まれるずっと前から設計の常識となり、今では有限要素ソルバーの検証例題の定番になっている。板の幅が有限になると、穴を除いた断面の応力を基準にした係数は3より小さくなり、穴が板幅の2割なら約2.51になる。

計算例

単位はすべてN・mm・MPa・t(トン)で揃える。鋼はE = 210,000 MPa、ポアソン比0.3、密度7.85×10⁻⁹ t/mm³とする。

例1:穴あき板の引張(静解析)

幅100 mm・板厚10 mmの板の中央に直径20 mmの穴があり、20 kNで引っ張る。

量値
遠方の応力 P/(W t)20.0 MPa
穴の断面の平均応力 P/((W−d) t)25.0 MPa
応力集中係数 Kt,net2.512
穴の縁の最大応力62.8 MPa
参考:無限板(キルシュ)の3倍60.0 MPa

Code_Asterでは、対称性を使って1/4のモデルにし、対称面にDDL_IMPOで法線方向の変位を0に、端面にFORCE_FACEまたはPRES_REPで荷重を与える。応力はCALC_CHAMPで節点の値(SIGM_NOEU)を求め、穴の縁の値を読む。要素を細かくしていき、62.8 MPaに近づいて変わらなくなれば検証は完了だ。

幅50 mm・高さ100 mm・長さ1,000 mmの鋼の片持ち梁。曲がる向きが2つあるので、それぞれの曲げの固有振動数を求める。

曲がる向き1次2次3次
薄い方向(幅50 mmの向き)41.8 Hz261.8 Hz733.1 Hz
厚い方向(高さ100 mmの向き)83.6 Hz523.6 Hz1,466.1 Hz

CALC_MODESが小さい順に出す結果は、41.8 Hz(薄い方向の1次)、83.6 Hz(厚い方向の1次)、261.8 Hz、…と2つの向きが混ざって並ぶ。3次元の要素ではねじりや軸方向のモードも加わる。モード形状を見て、どの向きの何次かを確かめてから理論値と比べる。

🙋

固有値解析をしたら、1次が0.0000418 Hzと出ました。静解析のたわみは合っていたのに、なぜですか?

🎓

密度の単位を疑う。N・mm・MPaの単位系では、力=質量×加速度の関係から質量の単位がトンになる。だから鋼の密度は7850(kg/m³の数値)ではなく、7.85×10⁻⁹(t/mm³)で入れなければならない。7850のまま入れると質量が10¹²倍になり、振動数はその平方根の逆数、つまり100万分の1になる。41.8 Hzが0.0000418 Hzになったのはこのためだ。静解析は密度を使わないので、この誤りは静解析では見えない。固有値解析を始めるときは、まず片持ち梁の1次が41.8 Hzになるかを確かめると、単位の誤りをすぐに見つけられる。

例3:柱の座屈荷重(線形座屈)

20 mm角・長さ1,000 mmの鋼の柱。下端を固定し、上端に基準荷重1 kNの圧縮をかける。

量値
断面二次モーメント13,333 mm⁴
座屈荷重(下端固定・上端自由)6,909 N
Code_Asterの固有値(座屈倍率)約6.91
座屈時の圧縮応力17.3 MPa
参考:両端ピンの座屈荷重27,635 N
参考:降伏荷重(235 MPa)94,000 N

座屈時の応力17.3 MPaは降伏応力235 MPaよりずっと小さいので、この柱は降伏する前に弾性座屈で壊れる。Code_Asterでは、1 kNで静解析(MECA_STATIQUE)を行い、その応力からCALC_MATR_ELEMのOPTION='RIGI_GEOM'で幾何剛性を作り、CALC_MODESのTYPE_RESU='MODE_FLAMB'で解く。得られた固有値(CHAR_CRIT)に基準荷重を掛けたものが座屈荷重だ。

🙋

座屈の固有値がマイナスの値で出てきました。計算が失敗したのですか?

🎓

失敗ではなく、荷重の向きが逆という意味だ。座屈の固有値は「基準荷重を何倍すると座屈するか」なので、マイナスの値は「荷重を反対向きにすると、その倍率で座屈する」ことを示す。圧縮のつもりが引張になっていないか、荷重の符号と向きを確かめる。また、線形座屈は形の初期不整や材料の降伏を考えない理想的な値なので、実際の設計では初期たわみを入れた非線形解析(STAT_NON_LINE)で、座屈荷重がどれだけ下がるかも確かめる。特に座屈応力が降伏応力に近い太く短い柱では、線形座屈の値をそのまま使うと危険側になる。

実務での進め方

  1. 単位系を決め、入力ファイルの冒頭にコメントで書く(N・mm・MPa・t など)。
  2. 静解析で、たわみや反力が手計算と合うことを確かめる。反力の合計が荷重と一致するかも見る。
  3. 固有値解析では密度を加え、片持ち梁などで振動数の桁を確かめてから本番のモデルに進む。
  4. 座屈では基準荷重を単位荷重にし、固有値をそのまま倍率として読めるようにする。
  5. 要素を細かくしても結果が変わらなくなることを確かめ、同じ機能を使う公式の検証例題も参照する。
Coffee Break よもやま話

「1次モードが思った向きではない」

例2の梁は、荷重をかける向き(高さ100 mmの向き)だけを考えていると、1次は83.6 Hzだと思い込みやすい。しかし実際の1次は、薄い方向に揺れる41.8 Hzだ。振動を避けたい設備の回転数が45 Hzだとすると、厚い方向だけを見れば「1次の83.6 Hzより十分低い」と判断してしまうが、実際は41.8 Hzの近くで揺れが大きくなる。固有値解析は、想定していない向きのモードを見つけるためにこそ行う。結果は振動数の数字だけでなく、必ずモード形状と一緒に確認したい。

よくある間違い

間違いと対策

間違い起きること対策
mm単位で密度を7850振動数が100万分の17.85×10⁻⁹ t/mm³
RHOの入れ忘れ質量行列が組めずエラーELASにRHOを追加
モード形状を見ない想定外の向きの1次を見落とす形状と振動数を対で確認
座屈の荷重の向きが逆固有値がマイナス圧縮になっているか確認
線形座屈の値をそのまま使う危険側の評価初期不整つき非線形解析
🙋

関連する内容も知りたいです。

この記事の評価
ご回答ありがとうございます!
参考に
なった
もっと
詳しく
誤りを
報告
参考になった
0
もっと詳しく
0
誤りを報告
0