Ni 基超合金における結晶異方性を考慮したクリープ損傷シミュレーション

(1)

Ni 基超合金における結晶異方性を考慮した

クリープ損傷シミュレーション

大見敏仁

*

_{, 横堀壽光}

**

Simulation of creep damage formation taking into consideration crystal anisotropy for

Ni-base super alloy

Toshihito Ohmi, A. Toshimitsu Yokobori Jr.

Abstract:

Polycrystalline Ni-base super alloys are used for high temperature components, such as turbine blades and turbine disks under high tempreture condition. In order to use Ni-base super alloys safely and efficiently, understanding behavior and mechanism of creep damage formation is necessary. In this study, the analysis of creep damage formation for polycrystalline Ni-base super alloy was conducted. On the basis of results, grain boundaries of the creep crack initiation were obtained.

KEY WORDS : Ni based super alloy, creep damage, numerical analysis, EBSD 要旨: 発電機器に用いられる多結晶Ni基超合金の損傷形成挙動を明らかにすることで，実機の安全性や経済性の向上が期待できる。実験的にクリープ損傷の形成過程を明らかにする研究は多いが，数値解析によりき裂の発生個所を特定できる研究は少ない。本研究では，結晶粒の形状および弾塑性異方性を数値解析に取り入れ，弾塑性クリープ解析手法の開発を試みた。異方性は，EBSD観察の出力データであるオイラー角を利用できることを想定し，弾性異方性と塑性異方性に関するパラメータを算出した。その結果，実験結果と比較してき裂の発生個所を明確に特定できる解析手法を開発することができた。キーワード：Ni基超合金, クリープ損傷, 数値解析, EBSD １．研究背景発電機器などに使用されるガスタービン動翼には，精密鋳造ニッケル基超合金製のものが用いられる。これらの中でも比較的安価な多結晶材の動翼は，一般的な材料に比べて結晶が粗大で，同様の使用履歴であってもき裂の発生する結晶粒界とそうでない結晶粒界が存在する。発電機器の安全性と経済性を両立するために，き裂の発生しやすい粒界の特定や損傷度の評価を行い，ガスタービン動翼の適切な交換時期の決定が必要になる。

一方，電子線後方散乱（electron back scatter

diffraction : EBSD）法の発達により，金属の結晶方位分布の詳細な計測が可能となっている。この手法を用いて使用前後の結晶方位分布を比較し，高温クリープ損傷と結び付ける研究がなされている。本研究では，これらの技術を応用し，ミクロ組織分布，結晶方位等の影響を模擬した力学モデルの構築と損傷シミュレーション技術の開発を行った。 2．異方性のモデル化本研究では，Ni 基超合金のき裂発生前後の EBSD 解析結果を用い，数値解析に用いるミクロ組織分布，結晶方位等の影響を模擬した有限要素モデルを構築する。その後，このモデルに対して弾塑性クリープ有限要素解析を行い，粒界局所応力によるひずみの変化とEBDS 解析結果を比較し，き裂が発生しやす * 湘南工科大学工学部機械工学科講師 **帝京大学工学部戦略的イノベーション研究拠点客員教授

(2)

い粒界の特徴を考察した。一般的に結晶粒界は複雑な形をしており局所的には高い多軸性を持つと考えられる。このため，より局所的な荷重（ひずみ）の負荷方向とそれに依存した弾塑性異方性を取り入れることで，割れ発生粒界をより高精度に予測することを試みた。 2.1 弾性異方性（ヤング率）の定義塑性多結晶体の弾塑性異方性を有限要素法でモデリングするために，本解析では結晶粒ごとに弾性異方性と塑性異方性を考慮した。弾性異方性については，EBSD 解析により得られる結晶方位を示すオイラー角と単結晶体の弾性剛性マトリクス（ヤング率に相当）を用いて座標変換により求める。一般的に有限要素法において用いられている弾性応力-ひずみの関係式は式(1), (2)により表される． (1) (2) ここで，は弾性応力増分，は弾性ひずみ増分，は弾性応力-ひずみ行列，は弾性コンプライアンス行列である。3 次元の異方性を表す場合，独立変数は９個となり，商用ソフトなどではヤング率(３方向)，ポアソン比(３方向)，横弾性係数(３方向) などを用いて記述している。本解析においては，弾性直交異方性の場合を考慮し，結晶構造が立方晶であれば，およびは次式のようになり，独立なパラメータが3 つで表される。 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (4) ここで， , , は弾性応力-ひずみ行列成分， , , は弾性コンプライアンス行列成分である。弾性応力-ひずみ行列の結晶方位座標と全体座標との間の変換，すなはち，剛性マトリクスの回転を考えるためには，を4 階のテンソルである弾性定数テンソルに戻す必要がある。そのときの弾性応力-ひずみの関係式は式(5)(6)のように示される。 (5) ここで，は弾性定数テンソル(４階のテンソル)である。次に，テンソルの回転を考える。4 階のテンソルである弾性定数テンソルの回転は，回転マトリクス[G] を用いた座標変換として指数表記すると，次式で与えられる。 ′ (7) (i, j, k, l, m, n, o, p) =1, 2, 3 ここで， : 全体座標系における弾性定数テンソル成分 ′ : 結晶座標系における弾性定数テンソル成分 G：回転マトリクスである。一方，回転マトリクス[G]はオイラー角(1，，2 ) を用いて式(8)と表せる。このオイラー角(1，，2 ) はEBSD のデータとして，計測点毎の数字が実験的に得られる。オイラー角の定義[1, 2]を Fig.1 に示す。 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (6)

 

_                                            cos sin 1 cos sin 1 sin sin 2 cos cos 2 cos 1 cos 2 sin 1 sin cos 2 cos 1 sin 2 sin 1 cos sin 2 sin cos 2 sin 1 cos 2 cos 1 sin cos 2 sin 1 sin 2 cos 1 cos ij G

(8)

(3)

Fig.1 Bunge の定義によるオイラー角[1, 2] また，座標変換元（結晶座標系）として用いる弾性係数テンソルは以下のようになる。 ′ ′ ′ (9) ここで， ′ ：結晶座標系弾性応力 ′ ：結晶座標系弾性ひずみ ′: 結晶座標系における弾性定数テンソル成分である。これと回転マトリクス[G]を用いて，座標変換を行い，回転後の全体座標系弾性定数テンソルを求める。回転後の全体座標系弾性定数テンソル・結晶系弾性定数テンソル・全体座標系弾性応力・全体座標系弾性ひずみの関係を以下に示す。 (10) ′ (11) ここで，：全体座標系弾性応力：全体座標系弾性ひずみ E : 回転後の全体座標系弾性定数テンソル ′: 回転前の結晶系弾性定数テンソル G：回転マトリクスである。上式を成分表示すると，式(12)のように示される。本解析では2 次元の平面ひずみ場を考慮するので， 0とすれば，結局，次式により表される。 ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ (13) 式(3)，(11) ，(12)の成分を比較すると，立方晶の弾性直交異方性では，d11，d12，d44の3 つのパラメータを用いれば剛性マトリクスを記述できることがわかる。弾性剛性行列の成分は様々な材料で調べられており。Ni 基超合金γ/γ’（870℃）での値を以下に示す。本研究ではこの値を結晶系弾性定数テンソルの成分として用いた。

Table 1 Elastic stuffiness coefficients of Ni base super alloy [3] Elastic stuffiness coefficients d11=201.5 GPa d12= 137.1 GPa d44=98.5 GPa 2.2 塑性異方性（Lankford の r 値）の定義塑性異方性をモデリングには，Hill の 2 次形式の降伏関数[4]を用いて，降伏条件 f を変化させた。用いた降伏条件fを式(14)に示す。ここで， (i,j=x,y,z)は真応力，は相当応力，F，G， H，L，M，Nは材料の異方性を定義する定数である。式(13)において，3F=3G=3H=L=M=Nのとき等方性材料となる。一般的に，材料の塑性異方性を表す特性値としてはLankford の r 値[5]と呼ばれる実験値が用いられる。このLankford の r 値と式(14)との関係は式(15)で表される。 ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ (12)





















3 2 1 2 2 2 2 2 2 2 2 2 2                       H G F xy N zx M yz L yy xx H xx zz G zz yy F f (14)

























_ _    90 1 0 1 2 1 45 , 90 1 , 0 1 r r r N r F r G 　　 (15)

(4)

式(15)中の，r0,r45, r90がLankford の r 値である。これらの値は，実験的に異方性の主軸に対して傾角 ＝0°,45°,90°を有する 3 本の試験片を用いて決定される。本解析では，荷重負荷方向に対して 0°,45°,90° 方向のシュミット因子を各結晶粒について計算し，その比を用いてLankford の r 値と仮定した。シュミット因子の計算は結晶座標系で行う。荷重負荷方向Floadを結晶座標系に変換するには，前節と同様の回転マトリクス[G]を用いる。具体的な計算方法を以下に示す。すべり面を(h1 k1 l1)，すべり方向を[u1 v1 w1]で表せるとする。荷重負荷方向 Floadが試料座標系のｘ軸[100]方向であるとき，Floadは結晶座標系では列ベクトル[1 0 0]に対して(16)式のように計算できる。

  

G 100



h2 k2 l2





u2 v2 w2



load F  　　 　　  　　 (16) このとき，Floadに対するシュミット因子は式(17)で表される。 FCC 構造では 12 のすべり系(すべり面 4 個・すべり方向3 個)があるので，12 個のすべり系のうち最も大きい値をその結晶でのシュミット因子とする。同様に，荷重負荷方向に対して0°,45°,90°のシュミット因子F0，F45，F90はそれぞれ次式で定義して計算を行った。

  

100, ₄₅

  

110, ₉₀

  

010 0 G 　　　F G 　　　F G 　　 F    (18) 最終的に用いるLankford の r 値は，各結晶粒で次式のように表される。 : : 1: : (19) 上記のようにして，結晶粒ごとに塑性異方性を考慮した。Fig.2 に試料座標系から結晶座標系への回転の概念図を示す。また，Fig.3 に荷重方向 F に対するシュミット因子の概念図を示す。 2.3 有限要素モデル解析モデルは平滑クリープ試験片を想定した。実験データ（平滑材き裂発生前後の結晶方位測定結果）を再現するため，12mm×4mm の解析範囲の境界以外を全て正三角形となるような離散化を行った。節点数は7254，要素数は 14260 である。解析に用いたメッシュ図をFig. 4 に示す。 Fig. 4 FEM モデルと弾塑性クリープ解析境界条件また,実験データを Fig.5 に示す。本データは，研究の都合上詳細は延べられないため，以後本論文では A 材と呼ぶこととする。Fig.5(a)は，試験前の EBSD 観察結果，Fig.5(b)は，き裂発生後の EBSD 観察結果である。Fig.5(b)では，試験片中央部に 3 か所の結晶粒界に沿ってき裂が発生しているのがわかる。Fig.5(b)は，破断直前（寿命比 t/tf=0.99）の観察結果である。EBSD の計測間隔は 5m とし,x 軸方向に2471 点，y 軸方向に 860 点の計測を行った。この結果を基に，数値解析に用いる物性値を定めた。 ) 2 2 2 2 2 2 )( 2 1 2 1 2 1 ( 2 1 2 1 2 1 ) 2 2 2 2 2 2 )( 2 1 2 1 2 1 ( 2 1 2 1 2 1 cos w v u w v u w w v v u u l k h l k h l l k k h h coa F S                 ₍₁₇₎ Fig.2 試料座標系から結晶座標系への回転の概念図 Fig.3 荷重方向F に対するシュミット因子の概念図

(5)

具体的には，各結晶粒の中心方位をオイラー角(1， ，2)として取得し，これを用いて式(3)の回転行列を計算して物性を決定し，有限要素法に用いた。解析に用いた多結晶モデルをFig.6 に示す。本解析のモデルは40 の結晶粒から構成され，それぞれの結晶粒にはFig.5(a)を基に，オイラーを与えた。

(a)試験前 IPF map

(b) t/tf=0.99 IPF map Fig.5 A 材の EBSD 観察結果（IPF Map）

Fig.6 Model A（A 材を模擬）

Fig. 7 弾塑性クリープ解析のフローチャート 2.4 解析の流れとクリープ解析構成測解析のフローチャートを Fig. 7 に示す。解析は EPIC-Ｉ[5]を改良した，弾塑性クリープ FEM プログラムを用いた[6]。このプログラムでは，クリープ損傷領域を明確にするため，塑性変形開始後にクリープ変形則が適用される。塑性変形における加工硬化係数Hp(MPa)は，n乗硬化則を相当塑性ひずみで偏微分を行うことで得られる次式を用いた[5] 。塑性領域での構成則は次式で与えられる。 1 1 1

(

)

1 





n p p

n

c

H





(20) ここで，Hp(MPa)は加工硬化係数n1, c, αは定数，



pは相当塑性ひずみである。降伏条件はVon-Mises の降伏条件に従うものとした。次にクリープ変形におけるクリープ変形抵抗係数 Hc_{(MPa)および時間増分}_t_{(hr)は，それぞれ Norton} 則とPrandle-Reusse 則を用いて得られる式(21)および式(22)を用いた。       _          1 1 1 2 2 2 2 1 1 n p n p c ) C t ( A n d d H    ₍₂₁₎ 2 n p A Δ Δt    (22) ここで，C2, n2およびA(MPa-n_hr-1_{)は定数である。} 本解析では時間の関数であるクリープ変形則から導かれるクリープ軟化則を剛性マトリクスに組み込んでいる。一般の汎用有限要素法解析ソフトでは，クリープ変形の効果はクリープひずみによる見かけの節点荷重として剛性方程式に与えられている[7]。この手法では，クリープ変形則は解析領域全体に入 START 塑性要素降伏条件終了条件クリープ解析塑性解析弾性解析クリープ変形則加工硬化則 yes yes no no STOP yes no

(6)

力条件として適用される。一方で，形状等の要因により応力集中部が存在した場合，クリープ損傷域は部材全面に一様に生じず，最大せん断応力方向や結晶粒界等に局所的に生じることが確認されている[8]。したがって，本解析では塑性変形開始後にクリープ変形則を剛性マトリクスに適用することでクリープ損傷により引き起こされる剛性低下を表現し，局所的なクリープ損傷を再現させている。 2.5 塑性仕事塑性仕事は次式で定義される。 pc

WP







(23) c p pc







(24) ここで (MPa)は相当応力，は相当塑性ひずみ，は相当クリープひずみである。この塑性仕事の分布は，EBSD 解析により得られる Grain Reference Orientation Deviation (GROD)値と定性的に一致しているため，本研究でも解析結果の評価指標として用いた。ここでGROD 値とは，測定点の方位と計測点と同一粒内の参照方位（平均方位など）との方位差の値である。 2.6 物性定数解析にはIN738LC(800℃)及び CMSX-3 (850℃) の物性を用いた。2.1 節で示した以外の物性定数を以下に示す。また材料が脆性材料であることから，解析は平面ひずみの2 次元問題として解析を行った。

Table 2 Mechanical property of nickel-based super alloy for numerical analysis Yield stress 500 MPa

work hardening coefficient c

6.5 GPa Work hardening exponent n1 1.0

Creep coefficient A 4.17×10-19_MPa-4.8_hr-1 Creep exponent n2 4.8 3．解析結果クリープ領域が十分拡大した状態での弾塑性クリープ要素の分布および塑性ひずみエネルギーWP の分布をFig.10 に示す。Fig.10(a)では，灰色で示された領域（要素）は，弾性変形範囲である。解析結果から，数個の結晶粒を残してクリープ変形が進んでいる様子がわかる。 Fig.10(b)からは，ほぼすべての結晶粒界で，塑性ひずみエネルギーWP の分布が不連続となっており， WP の分布は粒内でも均一でないことがわかる。なお，弾性要素に関してはWP=0 として表示している。さらに実験（Fig.5（b））と比較すると，弾性要素とクリープ要素の境界がき裂発生粒界にほぼ一致していることが分かる。弾性異方性を考慮しない解析モデルでは見られない，本手法の特徴である有意な結果が得られた。平滑材の高温クリープ変形を考えた場合，試験片内に異方性がなければ，平行部では一様にクリープ変形が進み，特定の粒や粒界にひずみが集中するとは考えにくい。しかし実際には，解析結果と同様に，結晶粒ごとの異方性に加えて，結晶粒や結晶粒界の形状に依存して，塑性およびクリープ変形が不均一に進行すると考えられる。本研究では，数値解析的にこれらを予測可能な解析手法を提示できたものと考えられる。本解析で行ったように，降伏条件やクリープ発生条件に関わる多様な応力状態，特に多軸応力状態を評価することで，き裂発生が容易な結晶粒界のより正確な予測が可能と考えられる。ところで，Fig.5(b)では，割れ発生粒界に隣接する粒内に滑り変形によると思われる結晶方位変化が見て取れる。この方位変化が，割れ発生の後に生じたのであれば，クリープ変形ではなくき裂の成長による塑性変形が方位変化の主たる要因と考えることが妥当である。今後，このような結晶内の方位変化の発生時期やその量を確認することで，解析手法の改良に役立てることができるであろう。また，本研究ではクリープによる劣化，すなはちの材料の物性変化は式(20)及び式(21)として取り入 □：クリープ要素 ■：塑性要素 ■：弾性要素 (a) 弾塑性クリープ要素分布 (b) 塑性ひずみエネルギーWP 分布 Fig.10 Model A の解析結果

(7)

れているが，その異方性が明確でないため数値解析では考慮できていない。このためクリープ変形中の方位変化と数値解析結果との比較は控えた。き裂発生以前のEBSD解析結果と比較することでクリープ変形則の異方性を定量的に予測することも可能と考えられる。 3．結言本解析では，弾性異方性を考慮することによって結晶学的な方位差のみならず結晶粒界の形状に依存した局所応力場をより正確に計算することが出来た。この数値解析によりき裂が発生しやすい粒界は，塑性ひずみWP の分布および弾塑性クリープ変形が不連続となっている結晶粒界であることが予想され，実験とよく一致した。予め試験片の微細組織に関する情報を得ておけば，き裂発生や損傷量の大まかな予測が可能であり，実験での損傷観察結果と併用することでクリープ損傷形成過程をより詳細に明らかにすることが可能であると考えられる。参考文献

[1] W. F. Hosford, The mechanics of crystals and textured polycrystals, Oxford University Press, New York, 1993.

[2] 鈴木清一，EBSD 読本，株式会社 TSL ソリューションズ．

[3] T. M. Pollock and A. S. Argon., Acta Metall Mater 42 (1994)1859-1874.

[4] W.T. Lankford, S.C. Snyder and J.A. Bausscher, Trans. ASM, 42,(1950)1197-1205.

[5] Y.Yamada, “Plasticity and Viscoelasticity”, Baifukan, (1972)pp180(in Japanese).

[6] H. Takeuchi, A. T. Yokobori Jr., S. Hosono, D. Kobayashi and K. Sato, Journal of the Japan Institute of Metals, 71(5),(2007)452-457 [7] 矢川元基，宮﨑則幸，有限要素法による熱応力・

クリープ・熱伝導解析 (1985), サイエンス社 [8] A.T. Yokobori,Jr., R. Sugiura, H. Takeuchi, G.

Ozeki, A. Ishida, D. Kobayashi and S. Hosono, Strength, Fracture and Complexity, Vol. 8, No. 1 (2013), pp. 25-44.

Ni 基超合金における結晶異方性を考慮したクリープ損傷シミュレーション