Mechanical Stability Analysis on Ideal Body-Centered Cubic Crystals under Finite Deformation

Synopsis:

This study investigates mechanical stability of ideal b.c.c. crystals (α-Fe, Mo and Nb) under finite deformation. Strain energy density W is introduced using a lattice dynamics model with the second- and third-order elastic constants measured by experiments. Critical stress and strains are determined from the failure of strong ellipticity condition for stress equilibrium equation. Present numerical analysis for (001) in-plane deformation revealed that mechanical stability of the crystals depends on the deformation mode including the sign. It also revealed that the number and direction of emerging characteristic curves vary with the critical point reflecting on the nature of inter-atomic bonding. We also conducted the stability analyses for simple shear deformations on the principal planes: (100), (110) and (111). The critical shear stress show semi-quantitative agreement with those reported by first-principle calculation. Average of the shear stress over shear modulus ratio turns out to be 11%, which is close to the displacement burst condition measured by nano-indentation experiments.

1. 緒言

一般的な結晶性金属材料に外力を負荷すると，変形初期には線形弾性理論で近似可能な微小変形が進行するが，負荷応力が初期欠陥を持つ材料の降伏点に到達すると，転位等の格子欠陥の運動を介した塑性変形へと移行する¹⁾。これは通常の金属材料が示す典型的な弾・塑性変形過程を示しているが，ここでもし外力を受けた結晶内部に転位等の格子欠陥が存在しない場合には(いわゆる完全結晶状態を仮定した場合には)，その変形機構は大きく様変わりすると考えられる。すなわち，負荷応力が降伏点に到達しても塑性変形へは移行できず，更に負荷を続けると非線形弾性理論によって記述される有限変形が生じる。更に変形が進行して応力がある臨界点に到達すると，結晶は力学的な不安定化を引き起こし，転位等の格子欠陥の発生・伝播を伴ったバースト的な塑性変形が進行する。完全結晶が示すこうした力学的な振る舞いは，これまでナノインデンテーションを用いた実験や，転位論に基づく理論解析，および分子動力学法や離散転位動力学法を用いたコンピューターシミュレーションによって研究が進められてきた^2,3,4,5,6,7)。これらの研究は完全結晶の持つ臨界点(臨界応力・臨界ひずみ)の予測や，臨界後の不安定変形機構について多くの有益な知見をもたらすものの，そこから普遍的な知見を導き出して体系化するためには，力学原理に基づいた数理解析が必要と考えられる。

完全結晶の示す不安定化現象の数理解析には，古典的にはBorn条件が適用されてきた。Born条件は，金属材料を含む固体物質全般の力学的不安定性を連続体力学の枠組みで記述した臨界条件であり，具体的には，完全結晶の力学的な安定性は結晶中を伝播する音速の正値性と等価と考えるものである⁸⁾。同様の臨界条件としては，Coleman-Nollの凸条件，Baker-Ericksen不等式，強楕円条件等が知られている^9,10)。後述するように，強楕円条件は応力の平衡方程式(力の釣り合い方程式)が持つ数学的な性質から導かれる臨界条件であり，解の正則性やひずみエネルギー密度のrank-one凸性に直結するという明瞭な特徴を持つだけでなく^9,10,11)，線形弾性体に対するBorn条件の非線形弾性体への自然な拡張と捉えることができるため，物理的にも重要な意味を持っている。実際，強楕円条件を臨界条件とした固体材料の力学的安定性解析はこれまでにも多数報告されており，その大半は等方性を有する非線形弾性体に対する解析結果と^{12,13,14,15,16,17)}，異方性を有する線形弾性体に対する解析に分類することができる^18,19,20)。完全結晶の力学的不安定化機構を解析するためには，異方性を有する非線形弾性体への解析が必要となるが，これまでのところ両者を同時に取り入れた強楕円条件の解析は行われていない。そこで本研究では，この解析を可能とする新しい数値計算方法を提案するとともに，それを用いてB.C.C.系完全結晶の力学的安定性解析を行うことを目的とする。異方性を有する非線形弾性体の変形は数学的に複雑で，それを解析な形で表現することは一般に困難であるものの，適切な数値計算方法を考案すればその解析は可能であり，また上記の臨界点の導出や臨界後の不安定変形機構についても一定の知見を導くことができると考えられる。

2. 解析方法

2・1　非線形弾性理論

解析対象は参照状態において単連結で一様，かつ客観性公理を満足した立方晶系の三次元弾性体Ωとする。いまΩ内の位置ベクトルをx_i∈Ω，ラグランジュ表記された変形ベクトルをφ_i=φ_i(x_i)，変位ベクトルをu_i=φ_i−x_iとし，参照状態における弾性体の密度をρ=const.とする。また，この弾性体に働く体積力の効果は無視できるものとする。

変形φ_iによって弾性体Ω内に生じる変形勾配テンソルF_ijは，変位ベクトルの勾配u_i,jを用いて次のように表すことができる。   

$$F_{ij}=\frac{\partial {\phi }_i}{\partial x_j}=\frac{\partial \left(u_i+x_i\right)}{\partial x_j}=\frac{\partial u_i}{\partial x_j}+{\delta }_{ij}=u_{i,j}+{\delta }_{ij}$$(1)

ここでδ_ijはクロネッカーのデルタ記号を表している。通常，弾性体のひずみエネルギー密度は変形勾配テンソルF_ijの関数としてW=W(F_ij)と表すことができるが，ここで(1)式よりF_ijがu_iの関数であることを考慮すると，これはW=W(u_i,j)と考えることができる。そのため，ひずみエネルギー密度の体積積分と外力のポテンシャルエネルギーの差として定義される系の全エネルギー(ギブスエネルギー)は，変位u_iの汎関数として次のように表すことができる。   

$$I\left(u_i\right)={\int }_{\Omega }W\left(u_{i,j}\right)dV-{\int }_{\partial \Omega }{g}_i{u}_{i}d\Gamma$$(2)

ここで上式右辺第一項は変形に伴って弾性体に蓄積されるひずみエネルギー項(ヘルムホルツエネルギー)を，同第二項は変形に伴って減少する外力g_iのポテンシャルエネルギー項を表している²¹⁾。弾性体の変分原理から，Ω内に生じる変位u_iはこの全エネルギー汎関数を最小化することが保証されている。このための必要条件は，式(2)に対する次のオイラー・ラグランジュ方程式が成り立つことである。   

$$\frac{\partial }{\partial x_j}\frac{\partial W}{\partial F_{ij}}=\frac{\partial T_{ij}}{\partial x_j}=\text{div}{T}_{ij}=0\left(in\Omega \right)$$(3)

ここでT_ij=∂W/∂F_ijはピオラ・キルヒホッフの第一応力テンソルを表している。通常，(3)式は応力の平衡方程式と呼ばれており，力学的には参照状態上で記述された力の釣り合い方程式と理解することができる。

2・2　解析モデル

非線形弾性体のひずみエネルギー密度Wには様々な形式が存在するが，本研究では格子力学研究で多用されている次の形式を用いてB.C.C.結晶の構成式を定義する。   

$$W=\frac{1}{2!}{C}_{ijkl}{E}_{ij}{E}_{kl}+\frac{1}{3!}{C}_{ijklmn}{E}_{ij}{E}_{kl}{E}_{mn}+O\left(E_{ij}^4\right).$$(4)

ここでE_ijは次式で定義されるグリーン・ラグランジュひずみテンソルである。   

$$E_{ij}=\frac{1}{2}\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)+\frac{1}{2}\frac{\partial u_k}{\partial x_i}\frac{\partial u_k}{\partial x_j}$$(5)

(4)式で表されたひずみエネルギー密度を採用するのは次の理由による。まず(5)式から明らかなように，ひずみテンソルE_ijはその右辺第二項に幾何学的な非線形項u_k,iu_k,j/2を備えている¹⁰⁾。これは線形弾性理論で用いられるコーシーひずみテンソルϵ_ij(同右辺第一項)の修正項であり，これが臨界点の解析に必要不可欠な有限変形(非無限小変形)の取り扱いを可能としている。またこの項の存在によって，(4)式で定義されるひずみエネルギー密度Wは客観性公理¹⁰⁾を満足する。すなわち，このひずみエネルギー密度を構成式とした場合，そこから計算されるコーシーの応力テンソルは観測座標の取り方によらず一定である。さらにこのひずみエネルギー密度は，(4)式の右辺第一項に表された力学的線形項に加えて，同第二項に力学的非線形項を備えている。線形弾性理論では力学的線形項しか考慮されておらず，従って原子間の相互作用力は調和近似されているが，この力学的非線形項を考慮することによって，(4)式で表されたひずみエネルギー密度は非調和(非対称)な原子間相互作用を表現することが可能となる。後述するように，この非調和効果は臨界点の符号依存性として本解析結果に現れる。(4)式におけるひずみの二次係数C_ijklと三次係数C_ijklmnは，それぞれ二次および三次の弾性定数と呼ばれており，本研究で解析の対象とする3種類のB.C.C.結晶(α-Fe, Mo, Nb)については実測値が報告されている(Table 1)²²⁾。したがって，実験で求められたこれらの弾性定数を(4)式へ代入し，得られたひずみエネルギー密度を用いて応力の平衡方程式(3)を導出すれば，個々のB.C.C.結晶の力学的性質を反映した非線形弾性変形の支配方程式が定義される。

Table 1. The second- and third-order elastic constants of B.C.C. and F.C.C. crystals. The units are GPa.

	C11	C12	C44	C111	C112	C123	C144	C155	C456
α-Fe	230	135	117	–2705	–626	–575	–836	–531	–721
Mo	465	163	109	–3557	–1333	–617	–269	–893	–555
Nb	245	132	28	–2560	–1140	–467	–343	–168	137

2・3　強楕円条件

非線形弾性理論では，弾性体のひずみエネルギー密度Wとピオラ・キルヒホッフの第一応力テンソルの間には関係式T_ij=∂W/∂F_ijが成立するが，この両辺を再度F_ijで微分すると次の4階のテンソルが定義される。   

$$\frac{{\partial }^{2}W}{\partial F_{ij}\partial F_{kl}}=\frac{\partial T_{ij}}{\partial F_{kl}}={\hat{C}}_{ijkl}$$(6)

この4階テンソル${\hat{C}}_{ijkl}={\hat{C}}_{ijkl}\left(F_{ij}\right)$は，変形の進行に伴ってその値が連続的に変化する性質を有しており，いわば変形状態において弾性体が持つ実効的な変形抵抗率と考えられる。この4階テンソル${\hat{C}}_{ijkl}$を用いると，先に導出した応力の平衡方程式(3)は次のように表すことができる。   

$${\hat{C}}_{ijkl}\frac{{\partial }^2{u}_k}{\partial x_j\partial x_l}=0$$(7)

この連立非線形偏微分方程式は，任意の単位ベクトルm_iおよびn_iについて，4階テンソル${\hat{C}}_{ijkl}$が次の不等式を満足するときに楕円型となる^9,10,11)。   

$${\hat{C}}_{ijkl}{m}_i{n}_j{m}_k{n}_l>0$$(8)

この微分方程式の持つ性質から不等式(8)は強楕円条件と呼ばれており，この条件が成立するとき応力の平衡方程式(7)の解にはC²級の正則性(すなわち変位の二階微分が存在しかつ連続である)が保証される。線形弾性体の場合，4階テンソル${\hat{C}}_{ijkl}$は変形勾配F_ijによらず常に二次の弾性定数C_ijklと等しく，変形の進行に伴ってその値が変化することはない。そこでいま，等方な線形弾性体を仮定して弾性定数をC_ijkl=λδ_ijδ_kl+μδ_ikδ_jl+μδ_ilδ_jkとおけば(ここでλとμはラメ定数)，(8)式が成立するための必要十分条件はλ+2μ>0かつμ>0となる。これは弾性体内を伝播する縦波と横波の音速が正であることを要請するBorn条件に一致する。このため，強楕円条件はBorn条件の非線形弾性体に対する自然な拡張として理解することができる。

通常，線形弾性体の弾性定数はBorn条件を満足しており，従って強楕円条件は常に満足されている。非線形弾性体の場合，変形初期には${\hat{C}}_{ijkl}$はC_ijklとほぼ等しく，したがって強楕円条件が成立するものの，変形の進行に伴って${\hat{C}}_{ijkl}$値が連続的に変化してある臨界点に到達すると，不等式(8)が不成立となって強楕円条件が崩壊する。その結果，平衡方程式(7)の解の持つC²級の正則性が損なわれ，その変形勾配F_ijには区分的な不連続性が発生する。この不連続性は，弾性体の力学的な不安定化に伴って引き起こされる格子欠陥の生成を表すものと考えられている⁹⁾。本研究では，強楕円条件(8)をB.C.C.結晶の力学的安定性を司る臨界条件として設定し，これ以降の解析を進める。

2・4　数値解析法の構築

先述の通り，異方性を有する非線形弾性体を解析対象とした場合には，強楕円条件(8)式の成否を解析的に求めることは一般に極めて困難である。そこで本研究では，数値計算を用いてこの条件を評価する新しい解析方法を導入する。まず初めに，解析対象とする変形モード(変形勾配F_ij)とその範囲を設定し，次いでその変形勾配空間をΔF_ij=10⁻³間隔の単純格子点で離散分割する。このようにして作成した離散格子点を順にα=1, 2, …とすれば，各格子点上ではそれぞれ異なる変形勾配値F^α_ijが割り当てられることになる。一方，(4)式，(6)式およびTable 1を用いれば，個々のB.C.C.結晶に対する4階テンソル${\hat{C}}_{ijkl}={\hat{C}}_{ijkl}\left(F_{ij}\right)$の値をF_ijの関数として解析的に求めることができる。こうして得られた${\hat{C}}_{ijkl}$に，先に準備した変形勾配F^α_ijを代入すると，全ての格子点上で4階テンソル${\hat{C}}^{\alpha }{}_{ijkl}={\hat{C}}_{ijkl}\left(F_{ij}^{\alpha }\right)$を決定することができる。最後に，乱数を用いて10⁴~10⁶通りの単位ベクトルm_i，n_iの組を作成し，これを格子点上の4階テンソル${\hat{C}}^{\alpha }{}_{ijkl}$へ作用させて(8)式の左辺を計算する。その結果，作成した全ての単位ベクトルの組み合わせに対して計算結果が正となれば，その格子点上では(すなわちその格子点が表す変形では)強楕円条件が成立して弾性体は力学的に安定と判定し，逆にもし一つでも値が負となる単位ベクトル組み合わせが存在すれば，強楕円条件は崩壊して弾性体は力学的に不安定と判定する。この解析を全ての格子点について行えば，強楕円条件を満足する格子点の集合が，B.C.C.結晶が力学的に安定となる変形勾配の範囲を与える。なお，本数値計算法の妥当性について検証するため，既に強楕円性の崩壊条件が解析的に与えられている等方性の二次元非線形弾性体を用いた検証を行い，本数値計算結果が既存の解析解と極めて良好な一致を示すことを確認している²³⁾。本数値解析法は弾性異方性の有無によらず臨界点の解析が可能であり，またその精度は変形勾配空間の分割間隔と(本研究ではΔF_ij=10⁻³)，(8)式の評価に用いる単位ベクトルの組み合わせ数(本研究では最大で10⁶通り)にのみ依存する。等方体に対する既存の解析結果との定量的な一致は，本数値計算法が異方性弾性体へも十分適用可能であることを意味している。

3. 解析結果と考察

3・1　(001)平面内における平面ひずみ変形

まず初めに(001)平面内における膨張・収縮変形に対するB.C.C.結晶の力学的な安定性解析を行う。いま，参照状態におけるB.C.C.結晶の3つの〈100〉方向を空間の基底と一致させれば，この変形モードの変形勾配テンソルF_ijは次の形式で表すことができる。   

$$F_{ij}=\left(\begin{array}{ccc}{F}_{11}& 0& 0\\ 0& F_{22}& 0\\ 0& 0& 1\end{array}\right)$$(9)

ここで，変形量を表すパラメータはF₁₁とF₂₂の二つのみである。そこで，これら二つのパラメータで張られた二次元変形勾配空間を定義し，これを間隔ΔF₁₁=ΔF₂₂=0.001の単純正方格子点へ離散分割後，各格子点に対する強楕円条件(8)の成否を前章で説明した数値計算法を用いて判定する。

得られた解析結果をFig.1に示す。図はα-Fe，Mo，Nbにより得られた数値解析結果を表しており，図の横軸はx₁軸方向の垂直変形を表す変形勾配F₁₁を，縦軸はx₂軸方向の垂直変形を表す変形勾配F₂₂を表している。無ひずみ状態を表す参照状態はF₁₁=F₂₂=1であり，この点を含む着色された領域(青色)が強楕円条件の成立範囲である。その周囲に存在する無着色の領域は強楕円条件が崩壊した不安定領域を表しており，したがって両者の境界がB.C.C.結晶の力学的な臨界点となる。なお臨界点の中でも特徴的な点についてはA~Fまでの記号を付けている。ここではまずα-Feに対する解析結果に注目する。図中のA点はα-Fe結晶を(001)面内で一様に膨張させた際の臨界点を表しており，その臨界変形勾配はF₁₁=F₂₂=1.06となった。一方，B点は(001)面内で一様に圧縮させた場合の臨界点を表しているが，その臨界変形勾配はF₁₁=F₂₂=0.79である。これらの結果は，α-Fe結晶を(001)面内で約6%膨張させると力学的に不安定化するのに対して，一様圧縮時の不安定化は圧縮率で約21%まで生じないことを意味している。これは(001)面内での体積変化という同一の変形モードであっても，膨張変形と収縮変形というひずみの符号に応じて臨界点が異なることを意味している。こうした力学的性質の符号依存性は，(4)式で定義されたひずみエネルギー密度に含まれる力学的非線形効果に起因しており，従来の結晶異方性とは本質的に異なる非線形弾性体特有の性質である。同様に，x₁軸方向の垂直変形を拘束し(F₁₁=1として固定し)，x₂軸方向の垂直変形量F₂₂のみを変化させた際に得られた臨界点がC点とD点である。これらの臨界変形勾配はそれぞれF₂₂=1.07およびF₂₂=0.67である。この結果は，単純引張り変形では力学的不安定化が伸長率7%で生じるのに対して，単純圧縮変形時の不安定化は収縮率で約33%となることを意味している。やはり同一の変形モードであっても，引張り変形と圧縮変形というひずみの符号の相違によって顕著な異方性が発現する。Fig.1(b)，(c)から明らかなように，同様の傾向はMoおよびNbについても認められた。

Fig. 1.

 Strongly elliptic regions for (a) α-Fe, (b) Mo and (c) Nb calculated from numerical analysis on Eq.(10). (Online version in color.)

次にFig.1における強楕円領域の形状について考察する。α-Fe結晶の強楕円領域は滑らかな境界を持った三角形状の領域として表されるのに対して，Moのそれはやや扁平した楕円形状を有しており，NbはMoがより先鋭化された形状を有している。MoとNbの解析結果に共通して現れた特筆事項はE点およびF点の存在である。これらはF₁₁=F₂₂の直線に対して対称な位置に現れる臨界点上の特異点で，Moの結果からは識別が容易でないが，Nbの結果からは容易に確認できるように，強楕円領域はこれら両点の方向に張り出す傾向が認められた。このような特異点はα-Feの解析結果には存在しておらず，したがって同一のB.C.C.構造を持つ結晶であっても，その臨界点の現れ方は各元素の持つ原子間結合力の性質に依存して定性的に異なることを示唆している。臨界点における各結晶の力学的な振る舞いについてより詳細に検討するため，次節では音響テンソルを用いた特性曲線解析を行う。

3・2　臨界点と特性曲線解析

前節では(001)平面内での膨張・収縮変形に対する力学的安定性解析を行ったが，この解析で明らかにできるのは臨界点のみであり，臨界後に生じる不安定変形のメカニズムは不明である。先述の通り，平衡方程式の楕円性が崩れると変形勾配には不連続点が発生するが，この不連続点を伝播させるキャリアは特性曲線と呼ばれており，その性質は音響テンソルを用いて解析することができる。

いま(6)式で定義されたB.C.C.結晶の実効的な弾性定数を${\hat{C}}_{ijkl}$とし，n_iをある単位ベクトルとすると，このベクトルベクトルn_iに対する音響テンソルQ_ijは次のように定義される。   

$$Q_{ij}={\hat{C}}_{ijkl}{n}_j{n}_l$$(10)

ここでQ_ijは実対称テンソルである。音響テンソルQ_ijを用いると，(10)式で表された強楕円条件は次式と同値となる。   

$$\det \left|Q_{ij}\right|>0$$(11)

これよりB.C.C.結晶の力学的な臨界点はdet|Q_ij|=0と表されるが，この条件を満足する単位ベクトルn_iが力学的な不安定化後に発生する特性曲線の法線方向を表している。簡単のため，この単位ベクトルn_iを解析対象である(001)平面内に限定すると，その二つの成分はx₃軸回りの回転角θを用いて次のように表すことができる。   

$$n_i=\left(\cos \theta ,\sin \theta \right)$$(12)

通常，det|Q_ij|=0で表される臨界点は二つの変形勾配値(F₁₁とF₂₂)と回転角θに依存するが，ここで臨界点を指定して変形勾配を固定すると(例えばFig.1(a)のA点を指定するとF₁₁=F₂₂=1.06と固定される)，臨界条件det|Q_ij|=0は回転角θにのみ依存する。この臨界条件を満足するθが指定された臨界点で生じる特性曲線の数とその伝播方向を表している。

線形代数によると，臨界条件det|Q_ij|=0は実対称テンソルQ_ijの固有値が少なくとも一つゼロになることと同値である。この数学的な性質を利用して，本研究ではFig.1(a)~(c)中に付記した臨界点に対する音響テンソルQ_ijの解析を行った。得られた結果をFig.2~4に示す。Fig.2はFig.1(a)に示されたα-Feの臨界点のうち，A~D点に対する固有値の解析結果を表しており，音響テンソルQ_ijの持つ二つの固有値が回転角θの関数として極座標形式でプロットされている。まず臨界点Aの解析結果を見ると，図中の青い曲線で表された固有値は全ての回転角θにおいて非ゼロの値を有しており，従ってα-Fe結晶の力学的不安定化には関与していない。一方，赤い曲線で表された固有値は，θ=π/4，3π/4，5π/4，7π/4の4つの回転角でゼロを取っており，これが不安定化の原因であることがわかる。次に一様圧縮変形の臨界点Bに対する解析結果を見ると，やはり一方の固有値のみがA点と同じ回転角でゼロとなっている。また，(001)面内でx₂軸方向のみの膨張・収縮変形を行った際に得られた臨界点C，Dでも，本質的にはA，B点のそれと同様の傾向が認められた。これらの結果は，α-Feを(001)面内で膨張・収縮変形させた場合には，結晶が力学的に不安定化する臨界点では4つの特異点が形成され，それぞれ特性曲線によって〈110〉系統の結晶方位へと伝播することを意味している。

Fig. 2.

 Orientation dependence of eigenvalues of acoustic tensorQ_ij obtained from critical points A to D denoted in Fig.1(a). (Online version in color.)

同様の解析をMoについて行った結果をFig.3に示す。図を見ると，(001)面内の一様膨張・収縮変形を示すA，B点では4つのθ値で固有値がゼロとなっているが，その角度は両者で互いにπ/4だけずれている。この結果は，何れの臨界点も不安定化後には4つの特性曲線が形成されるものの，その伝播方向はA点では〈100〉方向となるのに対して，B点では〈110〉方向となって両者は異なることを意味している。次にC点とE点における解析結果に着目すると，固有値が消失する際のθ値は2つしかなく，そのため不安定化後に形成される特性曲線は2本のみとなる。またFig.1(b)より明らかなように，C，E点は互いに隣接している(類似した変形勾配値で不安定化している)にも関わらず，特性曲線の伝播方向がC点では[010]方向となるのに対して，E点ではそれに直交する[100]方向となっている。この結果は，臨界点における僅かな変形勾配の差が，不安定化後の変形挙動を大きく左右することを意味している。こうした傾向は先に解析を行ったα-Fe結晶では観察されておらず，従って同一のB.C.C.構造を有する結晶であっても，力学的な不安定化後の特性曲線の発生数とその伝播方向は異なる。この結果は，B.C.C.結晶を含む金属材料の塑性変形機構が，最密面・最密方向をすべり系とした転位の運動によって整理されることと対照的であり，初期欠陥を持たない完全結晶の変形機構解析を困難とする一つの本質的な要因を創り出していると考えられる。Fig.4は同様の解析をNbに対して行った結果である。Nbはα-FeよりもMoに近い挙動を示しているが，Fig.3とFig.4の細部を見比べると固有値がゼロとなる数や方向についてはいくつかの相違点が存在していることから，やはりNbもMoとは異なる変形挙動を持つと考えるのが妥当である。

Fig. 3.

 Orientation dependence of eigenvalues of acoustic tensor Q_ij obtained from critical points A to F denoted in Fig.1(b). (Online version in color.)

Fig. 4.

 Orientation dependence of eigenvalues of acoustic tensor Q_ij obtained from critical points A to F denoted in Fig.1(c). (Online version in color.)

3・3　低指数面上の単純せん断変形

次に，(100)，(110)および(111)面上での単純せん断変形モードに対するB.C.C.結晶の力学的な安定性を解析する。単純せん断変形モードを解析対象とする場合，Fig.5に模式的に示すように，変形勾配空間はせん断ひずみの大きさと方向を表す二つの変数γとφのみによって表現することができる。実際，(100)面内での単純せん断変形による変形勾配F_ij⁽¹⁰⁰⁾は次の形式で表される。   

$$F_{ij}^{\left(100\right)}=\left(\begin{array}{ccc}1& 0& \cos \phi \tan \gamma \\ 0& 1& \sin \phi \tan \gamma \\ 0& 0& 1\end{array}\right)$$(13)

Fig. 5.

 Schematic illustrations on simple shear deformation for (a)(100) plane, (b)(110) plane and (c)(111) plane. Parameters γ and φ denote the magnitude and direction of the shear deformation. (Online version in color.)

同様に，Fig.5(b)で表された(110)面内での単純せん断変形勾配F_ij⁽¹¹⁰⁾は次の形式で表される。   

$$F_{ij}^{\left(110\right)}=\left(\begin{array}{ccc}1-\frac{1}{2}\cos \phi \tan \gamma & -\frac{1}{2}\cos \phi \tan \gamma & 0\\ \frac{1}{2}\cos \phi \tan \gamma & 1+\frac{1}{2}\cos \phi \tan \gamma & 0\\ \frac{1}{\sqrt{2}}\sin \phi \tan \gamma & \frac{1}{\sqrt{2}}\sin \phi \tan \gamma & 1\end{array}\right)$$(14)

全く同様にして，Fig.5(c)で表された(111)面内での単純せん断変形勾配F_ij⁽¹¹¹⁾は次の形式で表される。   

$$F_{ij}^{\left(111\right)}=\left(F_{i1}^{\left(111\right)},F_{i2}^{\left(111\right)},F_{i3}^{\left(111\right)}\right)$$(15)

ここでF_i1⁽¹¹¹⁾，F_i2⁽¹¹¹⁾，F_i3⁽¹¹¹⁾は(111)面上のせん断変形勾配F_ij⁽¹¹¹⁾を構成する次の列ベクトルである。   

$$F_{i1}^{\left(111\right)}=\left(\begin{array}{c}1-\frac{\sqrt{2}}{6}\tan \theta \left(\sin \theta +\sqrt{3}\cos \phi \right)\\ \frac{1}{3\sqrt{2}}\tan \gamma \left(\sqrt{3}\cos \phi -\sin \phi \right)\\ \frac{\sqrt{2}}{3}\sin \phi \tan \gamma \end{array}\right)$$(16)

  

$$F_{i2}^{\left(111\right)}=\left(\begin{array}{c}-\frac{1}{3\sqrt{2}}\tan \gamma \left(\sqrt{3}\cos \phi +\sin \phi \right)\\ 1+\frac{\sqrt{2}}{6}\tan \theta \left(\sqrt{3}\cos \phi -\sin \phi \right)\\ \frac{\sqrt{2}}{3}\sin \phi \tan \gamma \end{array}\right)$$(17)

  

$$F_{i3}^{\left(111\right)}=\left(\begin{array}{c}-\frac{1}{3\sqrt{2}}\tan \gamma \left(\sqrt{3}\cos \phi +\sin \phi \right)\\ \frac{1}{3\sqrt{2}}\tan \gamma \left(\sqrt{3}\cos \phi -\sin \phi \right)\\ 1+\frac{\sqrt{2}}{3}\sin \phi \tan \gamma \end{array}\right)$$(18)

変形勾配テンソルが定まると，これ以降の解析は3・1節と同様の手順で進めることができる。すなわち，変形勾配空間を二つの変数γとφによって表された極座標系上の離散格子点に分割し(Δγ=0.001，Δφ=π/360)，次に格子点α=1, 2, …に割り振られたγとφに対する変形勾配F^α_ijを式(15)~(20)を用いて計算後，これを用いて4階テンソル${\hat{C}}^{\alpha }{}_{ijkl}={\hat{C}}_{ijkl}\left(F_{ij}^{\alpha }\right)$の全成分を計算する。最後に乱数を用いて作成した10⁴~10⁶通りの単位ベクトルm_i，n_iの組を${\hat{C}}_{ijkl}$へ作用させて(8)式の左辺を計算し，その正負によって強楕円性の成否を判定する。

(100)面上での単純せん断変形に対する解析結果をFig.6に示す。これらの図では，半径方向にせん断ひずみ量γを，円周方向に回転角φを取っており，無ひずみ状態はγ=0となる図の中心にあたる。強楕円条件の崩壊を表す臨界点は図中の赤い曲線として表されており，この曲線の内側が強楕円領域内，すなわちB.C.C.結晶が力学的に安定であり，曲線の外部が不安定領域を表している。各図の背景は参考までに示したひずみエネルギー密度を表しており，せん断ひずみ量(半径方向)の増加に伴う低エネルギー値(青色)から高エネルギー値(赤色)へのエネルギー変化と，その結晶方位依存性が見て取れる。図を見ると，(100)面内での強楕円領域は解析を行ったα-Fe，Mo，Nbの何れの場合でも4回回転対称性を有しているが，これに対してひずみエネルギー密度は回転角φに依存しない等方性を示している。よく知られているように，B.C.C.結晶の(100)面は4回回転対称性を有していることから，臨界点もこれと同一の回転対称性を持つことが要請される。この結果は，(8)式に基づく強楕円性の判定はその解析にひずみエネルギー密度を利用しながらも，ひずみエネルギー密度そのものよりも妥当な解析結果を与えることを意味している。各B.C.C.結晶の臨界点を比べると，α-FeとMoが[100]方向のせん断変形に対して臨界点の最大値を持つのに対して(臨界値はそれぞれγ=0.15とγ=0.18)，Nbの場合には[110]方向で最大の臨界点(臨界値はγ=0.12)を迎えている。したがって，力学的に安定なせん断変形方向は元素の結合様式に応じて異なっている。同様の解析を(110)面，および(111)面に対して行った結果をFig.7，8に示す。これらの解析結果を見ると，強楕円領域とひずみエネルギー密度はともに結晶面の持つ回転対称性((110)面では二回回転対称性，(111)面では三回回転対称性)と一致しており，この点は(100)面に対する解析結果と異なる。一方，臨界点が最大となる方向は元素によって異なっており，これについては(100)面上での解析結果と同様である。

Fig. 6.

 Mechanical instability of B.C.C. crystals under (100) simple shear deformation. Strongly elliptic regions are defined within the red curves in the figures. The background contour gradations show strain energy density. (Online version in color.)

Fig. 7.

 Mechanical instability of B.C.C. crystals under (110) simple shear deformation. Strongly elliptic regions are defined within the red curves in the figures. The background contour gradations show strain energy density. (Online version in color.)

Fig. 8.

 Mechanical instability of B.C.C. crystals under (111) simple shear deformation. Strongly elliptic regions are defined within the red curves in the figures. The background contour gradations show strain energy density. (Online version in color.)

3・4　検証と考察

本研究で得られた解析結果の精度について定量的な検証を行うため，前節までに得られた臨界点を第一原理計算によって求められた報告値²³⁾と比較・検討する。完全結晶の理想強度は実験的にその値が確定しておらず，したがって第一原理計算の結果は定量的な実験検証が行われたものではない。またこの計算では熱振動効果も含まれていないが，現時点では最も客観的な検証方法と考えられる。

先述の通り，ひずみエネルギー密度Wとピオラ・キルヒホッフの第一応力テンソルの間には関係式T_ij=∂W/∂F_ijが成り立つが，W=W(F_ij)は(4)式によって変形勾配F_ijの関数として解析的に表されており，したがってその一次微分であるT_ij=T_ij(F_ij)もまた変形勾配F_ijにのみ依存する。この関数T_ijにFig.6-8で求めた臨界点における変形勾配値を代入すれば，その臨界点における応力テンソルT_ijを求めることができる。しかしながら，この応力テンソルは参照座標と現座標の双方を基底に持ったツーポイントテンソルであり，弾性理論で一般的に使用されるコーシー応力テンソルとは性質が異なる。そこで本研究では，まず初めに臨界点におけるピオラ・キルヒホッフの第一応力テンソルT_ijを計算し，次いで次式で表されるピオラ変換によってこれをコーシー応力テンソルT^φ_ijへと変換した¹⁰⁾。   

$$T_{ij}^{\phi }=\frac{1}{J}{T}_{ik}{F}_{kj}^T$$(19)

臨界応力の計算は，本研究で得られた数値解析結果に対する定量的な検証を目的としている。したがって，その解析対象は必ずしもB.C.C.結晶に限る必要はなく，むしろ他の結晶構造に対する解析結果を含めることで検証効果の向上が期待できる。そこで本研究では，上記で解析を行った3つのB.C.C.結晶に加えて，F.C.C.構造を有するCu，Au，Ag，Al，Ni，およびダイヤモンド構造を有するSiについても同様の解析を行うことにした。これらの立方晶系結晶は，ひずみエネルギー密度Wの形式がB.C.C.結晶のそれと一致するため，Wの形式は固定したまま弾性定数のみ変更すれば直ちに同様の解析が可能である上に，何れの結晶も第一原理計算による理想強度の解析結果が報告されているという利点がある。

臨界応力を比較した結果をTable 2に示す²⁴⁾。α-Feの(110)・〈111〉単純せん断変形に対する臨界応力は，強楕円条件に基づく本研究の解析結果で7.62 GPa，第一原理計算による解析結果で11.43 GPaとなり，本解析結果は第一原理計算結果と比べて33%ほど低く見積もられている。一方，Moに対する解析結果を見ると，本解析結果は第一原理計算と比べて32%ほど値が大きい。F.C.C.結晶に対する解析結果には−4%から−62%程度の開きが見られ，ダイヤモンド構造を有するSiについては解析結果に約10%程度の相違が認められる。これらの結果を見ると，強楕円条件を用いた本解析結果は第一原理計算と定量的には一致しないことが解る。

Table 2. Critical shear stress of ideal crystal obtained from present numerical analysis and first principle calculation reported in Ref.23. T^φ_ij/C₄₄ and T^φ_ij/C’₄₄ show normalized shear stress by shear moduli.

	Plane	Direction	Present Study Tφij (GPa)	First Principle (GPa)	Tφij/C44 (%)	Tφij/C’44 (%)
α-Fe	(110)	<111>	7.62	11.43	6.5	16.0
Mo	(110)	<111>	21.74	16.52	19.9	14.4
Nb	(110)	<111>	6.81	–	24.3	12.1
Cu	(111)	<112>	2.11	3.45	2.8	9.0
Au	(111)	<112>	0.99	1.42	2.3	6.8
Ag	(111)	<112>	0.96	2.57	2.1	6.4
Al	(111)	<112>	1.90	3.73	6.7	8.3
Ni	(111)	<112>	6.03	6.29	4.9	12.8
Si	(111)	<110>	12.28	11.13	15.5	24.1

両者の差が生じた原因としては，先に挙げたように比較対象とした第一原理計算結果の精度問題が考えられるが，これ以外にも，解析に用いた三次弾性定数の計測精度に関する問題が考えられる。一般に，三次弾性定数の計測には結晶中を伝播する音速の外部応力依存性が用いられる²⁵⁾。単結晶試料に対するこの計測は実験が極めて困難な上に，三次弾性定数の報告値に関する精度検証は現在でも十分行われていない。こうした問題点を考えると，Table 2に見られる本研究と第一原理計算の差−62~+33%は必ずしも大きいとは言い切れず，むしろ臨界応力の推定値が同一オーダーとなることは，本解析結果が半定量的な解析精度を有していることを示唆していると考えられる。

最後に，本研究で得られた臨界せん断応力を剛性率で規格化することを考える。立方晶系の線形弾性体は結晶方位によって剛性率が異なるが，その値はC₄₄とC'₄₄=(C₁₁−C₁₂)/2という二つの剛性率の間で変化する。そこで，本研究で得られた臨界点におけるせん断応力をこれら二つの剛性率で規格化した。得られた結果をTable 2に示す。この結果を見ると，臨界せん断応力T^φ_ijは剛性率に対して2.1~19.9%の範囲で変化しているが，C₄₄で規格化した場合の平均値が9.5%，C'₄₄で規格化した場合の平均値が12.2%となっており，これより臨界せん断応力の平均値は剛性率の約11%と見積もることができる。ここでナノインデンテーションを用いた実験結果を参照したい。前述のように，格子欠陥が乏しい微小領域の力学特性をナノインデンテーション法を用いて評価すると，その荷重−変位曲線にはしばしばポップインと呼ばれる現象が観察されている。Ohmura等は，ポップイン現象が生じる臨界応力を様々な結晶について系統的に計測しているが，本研究で解析を行ったAl，Ag，Nb，Cu，Ni，FeおよびMoを含むOhmura等の実験結果を見ると，最大せん断応力の平均推定値は，元素に依らず剛性率のほぼ1割前後であることが指摘されている⁶⁾。これは本研究で見積もられた臨界せん断応力の規格化平均値と概ね一致している。この結果は，ナノインデンテーション実験において荷重・変位曲線上のポップイン現象(変位のバースト現象)は，圧子下で有限変形を受けた結晶が力学的に不安定化した結果生じたものと解釈できることを示唆している。

4. 結言

本研究では，非線形弾性理論を用いてB.C.C.完全結晶の有限変形下における力学的安定性解析を行った。得られた結果は次のようにまとめることができる。

(1)幾何学的および力学的非線形項を取り入れた格子力学型の非線形ひずみエネルギー密度Wを導入し，ここへ実験によって求めらた二次と三次の弾性定数を代入することで，α-Fe，NbおよびWに対する非線形構成関係式を決定した。

(2)得られたひずみエネルギー密度Wを用いて応力の平衡方程式を導出し，その係数テンソル${\hat{C}}_{ijkl}$を用いて平衡方程式の楕円性を評価した。この際，強楕円条件の崩壊点を求める数値計算法を新しく考案し，これを用いてB.C.C.結晶の力学的な安定性解析を行った。

(3)(001)面内変形に対する解析の結果，解析を行った全てのB.C.C.結晶において力学的安定性の符号依存性が認められた。一例として，正符号のひずみテンソルを持つ一様膨張変形は，これと同一の変形モードで負の符号を持つ一様収縮変形と比べて臨界ひずみ値が小さくなることが挙げられ，この意味で，一様膨張変形は一様圧縮変形と比べて力学的に不安定であると結論された。また，音響テンソルを導入してその固有値解析を行うことで，臨界後に形成される特性曲線の数とその方向について解析を行った。その結果，同じB.C.C.構造を有する結晶であっても，臨界後の特性曲線の発生数とその伝播方向は本質的に異なり得ることが明らかとなった。

(4)B.C.C.結晶の3つの低指数面((100)面，(110)面および(111)面)上での単純せん断変形に対する力学的安定性解析を行い，得られた臨界応力を第一原理計算による報告値と比較した結果，少なくとも両者は同一オーダーとなることが明らかとなった。この結果は，本解析手法が半定量的な計算精度を有していることを示唆している。また，ナノインデンテーション実験によって観察されたポップイン時の臨界応力と剛性率の比を，本研究によって得られた臨界せん断応力の平均値と比較した結果，両者には定量的な一致が認められた。この結果は，ポップインを引き起こす変位バースト現象が，完全結晶の力学的な不安定化によって生じたこと，またその理解には本研究で用いた数理的な解析方法が有効であることを示唆している。

謝辞

本研究の一部は大阪大学大学院の宮部菜苗氏(現・特許庁)の協力の下で行われた。記して感謝致します。

文献

• 1)  幸田成康:金属物理学序論，コロナ社，東京，(1973).
• 2)   W.W.  Gerberich,  J.C.  Nelson,  E.T.  Lilleodden,  P.  Anderson and  J.T.  Wyrobek: Acta Mater., 44(1996), 3585.
• 3)   A.  Gouldstone,  H.-J.  Koh,  K.-Y.  Zeng,  A.E.  Giannakopoulos and  S.  Suresh: Acta Mater., 48(2000), 2277.
• 4)   Y.  Shibutani,  T.  Tsuru and  A.  Koyama: Acta Mater., 55(2007), 1813.
• 5)   T.  Tsuru,  Y.  Shibutani and  Y.  Kaji: Acta Mater., 58(2010), 3096.
• 6)   T.  Ohmura,  L.  Zhang,  K.  Sekido and  K.  Tsuzaki: J. Mater. Res., 27(2012), 1742.
• 7)   T.  Ohmura: Materia Jpn., 53(2014), 312.
• 8)   M.  Born and  K.  Huang: Dynamical Theory of Crystal Lattice, Oxford University Press, Oxford, (2002).
• 9)   J.K.  Knowles and  E.  Sternberg: Arch. Rat. Mech. Anal., 63(1976), 321.
• 10)   P.G.  Ciarlet: Mathematical Elasticity, Vol.1, North-Holland, Amsterdam, (2004).
• 11)   B.  Dacologna: Direct Methods in the Calculus of Variations, Springer, New York, (2008).
• 12)   R.C.  Abeyaratne: J. Elasticity, 10(1980), 255.
• 13)   P.J.  Davies: J. Elasticity, 26(1991), 291.
• 14)   H.C.  Simpson and  S.J.  Spector: Arch. Rat. Mech. Anal., 84(1983), 55.
• 15)   L.  Zee and  E.  Sternberg: Arch. Rat. Mech. Anal., 83(1983), 53.
• 16)   J.  Wang,  J.  Li,  S.  Yip,  S.  Phillpot and  D.  Wold: Phys. Rev. B, 52(1995), 12627.
• 17)   C.  Padovani: Meccanica, 37(2002), 515.
• 18)   S.  Chirita,  A.  Danescu and  M.  Ciarletta: J. Elasticity, 87(2007), 1.
• 19)   J.R.  Walton and  J.P.  Wilber: Int. J. Nonl. Mech., 38(2003), 441.
• 20)   D.  Han,  H.  Dai and  L.  Qi: J. Elasticity, 97(2007), 1.
• 21)  加藤雅治:入門転位論，裳華房，東京，(1999).
• 22)   A.G.  Every and  A.K.  McCurdy: Low Frequency Properties of Dielectric Crystals, Springer, Berlin, (1992).
• 23)  宮部菜苗:大阪大学修士論文，(2014年度).
• 24)   S.  Ogata,  J.  Li,  N.  Hirosaki,  Y.  Shibutani and  S.  Yip: Phys. Rev. B, 70(2004), 104104.
• 25)   R.N.  Thurston and  K.  Brugger: Phys. Rev., 133(1964), A1604.