2015/4/8 版有限要素法による応力解析１．有限要素法の概要代表的な構造解析手法である有限要素法は、構造物を仮想的な有限な大きさの要素に分割し、要素の集合体としてモデル化し解析する方法である．単純な形状を有する要素の特性は、単純な代数方程式により表現されるので、それらを重ね合わせて構造物の全体的な力学的挙動を近似することができる．要素の頂点には節点が配置される．解析モデルにおいて、構造物を構成する材料特性は要素の特性として考慮され、構造物に加わる荷重は節点に対して与えられる．重ねあわせによって得られたマトリクス方程式は、コンピュータを利用した数値計算によって解くことができる．静解析、固有振動数や座屈荷重を求めるための固有値解析、構造物の過渡的な応答を求めるための動的解析、大変形を扱う幾何学的非線形解析、弾塑性材料を扱うための材料非線形解析などが実施可能である．以下、変位法による線形弾性静解析を例にとり、有限要素解析の流れを示す． 1.1 有限要素解析の流れ ① 要素分割構造解析で用いられる有限要素は、図１に示すように連続体要素（continuum element）と構造要素（structural element）に分類される[1]．前者は、節点において連続体力学で用いる変位ベクトルをそのまま並進自由度として用いる．これに対して、後者は構造上の特徴から剛体変形を仮定する箇所（はりやシェルの断面）を有する要素で、並進自由度に加えて回転量も自由度として有している．図１に示すようにさまざまな形状を有する有限要素が開発され有限要素解析に用いられている [1][2][3][4][5]．原理的には三次元連続体要素を用いればあらゆるものをモデル化することができるが、骨組構造、薄肉構造をモデル化するためには、はり要素やシェル要素を用いた方が効率的である．また平面問題、軸対称問題に対しては二次元要素が用いられる．なお、どのような要素をどの程度用いるかについては、解析者の経験や、ガイドラインに基づく解析者の工学的判断にゆだねられている．図１さまざまな有限要素（連続体要素と構造要素） 1 ② 要素特性の計算個々の要素の挙動を表すマトリクス方程式は、直接法またはエネルギー原理を用いて導出することができる．直接法は、要素特性を力の釣り合い式を用いて求める方法であり、トラス、ラーメン、はりなどの単純な構造要素に対しての導出は容易であるが、複雑な問題や非線形問題に対しては適用困難となる．エネルギー原理が汎用的な導出手順を与える. ③ 方程式の組み立て個々の要素が共有する節点で結合していることを利用して、要素毎に得られる要素剛性方程式から全体剛性方程式を作成する．変位法による線形静弾性問題の解析においては、全体剛性方程式は、正定値対称の定数係数を有するマトリクス方程式となる．一般に節点に結合する要素数は高々数個であるので、係数マトリクスの全行列成分における非零項の割合は一般に小さい．このような行列は、疎行列またはスパース行列と呼ばれる． ④ 境界条件の付与境界条件として変位の拘束条件を付与する場合、拘束する節点の自由度に関連する係数マトリクスおよび右辺ベクトルを修正する．また、荷重は節点自由度に対する等価節点力に置き換えたうえで、対応する節点の自由度に関連する右辺ベクトルを修正する． ⑤ 方程式の求解節点変位についての連立一次方程式の求解方法としては、記憶容量および計算速度などの計算効率の観点から係数マトリクスのスパース性を利用した方法が用いられる． ⑥ 要素内のひずみと応力得られた節点変位を用いて、要素ごとに要素内のひずみ、さらには応力が計算される．変位法による定式化を用いた場合、一般に要素境界でひずみおよび応力の連続性は保証されない．したがって通常は要素剛性マトリクスを算出するために用いられるガウス積分における積分評価点において応力、ひずみを評価する． 1.2 FEM による微小変形弾性解析の定式化三次元固体について成立する仮想仕事の原理に基づき、微小変形弾性問題についての有限要素法による定式化を示す． ① 仮想仕事の原理微小変形線形弾性問題についての仮想仕事の原理は次式のように表わされる． ∫ σ dε ij V ui = ui ij dV = ∫ tidui dΓ + ∫ bidui dV Γt (1.1) V on Γu (1.2) 2 ここにσ ij , ε ij , u i は、応力テンソル、ひずみテンソル、変位ベクトルであり、 t i , bi は変位境界上で既定される表面力ベクトル（圧力など）、領域内部で既定される物体力ベクトル（自重、遠心力など）である．また、δu i は領域内部での任意の仮想変位であり、式(1.2)により、基本境界上でその値は零である．V は解析領域、Γu は変位境界、Γt は力学的境界を表す． ② ひずみ・変位関係式微小ひずみテンソルと変位ベクトルとの関係は次式で与えられる． ε ij = 1 (u i, j + u j ,i ) 2 (2) ③ 応力ひずみ関係式（構成方程式）線形弾性材料においては、応力とひずみとの関係式は４階の弾性テンソル C ijkl を用いて、次式で与えられる． σ ij = C ijkl ε kl (3) ④ 離散化解析領域を有限要素分割し、仮想仕事の原理式(1.1)を有限要素ごとに離散化して次式を得る． ∑ ∫ σ de dV = ∑ ∫ t du dS + ∑ ∫ b du dV ij e ij i e Γe t Ve i i i (4) e Ve ここに添え字 e は、要素毎に定義される諸量であることを示す．要素を構成する節点 I における変位量 uiI を用いて、要素内の変位場 ui を、内挿関数φI を用いて次式のように近似する． u i = ∑ φ I u iI (e) (5) I 式(5)は、マトリクス形式で次式のように表わされる． {u} = [ N ( e ) ]{U } (6) このときひずみベクトルは、式(2)に基づき次式のように与えられる． {e } = [ B ( e ) ]{U } (7) また、式(3)は、マトリクス形式で次式のように書くことができる． {σ } = [ D ( ε ) ]{ε } (8) 式(6) (7) (8)を式(4)に代入した式が、任意の仮想変位δ {U } について成立することにより次式を得る． 3    ∑ [ K ( e ) ] {U } = ∑{F ( e ) } e  e  (9) ここに [ K ( e ) ] = ∫ [ B ( e ) ]T [ D ( e ) ][ B ( e ) ]dV V (10.1) e {F ( e ) } = ∫ [N Γt e ] {t }dΓ + ∫ [ N ( e ) ]T {b }dΓ (e) T Γt (10.2) e 一方、変位境界条件式(1.2)は、変位境界上の節点自由度について次式のように表される {U } = {U } on Γu (11) 式(9)を拘束条件式(11)のもとで解くことにより、節点変位{U}を決定することができる. 2. ABAQUS/CAE によるＦＥＭ解析 2.1 有限要素解析の一般的手順 ① 解析領域及び解析条件の決定実際に解析すべき問題の多くは何らかの幾何学的及び力学的対称性があり、そのような問題において全領域を解析するのは計算コストの面で得策ではない．対称性を利用することは解析解の一部を利用することと等価であり、精度の向上に役立つ．また応力集中問題のように領域の一部分の応力分布を解析すればよい場合には適当な解析領域を決定し、領域の境界での既知応力またはその近似値を外力として解析することも可能である．適切な解析領域を決定することにより解析コストを軽減でき、また同じ計算コストでより高精度な解析を行うことができる． ② 有限要素モデル(入力データ)の作成解析者が決定した解析領域に対して有限要素メッシュを生成し、それに基づいて有限要素解析のための入力データを作成する．入力データには節点座標リスト、要素節点番号リスト、材料特性、節点に与えられる境界(拘束)条件等が含まれる．これらの処理を「前処理」あるいは「プリプロセッシング」と呼ぶ． ③ 有限要素解析の実行解析プログラムはプリプロセッシングにより作成された入力データを読み込み剛性方程式(剛性マトリックス、荷重ベクトル)の作成と解析を行う．剛性方程式は有限要素メッシュが密になると大規模な連立一次方程式となり、解析プログラム及びコンピュータの性能が要求される．剛性方程式の解析結果を用いて各要素の応力、ひずみ等が計算される．このような処理を「メイン処理」と呼ぶことがある． 4 ④ 解析結果の表示有限要素解析の結果は節点変位、要素応力及びひずみ等であり、それらは膨大なデータ量となり、解析結果を適切に解釈するためには何らかの可視化表示が不可欠である．それらの処理を「後処理」あるいは「ポストプロセッシング」と呼ぶ． 2.2 利用プログラム本実験では ABAQUS/CAE を用いた有限要素解析を行う．ABAQUS/CAE は、SIMULIA 社が開発した汎用 FEM 解析システムである．世界中の多くの製造業、研究開発機関で利用されている．ただし、本実験では ABAQUS/CAE の Student Edition を利用するものとする．ABAQUS/CAE の Student Edition は、ABAQUS が有するほぼ同じ機能を有しているが、解析できる FEM の節点数は 1,000 節点に制限されている．ABAQUS/CAE は、前処理、メイン処理、後処理を行うプログラムが統合化された計算環境を提供する．GUI（Graphical User Interface)を有しているので、対話型処理により比較的容易に FEM 解析を実行することができる． 2.3 例題による解析手順の説明有限要素解析を行うための具体的な手順を例題を用いて説明する．ここで取り上げる例題は図 2 に示すような一様引張り荷重を受ける円孔つき平板問題であり、その形状は荷重試験においてひずみ計測する試験片と同じである．ここでは、ヤング率 72000N/mm2 ポアソン比 0.3 とし、荷重 100kgf (980N) で一様に引張るものとする． 60 mm 200 mm D 図 2 円孔つき平板 ① 解析領域及び解析モデルの決定この問題は幾何学的及び力学的対称性から 1/4 領域をモデル化すればよい． 5 ② ABAQUS/CAE による手順 ABAQUS/CAE の起動「モデルデータベースの作成, Standard/Explicit モデル」を押す．入力画面 ABAQUS（や多くの CAE ソフトウエア）では、回転と角度を除いて、決められた単位系はない．回転自由度はラディアンで表現され、それ以外のすべての角度の単位は「度」で表現される．選択した単位系で矛盾がないことの保証は、ユーザの責任である。 SI 単位：質量[M]：kg、長さ[L]：m 、時間[T]：sec 、温度[θ]：K→応力[Pa] 推奨単位：質量[M]：ton、長さ[L]：mm 、時間[T]：sec 、温度[θ]：K →応力[MPa] 6 ワークディレクトリの設定ファイルメニュー：ワークディレクトリの設定自分で管理できるフォルダー（自分で読み書き可能なフォルダー）を指定する．パートの作成モデルツリー：パートをダブルクリックパートの作成ダイアログ：名前→plate_with_hole ←任意モデリング空間→2 次元平面タイプ→変形体ベースフィーチャー →シェルサイズ→200 程度 7 スケッチ解析領域の 1/4 部分（右上半分）を直線で結び矩形形状を作成する．スケッチャーのメニューの矩形形状作成機能を使い適当な大きさで描く．さらに円孔形状作成機能を使い円孔形状を重ね書きする．（ここでは直径 20mm（半径 10mm）とする．）トリム機能を用いて余計な線を削除する． 8 解析対象領域が完成する．材料特性の作成モデルツリー：材料特性をダブルクリック材料特性の編集ダイアログ：名前→Material-1（デフォルト）タブ機械的を選択材料挙動→弾性（デフォルト）タイプ→等方性（デフォルト）ヤング率→ 72000MPa（入力）ポアソン比→0.3（入力） →OK 9 要素特性の作成モデルツリー：要素特性をダブルクリック要素特性の作成ダイアログ：名前→Section-1（デフォルト）カテゴリー →ソリッド（デフォルト）タイプ→均質（デフォルト） →続ける要素特性の編集ダイアログ：材料特性→Material-1 平面応力／ひずみ要素厚 →1 mm（設定） →OK 10 要素特性の割り当てモデルツリー：パート、plate_with_hole 要素特性の割り当て →クリック作成した形状をハイライトさせ →完了要素特性の割り当ての編集ダイアログ要素特性→Section-1 →OK 11 モデルのアセンブリモデルツリー：アセンブリインスタンス →ダブルクリックインスタンスの作成のダイアログインディペンデントをクリック →OK ステップの作成モデルツリー→ステップ（ダブルクリック）ステップの作成ダイアログ線形摂動→続けるステップの編集→OK 12 境界条件の設定１(対称線）モデルツリー→境界条件（ダブルクリック）境界条件作成ダイアログ対称／反対称／完全固定→続ける 13 対称線を選択→完了境界条件の編集ダイアログ YSYMM を選択→OK 境界条件の設定２(対称線）モデルツリー→境界条件（ダブルクリック）境界条件作成ダイアログ対称／反対称／完全固定→続ける 14 対称線を選択→完了境界条件の編集ダイアログ XSYMM を選択→OK 荷重の作成モデルツリー→Step-1 →荷重条件（ダブルクリック）荷重の作成ダイアログ圧力→続ける 15 荷重のためのサーフェスを指定 →完了荷重の編集ダイアログ →Y 方向（第２方向）の負の方向に-16.3333 N/mm2 を設定メッシュの作成モデルツリー →インスタンス →plate_with_hole ←パートの定義名 →メッシュ（ダブルクリック） 16 メッシュメニュー →コントロール →メッシュコントロールダイアログ →OK メッシュメニュー →要素タイプダイアログ要素ライブラリ→Standard（デフォルト）ジオメトリ次数→線形（デフォルト）要素コントロール→低減積分のチェックをはずす →OK 17 シードメニュー →インスタンスメッシュメニュー →インスタンス →インスタンスをメッシュ分割する 18 注意 ABAQUS/CAE の SE 版は最大 1000 節点までしかモデルを作成できない 19 解析の実施モデルツリー→ 解析→ジョブ（ダブルクリック）ジョブの作成ダイアログ →続けるジョブの編集ダイアログ（デフォルト） →OK 20 モデルツリー→ 解析→ジョブ（右クリック）マネージャー →ジョブ投入正常終了後：応力解析結果の評価引張り方向（鉛直方向）の垂直ひずみ E22, S22 分布 21 コンターレベルの数値の大きさを拡大する方法ビューポート →ビューポート注釈オプション →凡例 →フォントの設定解析結果の出力解析結果のプリントアウトも可能であるが、報告書に載せる図面を得ることが目的であるのであれば、 TIFF ファイルなど電子ファイルに出力可能． ABAQUS/CAE の終了ファイルメニュー →保存終了 22 さらに詳細なメッシュを作成するためのテクニック（１）パーティション（分割）モデルツリーでパートを選択フェースのパーティション:スケッチを選択スケッチャーが起動スケッチャーで、パーティション（分割）ラインを描画 23 パーティション完了（２）モデルツリーでインスタンスを選択 24 インスタンスを移動応力が集中する円孔の縁を XY 座標の原点に移動する． 25 （３）境界条件、荷重条件の修正領域を分割（パーティション）することによって、もともとの形状で設定した境界条件や荷重条件が不適切になってしまう場合がある．そのために一旦境界条件や荷重条件を削除して再度設定する必要がある．設定方法は前述のとおり．（pp13-16 参照）（４）メッシュ再分割モデルツリーでメッシュを選択 26 メッシュメニュー →コントロール →３つの部分領域を選択し、メッシュコントロールダイアログで構造メッシュを選択シードメニューを用いて、分割された領域のエッジ（領域境界）の分割数を指定し、四辺形要素で分割． 27 （5）再解析引張り方向（鉛直方向）の垂直ひずみ E22, S22 分布文献 [1] 久田・野口 (1995)、非線形有限要素法、丸善. [2] 鷲津ほか(1980)、有限要素法ハンドブック, 培風館. [3] Hughes TJR (1987), The Finite Element Method Linear Static and Dynamic Finite Element Analysis, Dover . [4] Belytschko T, Liu, WK, Moran B (2000) Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, Ltd. [5] Bathe KJ (1982) Finite Element Procedures in Engineering Analysis. Prentice-Hall. 28 ＦＡＱ＃質問参考 0 解析の全体の流れは？ A0 1 スケッチツールの機能は？ A1 2 モデルの表示方法は？ A2 3 出力対象とする変数の選択方法は？ A3 4 出力変数の最大、最小の位置を表示するには？ A4 5 解析領域の表示範囲を変更するには？ 6 表示されている物理量の数値を表示するには？ A6 7 画像データへ出力するには？ A7 8 メッシュ数や要素数の情報を知るには？ A8 9 コンター表示の設定を変更するには？ A9 29 A5 ABAQUS/CAEの起動プリ処理パートの作成メッシュの作成材料特性の作成解析の実施要素特性の作成応力解析の評価メイン処理ポスト処理モデルのアセンブリ解析結果の出力要素特性の割り当て ABAQUS/CAEの終了境界条件の設定ステップの作成荷重の作成 A0 ABAQUS/CAE における解析の流れ直線の作成矩形円拘束寸法の追加削除 A1 スケッチツール 30 スケーリング拡大ズーム縮小ビュー表示を標準設定に戻したい場合ビューメニューからビューツールボックスを選択移動ビュー回転ビュー A2 モデルの表示方法結果メニュー → フィールド出力出力する変数を指定 A３出力変数の設定 31 オプションメニュー → コンター →範囲タブ場所の表示をチェック A４出力変数の最大値・最小値の場所の表示ツールメニュー → 表示グループ →マネージャーまたは作成 →選択の編集機能を使う A5 解析領域の表示範囲の設定 32 A6 A7 クエリ解析結果の画像ファイル出力 33 A8 クエリのメッシュを選択してダブルクリック A9 コンターオプション 34