(57)【特許請求の範囲】
【請求項１】
表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析シミュレータであって、
前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、
この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析するとともに、前記固体相におけるエネルギー散逸には、表面側の原子配列に対し深さ方向に位置する複数層が無限に続くと仮定することによって、表面側から入力されたエネルギーの反射をなくした計算であって、エネルギー散逸を表現する積分表示Ｇｒｅｅｎ関数の漸近解を用いた計算を用いる固体挙動解析シミュレータ。
【請求項２】
請求項１に記載の固体挙動解析シミュレータであって、
さらに、積分表示のＧｒｅｅｎ関数についての数値解を計算する手段を有し、Ｇｒｅｅｎ関数で取り扱う力の履歴時間に応じて、履歴時間が短い間は漸近解を用い、履歴時間が長くなった場合に数値解を用いる固体挙動解析シミュレータ。
【請求項３】
コンピュータにより、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析プログラムであって、
コンピュータに、
前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、
この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析させるとともに、前記固体相におけるエネルギー散逸には、表面側の原子配列に対し深さ方向に位置する複数層が無限に続くと仮定することによって、表面側から入力されたエネルギーの反射をなくした計算であって、エネルギー散逸を表現する積分表示Ｇｒｅｅｎ関数の漸近解を用いた計算を用いる固体挙動解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析シミュレーションに関する。
【背景技術】
【０００２】
機械的な動力の伝達などの解析において、固体間摩擦の解析が必要となる。原子レベルから固体間摩擦を解析する理論手法に関しては、いくつかの簡易的なモデルを用いたものが提案されている。非特許文献１は原子スケールにおける固体間摩擦モデルとして、図１（ａ）に示すような剛体の壁に原子がバネで接続されたもの（Ｔｏｍｌｉｎｓｏｎモデル）を考案した。この一次元原子列をある周期ポテンシャルの上で直線運動させることで、摩擦の解析が行われた。非特許文献２，３では、結晶転移モデルなどでよく使われるＦｒｅｎｋｅｌ−Ｋｏｎｔｒｏｖａモデル（ＦＫモデル）を原子レベルの固体間摩擦解析に用いた。このモデルは図１（ｂ）のように原子同士がバネで結合し一次元鎖を形成したものである。また、非特許文献４では、より現実的な固体表面モデルとして、図１（ｃ）に示すように表面原子は互いに結合し、さらに剛体壁にも接続されているモデルを用いている。このモデルはＦｒｅｎｋｅｌ−Ｋｏｎｔｒｏｖａ−Ｔｏｍｌｉｎｓｏｎモデル（ＦＫＴモデル）、もしくはｇｅｎｅｒａｌｉｚｅｄ−ＦＫモデルなどと呼ばれる。Ｔｏｍｌｉｎｓｏｎモデル（図１（ａ））、ＦＫモデル（図１（ｂ））およびＦＫＴモデル（図１（ｃ））により固体間摩擦のメカニズムが論じられている。
【０００３】
また、動的分子シミュレーションによって、固体表面における動的物理現象を解析することについては、特許文献１〜３などに記載がある。
【０００４】
上述した非特許文献１〜４に示される既存モデルは全て表面原子のみが陽に扱われ、固体バルク原子の存在は全く考慮されていない。現実の固体表面は固体バルク原子が存在するため、既存モデルのように簡易化された表面では、原子スケールにおける摩擦機構を十分には捕らえられない。
【０００５】
そこで、本発明者らは、固体バルク内におけるエネルギー散逸について、本出願人は特開２００８−３０３１８６号において、表面側の原子配列に対し深さ方向に位置する複数層が無限に続くと仮定することによって、表面側から入力されたエネルギーの反射をなくした計算手法について提案した。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００２−１５０２０３号公報
【特許文献２】特開平８−３５８３１号公報
【特許文献３】特開平１０−２３２８６０号公報
【非特許文献】
【０００７】
【非特許文献１】Ｇ．Ａ．Ｔｏｍｌｉｎｓｏｎ，  Ｐｈｉｌｏｓ．  Ｍａｇ．  ７（１９２９）  ９０５
【非特許文献２】Ｙ．Ｉ．Ｆｒｅｎｋｅｌ  ａｎｄ  Ｔ．Ｋｏｎｔｏｒｏｖａ，  Ｚｈ．Ｅｋｓｐ．  Ｔｅｏｒ．  Ｆｉｚ．  ８，  １３４０，  １９３８．
【非特許文献３】Ｋ．Ｓｈｉｎｊｏ  ａｎｄ  Ｍ．Ｈｉｒａｎｏ，  Ｓｕｒｆ．  Ｓｃｉ．  ２８３，  ４７３，  １９９３．
【非特許文献４】Ｈ．Ｍａｔｓｕｋａｗａ  ａｎｄ  Ｈ．Ｆｕｋｕｙａｍａ，  Ｐｈｙｓ．  Ｒｅｖ．  Ｂ４９，  １７２８６，  １９９４．
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、先に提案の手法おいて、エネルギー散逸を記述するＧｒｅｅｎ関数は、上下面の格子間隔が同一（１対１コメンシュレート）の場合を前提としている。すなわち、上下面の格子間隔が異なる（非１対１コメンシュレート）場合には積分表示されたＧｒｅｅｎ関数しか得られていない。このため、非１対１コメンシュレート系を扱う場合、積分表示のＧｒｅｅｎ関数を数値的に求めることになり、計算コストが非常に大きく（１万原子の場合は１ヶ月以上）なる。
【課題を解決するための手段】
【０００９】
本発明は、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析シミュレータであって、前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析するとともに、前記固体相におけるエネルギー散逸には、表面側の原子配列に対し深さ方向に位置する複数層が無限に続くと仮定することによって、表面側から入力されたエネルギーの反射をなくした計算であって、エネルギー散逸を表現する積分表示Ｇｒｅｅｎ関数の漸近解を用いた計算を用いる。
【００１０】
また、積分表示のＧｒｅｅｎ関数についての数値解を計算する手段を有し、Ｇｒｅｅｎ関数で取り扱う力の履歴時間に応じて、履歴時間が短い間は漸近解を用い、履歴時間が長くなった場合に数値解を用いることが好適である。
【００１１】
本発明は、コンピュータにより、表面側から与えられたエネルギーに対する固体の挙動を解析する固体挙動解析プログラムであって、コンピュータに、前記固体の一部を、表面側に面状に配列された原子と、この表面側の原子配列に対し深さ方向に位置する複数層の原子配列と、から形成された３次元の固体相で構成し、この固体相の表面側から与えられるエネルギーによる固体相の挙動を、表面側からのエネルギーの深さ方向へのエネルギー散逸を伴う深さ方向に配列される原子配列の動的な物理現象としてその運動方程式を解いて解析させるとともに、前記固体相におけるエネルギー散逸には、表面側の原子配列に対し深さ方向に位置する複数層が無限に続くと仮定することによって、表面側から入力されたエネルギーの反射をなくした計算であって、エネルギー散逸を表現する積分表示Ｇｒｅｅｎ関数の漸近解を用いた計算を用いる。
【発明の効果】
【００１２】
本発明によれば、非１対１コメンシュレートの場合にも、比較的は短い時間で固体における深さ方向の伝達によるエネルギーの散逸を考慮して、表面側から与えられたエネルギーに対する固体の挙動の計算を完了することができる。
【図面の簡単な説明】
【００１３】
【図１】（ａ）Ｔｏｍｌｉｎｓｏｎモデル、（ｂ）ＦＫモデル、（ｃ）ＦＫＴモデルの模式図である。
【図２】（左）は、連成振動子モデルの模式図であり、表面垂直方向にＮ個、表面平行方向にはＮ’個の原子が連成し、表面平行方向に周期的境界条件が課されている。（右）は、表面平行方向および垂直方向に伸びた原子間の結合とそのバネ定数を示す。
【図３】本摺動計算系の模式図を示す図である。
【図４】ｓ関数について、数値計算による計算と、漸近解を使う計算による計算結果を示した図である。
【図５】ｓ関数について、数値計算による計算と、漸近解を使う計算による計算結果を示した図である。
【発明を実施するための最良の形態】
【００１４】
以下、本発明の実施形態について、図面に基づいて説明する。
【００１５】
材料開発において、成膜におけるスパッタリングや蒸着およびポリッシング、分析におけるＨｅ散乱、物性評価においてはナノインデンテーションや摩擦試験などが行われている。しかし、これらの実験のみでは起こっている現象を詳しく解析できない。このため、これら固体表面における動的物理現象に対して原子、分子スケールでのシミュレーションをする必要がある。
【００１６】
このシミュレーションは、汎用のコンピュータを用いて行われ、予め求められたモデルに対し、各種のデータを入力して動作を模擬することによって行われる。従って、モデルとして、物理現象を正しく表現したものが好ましい。
【００１７】
本実施形態は、このようなモデルを用いたシミュレーションを行うシミュレータおよびそのプログラムに関し、ＣＤ、ＤＶＤ、通信などによって、コンピュータの記憶部に前記モデルを利用したシミュレーションを行うプログラムをインストールし、プログラムを実行し、各種条件、データを入力して、摩擦現象や、摩擦係数などの物性値を求めるためのシミュレーションを行う。なお、シミュレーションプログラム自体は、すでに販売されているものを用い、そのシミュレーションプログラム内のモデルを本実施形態のプログラムによって設定することも好適である。
【００１８】
そして、シミュレーションは、現実の固体材料について行われ、材料開発などに用いられる。例えば、各種の原子、原子配列の固体についての、摩擦係数を求めることで、所望の摩擦係数を得るために、どのような原子、原子配列の材料がよいかを求めることができ、現実の材料を特定することができる。例えば、自動車の動力伝達機構であるクラッチや、ギアなどにおける材料にどのような原子構造のものを採用するのがよいかなどを知ることができる。
【００１９】
「Ｃ  無限連成振動子方程式（Ｎ→∞モデル）」
固体表面を図２のように結晶バルク方向に原子をＮ個連ねた連成振動子モデルで表し、無限系（Ｎ→∞）の運動方程式を解く。このＮ→∞により表面から与えられた摺動エネルギーは無限の彼方まで伝わっていき、再び表面に戻ってこないので現実的なエネルギー散逸を再現できる。これは無限に大きなバルクの影響が熱浴として振る舞うと考えても良い。さらに、現実の固体においてＮは１０１０オーダーと非常に大きいため、このＮ→∞操作は実際の摩擦現象をよく再現すると考えられる。
【００２０】
以下に、無限連成振動子モデルの運動方程式の構築およびＧｒｅｅｎ関数の漸近解を用いた解法について説明する。
【００２１】
図２（左）に連成振動子モデルの模式図を示す。ユニットセルは表面垂直方向にＮ個、表面平行方向にはＮ’個の原子（合計Ｎ×Ｎ’個）が連成した正方格子である。そのユニットセルに対し表面平行方向に周期的境界条件を課すことで、表面を形成する。整数η，ξは原子指標で、それぞれη＝１，２，・・・Ｎ’，ξ＝１，２，・・・Ｎの範囲をとる。ξ＝１は表面原子を表わし、結晶バルク側の原子に進む方向に沿ってξを番号付けしている。ξ＝Ｎの原子は剛体壁に接続されている。原子間の結合は線形近似で表す。表面平行方向および垂直方向に伸びたバネ定数をそれぞれ（Ｋ’ｘ，Ｋ’ｚ），（Ｋｘ，Ｋｚ）とする（図２（右）参照）。ここでＫｘ，Ｋ’ｘ、Ｋｚ，Ｋ’ｚは、表面平行方向および垂直方向の平衡原子間距離からの変位に対するバネ定数である。ｘ，ｚを表面平行方向および垂直方向に対する原子の平衡位置からの変位とする。すなわち、変位ｕ＝ｘ  ｏｒ  ｚであり、ξ番目の層の原子変位ｕξ＝（ｕ１，ξ，・・・ｕ１，ξ，・・・ｕ１，ξ）と表される。また、外部から作用する力は表面原子にのみ作用すると仮定する。その時、原子変位ｕη，ξ＝ｘη，ξ，ｚη，ξの運動方程式は次のようになる。
【００２２】
【数１】
0005428694.tif 000002
【００２３】
ここで、Ｆｘ，Ｆｚはそれぞれ表面原子ξ＝１にのみ作用する力の表面平行方向および表面垂直方向成分である。本モデルを下面と上面に配置すれば、Ｆｕは二固体間の相互作用力を表す。ｍは原子の質量を表す。表面平行方向に対する周期的境界条件は、
【数２】
0005428694.tif 000003
である。
【００２４】
式４（数８）をベクトル表示すると、次のようになる。
【数３】
0005428694.tif 000004
0005428694.tif 000005
【００２５】
ここで∧Ｂ（Ｂハット），∧Ｃ（Ｃハット）は、Ｎ×Ｎ×Ｎ’×Ｎ’の４階テンソル、∧ＣはＮ’×Ｎ’の２階テンソル（行列）である。∧Ｂ，∧Ｃは表面垂直方向および面内の質点間のバネの相互作用力項を表している。∧１はＮ’×Ｎ’の２階テンソルの単位元である。面間の相互作用力Ｆｕは表面原子ξ＝１に働き、
【数４】
0005428694.tif 000006
である。
【００２６】
方程式６（数３）の解は斉次解ｕ０ξと４階テンソル形式のＧｒｅｅｎ関数∧Ｇ（Ｇハット）を用いて次式８のように書き下せる。
【数５】
0005428694.tif 000007
【００２７】
後に論ずる式展開と数値解の容易さのため、上式を時間で微分した速度に対する解が有用となる（式９）。
【数６】
0005428694.tif 000008
【００２８】
ここで∧ｇ＝ｄ／ｄｔ∧Ｇである。Ｇｒｅｅｎ関数∧Ｇと斉次解ｕ０ξはそれぞれ以下の式を満たす。
【数７】
0005428694.tif 000009
【００２９】
つまり、運動方程式６を解くには、式１０および式１１の２つの方程式により∧ｇとｄ／ｄｔ  ｕ０を求め、方程式９に代入して解けば良い。以下∧ｇおよびｄ／ｄｔｕ０を求める。
【００３０】
Ｃ−１∧ｇの導出
方程式１０にｅｘｐ（−ｓｔ）をかけてラプラス変換を施すと、式１２となる。
【数８】
0005428694.tif 000010
【００３１】
ここで、∧Ｂ，∧Ｃを対角化する行列を∧ＵＢおよび∧Ｕｃとする。すなわち、次式となる。
【数９】
0005428694.tif 000011
【００３２】
また、∧ＵＣを次のように定義する。
【数１０】
0005428694.tif 000012
【００３３】
方程式１２（数８）の両辺に∧ＵＣ∧ＵＢをかける。∧ＵＢと∧ＵＣは可換であることを使えば（∧ＵＢのξ，ｊ成分はＵξ，ｊ∧１であるから）、∧ＵＣ∧ＵＢ∧Ｇ（ｓ）≡〜Ｇ（ｓ）のη，ｌ，ξ，ｊ成分は以下の式１３を満たす。
【００３４】
【数１１】
0005428694.tif 000013
ここで（ωＢｊ）２＝（Ｋｕ／ｍ）λＢｊおよび（ωＣｌ）２＝（Ｋ’ｕ／ｍ）λＣｌである。∧Ｂ，∧ＣはエルミートであるからＵＢ，ＵＣはユニタリ行列である。よって両辺に（ＵＢ）−１（ＵＣ）−１＝（ＵＢ）ｔ（ＵＣ）ｔを作用させ、さらにラプラス逆変換を施せば∧Ｇ（ｔ）のη，ｌ，ξ，ｊ成分は、次式となる。
【００３５】
【数１２】
0005428694.tif 000014
【００３６】
よって∧ｇ（ｔ）は、式１４となる。
【数１３】
0005428694.tif 000015
【００３７】
相互作用力Ｆｕは、表面原子層ξ＝１にのみ作用すると仮定したため、∧ｇの行列要素として必要なのはξ’＝１成分のみである。そのためｇη，η’，ξ，ξ’＝１≡ｇη，η’，ξと表記し直す。
【００３８】
ωＢｊ，ＵＢξ，ｊは次のように表される。
【数１４】
0005428694.tif 000016
【００３９】
ここで、θＢ＝π（ｊ−１／２）／（２Ｎ＋１），α＝２√（Ｋｕ／ｍ）である。ｊの取りうる範囲はｊ＝１，２，・・・，Ｎである。式１４に式１５，１６を代入すれば、次式が得られる。
【数１５】
0005428694.tif 000017
【００４０】
ここで、バルク方向の原子数Ｎの無限大極限をとる。Ｎが十分大きければ、次式となる。
【数１６】
0005428694.tif 000018
【００４１】
Ｎ→∞によりｊに関する和を積分に変換し、θＢ→θ、ωＢｊ→ω（θ）＝αｓｉｎθと書き換えると、式１７が得られる。
【数１７】
0005428694.tif 000019
【００４２】
ωＣｌ，ＵＣη，ｌは次のように表される。
【数１８】
0005428694.tif 000020
【００４３】
ここで、α’＝２√（Ｋ’ｕ／ｍ）である。上式を式１７に代入し、ｓξ，ｌ（ｔ）＝ｓξ，−ｌ（ｔ）であることを使えば、式１９（数１９）が得られる。
【数１９】
0005428694.tif 000021
【００４４】
特に、ｌ＝０の時、ｓξ，ｌ（ｔ）の式１８の積分はｎ次の第一種Ｂｅｓｓｅｌ関数Ｊｎの積分表示に帰結し、式２０で表現できる。
【数２０】
0005428694.tif 000022
【００４５】
ここで、第一種Ｂｅｓｓｅｌ関数の積分表示πＪ２ξ（ｚ）＝∫ｄθｃｏｓ（ｚｓｉｎθ）ｃｏｓ２ξθ［θ＝０：π］＝２（−１）ξ∫ｄθｃｏｓ（ｚｃｏｓθ）ｃｏｓ２ξθ［θ＝０：π／２］を用いた。
【００４６】
ここで、ｇη，η’，ξの物理的意味を議論しておく。原子（η’，１）に作用する力Ｆ（ｔ）が原子（η，ξ）の速度に与える影響は、次式で表される。
【数２１】
0005428694.tif 000023
【００４７】
力Ｆは、関数ｓとの畳み込み積分により、バルク方向に伝わる波と面内方向に伝わる波に分解される。そして、ｅｘｐ（２πｉｌ（η−η’）／Ｎ’）は、面内方向へ伝搬する波の位相を表現している。その伝搬していく様々な周波数の波を足し合わせたものが∫ｇη，η’，ξ（ｔ−ｔ’）Ｆ（ｔ’）ｄｔ’［ｔ’＝０：ｔ］となる。つまり、周囲の原子への力の伝搬を表すのがｇη，η’，ξである。
【００４８】
Ｃ−２（ｄ／ｄｔ）ｕ０の導出
有限Ｎ系であれば摺動計算において時間発展させれば、ある一定時間後には初期条件によらず定常状態を得られ、その時の動摩擦係数を求めることができる。よって、最も簡単な初期条件ｕ０ξ（ｔ＝０）＝０，（ｄ／ｄｔ）ｕ０ξ（ｔ＝０）＝０における斉次解、つまりｕ０ξ（ｔ）＝０を用いて差し支えない。しかしながら、この初期条件では無限Ｎ系で定常状態を得ることができない。このことを簡単な例を用いて示す。
【００４９】
初期条件ｕ０ξ（ｔ＝０）＝０，（ｄ／ｄｔ）ｕ０ξ（ｔ＝０）＝０にある表面原子ξ＝１に一定の力Ｆが作用したとする。十分長い時間ｔ＊経過すれば、力の釣り合いからｕξ＝１（ｔ＊）＝ＸｕＮ，Ｘｕ＝Ｆ／Ｋｕの定常状態になる。有限Ｎの系ではｕξ＝１（ｔ＊）は有限の値であり、時間発展によって定常状態に到達できる。しかしながらＮ→∞の系ではｕξ＝１（ｔ＊）は発散してしまう。これは一つの原子間の結合のひずみが小さくとも、Ｎが無限であるので発散するためであり、無限系特有の症状といえる。つまり無限系では、初期条件ｕ０ξ（ｔ＝０）＝０，（ｄ／ｄｔ）ｕ０ξ（ｔ＝０）＝０から時間発展させもｕξ＝１は定常状態に到達することができない。
【００５０】
無限系で定常状態を達成するには、あらかじめ初期条件として原子間に一定のひずみＸｕを与えておき、それが摺動による相互作用力の時間平均と釣り合うように決定する手法を用いる。その手法を以下に詳しく述べる。
【００５１】
まず式１１の初期条件を設定する。初期時間ｔ＝０における初期条件を、表面垂直方向に配列した原子間はある一定ひずみＸｕ＝ｕ０η，ξ−ｕ０η，ξ＋１を持ち、面内原子間のひずみは起きずｕ０η，ξ＝ｕ０η’，ξとなって静止している、と設定する。これは結晶格子が表面原子にある一定かつ一様な力を受けて平衡状態になった配置で、ひずみを持ったまま静止している状態である。このとき面内の原子間ひずみは起きないので∧Ｃは消え、原子の運動にη依存性はなくなる。よってｕ０ξ＝ｕ０ξ（左側の項はベクトル（太字）ｕ）と表示し直すと、式１１は以下のようになる。
【００５２】
【数２２】
0005428694.tif 000024
【００５３】
以下、式２１を解いていく。式２１の両辺に∧ＵＢをかけ、〜ｕ０ｊ＝ΣξＵＢξ，ｊｕ０ξとおくと、次式となる。
【数２３】
0005428694.tif 000025
【００５４】
上式は簡単に解くことができ〜ｕ０ｊ（ｔ）＝ｃｊｃｏｓωＢｊｔを得る。ただし、（ｄ／ｄｔ）〜ｕ０ξ（ｔ＝０）＝０を用いた。
【００５５】
未定係数ｃｊは初期変位に関する条件ｃｊ＝ΣＵＢξ，ｊＸｕ（Ｎ−ξ＋１）［ξ＝１：Ｎ］から決定される。
【００５６】
式１６を用いれば、次式が得られる。
【数２４】
0005428694.tif 000026
【００５７】
公式Σｓｉｎ２ξθＢ［ξ＝１：Ｎ］＝ｓｉｎ（Ｎ＋１）θＢｓｉｎＮθＢ／ｓｉｎθＢにより、次式を得る。
【数２５】
0005428694.tif 000027
【００５８】
よってｕ０ξ（ｔ）は、次式となる。
【数２６】
0005428694.tif 000028
【００５９】
ここで、（ｄ／ｄｔ）ｕ０ξ（ｔ）は、次のように変形できる。
【数２７】
0005428694.tif 000029
【００６０】
最後は、Ｎが非常に大きいことを使った。Ｎ→∞を行いｊに関する和を積分に変換し、ωＢｊ→ω（θ）＝αｓｉｎθと書き換えると、次式を得ることができる。
【数２８】
0005428694.tif 000030
【００６１】
ここで、最後の等号は、第一種Ｂｅｓｓｅｌ関数の積分表示を用いた。この（ｄ／ｄｔ）ｕ０ξ（ｔ）を用いて定常状態を得る方法に関しては次の項で述べる。
【００６２】
Ｃ−３半無限連成振動子方程式と自己無撞着計算
上述の結果により、（η，ξ）原子の運動方程式の形式解は、式２３のように記述することができた。
【数２９】
0005428694.tif 000031
【００６３】
ここで、α＝２√（Ｋｕ／ｍ）、α’＝２√（Ｋ’ｕ／ｍ）、ω（θ）＝αｓｉｎθ、ωＣｌ＝α’ｓｉｎπｌ／Ｎ’である。この新しく得られた微積分方程式を数値的に解けば、半無限連成振動子系の原子の運動を求めることができる。また今回は２次元に配列した結晶格子で定式化を行ったが、同様の式展開を行えば３次元に配列した結晶格子においても方程式を得ることができる。
【００６４】
ところで、動摩擦係数μ、単位面積あたりの動摩擦力Ｐ／／、垂直圧力Ｐ⊥は、次のように求めることができる。
【数３０】
0005428694.tif 000032
【００６５】
ここで＜＞ｔは定常状態での時間平均操作、Ｓは表面原子η＝１，２，・・・Ｎ’が構成する面積を表す。
【００６６】
外部との相互作用力は表面原子層ξ＝１のみに作用するため、摩擦解析に必要な物理量を得るには、表面原子ξ＝１の運動が分かれば良い。もしも、運動方程式６自体により表面原子ξ＝１の運動を解こうとすれば、Ｎ→∞漸近解が扱えないのは当然であるが、バルク原子ξ≠１も同時に解かねばならない。一方、式２３のように形式的に運動方程式を解いたものを用いれば、ξ≠１を解かずに原子ξ＝１の運動を得ることができる。よって解くべきξ＝１の微積分方程式は、次のようになる。
【００６７】
【数３１】
0005428694.tif 000033
【００６８】
ここで、表面原子を表すξ＝１を省略し、（ｄ／ｄｔ）ｕ０ξ（ｔ）を用いて書き直すと、次のようになる。
【数３２】
0005428694.tif 000034
【００６９】
ここで、上述と同様に、ｕは表面原子のｘもしくはｚ方向の変位、ηは表面平行方向の原子ラベル、Ｎ’はユニツトセル内の表面平行方向に並んでいる原子数、Ｆは表面原子に加わる力、ｍは原子質量を表す。
【００７０】
また、Ａ＝２π／Ｎ’，ω(θ)=αｓｉｎθ、ωｌ＝α’ｓｉｎπｌ，α＝２√｛Ｋｕ／ｍ｝、α’＝２√｛Ｋ’ｕ／ｍ｝，Ｋｕ，Ｋ’ｕバネ定数である。（ｄ／ｄｔ）ｕ０は初期変位に起因する斉次解である。
【００７１】
式３２を式３１に代入すると、次のようになる。
【数３３】
0005428694.tif 000035
【００７２】
従って、関数ｓｌが得られれば無限連成振動方程式が計算できる。
【００７３】
ここで、図３に本摺動計算系の模式図を示す。下面は表面垂直方向に無限に原子が配列し、表面平行方向にも原子が周期的に配列した２次元正方格子である。原子間の結合は線形近似で表し、表面垂直方向に伸びたバネ定数を（Ｋｘ，Ｋｚ）とする。ここでＫｘ，Ｋｚは表面平行方向および垂直方向の平衡原子間距離からの変位に対するバネ定数である。上面を構成する原子は伸縮しない（無限に固い）結合で結ばれて配列している。
【００７４】
上面と下面は、同じ格子定数を持つ物質（１対１コメンシュレート系）と仮定すれば、表面平行方向の周期性により面内の原子は全て同じ運動をするため、面内の振動は起きず表面平行方向に伸びた結合の寄与は無い。この場合には、上述の無限連成振動子方程式を簡略した形式に書き換えて計算が可能となる。
【００７５】
一方、このような仮定をしない場合には、上述のＧｒｅｅｎ関数の計算が必要であり、すなわち、式３３のように積分表示で示される、ｓｌを計算する必要がある。ここで、この積分を数値的に実行すると、結局は有限Ｎの計算に帰結してしまう。なぜなら、固体原子層数Ｎ→∞の操作により、Σ→∫ｄθへと変換を行った帰結が式３３であるため、数値的にこの積分を差分化すると積分分割数と等しい原子層数Ｎのｓ関数を生成することと等価になるからである。実際に、数値的にｓｌを計算した結果を図４に示す。このように、ｓｌはｔが大きくなるにつれ振動しながら減衰していく。ｓｌは過去に受けた力が原子速度へ及ぼす影響であり、これは現在（ｔ＝０）から過去に向かうにつれ力の影響が減衰していくことを示している。すなわち、グリーン関数が過去に受けた力の影響をメモリーするメモリー項を含んでおり、遠い過去の影響ほど小さくなっていることを示している。
【００７６】
そして、図４に示すとおり、ｓｌは、ｔ＝５０付近で再び大きな振幅を持ってしまう。これは、ある時刻で以前受けた力の影響が有限原子層の端で振動が跳ね返って再帰したことを示している。
【００７７】
この問題を解決する最も一般的な方法は、摺動時間Ｔの間にエネルギーの反射が起きないほどθ積分の分割数を増やすことである。しかし、定常状態を得るまでのＴは非常に長く、反射が起きない分割数（およそ１万分割以上）にしてしまうとｓｌ作成に膨大な計算時間がかかってしまう。さらに、結局は有限Ｎの計算になってしまうため、Ｎ→∞の理論の枠組みとして不完全である。
【００７８】
そこで、積分を数値的ではなく、解析的に実行する必要がある。式３３の積分は解析的に実行できなかったが、ｔ→大の領域で積分が実行でき漸近解が得られることを以下に示す。式３３の変数を書き換えると、次のようになる。
【数３４】
0005428694.tif 000036
【００７９】
ここで、β２＝（Ｋ′ｕ／Ｋｕ）ｓｉｎ２（πｌ／Ｎ′）である。式３５をＴａｙｌｏｒ展開しθ積分をまず実行する。
【数３５】
0005428694.tif 000037
【００８０】
ここで、ｎＣｍはｎ個からｍ個を選ぶ組み合わせ総数（注：ｎＣｍ＝ｎ！／ｍ！（ｎ−ｍ）！）を表す。積分公式∫ｃｏｓ２θｓｉｎ２ｍθｄθ［０〜π／２］＝（π／２）（２ｍ−１）！！／（２ｍ＋２）！！を上式に用いて積分を実行すると（注：ｎ奇数：ｎ！！＝ｎ（ｎ−２）・・・３・１、ｎ偶数：ｎ！！＝ｎ（ｎ−２）・・・４・２）、次のようになる。
【００８１】
【数３６】
0005428694.tif 000038
【００８２】
（２ｍ−１）！！＝２ｍ！／２ｍｍ！と、２ｍ！！＝２ｍｍ！を用いた。
【００８３】
ｎ−ｍ＝ｋとして和の取り方を変更すれば、次のようになる。
【数３７】
0005428694.tif 000039
【００８４】
ｎ′＝ｎ−ｋとおけばｓｋｌは、次のように表される。
【数３８】
0005428694.tif 000040
【００８５】
ここでｚ＝αｔ／２とおいた。一方、ｓｋ−１ｌは、次の通りである。
【数３９】
0005428694.tif 000041
【００８６】
式３８を式３７に代入すればｓｋｌは次の漸化式を持つことが分かる。
【数４０】
0005428694.tif 000042
【００８７】
Ｂｅｓｓｅｌ関数の公式Ｊ１（２ｚ）／ｚ＝Σ｛（−１）ｎ’ｚ２ｎ／ｎ！（ｎ＋１）！［ｎ’＝０〜∞］を用いれば初項ｓ０ｌは、次の通りになる。
【数４１】
0005428694.tif 000043
【００８８】
漸化式３９に、式４０の初項を代入していき、ｋについて無限大の和を取ればｓｌが求まるが、解析的に和を実行することができない。そこで、ｔが大きい領域のみを考える。ｓ１ｌは、次のように表される。
【数４２】
0005428694.tif 000044
【００８９】
ｚが十分に大きければ∫ｄｚＪ１（ｚ）／ｚ→１［０〜ｚ］となるため、ｔが大きい時の漸近解（以降、漸近解を表す際には上付きの〜を付ける）は、下記のようになる。
【数４３】
0005428694.tif 000045
【００９０】
よって〜ｓｋｌは、次のようになる。
【数４４】
0005428694.tif 000046
【００９１】
よって、ｓｌの漸近解は、Ｂｅｓｓｅｌ関数公式を用いると、次のように表される。
【数４５】
0005428694.tif 000047
【００９２】
このようにして、ｓｌについての漸近解〜ｓｌが得られ式３２のグリーン関数ｇη，η’（ｔ）を求めることができる。そこで、積分表示のＧｒｅｅｎ関数を時間の大きい領域で解析的に積分を実行し、得られた漸近解により無限個の固体原子層をＧｒｅｅｎ関数に組み込むことができる。特に、漸近解を用いる手法は、β＝０としておらず、非１対１コメンシュレート系に適用できる。従って、非１対１コメンシュレート系における計算の高速化が実現することができる。
【００９３】
「漸近解についての検証」
次に、この数学手続きで得られた漸近解がどの程度のｔで一致するのかを検証する。ここで、ｔは、Ｇｒｅｅｎ関数（またはｓ関数）で取り扱う力の履歴時間である。
【００９４】
図４と５に〜ｓｌと数値積分によって求めたｓｌを示す。漸近解は、ｔが小さいところでは数値積分のものとズレがあるが、ｔが大きくなるにつれ差は小さくなりほぼ完全に一致している。また、ｔの増加によって〜ｓｌは速やかにｓｌに収束していることが分かり、この漸近解の性質の良さが確かめられた。
【００９５】
これは漸近解の条件∫ｄｚＪ１（ｚ）／ｚ→１［０〜ｚ］の被積分関数がｔ−３／２という時間に対して速やかに０に近づくため、ｔが比較的小さくとも漸近解の条件がよく満たされるからである。
【００９６】
最後に、ｓｌの作成方法を述べる。図４，５は、ｓｌ（ｔ）についての式３３について、数値積分（白丸）と漸近解の式３４（実線）を示している。なお、図４，５において数値積分の分割数は２０、α＝２を用い、図４においてβ＝ｓｉｎ（π／１０）、図５においてβ＝ｓｉｎ（π／２）を用いている。図４に示すようにβが小さい場合、漸近解はｔが小さい領域でも良い一致を示す。しかし、図５のようにβが大きい場合、漸近解はｔが小さい領域で大きくずれてしまう。
【００９７】
そこで、まず漸近解が数値積分による解とほぼ等しくなる時間をｔｃとおく（図４参照）。そして、ｔ≦ｔｃでは式３５の積分を数値的に行い、ｔ＞ｔｃでは式４４による漸近解を用いてｓｌを作成すればよい。ｔが小さければ、Ｎがそれほど大きくなくても、反射による影響がなく、比較的正確な計算をそれほど時間をかけずに行うことができる。すなわち、通常の計算において、数値積分において、反射の影響が出始める前に、数値積分と漸近解が一致する時点（ｔｃ）が得られる。そこで、そのｔｃまでを数値積分として、その後は漸近解を使用することで、反射の影響を排除して比較的正確な演算を行うことができる。例えば、数値積分を用いる計算と、漸近解を用いる計算を並列して行い、両者が一致した時点で、数値積分の計算を中止して、その後の履歴時間については漸近解を使用する演算のみを行えば良い。また、βの値に応じて、両者が一致する時間が異なるため、βの値に応じてｔｃを変更して、不要な数値積分の計算をなるべく少なくすることが好適である。
【００９８】
従って、このように２種類の演算を組み合わせることで、全体として正確な演算を比較的短時間で行うことができる。
【００９９】
ここで、定常状態に至るまでの時間とＴとした場合、数値積分でｓ１を計算しようとすると、Ｎ’×Ｔ２の計算時間が掛かる（Ｎ’は原子数）。一方、漸近解を使用すれば、Ｎ’×Ｔの時間で計算が可能である。例えば、Ｎ’＝１００００程度の場合、現状のコンピュータによる計算で、数値計算では１ヶ月程度の時間が掛かるのに対し、漸近解を用いることにより、１０分以下で計算できることがわかっている。
【０１００】
このようにして、本実施形態によれば、非１対１コメンシュレート系においても、Ｇｒｅｅｎ関数についての計算が可能となる。従って、上面側の固体の格子定数と、上面側の固体の格子定数が異なる系において、上面側と下面側の摩擦を両面の原子同士の衝突として、その運動の計算が可能となり、その際のエネルギーの散逸に応じた摩擦係数の計算が可能となる。特に、一方の面（上述の場合は上面）についてはその下方に無限に原子が連なると仮定して、エネルギーが反射してくることを防止した原子の運動方程式を解く計算が可能となる。
【図１】
0005428694.tif 000048
【図２】
0005428694.tif 000049
【図４】
0005428694.tif 000050
【図５】
0005428694.tif 000051
【図３】
0005428694.tif 000052