(57)【要約】      （修正有）【課題】極低サイクル疲労破壊をシミュレートできる有限要素解析装置、方法及びプログラムを提供する。【解決手段】有限要素解析は、有限要素モデルの剛性行列を作成し、剛性行列と前記解析条件とで規定される剛性方程式の解に基づいて、有限要素モデルの加工硬化特性を算出する特性計算を行う応力ひずみ算出手段と、有限要素モデルの内力と外力とが釣り合うまで、応力ひずみ算出手段に特性計算を繰り返し行わせる制御手段とを備える。応力ひずみ算出手段は、剛性方程式の解に基づく試行弾性歪みを用い、繰り返し負荷を受けている部分であって、破壊が発生していない部分の応力とひずみを算出する加工硬化モデルと、破壊が発生している部分の破壊進行度と応力を算出する材料損傷モデルとにより応力とひずみとを連成させた材料損傷方程式を作成し、材料損傷方程式の解を求め、求めた解に基づいてコンシステント接線剛性を作成し、剛性行列を作成する。【選択図】図４
2016004543.tif 000002
【特許請求の範囲】
【請求項１】
非線形有限要素解析を行うための有限要素解析装置であって、
解析対象の有限要素モデルと解析条件とを取得する取得手段と、
前記有限要素モデルの剛性行列を作成し、作成した剛性行列と前記解析条件とで規定される剛性方程式の解に基づいて、前記有限要素モデルの加工硬化特性を算出する特性計算を行う応力ひずみ算出手段と、
前記有限要素モデルの内力と外力とが釣り合うまで、前記応力ひずみ算出手段に前記特性計算を繰り返し行わせる制御手段とを備え、
前記応力ひずみ算出手段は、前回算出した前記剛性方程式の解に基づく試行弾性歪みを用い、繰り返し負荷を受けている部分であって、破壊が発生していない部分の応力とひずみを算出する加工硬化モデルと、破壊が発生している部分の破壊進行度と応力を算出する材料損傷モデルとにより応力とひずみとを連成させた材料損傷方程式を作成し、作成した材料損傷方程式の解を求め、求めた解に基づいてコンシステント接線剛性を求め、求めたコンシステント接線剛性から前記剛性行列を作成する
ことを特徴とする有限要素解析装置。
【請求項２】
前記応力ひずみ算出手段は、前記材料損傷方程式の解に基づいて、前記有限要素モデルの各要素の損傷値を求め、前記損傷値が所定の閾値を上回る要素は破壊されたと判断する
ことを特徴とする請求項１に記載の有限要素解析装置。
【請求項３】
前記加工硬化モデルは、分割背応力形Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｋ則による加工硬化モデルである
ことを特徴とする請求項１に記載の有限要素解析装置。
【請求項４】
前記コンシステント接線剛性は、完全陰解法により求める
ことを特徴とする請求項１ないし３のいずれか一項に記載の有限要素解析装置。
【請求項５】
前記制御手段は、前記有限要素モデルの内力と外力とが釣り合うと、前記解析条件のうちの強制変位又は外力を、所定の増分値だけ増加させて、前記応力ひずみ算出手段に前記特性計算を繰り返し行わせる
ことを特徴とする請求項１ないし４のいずれか一項に記載の有限要素解析装置。
【請求項６】
前記有限要素モデルを表示する表示部を、更に備え、
前記制御手段は、前記応力ひずみ算出手段が破壊されたと判断した要素を除いて、各要素をそれぞれの応力に応じた表示態様で、前記有限要素モデルを前記表示部に表示させる
ことを特徴とする請求項１ないし５のいずれか一項に記載の有限要素解析装置。
【請求項７】
前記剛性算出手段は、作成した材料損傷方程式の解を求める前に、前記有限要素モデルの各要素の応力に対して、共回転変換を行う
ことを特徴とする請求項１ないし６のいずれか一項に記載の有限要素解析装置。
【請求項８】
非線形有限要素解析を行うための有限要素解析装置で用いられる有限要素解析方法であって、
解析対象の有限要素モデルと解析条件とを取得する取得ステップと、
前記有限要素モデルの剛性行列を作成し、作成した剛性行列と前記解析条件とで規定される剛性方程式の解に基づいて、前記有限要素モデルの加工硬化特性を算出する特性計算を行う応力ひずみ算出ステップと、
前記有限要素モデルの内力と外力とが釣り合うまで、前記応力ひずみ算出ステップに前記特性計算を繰り返し行わせる制御ステップとを備え、
前記応力ひずみ算出ステップは、前回算出した前記剛性方程式の解に基づく試行弾性歪みを用い、繰り返し負荷を受けている部分であって、破壊が発生していない部分の応力とひずみを算出する加工硬化モデルと、破壊が発生している部分の破壊進行度と応力を算出する材料損傷モデルとにより応力とひずみとを連成させた材料損傷方程式を作成し、作成した材料損傷方程式の解を求め、求めた解に基づいてコンシステント接線剛性を求め、求めたコンシステント接線剛性から前記剛性行列を作成する
ことを特徴とする有限要素解析方法。
【請求項９】
非線形有限要素解析を行うための有限要素解析装置で用いられる有限要素解析プログラムであって、
解析対象の有限要素モデルと解析条件とを取得する取得手段と、
前記有限要素モデルの剛性行列を作成し、作成した剛性行列と前記解析条件とで規定される剛性方程式の解に基づいて、前記有限要素モデルの加工硬化特性を算出する特性計算を行う応力ひずみ算出手段と、
前記有限要素モデルの内力と外力とが釣り合うまで、前記応力ひずみ算出手段に前記特性計算を繰り返し行わせる制御手段として、コンピュータを機能させ、
前記応力ひずみ算出手段は、前回算出した前記剛性方程式の解に基づく試行弾性歪みを用い、繰り返し負荷を受けている部分であって、破壊が発生していない部分の応力とひずみを算出する加工硬化モデルと、破壊が発生している部分の破壊進行度と応力を算出する材料損傷モデルとにより応力とひずみとを連成させた材料損傷方程式を作成し、作成した材料損傷方程式の解を求め、求めた解に基づいてコンシステント接線剛性を求め、求めたコンシステント接線剛性から前記剛性行列を作成する
ことを特徴とする有限要素解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、加工（移動）硬化モデルと損傷モデルとを用いた有限要素解析方法に関する。
【背景技術】
【０００２】
近年，産業分野において機械部品の設計にＣＡＥ（Ｃｏｍｐｕｔｅｒ  Ａｉｄｅｄ  Ｅｎｇｉｎｅｅｒｉｎｇ）解析が行われている。この解析手法は、解析対象を要素に分割してモデル化し、物体に作用する外力により各要素に生ずるひずみや応力を逐次計算して、破壊に至る過程を計算する。
【０００３】
例えば、金属板から成形加工を経て製造させる金属部品のＣＡＥ解析用装置が、特許文献１に示されている。この装置は、加工部位の塑性ひずみをＣＡＥ解析の初期状態に反映させることで、成形加工時に生じる塑性変形による影響を考慮した非線形解析を行うＣＡＥ解析装置である。
【０００４】
ここで、機械部品の強度設計における問題の１つとして、極低サイクル疲労破壊問題がある。金属等の物体は、作用する外力がある大きさを超えると、構成材料に局部的な破壊（ボイド）が生じて、変形が徐々に大きくなる。このような物体の変形挙動（損傷）を記述する損傷モデル（力学構成式）として、例えば、Ｇｕｒｓｏｎモデルがある。また、繰り返し大負荷を受けた場合の物体の変形挙動を記述する移動（加工）硬化モデルとして、例えば、Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋ則といったモデルがある。このモデルは、弾塑性変形を伴う繰返し大負荷を受けると現れるＢａｕｓｃｈｉｎｇｅｒ効果を表現することができる。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００７−１７９３８５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
極低サイクル疲労破壊問題をシミュレートするためには、繰り返し負荷がかかる過程で生じる材料自体の損傷が、この繰り返しの過程における物体の変形挙動に及ぼす影響を考慮する必要がある。しかし、Ｇｕｒｓｏｎモデルでは、引張り負荷下で損傷が進展するが、逆に圧縮状態では損傷が回復してしまうため、繰り返し負荷を与えるような場合の損傷を表現するのに適切ではない。また、Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋ則といったモデルでは、Ｂａｕｓｃｈｉｎｇｅｒ効果を表現することはできるが、繰り返しかけられる負荷による材料自体の損傷の影響を物体の変形挙動として表現することはできない。
【０００７】
そこで、本発明は、解析対象物の有限要素解析モデルを用いて、極低サイクル疲労破壊をシミュレートできる有限要素解析装置、有限要素解析方法及び有限要素解析プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上述の目的を達成するために、本発明に係る有限要素解析装置は、非線形有限要素解析を行うための有限要素解析装置であって、解析対象の有限要素モデルと解析条件とを取得する取得手段と、前記有限要素モデルの剛性行列を作成し、作成した剛性行列と前記解析条件とで規定される剛性方程式の解に基づいて、前記有限要素モデルの加工硬化特性を算出する特性計算を行う応力ひずみ算出手段と、前記有限要素モデルの内力と外力とが釣り合うまで、前記応力ひずみ算出手段に前記特性計算を繰り返し行わせる制御手段とを備え、前記応力ひずみ算出手段は、前回算出した前記剛性方程式の解に基づく試行弾性歪みを用い、繰り返し負荷を受けている部分であって、破壊が発生していない部分の応力とひずみを算出する加工硬化モデルと、破壊が発生している部分の破壊進行度と応力を算出する材料損傷モデルとにより応力とひずみとを連成させた材料損傷方程式を作成し、作成した材料損傷方程式の解を求め、求めた解に基づいてコンシステント接線剛性を求め、求めたコンシステント接線剛性から前記剛性行列を作成することを特徴とする。
【０００９】
そして、本発明に係る有限要素解析方法は、非線形有限要素解析を行うための有限要素解析装置で用いられる有限要素解析方法であって、解析対象の有限要素モデルと解析条件とを取得する取得ステップと、前記有限要素モデルの剛性行列を作成し、作成した剛性行列と前記解析条件とで規定される剛性方程式の解に基づいて、前記有限要素モデルの加工硬化特性を算出する特性計算を行う応力ひずみ算出ステップと、前記有限要素モデルの内力と外力とが釣り合うまで、前記応力ひずみ算出ステップに前記特性計算を繰り返し行わせる制御ステップとを備え、前記応力ひずみ算出ステップは、前回算出した前記剛性方程式の解に基づく試行弾性歪みを用い、繰り返し負荷を受けている部分であって、破壊が発生していない部分の応力とひずみを算出する加工硬化モデルと、破壊が発生している部分の破壊進行度と応力を算出する材料損傷モデルとにより応力とひずみとを連成させた材料損傷方程式を作成し、作成した材料損傷方程式の解を求め、求めた解に基づいてコンシステント接線剛性を求め、求めたコンシステント接線剛性から前記剛性行列を作成することを特徴とする。
【００１０】
また、本発明に係る、コンピュータに実行させるための有限要素解析プログラムは、非線形有限要素解析を行うための有限要素解析装置で用いられる有限要素解析プログラムであって、解析対象の有限要素モデルと解析条件とを取得する取得手段と、前記有限要素モデルの剛性行列を作成し、作成した剛性行列と前記解析条件とで規定される剛性方程式の解に基づいて、前記有限要素モデルの加工硬化特性を算出する特性計算を行う応力ひずみ算出手段と、前記有限要素モデルの内力と外力とが釣り合うまで、前記応力ひずみ算出手段に前記特性計算を繰り返し行わせる制御手段として、コンピュータを機能させ、前記応力ひずみ算出手段は、前回算出した前記剛性方程式の解に基づく試行弾性歪みを用い、繰り返し負荷を受けている部分であって、破壊が発生していない部分の応力とひずみを算出する加工硬化モデルと、破壊が発生している部分の破壊進行度と応力を算出する材料損傷モデルとにより応力とひずみとを連成させた材料損傷方程式を作成し、作成した材料損傷方程式の解を求め、求めた解に基づいてコンシステント接線剛性を求め、求めたコンシステント接線剛性から前記剛性行列を作成することを特徴とする。
【００１１】
このような構成の有限要素解析装置、有限要素解析方法及び状有限要素解析プログラムによれば、加工硬化モデルと損傷モデルとを連成（強連成）したモデルを用いるので、相互に影響し合う非破壊部分と、破壊部分の剛性変化を考慮して、破壊を算出することができる。よって、サイクル負荷による破壊発生部位と破壊状況とを高精度にシミュレートすることが可能となる。
【００１２】
また、上述の有限要素解析装置において、前記応力ひずみ算出手段は、前記材料損傷方程式の解に基づいて、前記有限要素モデルの各要素の損傷値を求め、前記損傷値が所定の閾値を上回る要素は破壊されたと判断することが好ましい。
【００１３】
この構成によれば、相互に影響し合う非破壊部分と、破壊部分の剛性変化を考慮して、有限要素モデルの要素の破壊を算出することが可能となり、シミュレーションにおいて損傷値を求めることができるので、き裂が生じているかを判断することができ、破壊進展のシミュレーションを行うことが可能となる。
【００１４】
また、上述の有限要素解析装置において、前記加工硬化モデルは、分割背応力形Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｋ則による加工硬化モデルであることが好ましい。
【００１５】
この構成によれば、分割背応力モデルと損傷モデルとを連成（強連成）したモデルを用いるので、相互に影響し合う非破壊部分と、破壊部分の剛性変化を考慮して、破壊を算出することができる。
【００１６】
また、上述の有限要素解析装置において、前記コンシステント接線剛性は、完全陰解法により求めることが好ましい。
【００１７】
この構成によれば、完全陰解法による解析であるため、収束性良く、安定に計算を行うことが可能となる。
【００１８】
また、上述の有限要素解析装置において、前記制御手段は、前記有限要素モデルの内力と外力とが釣り合うと、前記解析条件のうちの強制変位又は外力を、所定の増分値だけ増加させて、前記応力ひずみ算出手段に前記特性計算を繰り返し行わせることが好ましい。
【００１９】
この構成によれば、強制変位又は外力を変えながら繰り返し前記特性計算を行わせることができるので、サイクル負荷による破壊発生部位と破壊状況を高精度に算出することが可能となる。
【００２０】
また、上述の有限要素解析装置において、前記有限要素モデルを表示する表示部を、更に備え、前記制御手段は、前記応力ひずみ算出手段が破壊されたと判断した要素を除いて、各要素をそれぞれの応力に応じた表示態様で、前記有限要素モデルを前記表示部に表示させることが好ましい。
【００２１】
この構成によれば、破壊シミュレーション結果を、表示することが可能となる。
【００２２】
また、上述の有限要素解析装置において、前記剛性算出手段は、作成した材料損傷方程式の解を求める前に、前記有限要素モデルの各要素の応力に対して、共回転変換を行うことが好ましい。
【００２３】
この構成によれば、Ｕｐｄａｔｅｄ−Ｌａｇｒａｎｇｅ形の有限変形モデルに対応することが可能となる。前記材料損傷方程式（材料損傷モデル）を有限変形モデルに対応させることで、繰り返し大負荷を受ける極低サイクル疲労破壊の中でも、変形が大きなものに対応することが可能になる。
【発明の効果】
【００２４】
このような構成の有限要素解析装置、有限要素解析方法及び有限要素解析プログラムでは、解析対象物の有限要素解析モデルを用いて、極低サイクル疲労破壊をシミュレートできる。
【図面の簡単な説明】
【００２５】
【図１】有限要素解析装置の構成を示すブロック図である。
【図２】図１の有限要素解析装置の解析処理のフローチャートである。
【図３】図１の有限要素解析装置の応力ひずみ解析処理のフローチャートである。
【図４】応力ひずみ線図の例を示す図である。
【図５】図５（ａ）は、有限要素モデルの例を示す図であり、図５（ｂ）は、材料定数の例を示す図である。
【図６】実測値によるヒステリシスループの例を示す図である。
【図７】比損傷モデルによるヒステリシスループの例を示す図である。
【図８】実施形態の材料損傷モデルによるヒステリシスループの例を示す図である。
【図９】サイクル数に伴う各ループの引張時のピーク応力を示した図である。
【図１０】実施形態の材料損傷モデルによるサイクル数に伴う損傷値の変化を示す図である。
【図１１】サイクル数の増加に伴うヤング率（剛性）の変化を示した図である。
【図１２】初期き裂進展挙動についてのシミュレーションを説明するための図である。
【図１３】図１３（ａ）は、板厚表面のき裂進展の様子を示す図であり、図１３（ｂ）は、板厚中央のき裂進展の様子を示す図である。
【発明を実施するための形態】
【００２６】
＜実施形態１＞
以下、本発明に係る実施形態を図面に基づいて説明する。
＜構成＞
図１は、有限要素解析装置の構成を示すブロック図である。図１において、有限要素解析装置１は、演算処理部１１、入力部１２、出力部１３、記憶部１４、及び、バス１９を備えて構成される。
【００２７】
演算処理部１１は、例えば、マイクロプロセッサ等を備えて構成され、機能的に、制御部１１１、初期設定部１１２、及び、応力ひずみ解析部１１３を備える。
【００２８】
制御部１１１は、初期設定部１１２、及び、応力ひずみ解析部１１３を制御して有限要素解析の演算を行わせ、また、入力部１２、出力部１３、及び、記憶部１４を制御する。
【００２９】
初期設定部１１２は、解析対象である有限要素モデルの解析処理（シミュレーション）を行う際の初期値を設定する。初期値は、例えば、材料定数、境界条件、ステップ増分値（変位、外力）などである。
【００３０】
応力ひずみ解析部１１３は、ステップ増分値に応じた、応力、ひずみ、損傷値等の物理量（加工硬化特性）を算出する。具体的には、損傷モデルと加工（移動）硬化モデルとを、強連成し完全陰的に解くための方程式（アルゴリズム）を作成し、応力ひずみ挙動を求める。応力ひずみ解析部１１３の処理の詳細は、＜動作＞の項で説明する。
【００３１】
入力部１２は、有限要素解析装置１の演算開始指示等の各種コマンドやシステムパラメータ等の各種データを有限要素解析装置１に入力する機器であり、例えば、キーボードやマウス等である。出力部１３は、入力部１２から入力されたコマンド及びデータや本有限要素解析装置１の演算結果等を出力する機器であり、例えばＬＣＤ（Ｌｉｑｕｉｄ  Ｃｒｙｓｔａｌ  Ｄｉｓｐｌａｙ）、有機ＥＬ（Ｅｌｅｃｔｒｏ  Ｌｕｍｉｎｅｓｅｎｃｅ）ディスプレイ等の表示装置である。
【００３２】
記憶部１４は、解析処理を実行する上で必要な各種プログラムや各種データを記憶する装置であり、例えば、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）やＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ  Ｅｒａｓａｂｌｅ  Ｐｒｏｇｒａｍｍａｂｌｅ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）等の不揮発性記憶素子、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等の揮発性記憶素子およびその周辺回路等を備えて構成される。
【００３３】
また、記憶部１４は、有限要素モデル記憶部１４１、解析条件記憶部１４２、及び、物理量記憶部１４３を備える。
【００３４】
有限要素モデル記憶部１４１は、解析対象の有限要素モデルを記憶する。
【００３５】
解析条件記憶部１４２は、有限要素モデル記憶部１４１に記憶されている有限要素モデルのシミュレーションに必要なデータ、例えば、材料定数、ステップ増分値等を記憶する。
【００３６】
物理量記憶部１４３は、応力ひずみ解析部１１３が算出した物理量を記憶しておく。この物理量は、ステップ増分値が加算される毎に算出され、記憶される。この物理量に基づいて、有限要素モデルの変形挙動が画像として表示される。
【００３７】
これら演算処理部１１、入力部１２、出力部１３、及び、記憶部１４は、データを相互に交換することができるようにバス１９にそれぞれ接続される。
＜動作＞
次に、本実施形態の動作について説明する。
【００３８】
図２は、本実施形態の有限要素解析装置１の解析処理のフローチャートである。尚、フローチャート内の括弧書きの番号は、後述の式の番号を示す。例えば、（１６）は、式（１６）を示す。
【００３９】
まず、ユーザは、入力部１２を介して、解析処理の開始を指示するコマンド等を入力する。入力部１２を介して、解析処理の開始を指示するコマンドが入力されたことを検知した演算処理部１１の制御部１１１は、初期設定部１１２に初期値の取得を依頼する。依頼を受けた初期設定部１１２は、出力部１３に初期設定の入力を要求するメッセージを表示する。
【００４０】
ユーザは、解析対象の有限要素モデルを作成し、入力部１２を介して、有限要素モデル記憶部１４１に記憶させる（ステップＳ１０）。そして、ユーザは、入力部１２を介して、境界条件として、解析対象モデルに対して、どの部位をどのように固定するか、あるいはどのくらい強制的な変位を与えるか、などの拘束条件（実施形態では、境界変位ｕＢＣ）、及び、どの部位にどのような方向でどのくらい力を加えるかなどの荷重条件（実施形態では、境界外力ｆＢＣ）を、入力する（ステップＳ１１）。また、ユーザは、入力部１２を介して、解析対象モデルの材料定数を入力し（ステップＳ１２）、変位ｕ、及び、外力ｆのステップ増分値（Δｕ、Δｆ）を入力する（ステップＳ１３）。尚、実施形態では、説明の便宜上、変位ｕ及び外力ｆの両方を変化させることとして説明するが、変位ｕ及び外力ｆのうちのいずれか一方を変化させることとしてもよい。つまり、変位ｕを変化させる方法、外力ｆを変化させる方法、変位ｕ及び外力ｆを変化させる方法の３つのうちのいずれかでシミュレーションすることとすればよい。いずれを変化させるかは、実施するシミュレーションに応じて決定される。
【００４１】
ステップＳ１１〜Ｓ１３において入力されたデータは、解析条件記憶部１４２に記憶される。尚、ステップＳ１０〜ステップＳ１３のデータは、予め有限要素モデル記憶部１４１及び解析条件記憶部１４２に記憶させておき、ユーザが選択することとしてもよい。
【００４２】
初期設定部１１２による初期設定が終了すると、制御部１１１は、現在のステップｎ（最初はｎ＝１）、つまり、応力ひずみ解析部１１３に渡す変位ｕｎ、及び、外力ｆｎを算出する（ステップＳ１４）。制御部１１１は、ステップごとにステップ増分値ずつ増加させて、予め定めたステップの数だけ応力ひずみ解析部１１３に物理量を繰り返し算出させる。従って、現在のステップｎにおける変位ｕｎ、外力ｆｎは、前回の変位ｕｎ−１、外力ｆｎ−１にそれぞれΔｕ、Δｆが加算されて求められる。
【００４３】
制御部１１１は、変位ｕｎ、及び、外力ｆｎを応力ひずみ解析部１１３に渡して、応力ひずみ解析処理を行わせて物理量、すなわち、応力σｎ＋１、等方硬化内部変数Ｒｎ＋１、全体背応力βｎ＋１、分割背応力βｋｎ＋１、損傷値Ｄｎ＋１を求め（ステップＳ１５）、求めた物理量を物理量記憶部１４３に記憶させる（ステップＳ１６）。
【００４４】
制御部１１１は、変位ｕｎ、外力ｆｎそれぞれが、境界変位ｕＢＣ、境界外力ｆＢＣとなるまで（ステップＳ１７：ＮＯ）、ステップＳ１４〜Ｓ１６までの処理を繰り返し実行する（ステップＳ１７：ＹＥＳ）。
【００４５】
尚、極低サイクル疲労問題等のシミュレーションを行う場合は、さらに、ステップＳ１７においてＹＥＳとなった後に、ステップ増分値（Δｕ、Δｆ）を変えて、ステップＳ１４〜Ｓ１７までの処理を繰り返す。例えば、負荷を解除した場合は、Δｕ及びΔｆをマイナスとして、ステップＳ１４〜Ｓ１７までの処理を繰り返す。そして、次のサイクルでは、Δｕ及びΔｆを、プラスとしてステップＳ１４〜Ｓ１７までの処理を繰り返し、Δｕ及びΔｆを、マイナスとしてステップＳ１４〜Ｓ１７までの処理を繰り返す。
【００４６】
シミュレーションが終了すると、制御部１１１は、ユーザの指示に応じて、物理量記憶部１４３に記憶されている物理量に基づいて、時系列に、有限要素モデルを出力部１３に表示する。例えば、有限要素モデルの各要素を、要素の応力値に応じた色で表示するなどである。また、破壊されたと判断された要素は、描画されないようにすることで、破壊の進展の様子が表示されることになる。
【００４７】
次に、図３を用いて、図２のステップＳ１５の応力ひずみ解析処理について説明する。図３は、応力ひずみ解析処理のフローチャートである。
【００４８】
応力ひずみ解析部１１３は、制御部１１１から渡されたステップｎの変位ｕｎ、及び、外力ｆｎを用いて、以下の式（Ｇ１）で示す剛性方程式を作成する。
【数Ｇ１】
2016004543.tif 000003
[Ｋ]は全体剛性マトリックス、｛Δｕ｝は節点変位ベクトル、｛Δｆ｝は節点力ベクトルを示す。応力とひずみの関係は、非線形であるため、剛性マトリックスは応力やひずみによって変化する。従って、最初のステップ１では、予め定めておいた剛性マトリックスが用いられ、ステップ２からは、後のステップＳ２８で作成された剛性マトリックスが用いられる。
【００４９】
次に、応力ひずみ解析部１１３は、作成した剛性方程式の解を算出し、算出した解から各要素のひずみを算出する（ステップＳ２０）。
【００５０】
次に、応力ひずみ解析部１１３は、各要素のＣａｕｃｈｙ応力に対して、式（１６）で示される共回転変換を行う（ステップＳ２１）。
【００５１】
【数１６】
2016004543.tif 000004
Ｌは増分内でのスピンが一定であると仮定した時の物質点の回転を示す直交テンソルである。
【００５２】
次に、応力ひずみ解析部１１３は、以下の式（Ａ２０）に示す材料損傷モデルである強連成方程式を作成する（ステップＳ２２）。具体的には、損傷モデルとして連続体損傷力学に基づく損傷モデルと、加工硬化モデルとして分割背応力型Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋモデルとの強連成方程式（以下、「材料損傷モデル」、「材料損傷方程式」という。）を作成する。材料損傷モデルの導出方法については、＜材料損傷モデル＞の項で説明する。
【数Ａ２０】
2016004543.tif 000005
【００５３】
次に、応力ひずみ解析部１１３は、作成した材料損傷方程式をＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法により、求解し、塑性常数Δγを求める（ステップＳ２３、Ｓ２４）。応力ひずみ解析部１１３は、求めた塑性常数Δγを、後述する式（１４）に代入して、ひずみに応じた物理量を算出する（ステップＳ２５）。
【００５４】
次に、応力ひずみ解析部１１３は、算出した物理量のうちの損傷値Ｄｎ＋１を基に、破壊判定を行う（ステップＳ２６）。損傷値Ｄｎ＋１が閾値ＤＦ（材料定数）よりも大きい場合、応力ひずみ解析部１１３は、その要素は破壊されたと判断して、除去する。つまり、応力ひずみ解析部１１３によって破壊されたと判断された要素は、以後の演算の対象から外されることになり、解析対象の有限要素モデルを描画する際にも表示されないことになる。尚、実施形態では、破壊されたと判断された要素は、除去することとしているが、その要素の剛性を極端に小さくし、実質的に除去されているものと近似して計算することとしてもよい。
【００５５】
次に、応力ひずみ解析部１１３は、Ｃｏｎｓｉｓｔｅｎｔ接線剛性を算出する（ステップＳ２７）。尚、Ｃｏｎｓｉｓｔｅｎｔ接線剛性を算出する方法は、＜Ｃｏｎｓｉｓｔｅｎｔ接線剛性＞の項で説明する。
【００５６】
ここで、図４に、材料に一定の負荷をかけ続けた場合の応力ひずみ線図を示す。縦軸は、応力を示し、横軸は、ひずみを示す。曲線のグラフは、実測値を示す。横軸の実線矢印は、シミュレーションにおける変形量（ひずみ増分）を示す。この変形量は、剛性方程式の解から算出された各要素のひずみの、合計である全ひずみ増分である。
【００５７】
シミュレーションにおける応力ひずみ線図が実測値と近似するためには、応力ひずみ線図の傾き、つまり、接線剛性が、より実測に近く求められる必要がある。シミュレーションでは、１ステップごとに加える負荷を増加させて計算を行うので、このステップ間（２点間）においては線形と考えて、Ｃｏｎｓｉｓｔｅｎｔ接線剛性を求めるものとする。より実測に近いＣｏｎｓｉｓｔｅｎｔ接線剛性を求めるためには、図４の応力ひずみ線図で示すように、より実測に近い積分応力を求める必要がある。実施形態の材料損傷モデルは、この応力を、より実測に近く求めることができる。尚、図４の横軸の実線矢印で示される全ひずみ増分Δεは、材料損傷モデルの式（Ａ２０）におけるＭとＹｎ＋１に含まれており、材料損傷モデルにより全ひずみ増分Δεに対応する積分応力を求めることができる。Ｍは、後述の式（Ａ１８）のようにｑＥを含み、ｑＥは、式（Ａ９）のように、ηＥを含み、ηＥは、式（Ａ７）のようにσｔｒｙｎ＋１を含み、σｔｒｙｎ＋１は、式（Ａ３）のようにεｅｔｒｙｎ＋１を含み、εｅｔｒｙｎ＋１は、式（１５）のようにΔεを含む。また、Ｙｎ＋１は、式（Ａ１２）のようにεｅｔｒｙｎ＋１を含み、εｅｔｒｙｎ＋１は、式（１５）のようにΔεを含む。
【００５８】
Ｃｏｎｓｉｓｔｅｎｔ接線剛性Ｄｅｐ（チルド）は、以下の式（Ｇ２）で示すように、剛性マトリクス[Ｋ]を求めるために用いられる。尚、「英文字記号（チルド）」は、式中で英文字記号（変数）に「〜」が付されているものを示す。
【数Ｇ２】
2016004543.tif 000006
[Ｂ]は、ひずみと変位の関係を示すＢマトリクス、ｖは要素の体積を示す。
【００５９】
応力ひずみ解析部１１３は、求めたＣｏｎｓｉｓｔｅｎｔ接線剛性Ｄｅｐ（チルド）を用いて、剛性マトリックス[Ｋ]を作成する（ステップＳ２８）。
【００６０】
応力ひずみ解析部１１３は、有限要素モデルの内力から外力を引いた値（残差力）を算出し、残差力がゼロ以下である場合（内力と外力とが釣り合った場合）は（ステップＳ２９：ＯＫ）、計算が収束したと判断して、処理を終了する。応力ひずみ解析部１１３は、この時の物理量を、現在のステップの物理量として、制御部１１１に出力し、処理を終了する。
【００６１】
一方、残差力がゼロより大きい場合は（ステップＳ２９：ＮＧ）、応力ひずみ解析部１１３は、ステップＳ２０からの処理を繰り返す。尚、ステップＳ２１〜Ｓ２７の処理は、有限要素モデルを構成する全要素について行われる。ステップＳ２０、Ｓ２８の処理は、ステップＳ２１〜Ｓ２７で行われた個々の要素に対する演算結果、つまり、全要素に対する処理結果を元にした、解析対象の有限要素モデル（例えば、部品）の形状全体についての振る舞いとしての処理である。
【００６２】
＜材料損傷モデル＞
次に、以下の式（Ａ２０）で示される方程式ｆ（Δγ）、実施形態の材料損傷モデルについて説明する。実施形態の材料損傷モデルは、分割背応力形Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｋ則による加工硬化モデルと損傷モデルとを強連成した静的、完全陰解法の計算手法である材料損傷方程式である。
【数Ａ２０】
2016004543.tif 000007
Ｇはせん断弾性係数、Ｍは後述の式（Ａ１８）第２式、Ｄｎは損傷値、Ｙｎ＋１は損傷エネルギ解放率、ｒ、Ｓは損傷発展速度を決める定数である。
【００６３】
応力ひずみ解析部１１３は、材料損傷方程式を作成し、作成した方程式をＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法で解いてΔγを求める。そして、塑性常数Δγを、後述する式（１４）に代入することで、５つの物理量、応力σｎ＋１、等方硬化内部変数Ｒｎ＋１、全体背応力βｎ＋１、分割背応力βｋｎ＋１、損傷値Ｄｎ＋１を求める。
【００６４】
まず、分割背応力型Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋモデルについて説明する。分割背応力形Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋ則は、Ｃｈａｂｏｃｈｅらによって提案されたモデルであり、以下の式（１）で示される。このモデルは、Ｎ組の分割背応力βｋを重ね合わせる事で、全体背応力を表現するものである。以下、表記に関してスカラ(０階テンソル)量は細字，２階テンソル量は太字，４階テンソル量は太字にチルド記号を付けて示す。尚、以下、式中で英文字記号（変数）に「・」又は「〜」を付すものは、「英文字記号（チルド）」と記載するものとする。
【００６５】
【数１】
2016004543.tif 000008
応力空間における降伏曲面の移動量は、この全体背応力を用いて表現される。また各分割背応力は、以下の式（２）のように、非線形の背応力挙動を表現するＡｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋ則によって表現される。
【００６６】
【数２】
2016004543.tif 000009
ここでａｋ、ｂｋは、各分割背応力に対する材料定数、εｐ（チルド）は塑性ひずみ、γ（チルド）は塑性乗数を示す。このように分割背応力形Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋ則は全体背応力を複数の非線形特性の組み合わせにより表現する事ができるため、単一の背応力を用いた場合のＡｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋ則よりも比較的複雑な背応力発展挙動を表現する事が可能となる。
【００６７】
また、等方硬化則として、以下の式（３）のように、Ｍ個の指数関数則を組み合わせた形式を用いる。δｙは等方硬化による相対相当応力、δｙ０は初期降伏応力、Ａｋ、Ｂｋは材料定数、Ｒは等方硬化の内部変数を示す。
【００６８】
【数３】
2016004543.tif 000010
分割背応力形Ａｒｍｓｔｒｏｎｇ−Ｆｒｅｄｅｒｉｃｋ則と連続体損傷力学に基づく材料損傷を考慮した散逸ポテンシャルΨは、以下の式（４）で示される。
【００６９】
【数４】
2016004543.tif 000011
ここでＦｐは降伏関数、Ｆβは背応力項、ＦＤは損傷項を示す。δはＣａｕｃｈｙ応力、Ｄは材料損傷値、ｒ、Ｓは損傷発展速度を決める定数である。
【００７０】
式（４）中のＹは損傷発展による内部エネルギ密度の変化量を示すもので、損傷エネルギ解放率と呼ばれ、以下の式（５）のように応力三軸度関数Ｒνと関与する量である。Ｅはヤング率、νはポアソン比、ｐは静水圧、σｅｑは相当応力(Ｍｉｓｅｓ応力)を示す。
【００７１】
【数５】
2016004543.tif 000012
更に、弾性材料剛性テンソルをＤｅ（チルド）、偏差射影テンソルをＩｄ（チルド）、対称射影テンソルをＩｓ（チルド）、４階の恒等テンソルをＩ（チルド）、２階の等方テンソルをＩとおくと、それぞれ以下の式（６）で表記される。ただしＩの成分は、クロネッカーのデルタδと等しくなる。λ、μはＬａｍｅの定数である。
【００７２】
【数６】
2016004543.tif 000013
式（４）のように散逸ポテンシャルに降伏関数以外の項が含まれるため、以下の式（７）のような非関連流れ則が導かれる。
【００７３】
【数７】
2016004543.tif 000014
ここで有効相対応力ηを、以下の式（８）に示すようにおくと、
【００７４】
【数８】
2016004543.tif 000015
式（７）中の各流れベクトルは、以下の式（９）となる。
【００７５】
【数９】
2016004543.tif 000016
ただし、式（７）第３式には、以下の式（１０）で示されるように、分割背応力βｋに関する自由エネルギΨβに関して熱力学的応力すなわち分割背応力βｋとその内部変数Ｘｋとの関係式を用いた。
【００７６】
【数１０】
2016004543.tif 000017
連続体損傷力学では、マイクロクラックやマイクロボイドなどの微視的な材料損傷によって内力を受け持つ有効面積の減少が、材料の巨視的な弾性材料剛性の低下を引き起こすものと考えており，構成則を式（１１）のように表現する。
【００７７】
【数１１】
2016004543.tif 000018
式（７）に示される分割背応力の発展則をそのまま使用して応力積分を行うと、全ての分割背応力βｋが未知数となり応力積分に多大の時間を要してしまう。これを避けるためにステップｎからステップｎ＋１に掛けて塑性ひずみ増分の方向が一定であると仮定して、式（７）第３式に対してこの時間増分間での解析解を求めると、最終的に式（１２）および式（１３）が得られる。ただし添字ｎ、ｎ＋１は、そのステップにおける量を示す。式（１２）および式（３）のように時間増分間における解析解を用いる事で、各分割背応力はステップｎの既知量となり全体背応力βｎ＋１のみが背応力の未知数となる。
【００７８】
【数１２】
2016004543.tif 000019
ただし、
【００７９】
【数１３】
2016004543.tif 000020
である。式（７）の流れ則を後退形Ｅｕｌｅｒ表記したものと、式（１２）および弾塑性状態における降伏関数の適合条件から、式（１４）で示すＲｅｔｕｒｎ−Ｍａｐｐｉｎｇ方程式が非線形連立方程式として得られる。
【００８０】
【数１４】
2016004543.tif 000021
ただし、εｅｔｒｙｎ＋１は試行弾性ひずみと呼ばれ、式（１５）のように前ステップの弾性に全ひずみ増分を加えたもので定義する。
【００８１】
【数１５】
2016004543.tif 000022
よって、式（１４）は、応力σｎ＋１、等方硬化内部変数Ｒｎ＋１、全体背応力βｎ＋１、損傷値Ｄｎ＋１、塑性常数Δγの合計１５個の独立な未知変数成分を含む方程式となる。
【００８２】
さらに、式（１４）第１式より、
【数Ａ１】
2016004543.tif 000023
となる。
また有効相対応力ηｎ＋１は、式（８）に式（Ａ１）、式（１２）を用いる事で、
【数Ａ２】
2016004543.tif 000024
となる。ただし試行応力を以下の式（Ａ３）のように定義した。
【数Ａ３】
2016004543.tif 000025
さらに降伏曲面に対する法線ベクトルは式（１３）より、
【数Ａ４】
2016004543.tif 000026
となる。ただし相対相当応力を、
【数Ａ５】
2016004543.tif 000027
と定義した。これらを式（Ａ２）に代入して整理すると、
【数Ａ６】
2016004543.tif 000028
が得られる。ここで、
【数Ａ７】
2016004543.tif 000029
とする。さらに式（Ａ６）の両辺に√３／２をかけて、両辺のノルムを取ると、
【数Ａ８】
2016004543.tif 000030
となる。さらに式（１４）第５式より、
【数Ａ９】
2016004543.tif 000031
とおくと、式（Ａ８）は、
【数Ａ１０】
2016004543.tif 000032
となる。式（Ａ１０）と、損傷発展式である式（１４）第４式とを連立させる事で、Ｒｅｔｕｒｎ−Ｍａｐｐｉｎｇ方程式は、
【数Ａ１１】
2016004543.tif 000033
となる。この段階で、Ｒｅｔｕｒｎ−Ｍａｐｐｉｎｇ塑性方程式は、塑性乗数Δγと損傷値Ｄｎ＋１の２個のスカラ未知数にまで縮約させた。
次に、まず損傷エネルギ解放率Ｙについて考える。式（５）より、
【数Ａ１２】
2016004543.tif 000034
となる。式（Ａ１２）の右辺第１項は、式（８）および式（１２）を用いて、
【数Ａ１３】
2016004543.tif 000035
と表記できる。さらに式（Ａ６）より、
【数Ａ１４】
2016004543.tif 000036
となる。これを式（Ａ１３）に代入すると、
【数Ａ１５】
2016004543.tif 000037
となり。式（Ａ１５）右辺をｇと置く。以上より、損傷エネルギ解放率Ｙを、式（Ａ１６）のように書き直すことができる。
【数Ａ１６】
2016004543.tif 000038
次に、式（Ａ１１）第1式より、
【数Ａ１７】
2016004543.tif 000039
ここで、
【数Ａ１８】
2016004543.tif 000040
とおくと、式（Ａ１７）から、
【数Ａ１９】
2016004543.tif 000041
となる。式（Ａ１８）と（Ａ１９）とを、式（Ａ１１）第２式に代入して整理すると、式（Ａ２０）に示すような、塑性乗数Δγのみのスカラ１変数に縮約されたＲｅｔｕｒｎ−Ｍａｐｐｉｎｇ方程式が得られる。
【数Ａ２０】
2016004543.tif 000042
この式（Ａ２０）で示す方程式ｆ（Δγ）が、実施形態の材料損傷方程式（モデル）である。
【００８３】
実施形態では、さらに、応力ひずみ解析部１１３は、Ｕｐｄａｔｅｄ−Ｌａｇｒａｎｇｅ形の有限変形モデルに対応できるよう、応力積分計算（式（Ａ２０）の求解）の前にＣａｕｃｈｙ応力に対して、式（１６）で示される共回転変換を行う。
【００８４】
【数１６】
2016004543.tif 000043
ただし、Ｌは増分内でのスピンが一定であると仮定した時の物質点の回転を示す直交テンソルである。実施形態の材料損傷方程式（材料損傷モデル）を有限変形モデルに対応させることで、繰り返し大負荷を受ける極低サイクル疲労破壊の中でも、変形が大きなものに対応することが可能になる。このような変形が大きな破壊とは、例えば、「針金を繰り返し何度も大きく曲げると数回目で折れる」、「ジュースの空き缶のフタを何度も繰り返し大きく曲げていると数回目で折れる」といった事例である。
【００８５】
＜Ｃｏｎｓｉｓｔｅｎｔ接線剛性＞
応力ひずみ解析部１１３が、Ｃｏｎｓｉｓｔｅｎｔ接線剛性を算出する方法について説明する。Ｃｏｎｓｉｓｔｅｎｔ接線剛性は、完全陰解法によるものである。
【００８６】
式（Ａ６）より、ηｎ＋１とηＥとは互いにスカラ倍されている関係に過ぎないので、両者の方向は一致する。よって、
【数Ｂ１】
2016004543.tif 000044
となる。式（Ｂ１）および式（Ａ１８）を用いて、式（Ａ１）で示される構成側を以下のように変形する。
【数Ｂ２】
2016004543.tif 000045
これよりσｎ＋１は内部変数を用いて、
【数Ｂ３】
2016004543.tif 000046
と表記できる。以上より、Ｃｏｎｓｉｓｔｅｎｔ接線剛性Ｄｅｐ（チルド）は、σｎ＋１の各内部変数を考慮し、以下の式（Ｂ４）に示すように、
【数Ｂ４】
2016004543.tif 000047
と表記できる。各微分について導出すると以下の式（Ｂ５）が得られる。
【数Ｂ５】
2016004543.tif 000048
また式（Ｂ４）の塑性乗数の試行弾性ひずみによる微分は、式（Ａ２０）で示した縮約されたＲｅｔｕｒｎ−Ｍａｐｐｉｎｇ方程式を、試行弾性ひずみで微分する事により以下の式（Ｂ６）のように得られる。
【数Ｂ６】
2016004543.tif 000049
右辺の分母の微分は、式（Ａ２２）から既に求まっている。そこで分子の微分について考えると、以下の式（Ｂ７）に示すように、
【数Ｂ７】
2016004543.tif 000050
となる。各微分について計算すると、以下の式（Ｂ８）を得る。
【数Ｂ８】
2016004543.tif 000051
以上より、陽的にＣｏｎｓｉｓｔｅｎｔ接線剛性を計算する事が可能となる。
＜シミュレーションの例＞
ここでは、図５に示す有限要素モデルを用い、上面に対してＺ方向の両振りの繰返し強制変位を与えた場合のシミュレーション結果を説明する。尚、有限要素モデルの上面以外の他の面については、単軸負荷の対称条件を満たす様に拘束するものとする。図５（ａ）は、有限要素モデルを示す図であり、図５（ｂ）は、シミュレーションで用いた材料定数の表を示す図である。シミュレーションで用いたヤング率、ポアソン比、初期降伏応力、等方硬化定数、移動硬化定数、損傷定数は、図５（ｂ）の表に示すとおりである。
【００８７】
図６は、実測値の応力ひずみ線図であり、ひずみ振幅一定で繰り返し負荷を与えた時の、応力、ひずみ挙動を示す。詳細には、ひずみは、−０．０４から＋０．０４まで変化させられ、「Ｓ５００ｓ，Δε＝８％，４５ｄａｙｓ」が示すように、入力負荷のひずみがΔ８％であり、モデルの材料はＳ５００ｓ鋼であって、４５日間腐食させた材料である。図６の応力ひずみ線図は、Ｂａｕｓｃｈｉｎｇｅｒ効果が表れたヒステリシスループとなっている。図７は、非損傷モデルによるシミュレーションの結果を示すヒステリシスループであり、図８は、実施形態の材料損傷モデルによるシミュレーションの結果を示すヒステリシスループである。
【００８８】
比較用の実測データの選定として、一般的な鋼材やアルミ材でも同様の傾向は見られるため差し支えないが、図６に示したように損傷に伴う変形挙動（矢印で示すような、サイクルの増加につれて生じるズレ）がより顕著に見られるため事前に腐食させたＳ５００ｓ鋼のデータを参考にした。ヤング率が一般的な鋼に対して低めとなっている理由は腐食によるものと考えられる。
【００８９】
また、ここで参考にした実材料の詳細な損傷定数が不明なため各材料定数として、図５（ｂ）の表に示すような現実的な金属材料の範囲で便宜的な値を入力し、定性的な変形挙動を表現できるかのみ検証した。
【００９０】
図７は、図５の有限要素モデルに公称ひずみ振幅４%の両振負荷を１０ｃｙｃｌｅ掛けた時のヒステリシスループであり、非損傷モデルを用いた場合の図である。図８は、実施形態の材料損傷モデルを用いた場合の図である。ヒステリシスループにおいて、横軸は強制変位量を元の長さで割った公称ひずみを示し，縦軸は引張面の荷重合計値を変形前の断面積で割った公称応力を示す。実測結果では図６のようにサイクル数の増加に伴う明瞭な加工軟化挙動が見られるが、非損傷モデルでは図７のようにヒステリシスループが塑性シェイクダウンを起こし、実測のような加工軟化挙動が見られない。一方、実施形態の材料損傷モデルでは図８のサイクル数の増加に伴う加工軟化挙動が確認できる（矢印参照）。
【００９１】
図９は、サイクル数に伴う各ループの引張時のピーク応力を示した図であり、縦軸は応力を示し、横軸はサイクル数を示す。白丸が非損傷モデルでの結果を示し、灰色の丸が実施形態の材料損傷モデルでの結果を示す。この図からも、実施形態の材料損傷モデルによる結果は、図６の実測と同様にサイクル数に伴う応力低下が見られ、実施形態の材料損傷モデルは、加工軟化挙動を定性的に表現できている事がわかる。
【００９２】
また、図１０は、実施形態の材料損傷モデルによるサイクル数に伴う損傷値の変化を示す図であり、縦軸が損傷値を示し、横軸がサイクル数を示す。図１０が示すように、実施形態の材料損傷モデルでは、サイクル数の増加に伴い損傷値が増加してゆく傾向が見られる。すなわち、実施形態の材料損傷モデルにおいて加工軟化傾向が表現されるのは、この損傷値の増加に対して、上記式（３）のように降伏関数の収縮による応力積載能力の低下が要因であると考えられる。
【００９３】
さらに、実測結果では図６のようにサイクル数の増加に伴い弾性領域のヤング率（剛性）が低下している様子が確認できる。図１１は、サイクル数の増加に伴うヤング率（剛性）の変化を示した図であり、縦軸がヤング率（剛性）を示し、横軸がサイクル数を示す。白丸が非損傷モデルでの結果を示し、灰色の丸が実施形態の材料損傷モデルでの結果を示す。非損傷モデルでは、弾塑性変形の履歴に関わらずヤング率（剛性）は一定値となっているが、実施形態の材料損傷モデルでは、実測結果と同様にサイクル数の増加に対してヤング率（剛性）が低下している。これは図１０のようにサイクル数の増加に伴って大きくなる損傷値に対して、式（９）で示されるように巨視的弾性材料剛性の低下が表現されるためである。
【００９４】
以上のように、実施形態の材料損傷モデルを用いる事で、繰返し大負荷を掛けた時の応力、ひずみ挙動を定性的に良く表現できている事がわかる。
【００９５】
次に、初期き裂進展挙動についてのシミュレーションを説明する。図１２に示すように、両端単純支持ばりを表す有限要素モデルに対し、矢印で示される集中負荷をかける場合をシミュレートする。有限要素モデルである梁は、Ａ５３３Ｂ鋼からなる高さＷ（＝２０ｍｍ）、幅Ｔ（＝４ｍｍ）、長さＬ（＝８０ｍｍ）の直方体であり、シミュレーションの計算は、梁の幅方向２ｍｍ分の２分の１モデルで行う。また、長さ方向の一端からＬ０（＝１２ｍｍ）の位置に、初期き裂が生じているものとする。この状態から、梁の長さ方向中央に５ｍｍの曲げ強制変位を生じさせる荷重をかけ続けた場合の結果を、図１３に示す。尚、幅方向の手前側の端面を板厚表面といい、この板厚表面から幅方向２ｍｍの位置における前記板厚表面と平行な面を板厚中央面というものとする。
【００９６】
次に、初期き裂進展挙動についてのシミュレーションを説明する。図１２に示すように、両端単純支持ばりを表す有限要素モデルに対し、矢印で示されるように集中荷重をかける場合をシミュレートする。
【００９７】
図１３（ａ）は、板厚表面のき裂進展の様子を示す図であり、画像１３１は、き裂が進展した実際の材料の写真（実験結果）であり、実線矢印は、き裂進展方向を示す。また、画像１３２は、実施形態の材料損傷モデルを用いたシミュレーション結果を示す画像である。また、図１３（ｂ）は、板厚中央面のき裂進展の様子を示す図であり、画像１３３は、き裂が進展した実際の材料の切断面の写真（実験結果）であり、破線は、き裂進展前の予き裂面を示し、実線矢印は、き裂進展方向を示す。また、画像１３４は、実施形態の材料損傷モデルを用いたシミュレーション結果を示す画像である。シミュレーション結果の画像１３２、１３４において要素（三角形）が記載されていない部分は、き裂部分を示す。
【００９８】
板厚表面では、図３（ａ）の実験結果の画像１３１に示すように、二手に分かれてき裂進展し、シミュレーション結果の画像１３２でも、実験結果の画像１３１と同様に、二手に分かれてき裂が進展していることが分かる。また、板厚中央面では、図３（ｂ）の実験結果の画像１３３に示すように、予き裂面から紙面右側に偏ってき裂が進展し、シミュレーション結果の画像１３４でも、実験結果の画像１３３と同様に、紙面右手に偏ってき裂が進展していることが分かる。
【００９９】
尚、図１３の画像１３１、１３３で示す破壊は、き裂先端近傍が大きく変形している破壊である。上述した式１６を用いて実施形態の材料損傷方程式（材料損傷モデル）を有限変形モデルに対応させたことにより、画像１３１、１３３で示すような、繰り返し大負荷を受ける極低サイクル疲労破壊の中でも変形が大きなものにも対応することが可能となり、画像１３２、１３４に示すようなシミュレーション結果を得ることが可能となっている。
【０１００】
実施形態の材料損傷モデルでは、損傷値が算出されるので、材料定数で示される所定の閾値を上回った場合に、き裂が生じたと判断することが可能となる。従って、損傷値が閾値を超えた要素を描画しないことで、図１３の画像１３２、画像１３４に示すようなき裂が進展した画像を表示することが可能となる。
【符号の説明】
【０１０１】
１  有限要素解析装置
１１  演算処理部
１２  入力部
１３  出力部
１４  記憶部
１１１  制御部
１１２  初期設定部
１１３  応力ひずみ算出部
１４１  有限要素モデル記憶部
１４２  解析条件記憶部
１４３  物理量記憶部
【図１】
2016004543.tif 000052
【図２】
2016004543.tif 000053
【図３】
2016004543.tif 000054
【図４】
2016004543.tif 000055
【図５】
2016004543.tif 000056
【図６】
2016004543.tif 000057
【図７】
2016004543.tif 000058
【図８】
2016004543.tif 000059
【図９】
2016004543.tif 000060
【図１０】
2016004543.tif 000061
【図１１】
2016004543.tif 000062
【図１２】
2016004543.tif 000063
【図１３】
2016004543.tif 000064