(57)【要約】【課題】各々指定された歪度及び尖度を有する粗面同士の接触状況を算出することができる。【解決手段】分布算出部３２が、二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する。接触算出部３４が、前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する。【選択図】図７
2018152033.tif 000002
【特許請求の範囲】
【請求項１】
二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置であって、
前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部と、
前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部と、
前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部と、
を含む演算装置。
【請求項２】
前記分布算出部は、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布を算出する請求項１に記載の演算装置。
【請求項３】
前記分布算出部は、前記突起頂点分布としての、基準面からの突起頂点の高さが各高さであり、かつ、突起頂点の曲率が各曲率である確率を規定する確率密度関数を表すGreenwood simplified ellipticモデルを、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布に変換することにより、前記突起頂点分布を算出する請求項２に記載の演算装置。
【請求項４】
前記分布算出部は、前記突起頂点分布の変換の際に、前記粗さ分布がJohnson分布によるときの前記突起頂点分布の変換に応じて、突起頂点の各高さに対する突起頂点の曲率を補正するように、前記突起頂点分布を算出する請求項３記載の演算装置。
【請求項５】
前記接触算出部は、前記突起頂点分布について予め計算された積分積算を表す演算子、及び二つの前記突起頂点分布の組み合わせについて予め計算された積分積算を表す演算子を用いて、前記真実接触面積又は垂直抗力を算出する請求項１〜請求項４の何れか１項に記載の演算装置。
【請求項６】
二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置を制御するためのプログラムであって、
コンピュータを、前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部、
前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部、及び
前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部、
として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、演算装置、及びプログラムに係り、特に、面同士の接触状況を算出するための演算装置、及びプログラムに関する。
【背景技術】
【０００２】
従来技術では、面同士の接触について、Greenwood and Williamson(1966)のモデル（ＧＷモデル）を用いて接触面積および垂直抗力を算出している（特許文献１参照）。特許文献１には、運動摩擦を最小限に抑える歪度値および尖度値を決定するための接触モデリングの確率密度関数を用いた、表面粗さ分布の設計方法が記載されている。このような従来技術の手法は、簡易計算としては有用な手法である。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特表２０００−５００６０４号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかし、特許文献１に記載されているような従来技術には、以下の問題点により、その計算精度と適用対象の範囲が不十分といえる。
【０００５】
一つ目の問題は、ＧＷモデルでは、図１１に示すように、粗面と平滑面との接触を仮定しているため、粗面同士の接触を直接取り扱うことができないことが挙げられる。両粗面の突起頂点分布が確率論や統計学で用いられる正規分布であると仮定した場合は、粗面同士の接触であっても合成粗面と平滑面の接触とみなすことができるが、突起頂点分布が非正規分布の場合は従来手法を用いることができない。よって、従来技術では粗面及び平滑面を対象とした限定された範囲の計算しか行えない。
【０００６】
二つ目の問題は、ＧＷモデルでは突起頂点分布を正規分布としているため、従来技術では突起頂点分布をPearson分布に置き換えて計算していることである。従来の解析（Nayak(1971)）から明らかなように、粗さ分布が正規分布であっても突起頂点分布は一般には正規分布とはならない。したがって、粗さ分布が正規分布の突起頂点分布を直接Pearson分布で変換とするのは誤りである。また、粗面の持つ歪度、及び尖度は粗さ分布について与えられるものであり、突起頂点分布を指定された歪度、及び尖度になるように変換するのは不適切である。
【０００７】
本発明は、上記問題点を解決するために成されたものであり、各々指定された歪度及び尖度を有する粗面同士の接触状況を算出することができる演算装置、方法及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するために、第１の発明に係る演算装置は、二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置であって、前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部と、前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部と、前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部と、を含んで構成されている。
【０００９】
また、第１の発明に係る演算装置において、前記分布算出部は、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布を算出するようにしてもよい。
【００１０】
また、第１の発明に係る演算装置において、前記分布算出部は、前記突起頂点分布としての、基準面からの突起頂点の高さが各高さであり、かつ、突起頂点の曲率が各曲率である確率を規定する確率密度関数を表すGreenwood simplified ellipticモデルを、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布に変換することにより、前記突起頂点分布を算出するようにしてもよい。
【００１１】
また、第１の発明に係る演算装置において、前記分布算出部は、前記突起頂点分布の変換の際に、前記粗さ分布がJohnson分布によるときの前記突起頂点分布の変換に応じて、突起頂点の各高さに対する突起頂点の曲率を補正するように、前記突起頂点分布を算出するようにしてもよい。
【００１２】
また、第１の発明に係る演算装置において、前記接触算出部は、前記突起頂点分布について予め計算された積分積算を表す演算子、及び二つの前記突起頂点分布の組み合わせについて予め計算された積分積算を表す演算子を用いて、前記真実接触面積又は垂直抗力を算出するようにしてもよい。
【００１３】
第２の発明に係るプログラムは、二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置を制御するためのプログラムであって、コンピュータを、前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部、前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部、及び前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部、として機能させるためのプログラムである。
【発明の効果】
【００１４】
本発明の演算装置、及びプログラムによれば、二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出し、前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出することにより、各々指定された歪度及び尖度を有する粗面同士の接触状況を算出することができる。
【図面の簡単な説明】
【００１５】
【図１】粗面と粗面との接触を仮定した本実施の形態の一例を示す図である。
【図２】歪度、及び尖度の組合せと、組み合わせに対応する分布の種類の一例を示す図である。
【図３】歪度、及び尖度による分布形状の一例を示す図である。
【図４】粗さ分布と突起頂点分布の一例を示す図である。
【図５】数値的に生成した粗面における平均突起形状の一例を示す図である。
【図６】変換による突起形状の変化の一例を示す図である。
【図７】本発明の実施の形態に係る演算装置の構成を示すブロック図である。
【図８】本発明の実施の形態に係る演算装置における演算処理ルーチンを示すフローチャートである。
【図９】実験における膜圧比に対する真実接触面積の変化を示す図である。
【図１０】実験における膜圧比に対する真実接触面積の変化を示す図である。
【図１１】粗面と平滑面との接触を仮定した従来技術の一例を示す図である。
【発明を実施するための形態】
【００１６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１７】
＜本発明の実施の形態に係る概要＞
【００１８】
まず、本発明の実施の形態における概要を説明する。
【００１９】
上記課題の一つ目の問題に対し、本発明の実施の形態では、粗さ突起接触モデルであるGreenwood and Trippモデル(ＧＴモデル)の考え方を用いることによって、図１に示すように、非正規分布を持つ粗面同士の接触を取り扱うことができる。これらの粗面は互いに歪度、及び尖度が異なっていても計算可能である。
【００２０】
上記課題の二つ目の問題に対し、本発明の実施の形態では、粗さ分布が指定された歪度、及び尖度に変換されるときに、歪度、及び尖度に従属する形で突起頂点分布を変換する。変換前の突起頂点分布はGreenwood(2006)のsimplified ellipticモデル(ＧＳＥモデル)で与えられる（非特許文献１参照）。この計算手順については、以下の原理の説明で詳細を示す。
【００２１】
［非特許文献１］：Greenwood, J. A. (2006). A simplified elliptic model of rough surface contact. Wear, 261(2), 191-200.
【００２２】
ＧＴモデルをベースに、以下の改良を施すことで精度を大幅に向上させる。ＧＴモデルでは突起頂点分布を正規分布とし、突起曲率半径を一定と仮定しており、このままでは精度が不十分なため、より実粗面に近い突起頂点分布および突起曲率半径分布を用いているＧＳＥモデルを用いてＧＴモデルの式を改良する。
【００２３】
また、本発明の実施の形態では、取り扱いやすさと精度を考慮し、非正規分布としてPearson分布に代わりJohnson分布を用いる。上記の改良では粗さ分布が正規分布である粗面しか扱えないため、歪度、尖度に特徴を持つ粗面について取り扱うことができない。そこで、粗さ分布にJohnson分布を用いて任意の歪度、尖度を持つ粗面を取り扱えるようにした。粗さ分布にJohnson分布を用いた場合の突起頂点分布の変換について厳密に計算を行い、曲率半径の変化についても補正を与えることによって精度を維持する。
【００２４】
また、既存の突起接触モデルではＨｅｒｔｚの接触理論に基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積、及び垂直抗力を計算している。このとき、突起を一定の曲率を持つ半球であると仮定しているが、実際の突起は先端から離れるにしたがって曲率が小さくなる。そこで、本発明の実施の形態では突起の曲率変化を考慮して、元の接触モデルにおける曲率に一定倍率の補正を与える。
【００２５】
以上の改良によって、検討可能な粗面の条件が大幅に拡張され、さらに直接数値計算と比較しても非常に良好な真実接触面積の予測が可能となる。
【００２６】
＜本発明の実施の形態に係る原理＞
【００２７】
次に、本発明の実施の形態における原理を説明する。
【００２８】
本発明の演算装置に入力するパラメータは、ηi：突起数密度、σki：断面曲率又はβi：平均突起曲率半径、σi：二乗平均平方根粗さ、Ｓski：歪度（スキューネス）、Ｓkui：尖度（クルトシス）、Ｅi：ヤング率（縦弾性係数）、ｖi：ポアソン比、Ａn：見かけの接触面積である。添字のｉは粗面１、及び粗面２それぞれの持つ物理量であることを示す。以下に、これらのパラメータについて概要を述べる。
【００２９】
以下の原理の説明では、等方粗さを持つ粗面を仮定する。以下（１）式の粗さ断面曲線のパワースペクトル密度Ｐ（k）（k：波数）のモーメントを用いて、以下（２）式で二乗平均平方根粗さσ、断面曲率σｋ、突起数密度ηは算出される。
【００３０】
【数１】
2018152033.tif 000003
・・・（１）
2018152033.tif 000004
・・・（２）
【００３１】
上記（１）式、（２）式のパラメータをパワースペクトル密度から算出する以外に、これらのパラメータに実測した値を直接用いても良い。本発明の実施の形態で援用するＧＳＥモデルの場合、実測で得られる平均突起曲率半径βと断面曲率σｋとの間には以下（３）式の関係がある。
【００３２】
【数２】
2018152033.tif 000005
・・・（３）
【００３３】
３次元の離散点で表現される粗面において、ある点の基準面からの高さがｚである確率が確率密度関数（以下、粗さ分布関数）φ(z)で規定されるとき、歪度Ｓsk、及び尖度Ｓkuは以下（４）式、及び（５）式となる。
【００３４】
【数３】
2018152033.tif 000006
・・・（４）
2018152033.tif 000007
・・・（５）
【００３５】
ここで、ｍは粗さ基準面高さであり、一般にデータ処理前に取り除かれるため、以降はｍ＝０とする。粗さ分布関数が正規分布である場合は、Ｓsk＝０、Ｓku＝３となる。これらのパラメータについても、光学顕微鏡などで自動的に算出された値をそのまま用いても良い。
【００３６】
ヤング率Ｅi、ポアソン比ｖiは材料固有の値であり、垂直抗力の算出には等価ヤング率を用いる。
【００３７】
【数４】
2018152033.tif 000008
・・・（６）
【００３８】
（突起頂点分布への変換）
【００３９】
まず、指定された歪度、及び尖度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を決定する。歪度、及び尖度を制御できる分布関数として、Pearson分布やJohnson分布関数などがあるが、本発明の実施の形態ではJohnson分布を用いる。Johnson分布は歪度、及び尖度の組合せに対応した以下（７）式の２つの異なる式で与えられる。
【００４０】
【数５】
2018152033.tif 000009
・・・（７）
【００４１】
ここで、ξは粗面上の点の無次元高さであり、高さｚを二乗平均平方根粗さで無次元化したものである。上記（７）式のφJSU(ξ)及びφJSB(ξ)はそれぞれJohnsonSU(JSU)分布、JohnsonSB(JSB)分布と呼ばれる。図２に歪度、及び尖度の組合せと、組み合わせに対応する分布の種類を示す。（７）式中のγ、δ、λ、ζは、Johnson分布パラメータと呼ばれ、分布の種類と同様に歪度、及び尖度の組合せによって決まる。数値計算を用いて歪度、及び尖度の組合せに対応する分布の種類とJohnson分布パラメータを決定し、得られた結果をデータベース化することによって、速やかに分布関数を求めることができる。表１に代表的な歪度、及び尖度に対応するパラメータの組合せを示す。また、図３に歪度、及び尖度による分布形状を示す。
【００４２】
【表１】
2018152033.tif 000010
【００４３】
Johnson分布のパラメータの導出方法は、以下（１）〜（５）の手順に従う。
【００４４】
（１）まず、必要なＳsk、Ｓkuの値から、上記図２に示す関係に応じて、ＪＳＢ、又はＪＳＵのいずれの式を用いるか決定する。
（２）γ,δ,λ,ζに適当な初期値を与える。
（３）γ,δ,λ,ζから分布Φを計算し、（４）式、及び（５）式に代入し、Ｓsk、Ｓkuを導出する。この際、（４）式、及び（５）式の無次元化が必要である。
（４）上記（３）で導出したＳsk、Ｓkuと必要なＳsk、Ｓkuの値が異なる場合、γ,δ,λ,ζの値を更新する。
（５）上記（２）〜（４）をＳsk、Ｓkuの値が必要なＳsk、Ｓkuの値に十分近づくまで繰り返す。このとき（２）に関しては適切な値を与えることで収束を早くさせることができる。また、（４）の更新の時は最適化手法に従って更新を行うと効率が良い。最適化手法としては、例えば、焼きなまし法やタブー探索法などが挙げられる。
【００４５】
上記のJohnson分布を用いれば歪度、及び尖度の値に対応する粗さ分布を算出できるが、突起接触モデルで計算に用いるのは突起頂点分布である。図４に示すように一般的には粗さ分布と突起頂点分布とは異なり、突起頂点分布は、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である。すなわち、粗さ分布をJohnson分布によって変換した場合、突起頂点分布はそれに対応して異なる分布形状に変換されることになる。
【００４６】
以下に突起頂点分布の変換方法を示す。Johnson分布関数の累積分布関数φJSU(ξ)及びφJSB(ξ)の正規分布の累積分布関数Φ(ξ)の間には以下（８）式の関係がある。
【００４７】
【数６】
2018152033.tif 000011
・・・（８）
【００４８】
このことは、正規分布に従う粗さ分布における離散点の無次元高さをξ０とし、Johnson変換後の離散点の無次元高さをξとしたときに、両者の間に以下（９）式の関係が成り立つことを意味している。
【００４９】
【数７】
2018152033.tif 000012
・・・（９）
【００５０】
突起頂点分布もJohnson分布と同様に変換されるため、これらの関係を用いて突起頂点分布の変換を求めることが可能である。変換前の突起頂点分布関数をφ０(ξ０)、変換後の突起頂点分布関数をφ(ξ)とする。変換前に高さξに位置していた突起が変換後に高さξ０に移ったとすると、次の関係式が成り立つ。
【００５１】
【数８】
2018152033.tif 000013
・・・（１０）
【００５２】
（１０）式を高さξについて微分すると以下（１１）式が得られる。
【００５３】
【数９】
2018152033.tif 000014
・・・（１１）
【００５４】
よって、変換前の突起頂点分布φ０(ξ)が既知であれば変換後の突起頂点分布関数を解析的に求めることができる。本発明の実施の形態では変換前の突起頂点分布関数として、突起の高さと曲率についての分布を表すＧＳＥモデルを用いる。ＧＳＥモデルの突起頂点分布関数は以下（１２）式により表される。
【００５５】
【数１０】
2018152033.tif 000015
・・・（１２）
【００５６】
ｇは断面曲率σｋで無次元化した突起曲率（無次元曲率）である。この式をｇの全範囲で積分すれば突起頂点高さについての分布となる。αはNayakのパラメータと呼ばれる無次元数で、入力パラメータη、σ、σｋを用いて以下のように表すことができる。
【００５７】
【数１１】
2018152033.tif 000016
・・・（１３）
【００５８】
（突起曲率の変換）
【００５９】
ここで、非特許文献１記載のＧＳＥモデルにおいて、粗面及び平滑面についての真実接触面積Ａｒは以下（１４）式で表される。
【００６０】
【数１２】
2018152033.tif 000017
・・・（１４）
【００６１】
また、垂直抗力Ｗは以下（１５）式で表される。
【００６２】
【数１３】
2018152033.tif 000018
・・・（１５）
【００６３】
〈・〉はＭａｃａｕｌｅｙのブラケットで、以下（１６）式のように定義される。
【００６４】
【数１４】
2018152033.tif 000019
・・・（１６）
【００６５】
ここで、異なる面性状を持つ粗面同士の接触における真実接触面積Ａｒおよび垂直抗力Ｗの導出を行うと、以下（１７）式、及び（１８）式が得られる。
【００６６】
【数１５】
2018152033.tif 000020
・・・（１７）
2018152033.tif 000021
・・・（１８）
【００６７】
ここで、ｄは両粗面の基準面間距離であり、これを合成粗さσ*で無次元化したものを膜厚比と呼ぶ。φ1が、粗面１の突起頂点分布関数であり、φ2が、粗面２の突起頂点分布関数である。
【００６８】
ＧＳＥモデルでは粗面の無次元高さξおよび無次元曲率ｇについての二重積分であるが、本発明の実施の形態では、上記（１７）式、及び（１８）式のように、粗面１、及び粗面２のそれぞれの無次元高さξｉおよび無次元曲率ｇｉについての四重積分である。
【００６９】
（１７）式、及び（１８）式に前工程で算出した突起頂点分布および入力パラメータを代入すれば、真実接触面積および垂直抗力を求めることができる。さらに、本発明では（１７）式、及び（１８）式に以下の２つの手法で曲率についての補正を行う。
【００７０】
第１に、従来の粗さ突起接触モデルはいずれも突起形状を一定の曲率半径を持つ半球と仮定している。しかし、実際の粗さ突起は真球ではなく、先端から徐々に曲率が変化していくものと考えるのが妥当である。そこで、数値的に生成した粗面における突起形状を調べ、それらを平均化した。その結果を図５に示す。突起は先端曲率円よりも広がっていることがわかる。厳密に突起形状の変化を考慮して計算に入れ込むことは困難であるため、ここでは簡易的に先端の曲率半径に一定倍率の補正を掛けるものとする。
【００７１】
第２に、Johnson変換における突起曲率の補正である。Johnson分布は粗さの高さ分布を変換するが、それに伴ってそれぞれの突起の曲率も変化する。例えば、図６に示す通り、変換によって突起が高さ方向に引き伸ばされたような形状になった場合、突起の曲率半径は小さくなる。逆に変換によって高さ方向に縮んだ場合は、突起の曲率半径は大きくなる。突起が引き伸ばされるか、縮められるかは突起の位置する高さによって決まる。突起先端の高さがξ０のとき、変換後の高さξにおける曲率の変化は上記（１１）式のｆ（ξ）に等しい。しかし、先端以外の突起曲率については同様にその位置の高さによって変化するため、突起全体の曲率変化を計算することは非常に困難である。そこで、本発明の実施の形態では突起頂点よりも一定距離下方での高さを用いて変換を行う。
【００７２】
以上の２つの突起曲率についての変換を一つの式で表すと、次のようになる。
【００７３】
【数１６】
2018152033.tif 000022
・・・（１９）
【００７４】
ここで、ｇは変換前の突起曲率、ｇ′は変換後の突起曲率であり、εは実粗面における突起形状の広がりを補正する係数、ｃはJohnson変換における突起曲率変化の基準距離である。
【００７５】
補正した突起曲率を用いて（１７）式、及び（１８）式を修正すると以下（２０）式、及び（２１）式となる。
【００７６】
【数１７】
2018152033.tif 000023
・・・（２０）
2018152033.tif 000024
・・・（２１）
【００７７】
これらを計算すれば、任意の面性状を持つ粗面同士の接触についての真実接触面積および垂直抗力が算出できる。
【００７８】
しかし、（２０）式、及び（２１）式はこのままでは四重積分となり膨大な計算コストが掛かってしまうため、以下（２２）式、及び（２３）式のように演算子を定義する。
【００７９】
【数１８】
2018152033.tif 000025
・・・（２２）
2018152033.tif 000026
・・・（２３）
【００８０】
定義したＦ及びＧの演算子を用いれば、真実接触面積Ａｒ及び垂直抗力Ｗは以下（２４）式、及び（２５）式となる。
【００８１】
【数１９】
2018152033.tif 000027
・・・（２４）
2018152033.tif 000028
・・・（２５）
【００８２】
実際の数値計算ではこれらの演算子をξについてあらかじめ計算しておき、それを代入することで二重積分と同等のコストで真実接触面積Ａｒ及び垂直抗力Ｗを算出することができる。
【００８３】
＜本発明の実施の形態に係る演算装置の構成＞
【００８４】
次に、本発明の実施の形態に係る演算装置の構成について説明する。図７に示すように、本発明の実施の形態に係る演算装置１００は、ＣＰＵと、ＲＡＭと、後述する演算処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この演算装置１００は、機能的には図７に示すように入力部１０と、演算部２０と、出力部５０とを備えている。
【００８５】
入力部１０は、二つの粗面の各々の接触時におけるパワースペクトル密度Ｐi（k）、ヤング率Ｅi、ポアソン比ｖi、及び見かけの接触面積Ａnを受け付ける。また、入力部１０は、二つの粗面の各々の歪度Ｓski及び尖度Ｓkuiを受け付ける。
【００８６】
演算部２０は、データベース２８と、設定部３０と、分布算出部３２と、接触算出部３４とを含んで構成されている。
【００８７】
データベース２８には、上記表１に示すような歪度、及び尖度に対応するJohnson分布のパラメータγ,δ,λ,ζの値が格納されている。
【００８８】
設定部３０は、入力部１０で受け付けた二つの粗面の各々の接触時におけるパワースペクトル密度Ｐi（k）、ヤング率Ｅi、ポアソン比ｖi、見かけの接触面積Ａn、歪度Ｓski、及び尖度Ｓkuiに基づいて、二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する。二乗平均平方根粗さσ、断面曲率σｋi、突起数密度ηiはパワースペクトル密度Ｐi（k）に基づいて上記（２）式で算出して、計算条件とする。また、等価ヤング率Ｅ′を上記（６）式で算出して、計算条件とする。
【００８９】
分布算出部３２は、設定部３０で設定した計算条件に基づいて、二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布φi(ξ)を算出する。ここでは、変換前の突起頂点分布関数φ０(ξ)として、上記（１２）式で表されるＧＳＥモデルを用いて、粗さ分布から、変換後の突起頂点分布φi(ξ)を算出する。ＧＳＥモデルは、突起頂点分布としての、基準面からの突起頂点の高さが各高さであり、かつ、突起頂点の曲率が各曲率である確率を規定する確率密度関数を表す。また、ＧＳＥモデルは、粗さ分布が正規分布によるときの突起頂点分布である。上記（１１）式に従って、ＧＳＥモデルを、粗さ分布が、尖度及び歪度の組み合わせによって定まるパラメータを用いて表される上記（８）式のJohnson分布によるときの突起頂点分布φi(ξ)に変換することにより、突起頂点分布φi(ξ)を算出する。また、上記（１９）式により、突起頂点分布の変換の際に、粗さ分布がJohnson分布によるときの突起頂点分布の変換に応じて、突起頂点の各高さに対する突起頂点の曲率を補正するように、突起頂点分布φi(ξ)を算出する。
【００９０】
接触算出部３４は、上記（２４）式、及び（２５）式に従って、設定部３０で計算条件とされた等価ヤング率Ｅ′及び見かけの接触面積Ａnと、二つの粗面の各々について算出された突起頂点分布φi(ξ)とに基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積Ａｒ及び垂直抗力Ｗを算出する。ここでは、上記（２２）式により突起頂点分布について予め計算された積分積算を表す演算子Ｆ、及び上記（２３）式により二つの突起頂点分布の組み合わせについて予め計算された積分積算を表す演算子Ｇを用いて、真実接触面積Ａｒ及び垂直抗力Ｗを算出する。
【００９１】
＜第１の実施の形態に係る演算装置の作用＞
【００９２】
次に、第１の実施の形態に係る演算装置１００の作用について説明する。入力部１０において二つの粗面の各々の接触時におけるパワースペクトル密度Ｐi（k）、ヤング率Ｅi、ポアソン比ｖi、見かけの接触面積Ａn、歪度Ｓski、及び尖度Ｓkuiを受け付けると、演算装置１００は、図８に示す演算処理ルーチンを実行する。
【００９３】
まず、ステップＳ１００では、入力部１０で受け付けた二つの粗面の各々の接触時におけるパワースペクトル密度Ｐi（k）、ヤング率Ｅi、ポアソン比ｖi、見かけの接触面積Ａn、歪度Ｓski、及び尖度Ｓkuiに基づいて、二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する。
【００９４】
次に、ステップＳ１０２では、入力部１０で受け付けた二つの粗面の各々の歪度Ｓski及び尖度Ｓkuiに基づいて、データベース２８からJohnson分布のパラメータγ,δ,λ,ζの値を読み出す。
【００９５】
ステップＳ１０４では、設定部３０で設定した計算条件と、ステップＳ１０２で読み出したJohnson分布のパラメータγ,δ,λ,ζとに基づいて、上記（８）式に従って、二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、上記（１１）式に従って、粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布φi(ξ)を算出する。
【００９６】
ステップＳ１０６では、上記（２２）式により突起頂点分布について計算された積分積算を表す演算子Ｆを算出し、上記（２３）式により二つの突起頂点分布の組み合わせについて計算された積分積算を表す演算子Ｇを算出する。
【００９７】
ステップＳ１０８では、上記（２４）式、及び（２５）式に従って、入力部１０で受け付けたヤング率Ｅi及び見かけの接触面積Ａnと、二つの粗面の各々について算出された突起頂点分布φi(ξ)と、上記ステップＳ１０６で算出された演算子Ｆ、演算子Ｇとに基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積Ａｒ及び垂直抗力Ｗを算出し、出力部５０に出力して処理を終了する。
【００９８】
（実験結果）
【００９９】
同じ歪度、及び尖度を持つ粗面同士を接触させたときの真実接触面積について、本発明の実施の形態の手法で算出した結果と、直接数値計算で算出した結果との比較を図９に示す。
【０１００】
粗面を数値的に生成する手法として、2D-Weierstrass-Mandelbrot関数（以下、ＷＭ関数）を用いる。ＷＭ関数についての詳細な手法は、例えば非特許文献２の手法を用いることができる。
【０１０１】
［非特許文献２］：Ciavarella, M., Greenwood, J. A., & Paggi, M. (2008). Inclusion of “interaction” in the Greenwood and Williamson contact theory. Wear, 265(5), 729-734.
【０１０２】
ここで、真実接触面積Ａ'ｒは変形を無視した場合での粗面同士のオーバーラップしている部分の投影面積を指しており、弾性変形を仮定した場合の真実接触面積Ａｒの２倍に等しい。図９のように本発明の実施の形態の手法で算出した真実接触面積と、直接数値計算で算出した真実接触面積は非常に良い一致を示す。一方、従来技術では同じ歪度、尖度を持つ粗面同士を接触させるような計算を行うことはできない。また、図１０に同じ歪度、尖度を持つ粗面同士を接触させたときの垂直抗力について示す。なお、図９及び図１０におけるＲｅｇは正規分布の場合を表す。
【０１０３】
以上説明したように、本発明の実施の形態に係る演算装置によれば、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出し、二つの粗面の各々について算出された突起頂点分布に基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出することにより、各々指定された歪度及び尖度を有する粗面同士の接触状況を算出することができる。
【０１０４】
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１０５】
例えば、上述した実施の形態では、（１９）式により、突起頂点の各高さに対する突起頂点の曲率を補正する場合を例に説明したが、これに限定されるものではない。突起頂点の曲率を補正しない場合には、上記ｇ′をｇに置き換えて計算して、真実接触面積Ａｒ及び垂直抗力Ｗを算出すればよい。
【０１０６】
例えば、上述した実施の形態では、真実接触面積Ａｒ及び垂直抗力Ｗを算出する場合を例に説明したが、これに限定されるものではなく、真実接触面積Ａｒ又は垂直抗力Ｗのいずれかを算出するようにしてもよい。
【符号の説明】
【０１０７】
１０   入力部
２０   演算部
２８   データベース
３０   設定部
３２   分布算出部
３４   接触算出部
５０   出力部
１００ 演算装置
【図１】
2018152033.tif 000029
【図２】
2018152033.tif 000030
【図３】
2018152033.tif 000031
【図４】
2018152033.tif 000032
【図５】
2018152033.tif 000033
【図６】
2018152033.tif 000034
【図７】
2018152033.tif 000035
【図８】
2018152033.tif 000036
【図９】
2018152033.tif 000037
【図１０】
2018152033.tif 000038
【図１１】
2018152033.tif 000039