(57)【要約】【課題】表面側から与えられたエネルギーに対する固体の挙動を正確に解析する。【解決手段】前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成する。そして、この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析する。【選択図】図１３
2010127783.tif 000002
【特許請求の範囲】
【請求項１】
表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析シミュレータであって、
前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、
この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析する固体挙動解析シミュレータ。
【請求項２】
請求項１に記載の固体挙動解析シミュレータにおいて、
前記固体層におけるエネルギー散逸は、表面側の原子配列に対し深さ方向に位置する複数層の原子配列における速度減衰項により計算する固体挙動解析シミュレータ。
【請求項３】
請求項２に記載の固体挙動解析シミュレータにおいて、
前記速度減衰項を表面から離れるに従って大きく設定する固体挙動解析シミュレータ。
【請求項４】
請求項１に記載の固体挙動解析シミュレータにおいて、
前記固体層におけるエネルギー散逸は、表面側の原子配列に対し深さ方向に位置する複数層が無限に続くと仮定することによって、表面側から入力されたエネルギーの反射をなくすことにより計算する固体挙動解析シミュレータ。
【請求項５】
請求項４に記載の固体挙動解析シミュレータにおいて、
前記固体層におけるエネルギー散逸は、原子間の結合を線形近似した上で計算する固体挙動解析シミュレータ。
【請求項６】
請求項４または５に記載の固体挙動解析シミュレータにおいて、
深さ方向における原子配列同士において、所定の距離ずれ、反力がある状態を初期値とし、この状態で力を印加して原子位置の時間平均の変化のなくなる定常状態を決定する固体挙動解析シミュレータ。
【請求項７】
請求項１〜６のいずれか１つに記載の固体挙動解析シミュレータにおいて、
さらに、
前記表面側から与えられたエネルギーは、固体の表面および分子についての、分子シミュレーションにより挙動を特定する分子シミュレータを有し、
この分子シミュレータによって特定されたエネルギーについて、固体相の挙動を解析する固体挙動解析シミュレータ。
【請求項８】
請求項１〜７のいずれか１つに記載の固体挙動解析シミュレータにおいて、
固体表面の摩擦現象を解析する固体挙動解析シミュレータ。
【請求項９】
コンピュータにより、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析プログラムであって、
コンピュータに、
前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、
この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析させる固体挙動解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析シミュレーションに関する。
【背景技術】
【０００２】
機械的な動力の伝達などの解析において、固体間摩擦の解析が必要となる。原子レベルから固体間摩擦を解析する理論手法に関しては、いくつかの簡易的なモデルを用いたものが提案されている。非特許文献１は原子スケールにおける固体間摩擦モデルとして、図１（ａ）に示すような剛体の壁に原子がバネで接続されたもの（Ｔｏｍｌｉｎｓｏｎモデル）を考案した。この一次元原子列をある周期ポテンシャルの上で直線運動させることで、摩擦の解析が行われた。非特許文献２，３では、結晶転移モデルなどでよく使われるＦｒｅｎｋｅｌ−Ｋｏｎｔｒｏｖａモデル（ＦＫモデル）を原子レベルの固体間摩擦解析に用いた。このモデルは図１（ｂ）のように原子同士がバネで結合し一次元鎖を形成したものである。また、非特許文献４では、より現実的な固体表面モデルとして、図１（ｃ）に示すように表面原子は互いに結合し、さらに剛体壁にも接続されているモデルを用いている。このモデルはＦｒｅｎｋｅｌ−Ｋｏｎｔｒｏｖａ−Ｔｏｍｌｉｎｓｏｎモデル（ＦＫＴモデル）、もしくはｇｅｎｅｒａｌｉｚｅｄ−ＦＫモデルなどと呼ばれる。Ｔｏｍｌｉｎｓｏｎモデル（図１（ａ））、ＦＫモデル（図１（ｂ））およびＦＫＴモデル（図１（ｃ））により固体間摩擦のメカニズムが論じられている。
【０００３】
また、動的分子シミュレーションによって、固体表面における動的物理現象を解析することについては、特許文献１〜３などに記載がある。
【０００４】
【特許文献１】特開２００２−１５０２０３号公報
【特許文献２】特開平８−３５８３１号公報
【特許文献３】特開平１０−２３２８６０号公報
【非特許文献１】Ｇ．Ａ．Ｔｏｍｌｉｎｓｏｎ，  Ｐｈｉｌｏｓ．  Ｍａｇ．  ７（１９２９）  ９０５
【非特許文献２】Ｙ．Ｉ．Ｆｒｅｎｋｅｌ  ａｎｄ  Ｔ．Ｋｏｎｔｏｒｏｖａ，  Ｚｈ．Ｅｋｓｐ．  Ｔｅｏｒ．  Ｆｉｚ．  ８，  １３４０，  １９３８．
【非特許文献３】Ｋ．Ｓｈｉｎｊｏ  ａｎｄ  Ｍ．Ｈｉｒａｎｏ，  Ｓｕｒｆ．  Ｓｃｉ．  ２８３，  ４７３，  １９９３．
【非特許文献４】Ｈ．Ｍａｔｓｕｋａｗａ  ａｎｄ  Ｈ．Ｆｕｋｕｙａｍａ，  Ｐｈｙｓ．  Ｒｅｖ．  Ｂ４９，  １７２８６，  １９９４．
【発明の開示】
【発明が解決しようとする課題】
【０００５】
上述した非特許文献１〜４に示される既存モデルは全て表面原子のみが陽に扱われ、固体バルク原子の存在は全く考慮されていない。現実の固体表面は固体バルク原子が存在するため、既存モデルのように簡易化された表面では、原子スケールにおける摩擦機構を十分には捕らえられない。
【０００６】
また、既存モデルのもう一つの欠点は、格子振動によるエネルギー散逸の記述の不完全さである。実際の摩擦現象において、摺動により表面原子に与えられたエネルギーは格子振動を誘起する。そして、その格子振動は固体（結晶）バルクへと伝搬しエネルギーは散逸されていく。散逸したエネルギーは外部と熱交換され、固体中の不純物を考慮しない限り表面に再帰することは無い。ＦＫＴモデルなどではこのエネルギー散逸の記述を原子速度に比例する抵抗力で表現するが、この速度依存の抵抗力には任意的なパラメータが用いられており、物理的な妥当性の無い項になっている。摩擦によるエネルギーロスは、固体のエネルギー散逸機構が担うため、エネルギー散逸の記述の不完全さは摩擦モデルの不完全さと直結すると考えられる。
【０００７】
このように既存モデルは、固体バルク原子の不考慮、格子振動によるエネルギー散逸の記述の不完全さ、という重大な欠点を内包している。また、特許文献１〜３においても、動的分子シミュレーション（分子動力学、第一原理計算、動的モンテカルロ法）によって固体表面における動的物理現象の解析が行われているが、動的過程によって発生するエネルギーの散逸項を、従来では気体や流体で用いられる速度スケーリングや熱統計平均による熱浴によって代用してきた。しかしながら、動的過程によって表面に与えられたエネルギーは固体の格子振動励起によって固体内部に散逸されるため、気体や流体で用いられる速度スケーリングで表面のエネルギー散逸を記述することに物理的正当性は無く、例えば、摩擦現象においては摩擦係数が正しく算出されないと考えられる。
【０００８】
ここで、従来技術で挙げた表面一原子のみを扱った既存モデルが、非常に大きな固体バルク原子数Ｎをもつ現実的固体表面とどの程度一致しているかを検証する。２つの表面の摺動により格子振動の縦波と横波が表面平行方向と表面垂直方向（バルク方向）それぞれに誘起される。今、表面が平坦で無限に広いとするならば明らかに表面平行方向に伝搬する縦波・横波はエネルギー散逸に寄与しない。なぜなら図２（ａ）に示すようにある表面原子に誘起された振動エネルギーが表面平行方向に伝搬しても、他の表面サイトで誘起された振動エネルギーが流入してくるからである。摺動により与えられた振動エネルギーは図２（ｂ）のようバルク方向に伝搬することでエネルギー散逸となる。なお、図２は、表面に与えられた振動エネルギーの拡散の概念図であり、矢印は振動エネルギーの伝搬の方向を表す。また、図２（ａ）は、振動エネルギーが表面平行方向へ伝搬する場合、図２（ｂ）振動エネルギーがバルク方向へ伝搬する場合を示している。
【０００９】
このエネルギー散逸による表面原子の振動の減衰挙動を線形近似によるバルク方向に伸
びた１次元Ｎ体振動子を用いて解析する。外力が無いときのこの振動子の運動方程式は
【数１】
2010127783.tif 000003
である。ここで、ｕはｘとｚを表し、Ｋｕはｘ、ｚ方向の変位に対するバネ定数でそれぞれ横波、縦波の周波数に関係する。この連成振動子の最表面原子（ｉ＝１）にのみ変位Ｘを与えた時、つまり初期条件
【数２】
2010127783.tif 000004
の振動解を導出する。振動解は以下のようになる。
【数３】
2010127783.tif 000005
【００１０】
有限Ｎの保存系では摺動により誘起される振動はある一定時間経過すると反射して表面原子まで戻ってくる。しかしながら、純粋な結晶系では摺動により与えられた振動エネルギーはバルクに拡散していき、外界との熱交換により失われ表面原子にエネルギーが戻ってくることは無い。また、現実的な固体表面のＮは１０１０オーダーと非常に大きい。この状況を表現するためにはバルク側を無限に大きな原子集団として、Ｎ→∞解を求めれば良いことになる。Ｎ→∞の振動解析解は次のように求めることができる。
【００１１】
【数４】
2010127783.tif 000006
【００１２】
ここで、Ｊｎはｎ次の第一種ベッセル関数である。
【００１３】
図３に、解析解とＦＫＴモデルとの比較を示す。図３（ａ）に表面原子の変位、図３（ｂ）にエネルギー時間微分の時間発展を示す。ここでは、Ｋｕ＝１，ｍ＝１，γ＝１．１を用いた。図３（ａ）に表面原子の変位ｕ１の解析解と、同様の初期条件１を与えたＦＫＴモデルでの解を示す。ここで、Ｘ＝１，Ｋｕ＝１，ｍ＝１を用い、ＦＫＴモデルではγ＝１．１を使用した。
【００１４】
また、表面原子のエネルギーを、Ｅ＝（１／２）（ｄｕ１／ｄｔ）２＋（１／２）（ｕ１−ｕ２）２で定義し、その時間微分でエネルギーの散逸を示したものを図３（ｂ）に示す。最表面原子のエネルギーがバルクに伝搬することにより減衰していく様子が分かる。ＦＫＴモデルは、変位、エネルギー時間微分の時間発展ともに解析解と明らかにずれている。この結果は、表面一原子のみを扱っている既存モデルを表面モデルとして扱うことは不適切であることを示している。また、γ＝１．１以外を使用してもこのずれは修正されない。
【課題を解決するための手段】
【００１５】
本発明は、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析シミュレータであって、前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析することを特徴とする。
【００１６】
また、前記固体層におけるエネルギー散逸は、表面側の原子配列に対し深さ方向に位置する複数層の原子配列における速度減衰項により計算することが好適である。
【００１７】
また、前記速度減衰項を表面から離れるに従って大きく設定することが好適である。
【００１８】
また、前記固体層におけるエネルギー散逸は、表面側の原子配列に対し深さ方向に位置する複数層が無限に続くと仮定することによって、表面側から入力されたエネルギーの反射をなくすことにより計算することが好適である。
【００１９】
また、前記固体層におけるエネルギー散逸は、原子間の結合を線形近似した上で計算することが好適である。
【００２０】
また、深さ方向における原子配列同士において、所定の距離ずれ、反力がある状態を初期値とし、この状態で力を印加して原子位置の時間平均の変化のなくなる定常状態を決定することが好適である。
【００２１】
さらに、前記表面側から与えられたエネルギーは、固体の表面についての、分子シミュレーションにより挙動を特定する分子シミュレータを有し、この分子シミュレータによって特定されたエネルギーについて、固体相の挙動を解析することが好適である。
【００２２】
また、固体表面の摩擦現象を解析する挙動とすることが好適である。
【００２３】
また、本発明は、コンピュータにより、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析プログラムであって、コンピュータに、前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析させることを特徴とする。
【発明の効果】
【００２４】
本発明によれば、固体における深さ方向の伝達によるエネルギーの散逸を考慮する。このため、表面側から与えられたエネルギーに対する固体の挙動をより正確に解析することができる。
【発明を実施するための最良の形態】
【００２５】
以下、本発明の実施形態について、図面に基づいて説明する。
【００２６】
材料開発において、成膜におけるスパッタリングや蒸着およびポリッシング、分析におけるＨｅ散乱、物性評価においてはナノインデンテーションや摩擦試験などが行われている。しかし、これらの実験のみでは起こっている現象を詳しく解析できない。このため、これら固体表面における動的物理現象に対して原子、分子スケールでのシミュレーションをする必要がある。
【００２７】
このシミュレーションは、汎用のコンピュータを用いて行われ、予め求められたモデルに対し、各種のデータを入力して動作を模擬することによって行われる。従って、モデルとして、物理現象を正しく表現したものが好ましい。
【００２８】
本実施形態は、このようなモデルを用いたシミュレーションを行うシミュレータおよびそのプログラムに関し、ＣＤ、ＤＶＤ、通信などによって、コンピュータの記憶部に前記モデルを利用したシミュレーションを行うプログラムをインストールし、プログラムを実行し、各種条件、データを入力して、摩擦現象や、摩擦係数などの物性値を求めるためのシミュレーションを行う。なお、シミュレーションプログラム自体は、すでに販売されているものを用い、そのシミュレーションプログラム内のモデルを本実施形態のプログラムによって設定することも好適である。
【００２９】
そして、シミュレーションは、現実の固体材料について行われ、材料開発などに用いられる。例えば、各種の原子、原子配列の固体についての、摩擦係数を求めることで、所望の摩擦係数を得るために、どのような原子、原子配列の材料がよいかを求めることができ、現実の材料を特定することができる。例えば、自動車の動力伝達機構であるクラッチや、ギアなどにおける材料にどのような原子構造のものを採用するのがよいかなどを知ることができる。
【００３０】
「Ａ  エネルギーバッファースラブモデル（Ｎ有限モデル）」
上述のように、既存モデルでは、連成振動子系の解析解と振動の挙動が大きく異なる。そこで、本実施形態では、より解析解と近くなる振動挙動を示すモデルを利用する。まず、バルク原子をＮ個連ねたスラブモデルにより表面を作る。減衰項は、ＦＫＴモデルなどと同様に速度に比例する形で入れるが、表面原子に近づくにつれ減衰項を小さくする。速度依存の減衰係数は上述の式３［数４］から分かるように解析振動解の時間依存性が√（Ｋｕ／ｍｔ）であることから、この時間依存性を持つように要請する。上記をふまえ以下の運動方程式で表されるモデルを用いる。
【００３１】
【数５】
2010127783.tif 000007
【００３２】
なお、以下において、このモデルをエネルギーバッファースラブモデル（ＥＢＳモデル）と呼ぶ。なお、このモデルは、Ｎ＝１の時はＦＫＴモデルと等価になる。γは減衰項の係数であり、Ｋｕやｍといった物性値依存性は無い。
【００３３】
図４に、Ｎ→∞における解析解、ＦＫＴモデルおよびＥＢＳモデルの表面原子に初期変位Ｘを与えた際の振幅（変位）（図４（ａ））、エネルギー時間微分の時間変化（図４（ｂ））を示す。ここでもＸ＝１，Ｋｕ＝１，ｍ＝１，γ＝１．１を使用した。ＥＢＳモデルのバルク原子数は、Ｎ＝２０を用いた。
【００３４】
ＥＢＳモデルは、ＦＫＴモデルと比べ振幅、エネルギー時間微分ともに解析解と近くなることが分かる。次に、有限Ｎにより振動の反射が表面原子まで戻ってくることがないかを検証した。結果を図５に示す。なお、図５も解析解、ＥＢＳモデル、連成振動子モデルとの比較であり、連成振動子モデルはＥＢＳモデルでγ＝０と置いたものとする。また、初期変位Ｘ＝１を与えたときの表面原子の変位（図５（ａ））と（エネルギー時間微分（図５（ｂ））を示している。なお、図４と同様に、Ｋｕ＝１，ｍ＝１，γ＝１．１、ＥＢＳモデルのバルク原子数はＮ＝２０を用いた。
【００３５】
図５より、γ＝０のＥＢＳモデル、すなわち単なる連成振動子モデルの場合、ｔ＝３５経過後に振動が反射している。一方、ＥＢＳモデルの場合、解析解と同様に振動の反射はほとんど起こっていないことが分かる。また、当然ながら、図４，５で示した結果は、ＥＢＳモデルのＮ数を大きくすればするほどより解析解に近くなる。
【００３６】
本実施形態のＥＢＳモデルは、ＦＫＴモデルと比べ実際の表面の格子振動によるエネルギー散逸を正確に模擬する。このことから、ＥＢＳモデルは摩擦によって引き起こされる格子ダイナミクスをよりリアリスティックに表現できるモデルであることが分かる。
【００３７】
「Ｂ  動摩擦係数の比較」
次に、ＥＢＳモデルで摩擦係数を導出する。周期的境界条件を課すことにより、図６に示すように、摺動系として下方に立方格子のＥＢＳモデル、上方には格子間隔を固定した立方格子を配置した。両者は、同じ原子間距離σを持つ（つまり、格子整合）。ＥＢＳモデルの物性値は、ダイヤモンドを想定して定めた。
【００３８】
すなわち、ｍ＝１２［質量数］、原子間距離σ＝１．５２［Å］。バネ定数はダイヤモ
ンド［１１１］表面で第一原理計算を行い、Ｋｘ＝０．０９，Ｋｚ＝０．１１［ａ．ｕ．］を得た。簡単のため、Ｋｘ＝Ｋｚ＝０．１［ａ．ｕ．］を計算で用いた。上面と下面との相互作用ＦＩは簡単に斥力のみを考えたソフトコア斥力を用いた。相互作用ＦＩは、次のように表される。
【００３９】
【数６】
2010127783.tif 000008
【００４０】
ここで、ｒは上面原子と下面の表面原子の距離を表す。斥力の係数Ｋは、第一原理計算によりダイヤモンドの結合Ｃ−Ｃボンドを縮ませたときの斥力によって定めた。運動方程式は、４次のｒｕｎｇｅ−ｋｕｔｔａ法を用いて解き、タイムステップは０．５７［ｆｓｅｃ］とした。上面は１０［ｍ／ｓｅｃ］で摺動方向に等速直線運動をさせて摺動を行った。面圧Ｐおよび動摩擦係数μは以下の式によって算出される。
【００４１】
【数７】
2010127783.tif 000009
ここで＜＞ｔは時間平均を表す。
【００４２】
ＥＢＳモデルでバルク原子数Ｎを１〜１０００まで変化させた時の摩擦係数と面圧の関係を図７に示す。Ｎの増加に比例して摩擦係数は上昇し、Ｎ＝３００で摩擦係数は収束する方向にある。Ｎ＝１の場合、つまりＦＫＴモデルの場合、摩擦係数はほとんど０となる。このように、摩擦係数はバルク原子の運動まで考慮することで大きく変化する。当然実際の表面で起こっていることはバルク原子まで関与した格子ダイナミクスであるため、ＦＫＴモデルは動摩擦係数の議論に適さず、本実施形態に係るＥＢＳモデルによる解析方法が有効であることが分かる。
【００４３】
以上のように、本実施形態におけるエネルギーバッファースラブ（ＥＢＳ）モデルは、バルク側原子のダイナミクスまで考慮しており、固体間摩擦における動摩擦係数の解析ができる。このＥＢＳモデルの表面原子に与えられた振幅とエネルギー減衰はバルク原子数Ｎ→∞の解析解とよい一致を示し、従来の表面一原子のみを扱うＦＫＴモデルよりもリアリティのあるモデルである。また、摺動計算により動摩擦係数を算出したところ、ＦＫＴモデルではほぼ０摩擦という不自然な結果を与えるのに対し、ＥＢＳモデルでは有限の動摩擦係数を与えられる。これらの結果は、ＥＢＳモデルのようにバルク原子の運動まで考慮した格子ダイナミクスを考えなければ摩擦係数の解析が行えないことを示している。
【００４４】
「Ｃ  半無限連成振動子方程式（Ｎ無限大モデル）」
上述のＥＢＳモデルでは、速度依存の抵抗力が与える表面ダイナミクスへの影響を少なくするため、抵抗力の係数を固体（結晶）バルク内部から表面に向かうにつれ小さく設定した。そして抵抗力の係数を、バルク原子数Ｎ→∞の系における表面原子のエネルギー散逸の時間変化とフィッティングすることで定めた。しかしながら、本来、速度依存の抵抗力は人工的な項であるので、それを除くために固体表面を図８のように結晶バルク方向に原子をＮ個連ねた連成振動子モデルで表し、新たに半無限系（Ｎ→∞）の運動方程式を解く方法を考案した。このＮ→∞により表面から与えられた摺動エネルギーは無限の彼方まで伝わっていき、再び表面に戻ってこないので現実的なエネルギー散逸を再現できる。これは無限に大きなバルクの影響が熱浴として振る舞うと考えても良い。さらに、現実の固体においてＮは１０１０オーダーと非常に大きいため、このＮ→∞操作は実際の摩擦現象をよく再現すると考えられる。
【００４５】
以下に、半無限連成振動子モデルの運動方程式が、形式的に解け、それにより表面原子の運動を記述できる、本実施形態に係る新しい微積分方程式が得られることを示す。
【００４６】
図８（左）に連成振動子モデルの模式図を示す。ユニットセルは表面垂直方向にＮ個、表面平行方向にはＮ’個の原子（合計Ｎ×Ｎ’個）が連成した正方格子である。そのユニットセルに対し表面平行方向に周期的境界条件を課すことで、表面を形成する。整数η，ξは原子指標で、それぞれη＝１，２，・・・Ｎ’，ξ＝１，２，・・・Ｎの範囲をとる。ξ＝１は表面原子を表わし、結晶バルク側の原子に進む方向に沿ってξを番号付けしている。ξ＝Ｎの原子は剛体壁に接続されている。原子間の結合は線形近似で表す。表面平行方向および垂直方向に伸びたバネ定数をそれぞれ（Ｋ’ｘ，Ｋ’ｚ），（Ｋｘ，Ｋｚ）とする（図８（右）参照）。ここでＫ（’）ｘ，Ｋ（’）ｚは、表面平行方向および垂直方向の平衡原子間距離からの変位に対するバネ定数である。ｘ，ｚを表面平行方向および垂直方向に対する原子の平衡位置からの変位とする。外部から作用する力は表面原子にのみ作用すると仮定する。その時、原子変位ｕη，ξ＝ｘη，ξ，ｚη，ξの運動方程式は次のようになる。
【００４７】
【数８】
2010127783.tif 000010
【００４８】
ここで、Ｆｘ，Ｆｚはそれぞれ表面原子ξ＝１にのみ作用する力の表面平行方向および表面垂直方向成分である。本モデルを下面と上面に配置すれば、Ｆｕは二固体間の相互作用力を表す。ｍは原子の質量を表す。表面平行方向に対する周期的境界条件は、
【数９】
2010127783.tif 000011
である。
【００４９】
式４（数８）をベクトル表示すると、次のようになる。
【数１０】
2010127783.tif 000012
2010127783.tif 000013
【００５０】
ここで∧Ｂ（Ｂハット），∧Ｃ（Ｃハット）は、Ｎ×Ｎ×Ｎ’×Ｎ’の４階テンソル、∧ＣはＮ’×Ｎ’の２階テンソル（行列）である。∧Ｂ，∧Ｃは表面垂直方向および面内の質点間のバネの相互作用力項を表している。∧１はＮ’×Ｎ’の２階テンソルの単位元である。面間の相互作用力Ｆｕは表面原子ξ＝１に働き、
【数１１】
2010127783.tif 000014
である。
【００５１】
方程式６（数９）の解は斉次解ｕ０ξと４階テンソル形式のＧｒｅｅｎ関数∧Ｇ（Ｇハット）を用いて次式（８）のように書き下せる。
【数１２】
2010127783.tif 000015
【００５２】
後に論ずる式展開と数値解の容易さのため、上式を時間で微分した速度に対する解が有用となる（式９）。
【数１３】
2010127783.tif 000016
【００５３】
ここで∧ｇ＝ｄ／ｄｔ∧Ｇである。Ｇｒｅｅｎ関数∧Ｇと斉次解ｕ０ξはそれぞれ以下の式を満たす。
【数１４】
2010127783.tif 000017
【００５４】
つまり、運動方程式６を解くには、式１０および式１１の２つの方程式により∧ｇとｄ／ｄｔ  ｕ０を求め、方程式９に代入して解けば良い。以下∧ｇおよびｄ／ｄｔｕ０を求める。
【００５５】
Ｃ−１∧ｇの導出
方程式１０にｅｘｐ（−ｓｔ）をかけてラプラス変換を施すと、式（１２）となる。
【数１５】
2010127783.tif 000018
【００５６】
ここで、∧Ｂ，∧Ｃを対角化する行列を∧ＵＢおよび∧Ｕｃとする。すなわち、次式となる。
【数１６】
2010127783.tif 000019
【００５７】
また、∧ＵＣを次のように定義する。
【数１７】
2010127783.tif 000020
【００５８】
方程式１２（数１２）の両辺に∧ＵＣ∧ＵＢをかける。∧ＵＢと∧ＵＣは可換であることを使えば（∧ＵＢのξ，ｊ成分はＵξ，ｊ∧１であるから）、∧ＵＣ∧ＵＢ∧Ｇ（ｓ）≡〜Ｇ（ｓ）のη，ｌ，ξ，ｊ成分は以下の式（１３）を満たす。
【００５９】
【数１８】
2010127783.tif 000021
ここで（ωＢｊ）２＝（Ｋｕ／ｍ）λＢｊおよび（ωＣｌ）２＝（Ｋ’ｕ／ｍ）λＣｌである。∧Ｂ，∧ＣはエルミートであるからＵＢ，ＵＣはユニタリ行列である。よって両辺に（ＵＢ）−１（ＵＣ）−１＝（ＵＢ）ｔ（ＵＣ）ｔを作用させ、さらにラプラス逆変換を施せば∧Ｇ（ｔ）のη，ｌ，ξ，ｊ成分は、次式となる。
【００６０】
【数１９】
2010127783.tif 000022
【００６１】
よって∧ｇ（ｔ）は、式（１４）となる。
【数２０】
2010127783.tif 000023
【００６２】
相互作用力Ｆｕは、表面原子層ξ＝１にのみ作用すると仮定したため、∧ｇの行列要素として必要なのはξ’＝１成分のみである。そのためｇη，η’，ξ，ξ’＝１≡ｇη，η’，ξと表記し直す。
【００６３】
ωＢｊ，ＵＢξ，ｊは次のように表される。
【数２１】
2010127783.tif 000024
【００６４】
ここで、θＢ＝π（ｊ−１／２）／（２Ｎ＋１），α＝２√（Ｋｕ／ｍ）である。ｊの取りうる範囲はｊ＝１，２，・・・，Ｎである。式１４に式１５，１６を代入すれば、次式が得られる。
【数２２】
2010127783.tif 000025
【００６５】
ここで、バルク方向の原子数Ｎ  の無限大極限をとる。Ｎが十分大きければ、次式となる。
【数２３】
2010127783.tif 000026
【００６６】
Ｎ→∞によりｊに関する和を積分に変換し、θＢ→θ、ωＢｊ→ω（θ）＝αｓｉｎθと書き換えると、式１７が得られる。
【数２４】
2010127783.tif 000027
【００６７】
ωＣｌ，ＵＣη，ｌは次のように表される。
【数２５】
2010127783.tif 000028
【００６８】
ここで、α’＝２√（Ｋ’ｕ／ｍ）である。上式を式１７に代入し、ｓξ，ｌ（ｔ）＝ｓξ，−ｌ（ｔ）であることを使えば、式１９（数２６）が得られる。
【数２６】
2010127783.tif 000029
【００６９】
特に、ｌ＝０の時、ｓξ，ｌ（ｔ）の式１８の積分はｎ次の第一種Ｂｅｓｓｅｌ関数Ｊｎの積分表示に帰結し、式２０で表現できる。
【数２７】
2010127783.tif 000030
【００７０】
ここで、第一種Ｂｅｓｓｅｌ関数の積分表示πＪ２ξ（ｚ）＝∫ｄθｃｏｓ（ｚｓｉｎθ）ｃｏｓ２ξθ［θ＝０：π］＝２（−１）ξ∫ｄθｃｏｓ（ｚｃｏｓθ）ｃｏｓ２ξθ［θ＝０：π／２］を用いた。
【００７１】
ここで、ｇη，η’，ξの物理的意味を議論しておく。原子（η’，１）に作用する力Ｆ（ｔ）が原子（η，ξ）の速度に与える影響は、次式で表される。
【数２８】
2010127783.tif 000031
【００７２】
力Ｆは、関数ｓとの畳み込み積分により、バルク方向に伝わる波と面内方向に伝わる波に分解される。そして、ｅｘｐ（２πｉｌ（η−η’）／Ｎ’）は、面内方向へ伝搬する波の位相を表現している。その伝搬していく様々な周波数の波を足し合わせたものが∫ｇη，η’，ξ（ｔ−ｔ’）Ｆ（ｔ’）ｄｔ’［ｔ’＝０：ｔ］となる。つまり、周囲の原子への力の伝搬を表すのがｇη，η’，ξである。
【００７３】
Ｃ−２（ｄ／ｄｔ）ｕ０の導出
有限Ｎ系であれば摺動計算において時間発展させれば、ある一定時間後には初期条件によらず定常状態を得られ、その時の動摩擦係数を求めることができる。よって、最も簡単な初期条件ｕ０ξ（ｔ＝０）＝０，（ｄ／ｄｔ）ｕ０ξ（ｔ＝０）＝０における斉次解、つまりｕ０ξ（ｔ）＝０を用いて差し支えない。しかしながら、この初期条件では無限Ｎ系で定常状態を得ることができない。このことを簡単な例を用いて示す。
【００７４】
初期条件ｕ０ξ（ｔ＝０）＝０，（ｄ／ｄｔ）ｕ０ξ（ｔ＝０）＝０にある表面原子ξ＝１に一定の力Ｆが作用したとする。十分長い時間ｔ＊経過すれば、力の釣り合いからｕξ＝１（ｔ＊）＝ＸｕＮ，Ｘｕ＝Ｆ／Ｋｕの定常状態になる。有限Ｎの系ではｕξ＝１（ｔ＊）は有限の値であり、時間発展によって定常状態に到達できる。しかしながらＮ→∞の系ではｕξ＝１（ｔ＊）は発散してしまう。これは一つの原子間の結合のひずみが小さくとも、Ｎが無限であるので発散するためであり、無限系特有の症状といえる。つまり無限系では、初期条件ｕ０ξ（ｔ＝０）＝０，（ｄ／ｄｔ）ｕ０ξ（ｔ＝０）＝０から時間発展させもｕξ＝１は定常状態に到達することができない。
【００７５】
無限系で定常状態を達成するには、あらかじめ初期条件として原子間に一定のひずみＸｕを与えておき、それが摺動による相互作用力の時間平均と釣り合うように決定する手法を用いる。その手法を以下に詳しく述べる。
【００７６】
まず式１１の初期条件を設定する。初期時間ｔ＝０における初期条件を、表面垂直方向に配列した原子間はある一定ひずみＸｕ＝ｕ０η，ξ−ｕ０η，ξ＋１を持ち、面内原子間のひずみは起きずｕ０η，ξ＝ｕ０η’，ξとなって静止している、と設定する。これは結晶格子が表面原子にある一定かつ一様な力を受けて平衡状態になった配置で、ひずみを持ったまま静止している状態である。このとき面内の原子間ひずみは起きないので∧Ｃは消え、原子の運動にη依存性はなくなる。よってｕ０ξ＝ｕ０ξ（左側の項はベクトル（太字）ｕ））と表示し直すと、式１１は以下のようになる。
【００７７】
【数２９】
2010127783.tif 000032
【００７８】
以下、式２１を解いていく。式２１の両辺に∧ＵＢをかけ、〜ｕ０ｊ＝ΣξＵＢξ，ｊｕ０ξとおくと、次式となる。
【数３０】
2010127783.tif 000033
【００７９】
上式は簡単に解くことができ〜ｕ０ｊ（ｔ）＝ｃｊｃｏｓωＢｊｔを得る。ただし、（ｄ／ｄｔ）〜ｕ０ξ（ｔ＝０）＝０を用いた。
【００８０】
未定係数ｃｊは初期変位に関する条件ｃｊ＝ΣＵＢξ，ｊＸｕ（Ｎ−ξ＋１）［ξ＝１：Ｎ］から決定される。
【００８１】
式１６を用いれば、次式が得られる。
【数３１】
2010127783.tif 000034
【００８２】
公式Σｓｉｎ２ξθＢ［ξ＝１：Ｎ］＝ｓｉｎ（Ｎ＋１）θＢｓｉｎＮθＢ／ｓｉｎθＢにより、次式を得る。
【数３２】
2010127783.tif 000035
【００８３】
よってｕ０ξ（ｔ）は、次式となる。
【数３３】
2010127783.tif 000036
【００８４】
ここで、（ｄ／ｄｔ）ｕ０ξ（ｔ）は、次のように変形できる。
【数３４】
2010127783.tif 000037
【００８５】
最後は、Ｎが非常に大きいことを使った。Ｎ→∞を行いｊに関する和を積分に変換し、ωＢｊ→ω（θ）＝αｓｉｎθと書き換えると、次式を得ることができる。
【数３５】
2010127783.tif 000038
【００８６】
ここで、最後の等号は、第一種Ｂｅｓｓｅｌ関数の積分表示を用いた。この（ｄ／ｄｔ）ｕ０ξ（ｔ）を用いて定常状態を得る方法に関しては次の項で述べる。
【００８７】
Ｃ−３半無限連成振動子方程式と自己無撞着計算
項Ｃ３−１，Ｃ３−２の結果により、（η，ξ）原子の運動方程式の形式解は、式２３のように記述することができた。
【数３６】
2010127783.tif 000039
【００８８】
ここで、α＝２√（Ｋｕ／ｍ）、α’＝２√（Ｋ’ｕ／ｍ）、ω（θ）＝αｓｉｎθ、ωＣｌ＝α’ｓｉｎπｌ／Ｎ’である。この新しく得られた微積分方程式を数値的に解けば、半無限連成振動子系の原子の運動を求めることができる。また今回は２次元に配列した結晶格子で定式化を行ったが、同様の式展開を行えば３次元に配列した結晶格子においても方程式を得ることができる。
【００８９】
ところで、動摩擦係数μ、単位面積あたりの動摩擦力Ｐ／／、垂直圧力Ｐ⊥は、次のように求めることができる。
【数３７】
2010127783.tif 000040
【００９０】
ここで＜＞ｔは定常状態での時間平均操作、Ｓは表面原子η＝１，２，・・・Ｎ’が構成する面積を表す。
【００９１】
外部との相互作用力は表面原子層ξ＝１のみに作用するため、摩擦解析に必要な物理量を得るには、表面原子ξ＝１の運動が分かれば良い。もしも、運動方程式６自体により表面原子ξ＝１の運動を解こうとすれば、Ｎ→∞漸近解が扱えないのは当然であるが、バルク原子ξ≠１も同時に解かねばならない。一方、式２３のように形式的に運動方程式を解いたものを用いれば、ξ≠１を解かずに原子ξ＝１の運動を得ることができる。よって解くべきξ＝１の微積分方程式は、次のようになる。
【００９２】
【数３８】
2010127783.tif 000041
【００９３】
さて、前項の最後に述べたように、半無限系において定常状態を実現するＸｕを求めなくてはならない。まず、表面垂直方向のゆがみＸｚは初期条件として適当な値に設定する。Ｘｚは垂直荷重の反作用力に起因するゆがみであり、定常状態における垂直圧力を決定するパラメータとなる。つまり問題は、あるＸｚ（＝垂直圧力）および摺動速度Ｖが与えられた時、どのようなＸｘが定常状態を実現するか、ということである。ｔが十分に長いとき、Ｊｎ（αｔ）→０，∫Ｊｎ（αｔ’）ｄｔ’［ｔ’＝０：ｔ］→１／αが成り立つため、式２５のｘに関する方程式は以下のようになる。
【００９４】
【数３９】
2010127783.tif 000042
【００９５】
定常状態であるとき、＜（ｄ／ｄｔ）ｘη，ξ＝１（ｔ）＞ｔ＝０となることを使えば、定常状態を実現するＸ＊ｘは、次の関係式を満たす。
【数４０】
2010127783.tif 000043
【００９６】
ここで、本計算法を用いた具体的な計算手順を示す。（１）摺動速度Ｖ，Ｘｚ，Ｘｘを設定する。（２）微積分方程式２５により表面原子の運動ｕη，ξ＝１（ｔ）を解く。（３）式２６の右辺を計算したものをＸ＊’ｘと置く。｜（Ｘ＊’ｘ−Ｘｘ）／Ｘｘ｜＝δＸｘ＜Δであれば定常状態を達成したとし、Ｘ＊ｘ＝Ｘ＊’ｘを得る。ここでΔは非常に小さい値とする。δＸｘ＞Δであれば定常状態に達していないとして、βＸ＊’ｘ＋（１−β）Ｘｘを新しいＸｘとして手順（１）に戻す。ここで０＜β＜１である。
【００９７】
上記の手順（１）〜（３）を繰り返していけば、定常状態を得ることができる。これは、定常状態を実現する初期変位Ｘ＊ｘを求める自己無撞着計算と言える（参考：例えば密度汎関数法による第一原理計算では、基底エネルギーを実現する電子密度を求める自己無撞着計算が行われる）。そして、そのＸ＊ｘにおける軌道ｕη，ξ＝１（ｔ）を式２４に用いれば、動摩擦係数μ、単位面積あたりの動摩擦力Ｐ／／、垂直圧力Ｐ⊥といった摩擦係数を得ることができる。
【００９８】
「Ｄ  計算方法」
本節では２固体が同一の格子定数を持つ摺動系に対し、考案した半無限連成振動子方程
式を用いて計算を行った。図９に本摺動計算系の模式図を示す。下面は表面垂直方向に無限に原子が配列し、表面平行方向にも原子が周期的に配列した２次元正方格子である。原子間の結合は線形近似で表し、表面垂直方向に伸びたバネ定数を（Ｋｘ，Ｋｚ）とする。ここでＫｘ，Ｋｚは表面平行方向および垂直方向の平衡原子間距離からの変位に対するバネ定数である。上面を構成する原子は伸縮しない（無限に固い）結合で結ばれて配列している。
【００９９】
上面と下面は同じ格子定数を持つ物質と仮定する。この場合、表面平行方向の周期性により面内の原子は全て同じ運動をするため、面内の振動は起きず表面平行方向に伸びた結合の寄与は無い。
【０１００】
時間ｔ＝０において、上面を一定速度Ｖで表面平行方向に等速直線運動させ摺動を開始する。下面のある表面原子に着目し、その表面原子のｔ＝０における位置を座標原点に置き、着目した表面原子の位置座標をｕ＝ｘ，ｚとする。ここでｘ，ｚはそれぞれ表面平行方向、表面垂直方向の座標を表す。位置座標ｕに作用する上面からの相互作用力を（Ｆｘ（ｘ，ｚ，ｔ），Ｆｚ（ｘ，ｚ，ｔ））とすれば、前節より半無限連成振動子方程式（Ｎ’＝１，ｌ＝０）は、次式となる。
【０１０１】
【数４１】
2010127783.tif 000044
【０１０２】
ここで、α＝２√（Ｋｕ／ｍ）、ｍは原子の質量である。Ｊｎはｎ次の第１種Ｂｅｓｓｅｌ関数である。Ｘｕはｔ＝０におけるｚ方向に配列した全ての原子間の一定ゆがみである。すなわち、原子間距離σとし、ｚ方向に隣り合わせたある２つの原子の位置座標をｕξおよびｕξ＋１とすれば、ｘξ−ｘξ＋１＝Ｘｘ、ｚξ−ｚξ＋１＝Ｘｚ＋σである。
【０１０３】
まず、Ｘｚをある負の値に設定し、下面全体を圧縮する。Ｘｚは垂直荷重に起因するゆがみであり、定常状態における垂直圧力を決定するパラメータとなる。次に、Ｘｚ（＝垂直圧力）および摺動速度Ｖを与えた時の、定常状態を実現する初期表面平行ゆがみＸ＊ｘを求めなくてはならない。なぜなら結晶バルク原子が無限に配列しているため、Ｘｘ＝０で摺動させても格子全体のゆがみは無限の時間歪み続け、永遠に定常状態が得られないからである。定常状態のとき、（ｄ／ｄｔ）ｘ（ｔ）の時間平均が０となることを使えば、定常状態を実現するＸ＊ｘは、次の関係式を満たす。
【０１０４】
【数４２】
2010127783.tif 000045
【０１０５】
ここで＜＞ｔは時間平均操作を表す。
【０１０６】
式２７、２８を用いた摺動計算手順は次の通りである。（ｉ）摺動速度Ｖ，Ｘｚ，Ｘｘを設定する。（ｉｉ）式２７により表面原子の運動ｕ（ｔ）をある一定時間解く。（ｉｉｉ）式２８の右辺を計算したものをＸ’ｘと置く。｜（Ｘ’ｘ−Ｘｘ）／Ｘｘ｜＝δＸｘ＜Δであれば定常状態を達成したとし、Ｘ＊ｘ＝Ｘ’ｘを得る。ここでΔは非常に小さい値とする。δＸｘ＞Δであれば定常状態に達していないとして、βＸ’ｘ＋（１−β）Ｘｘ→Ｘｘとして手順（ｉ）に戻す。ここで０＜β＜１である。上記の収束計算の手順（ｉ）〜（ｉｉｉ）を繰り返していけば定常状態が得られる。
【０１０７】
本計算では、摺動速度Ｖを１０ｍ／ｓｅｃに設定し、固体モデルのパラメータはダイヤモンドの物性値を用いた。すなはち、ｍは炭素の質量、σはダイヤモンドの格子定数２．８８［ａ．ｕ．］である。また第一原理計算によりダイヤモンド結晶のバネ定数を見積もり、およそＫｘ＝０．１，Ｋｚ＝０．１［ａ．ｕ．］を用いた。今回、上面原子と下面原子との相互作用力には着目しないため、簡単な斥力のみのソフトコアポテンシャルで表現した。ある上面原子ｉと下面原子間の距離をｒｉ、ｘ方向となす角度をθｉとすると、相互作用力は、次式である。
【０１０８】
【数４３】
2010127783.tif 000046
【０１０９】
ここで、和はｒｉ＜σを満たす上面原子全てに対して行う。斥力は非現実的なオーダーにならなければ良いので、ここではダイヤモンド原子間の斥力を第一原理計算で求め、それにフィッティングしてＫ＝０．１［ａ．ｕ．］、Ｍ＝５と定めた。
【０１１０】
定常状態を実現するＸ＊ｘの収束計算は動摩擦係数の誤差が１％以下になるまで繰り返した。βは０．７５程度を用いれば、広範囲の面圧で安定な収束が得られた。動摩擦係数μ、単位面積あたりの摩擦力Ｐ／／、垂直圧力Ｐ⊥は次のように求めることができる。
【０１１１】
【数４４】
2010127783.tif 000047
【０１１２】
ここで、Ｓは表面１原子が構成する面積を表し、Ｓ＝σ２である。
【０１１３】
「Ｅ  計算結果・考察」
図１０に面圧と摩擦係数の関係と結晶バルク原子層の数Ｎとの依存性を示す。縦軸は摩擦係数、横軸は面圧を表す。Ｎ→∞の摩擦係数は半無限連成振動子方程式を用い、Ｎ有限の摩擦係数はＥＢＳモデルによって算出した。Ｎ＝１、すなわち、ＦＫＴモデルにおいて摩擦係数は面圧に関わらずほぼ０摩擦を示す。Ｎを増加させるに従い摩擦係数は増大していき、Ｎ→∞で０．２５程度の有限の摩擦係数に収束する。Ｎによるこの著しい摩擦係数の差は、固体間摩擦現象において結晶バルク原子が本質的役割を演じていることを明らかに示している。
【０１１４】
次に、表面原子に印加された力および表面原子の変位を図１１に示す。Ｎ→∞における（ａ）摩擦力、（ｂ）表面原子の変位の時間変化およびＮ＝１（ＦＫＴモデル）における（ｃ）摩擦力、（ｄ）表面原子の変位の時間変化である。黒丸は相手面原子と接触した時間を表し、白丸は相手面原子と離れた時間を表す。Ｎ→∞の場合、表面原子が相手表面の原子から受ける摩擦力Ｐ／／（ｔ）は図１１（ａ）に示すようにノコギリ歯形状となる。また、表面原子の変位は図１１（ｂ）のように相手面原子と接触すると引きずられ、ある時点で相手面表面原子から離れ元の位置に戻っていく。摺動中において表面原子は常に運動している。Ｎ＝１の摩擦力および変位を図１１（ｃ）および（ｄ）に示す。摩擦力は正と負両方の値をとり、ｓｉｎカーブのようにほぼ対称的な挙動をとるため、ほぼ０摩擦となる。また変位も同様にｓｉｎカーブのような波形となる。また、Ｎ→∞の場合と対照的に、相手面原子と接触している時間のみ変位を起こし、接触していない時は静止していることが分かる。Ｎ→∞およびＮ＝１における摺動中の表面原子のダイナミクスの模式図を図１２に示す。このように、Ｎによって表面原子のダイナミクスは大きく変化しており、それが摩擦係数の違いに影響を与えている。図１２は、（Ａ１）〜（Ａ４）Ｎ→∞と、（Ｂ１）〜（Ｂ４）Ｎ＝１（ＦＫＴモデル）における摺動中の表面原子のダイナミクスの模式図である。
【０１１５】
Ｔｏｍｌｉｎｓｏｎは、摺動によって誘起される表面原子の非断熱運動によるエネルギー散逸が摩擦の起源であると説明した。ところが平野・新庄により、原子の固有振動数ωに対応する速度に比べ現実の摺動速度は十分遅く、非断熱運動によるエネルギー散逸は起きないとの主張により否定された。この主張は本計算結果においてＮ＝１でほぼ０の摩擦係数が得られたことからも裏づけられる。
【０１１６】
それでは、なぜＮ→∞の場合に摩擦が生じたのかは次のように説明できる。固体表面を結晶バルク原子まで考慮することにより、連成振動子を形成することになる。すると実際の表面原子の運動は単一の固有振動数ωのモードで規定されるのではなく、０〜２ωのモードを連続的にとれるのである。そのため摺動速度がどんなに遅くとも、それに対応する０付近の低振動数モードが誘起されて表面原子は運動する。そして表面原子に与えられた運動エネルギーが結晶バルクへと伝搬してエネルギー散逸となる。このように非断熱運動は低振動モードの運動により生じ、摩擦となるのである。
【０１１７】
以上の結果から示される重要な事実は、ＦＫＴモデルなどの表面を１原子のみで模したモデルは固体間摩擦現象の本質を捉えていない、ということである。そして理論の枠組みに結晶バルク原子の存在を組み込むということは、単なる摩擦係数の定量的な精度の向上という意味合いではなく、固体間摩擦においてエネルギー散逸機構を取り入れるという不可欠な操作である。
【０１１８】
半無限連成振動子方程式により、無限に結晶バルク原子が連なった表面での摩擦係数および表面原子のダイナミクスを解析した。摩擦係数は結晶バルク原子の数が多くなるほど増加し、表面原子のダイナミクスも大きく変化した。これは連成振動子系を成すことで、摺動によって表面原子の低振動モードが誘起され、それにより結晶バルク原子へエネルギーが散逸するためである
【０１１９】
「Ｆ  分子シミュレーションとの組み合わせ」
上述のように、材料開発において、成膜におけるスパッタリングや蒸着およびポリッシング、分析におけるＨｅ散乱、物性評価においてはナノインデンテーションや摩擦試験などが行われている。しかしながらこれらの実験のみでは起こっている現象を詳しく解析できないため、これら固体表面における動的物理現象に対して原子、分子スケールでのシミュレーションをする必要がある。
【０１２０】
動的分子シミュレーション（分子動力学、第一原理計算、動的モンテカルロ法）によって固体表面における動的物理現象の解析が行われているが、動的過程によって発生するエネルギーの散逸項を、従来では気体や流体で用いられる速度スケーリングや熱統計平均による熱浴によって代用してきた。しかしながら、動的過程によって表面に与えられたエネルギーは固体の格子振動励起によって固体内部に散逸されるため、気体や流体で用いられる速度スケーリングで表面のエネルギー散逸を記述することに物理的正当性は無く、例えば、摩擦現象においては摩擦係数が正しく算出されない。
【０１２１】
本実施形態では、固体の格子振動励起によるエネルギー散逸を現実的に取り扱うことを可能とする。これにより固体表面における動的物理現象、例えば結晶成長、ナノインデンテーションやＨｅ散乱などによる表面分析の有効な解析ツールが提供される。具体的には以下の２つの方法がある。
（ｉ）固体基板を表面平行方向には周期的境界条件、深さ方向に対する原子配列数Ｎを持つ原子層で構成する。そしてその構成原子の運動方程式に速度に依存する減衰項を付加することで固体のエネルギー散逸を記述する。このエネルギー散逸機能を持った固体基板上に動的分子シミュレーションによって表面層を作成したり分子の動的挙動を計算する方法。
（ｉｉ）固体基板を表面平行方向には周期的境界条件、深さ方向に対する原子配列数Ｎを持つ原子層で構成する。そして原子間の結合を線形近似した上で、原子配列数Ｎ→無限大操作（半無限連成振動子方程式）を行うことで、固体のエネルギー散逸を記述する。このエネルギー散逸機能を持った固体基板上に動的分子シミュレーションによって表面層を作成したり分子の動的挙動を計算する。より具体的には連成振動子方程式と動的分子シミュレーションの各タイムステップごとに座標データと相互作用力データをやり取りし、これらを連成して時間発展計算を行うことで実現される。例として、図１３に、連成振動子方程式と動的分子シミュレーションとの連成計算による、２固体表面間に分子が介在している系の摺動計算の模式図を示す。ここで、Ｐは垂直荷重、Ｖは滑り速度を表す。
【０１２２】
このように、液体などの分子が存在する部分と、これと接触する１層または数層については、動的分子シミュレーション（分子動力学、第一原理計算、動的モンテカルロ法）によってシミュレーションを行い、この結果として出てくる原子配列の運動を本実施形態による、エネルギー散逸を考慮した計算とする。これによって、液体も含めた、摺動系全体のシミュレーションを行うことができる。
【０１２３】
「シミュレーション」
上述のようにして、本実施形態によれば、深さ方向におけるエネルギー散逸を考慮したシミュレーションが可能となる。このためには、まず構成原子の質量、間隔、バネ定数などの結晶としての性状を決定する。そして、表面の１層に対し印加される力を決定することで、これに対する固体表面の挙動をシミュレーションすることができ、摩擦係数などを算出することができる。
【０１２４】
図１４に、実施形態における物理量（摩擦係数など）を算出するための処理（シミュレーション）についてのフローチャートを示す。
【０１２５】
まず、固体基板、表面層、介在分子の構造、原子種、格子定数および時間発展のタイムステップ、温度などを設定する（Ｓ１）。Ｓ１の入力データを利用して、動的分子シミュレーションによって表面層、介在分子の時間発展を計算する（Ｓ２）。また、Ｓ２の入力データを利用して、速度依存の減衰項を持つ基板原子層数Ｎを決定するか、もしくは半無限連成振動子方程式により時間発展を計算する（Ｓ３）。そして、Ｓ３により得られた時間発展の結果から、定常状態か否かを判定し（Ｓ４）、ＮＯであれば計算結果を利用する次のタイムステップの計算のためＳ３に戻る。
【０１２６】
ここで、時間発展の各タイムステップ毎にＳ２およびＳ３で得られた原子位置座標を用い、相互作用力を計算し（Ｓ５）、これをＳ２およびＳ３にフィードバックし、次のタイムステップの計算を行う。これによって、Ｓ２の動的分子シミュレーションによる計算結果と、Ｓ３の固体表面解析の計算結果の相互作用を考慮した計算が行える。
【０１２７】
そして、Ｓ５において定常状態と判定された場合に、目的の物理量（例えば摩擦係数）の時間平均を算出して（Ｓ６）、処理を終了する。
【図面の簡単な説明】
【０１２８】
【図１】（ａ）Ｔｏｍｌｉｎｓｏｎモデル、（ｂ）ＦＫモデル、（ｃ）ＦＫＴモデルの模式図である。
【図２】表面に与えられた振動エネルギーの拡散の概念図であり、矢印は振動エネルギーの伝搬の方向を表す。（ａ）振動エネルギーが表面平行方向へ伝搬する場合、（ｂ）振動エネルギーがバルク方向へ伝搬する場合を示す。
【図３】解析解とＦＫＴモデルとの比較を示す図である。初期変位Ｘ＝１を与えたときの表面原子の（ａ）変位、（ｂ）エネルギー時間微分の時間発展を示す。
【図４】Ｎ→∞における解析解、ＦＫＴモデルおよびＥＢＳモデルとの比較を示す図である。初期変位Ｘ＝１を与えたときの表面原子の（ａ）変位、（ｂ）エネルギー時間微分を示す。
【図５】解析解、ＥＢＳモデル、連成振動子モデルとの比較を示す図である。初期変位Ｘ＝１を与えたときの表面原子の（ａ）変位、（ｂ）エネルギー時間微分を示す。
【図６】ＥＢＳモデルを用いた摺動計算の概念図である。
【図７】ＥＳＢモデルで原子数Ｎを変化させたときの摩擦係数と面圧の関係を示す図である。
【図８】（左）は、連成振動子モデルの模式図であり、表面垂直方向にＮ個、表面平行方向にはＮ’個の原子が連成し、表面平行方向に周期的境界条件が課されている。（右）は、表面平行方向および垂直方向に伸びた原子間の結合とそのバネ定数を示す。
【図９】本摺動計算系の模式図を示す図である。
【図１０】結晶バルク原子層の数Ｎと摩擦係数との関係を示す図である。
【図１１】Ｎ→∞における（ａ）摩擦力、（ｂ）表面原子の変位の時間変化、Ｎ＝１（ＦＫＴモデル）における（ｃ）摩擦力、（ｄ）表面原子の変位の時間変化を示す。
【図１２】（Ａ１）〜（Ａ４）Ｎ→∞と、（Ｂ１）〜（Ｂ４）Ｎ＝１（ＦＫＴモデル）における摺動中の表面原子のダイナミクスの模式図である。
【図１３】連成振動子方程式と動的分子シミュレーションとの連成計算による、２固体表面間に分子が介在している系の摺動計算の模式図である。
【図１４】物理量の算出処理のフローチャートである。
【図１】
2010127783.tif 000048
【図２】
2010127783.tif 000049
【図３】
2010127783.tif 000050
【図４】
2010127783.tif 000051
【図５】
2010127783.tif 000052
【図６】
2010127783.tif 000053
【図７】
2010127783.tif 000054
【図８】
2010127783.tif 000055
【図１０】
2010127783.tif 000056
【図１１】
2010127783.tif 000057
【図１２】
2010127783.tif 000058
【図１４】
2010127783.tif 000059
【図９】
2010127783.tif 000060
【図１３】
2010127783.tif 000061
【手続補正書】
【提出日】平成22年3月25日(2010.3.25)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】００５１
【補正方法】変更
【補正の内容】
【００５１】
方程式６（数１０）の解は斉次解ｕ0ξと４階テンソル形式のＧｒｅｅｎ関数∧Ｇ（Ｇハット）を用いて次式（８）のように書き下せる。
【数１２】
2010127783.tif 000062