特許 6597716 演算装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6597716

(24)【登録日】2019年10月11日

(45)【発行日】2019年10月30日

(54)【発明の名称】演算装置

(51)【国際特許分類】

   G06F 17/10 20060101AFI20191021BHJP

   G16Z 99/00 20190101ALI20191021BHJP

【ＦＩ】

   G06F17/10 Z

   G16Z99/00

【請求項の数】6

【全頁数】18

(21)【出願番号】特願2017-119810(P2017-119810)

(22)【出願日】2017年6月19日

(65)【公開番号】特開2018-152033(P2018-152033A)

(43)【公開日】2018年9月27日

【審査請求日】2018年10月31日

(31)【優先権主張番号】特願2017-44992(P2017-44992)

(32)【優先日】2017年3月9日

(33)【優先権主張国】JP

(73)【特許権者】

【識別番号】000003609

【氏名又は名称】株式会社豊田中央研究所

(74)【代理人】

【識別番号】100079049

【弁理士】

【氏名又は名称】中島 淳

(74)【代理人】

【識別番号】100084995

【弁理士】

【氏名又は名称】加藤 和詳

(72)【発明者】

【氏名】友田 達規

(72)【発明者】

【氏名】近藤 靖裕

(72)【発明者】

【氏名】大森 俊英

【審査官】 加藤 優一

(56)【参考文献】

【文献】 特表２０００−５００６０４（ＪＰ，Ａ）

【文献】 特開２００６−０２１５５１（ＪＰ，Ａ）

【文献】 特開２０００−２０５２７５（ＪＰ，Ａ）

【文献】 中国特許出願公開第１０２６０９５６０（ＣＮ，Ａ）

【文献】 野々垣 稔，“種々の接触理論による接触近似式表面形状パラメータの考察”，奈良工業高等専門学校 研究紀要，日本，奈良工業高等専門学校，２０００年 ３月１５日，第35号，第1頁-第4頁，ISSN 0387-1150，ＵＲＬ，https://www.nara-k.ac.jp/nnct-library/publication/pdf/h11kiyo_all.pdf

【文献】 久門 輝正，“固体仕上面間の接触機構 （第５報，突起の弾性変形を考慮した場合）”，日本機械学會論文集，日本，日本機械学會，１９７２年１０月２５日，38巻, 314号，第2657頁-第2665頁，ISSN 2185-9485，ＵＲＬ，https://www.jstage.jst.go.jp/article/kikai1938/38/314/38_314_2657/_article/-char/ja/

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １７／１０

Ｇ１６Ｚ ９９／００

Ｇ０１Ｂ ２１／３０

Ｇ０１Ｂ ７／３４

Ｇ１１Ｂ ５／８２

Ｇ１１Ｂ ５／８４

(57)【特許請求の範囲】

【請求項1】

二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置であって、
前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部と、
前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部と、
前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部と、
を含む演算装置。

【請求項2】

前記分布算出部は、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布を算出する請求項１に記載の演算装置。

【請求項3】

前記分布算出部は、前記突起頂点分布としての、基準面からの突起頂点の高さが各高さであり、かつ、突起頂点の曲率が各曲率である確率を規定する確率密度関数を表すGreenwood simplified ellipticモデルを、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布に変換することにより、前記突起頂点分布を算出する請求項２に記載の演算装置。

【請求項4】

前記分布算出部は、前記突起頂点分布の変換の際に、前記粗さ分布がJohnson分布によるときの前記突起頂点分布の変換に応じて、突起頂点の各高さに対する突起頂点の曲率を補正するように、前記突起頂点分布を算出する請求項３記載の演算装置。

【請求項5】

前記接触算出部は、前記突起頂点分布について予め計算された積分積算を表す演算子、及び二つの前記突起頂点分布の組み合わせについて予め計算された積分積算を表す演算子を用いて、前記真実接触面積又は垂直抗力を算出する請求項１〜請求項４の何れか１項に記載の演算装置。

【請求項6】

二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置を制御するためのプログラムであって、
コンピュータを、前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部、
前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部、及び
前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部、
として機能させるためのプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、演算装置、及びプログラムに係り、特に、面同士の接触状況を算出するための演算装置、及びプログラムに関する。

【背景技術】

【0002】

従来技術では、面同士の接触について、Greenwood and Williamson(1966)のモデル（ＧＷモデル）を用いて接触面積および垂直抗力を算出している（特許文献１参照）。特許文献１には、運動摩擦を最小限に抑える歪度値および尖度値を決定するための接触モデリングの確率密度関数を用いた、表面粗さ分布の設計方法が記載されている。このような従来技術の手法は、簡易計算としては有用な手法である。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特表２０００−５００６０４号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

しかし、特許文献１に記載されているような従来技術には、以下の問題点により、その計算精度と適用対象の範囲が不十分といえる。

【0005】

一つ目の問題は、ＧＷモデルでは、図１１に示すように、粗面と平滑面との接触を仮定しているため、粗面同士の接触を直接取り扱うことができないことが挙げられる。両粗面の突起頂点分布が確率論や統計学で用いられる正規分布であると仮定した場合は、粗面同士の接触であっても合成粗面と平滑面の接触とみなすことができるが、突起頂点分布が非正規分布の場合は従来手法を用いることができない。よって、従来技術では粗面及び平滑面を対象とした限定された範囲の計算しか行えない。

【0006】

二つ目の問題は、ＧＷモデルでは突起頂点分布を正規分布としているため、従来技術では突起頂点分布をPearson分布に置き換えて計算していることである。従来の解析（Nayak(1971)）から明らかなように、粗さ分布が正規分布であっても突起頂点分布は一般には正規分布とはならない。したがって、粗さ分布が正規分布の突起頂点分布を直接Pearson分布で変換とするのは誤りである。また、粗面の持つ歪度、及び尖度は粗さ分布について与えられるものであり、突起頂点分布を指定された歪度、及び尖度になるように変換するのは不適切である。

【0007】

本発明は、上記問題点を解決するために成されたものであり、各々指定された歪度及び尖度を有する粗面同士の接触状況を算出することができる演算装置、方法及びプログラムを提供することを目的とする。

【課題を解決するための手段】

【0008】

上記目的を達成するために、第１の発明に係る演算装置は、二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置であって、前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部と、前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部と、前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部と、を含んで構成されている。

【0009】

また、第１の発明に係る演算装置において、前記分布算出部は、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布を算出するようにしてもよい。

【0010】

また、第１の発明に係る演算装置において、前記分布算出部は、前記突起頂点分布としての、基準面からの突起頂点の高さが各高さであり、かつ、突起頂点の曲率が各曲率である確率を規定する確率密度関数を表すGreenwood simplified ellipticモデルを、前記粗さ分布が、前記尖度及び歪度の組み合わせによって定まるパラメータを用いて表されるJohnson分布によるときの、前記突起頂点分布に変換することにより、前記突起頂点分布を算出するようにしてもよい。

【0011】

また、第１の発明に係る演算装置において、前記分布算出部は、前記突起頂点分布の変換の際に、前記粗さ分布がJohnson分布によるときの前記突起頂点分布の変換に応じて、突起頂点の各高さに対する突起頂点の曲率を補正するように、前記突起頂点分布を算出するようにしてもよい。

【0012】

また、第１の発明に係る演算装置において、前記接触算出部は、前記突起頂点分布について予め計算された積分積算を表す演算子、及び二つの前記突起頂点分布の組み合わせについて予め計算された積分積算を表す演算子を用いて、前記真実接触面積又は垂直抗力を算出するようにしてもよい。

【0013】

第２の発明に係るプログラムは、二つの粗面の接触時に発生する接触面積又は垂直抗力を予測する演算装置を制御するためのプログラムであって、コンピュータを、前記二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する設定部、前記計算条件に基づいて、前記二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出する分布算出部、及び前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出する接触算出部、として機能させるためのプログラムである。

【発明の効果】

【0014】

本発明の演算装置、及びプログラムによれば、二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、前記粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出し、前記二つの粗面の各々について算出された前記突起頂点分布に基づいて、前記二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出することにより、各々指定された歪度及び尖度を有する粗面同士の接触状況を算出することができる。

【図面の簡単な説明】

【0015】

【図1】粗面と粗面との接触を仮定した本実施の形態の一例を示す図である。

【図2】歪度、及び尖度の組合せと、組み合わせに対応する分布の種類の一例を示す図である。

【図3】歪度、及び尖度による分布形状の一例を示す図である。

【図4】粗さ分布と突起頂点分布の一例を示す図である。

【図5】数値的に生成した粗面における平均突起形状の一例を示す図である。

【図6】変換による突起形状の変化の一例を示す図である。

【図7】本発明の実施の形態に係る演算装置の構成を示すブロック図である。

【図8】本発明の実施の形態に係る演算装置における演算処理ルーチンを示すフローチャートである。

【図9】実験における膜圧比に対する真実接触面積の変化を示す図である。

【図10】実験における膜圧比に対する真実接触面積の変化を示す図である。

【図11】粗面と平滑面との接触を仮定した従来技術の一例を示す図である。

【発明を実施するための形態】

【0016】

以下、図面を参照して本発明の実施の形態を詳細に説明する。

【0017】

＜本発明の実施の形態に係る概要＞

【0018】

まず、本発明の実施の形態における概要を説明する。

【0019】

上記課題の一つ目の問題に対し、本発明の実施の形態では、粗さ突起接触モデルであるGreenwood and Trippモデル(ＧＴモデル)の考え方を用いることによって、図１に示すように、非正規分布を持つ粗面同士の接触を取り扱うことができる。これらの粗面は互いに歪度、及び尖度が異なっていても計算可能である。

【0020】

上記課題の二つ目の問題に対し、本発明の実施の形態では、粗さ分布が指定された歪度、及び尖度に変換されるときに、歪度、及び尖度に従属する形で突起頂点分布を変換する。変換前の突起頂点分布はGreenwood(2006)のsimplified ellipticモデル(ＧＳＥモデル)で与えられる（非特許文献１参照）。この計算手順については、以下の原理の説明で詳細を示す。

【0021】

［非特許文献１］：Greenwood, J. A. (2006). A simplified elliptic model of rough surface contact. Wear, 261(2), 191-200.

【0022】

ＧＴモデルをベースに、以下の改良を施すことで精度を大幅に向上させる。ＧＴモデルでは突起頂点分布を正規分布とし、突起曲率半径を一定と仮定しており、このままでは精度が不十分なため、より実粗面に近い突起頂点分布および突起曲率半径分布を用いているＧＳＥモデルを用いてＧＴモデルの式を改良する。

【0023】

また、本発明の実施の形態では、取り扱いやすさと精度を考慮し、非正規分布としてPearson分布に代わりJohnson分布を用いる。上記の改良では粗さ分布が正規分布である粗面しか扱えないため、歪度、尖度に特徴を持つ粗面について取り扱うことができない。そこで、粗さ分布にJohnson分布を用いて任意の歪度、尖度を持つ粗面を取り扱えるようにした。粗さ分布にJohnson分布を用いた場合の突起頂点分布の変換について厳密に計算を行い、曲率半径の変化についても補正を与えることによって精度を維持する。

【0024】

また、既存の突起接触モデルではＨｅｒｔｚの接触理論に基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積、及び垂直抗力を計算している。このとき、突起を一定の曲率を持つ半球であると仮定しているが、実際の突起は先端から離れるにしたがって曲率が小さくなる。そこで、本発明の実施の形態では突起の曲率変化を考慮して、元の接触モデルにおける曲率に一定倍率の補正を与える。

【0025】

以上の改良によって、検討可能な粗面の条件が大幅に拡張され、さらに直接数値計算と比較しても非常に良好な真実接触面積の予測が可能となる。

【0026】

＜本発明の実施の形態に係る原理＞

【0027】

次に、本発明の実施の形態における原理を説明する。

【0028】

本発明の演算装置に入力するパラメータは、η_i：突起数密度、σ_ki：断面曲率又はβ_i：平均突起曲率半径、σ_i：二乗平均平方根粗さ、Ｓ_ski：歪度（スキューネス）、Ｓ_kui：尖度（クルトシス）、Ｅ_i：ヤング率（縦弾性係数）、ｖ_i：ポアソン比、Ａ_n：見かけの接触面積である。添字のｉは粗面１、及び粗面２それぞれの持つ物理量であることを示す。以下に、これらのパラメータについて概要を述べる。

【0029】

以下の原理の説明では、等方粗さを持つ粗面を仮定する。以下（１）式の粗さ断面曲線のパワースペクトル密度Ｐ（k）（k：波数）のモーメントを用いて、以下（２）式で二乗平均平方根粗さσ、断面曲率σ_ｋ、突起数密度ηは算出される。

【0030】

【数1】

・・・（１）

・・・（２）

【0031】

上記（１）式、（２）式のパラメータをパワースペクトル密度から算出する以外に、これらのパラメータに実測した値を直接用いても良い。本発明の実施の形態で援用するＧＳＥモデルの場合、実測で得られる平均突起曲率半径βと断面曲率σ_ｋとの間には以下（３）式の関係がある。

【0032】

【数2】

・・・（３）

【0033】

３次元の離散点で表現される粗面において、ある点の基準面からの高さがｚである確率が確率密度関数（以下、粗さ分布関数）φ(z)で規定されるとき、歪度Ｓ_sk、及び尖度Ｓ_kuは以下（４）式、及び（５）式となる。

【0034】

【数3】

・・・（４）

・・・（５）

【0035】

ここで、ｍは粗さ基準面高さであり、一般にデータ処理前に取り除かれるため、以降はｍ＝０とする。粗さ分布関数が正規分布である場合は、Ｓ_sk＝０、Ｓ_ku＝３となる。これらのパラメータについても、光学顕微鏡などで自動的に算出された値をそのまま用いても良い。

【0036】

ヤング率Ｅ_i、ポアソン比ｖ_iは材料固有の値であり、垂直抗力の算出には等価ヤング率を用いる。

【0037】

【数4】

・・・（６）

【0038】

（突起頂点分布への変換）

【0039】

まず、指定された歪度、及び尖度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を決定する。歪度、及び尖度を制御できる分布関数として、Pearson分布やJohnson分布関数などがあるが、本発明の実施の形態ではJohnson分布を用いる。Johnson分布は歪度、及び尖度の組合せに対応した以下（７）式の２つの異なる式で与えられる。

【0040】

【数5】

・・・（７）

【0041】

ここで、ξは粗面上の点の無次元高さであり、高さｚを二乗平均平方根粗さで無次元化したものである。上記（７）式のφ_JSU(ξ)及びφ_JSB(ξ)はそれぞれJohnsonSU(JSU)分布、JohnsonSB(JSB)分布と呼ばれる。図２に歪度、及び尖度の組合せと、組み合わせに対応する分布の種類を示す。（７）式中のγ、δ、λ、ζは、Johnson分布パラメータと呼ばれ、分布の種類と同様に歪度、及び尖度の組合せによって決まる。数値計算を用いて歪度、及び尖度の組合せに対応する分布の種類とJohnson分布パラメータを決定し、得られた結果をデータベース化することによって、速やかに分布関数を求めることができる。表１に代表的な歪度、及び尖度に対応するパラメータの組合せを示す。また、図３に歪度、及び尖度による分布形状を示す。

【0042】

【表1】

【0043】

Johnson分布のパラメータの導出方法は、以下（１）〜（５）の手順に従う。

【0044】

（１）まず、必要なＳ_sk、Ｓ_kuの値から、上記図２に示す関係に応じて、ＪＳＢ、又はＪＳＵのいずれの式を用いるか決定する。
（２）γ,δ,λ,ζに適当な初期値を与える。
（３）γ,δ,λ,ζから分布Φを計算し、（４）式、及び（５）式に代入し、Ｓ_sk、Ｓ_kuを導出する。この際、（４）式、及び（５）式の無次元化が必要である。
（４）上記（３）で導出したＳ_sk、Ｓ_kuと必要なＳ_sk、Ｓ_kuの値が異なる場合、γ,δ,λ,ζの値を更新する。
（５）上記（２）〜（４）をＳ_sk、Ｓ_kuの値が必要なＳ_sk、Ｓ_kuの値に十分近づくまで繰り返す。このとき（２）に関しては適切な値を与えることで収束を早くさせることができる。また、（４）の更新の時は最適化手法に従って更新を行うと効率が良い。最適化手法としては、例えば、焼きなまし法やタブー探索法などが挙げられる。

【0045】

上記のJohnson分布を用いれば歪度、及び尖度の値に対応する粗さ分布を算出できるが、突起接触モデルで計算に用いるのは突起頂点分布である。図４に示すように一般的には粗さ分布と突起頂点分布とは異なり、突起頂点分布は、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である。すなわち、粗さ分布をJohnson分布によって変換した場合、突起頂点分布はそれに対応して異なる分布形状に変換されることになる。

【0046】

以下に突起頂点分布の変換方法を示す。Johnson分布関数の累積分布関数φ_JSU(ξ)及びφ_JSB(ξ)の正規分布の累積分布関数Φ(ξ)の間には以下（８）式の関係がある。

【0047】

【数6】

・・・（８）

【0048】

このことは、正規分布に従う粗さ分布における離散点の無次元高さをξ_０とし、Johnson変換後の離散点の無次元高さをξとしたときに、両者の間に以下（９）式の関係が成り立つことを意味している。

【0049】

【数7】

・・・（９）

【0050】

突起頂点分布もJohnson分布と同様に変換されるため、これらの関係を用いて突起頂点分布の変換を求めることが可能である。変換前の突起頂点分布関数をφ_０(ξ_０)、変換後の突起頂点分布関数をφ(ξ)とする。変換前に高さξに位置していた突起が変換後に高さξ_０に移ったとすると、次の関係式が成り立つ。

【0051】

【数8】

・・・（１０）

【0052】

（１０）式を高さξについて微分すると以下（１１）式が得られる。

【0053】

【数9】

・・・（１１）

【0054】

よって、変換前の突起頂点分布φ_０(ξ)が既知であれば変換後の突起頂点分布関数を解析的に求めることができる。本発明の実施の形態では変換前の突起頂点分布関数として、突起の高さと曲率についての分布を表すＧＳＥモデルを用いる。ＧＳＥモデルの突起頂点分布関数は以下（１２）式により表される。

【0055】

【数10】

・・・（１２）

【0056】

ｇは断面曲率σ_ｋで無次元化した突起曲率（無次元曲率）である。この式をｇの全範囲で積分すれば突起頂点高さについての分布となる。αはNayakのパラメータと呼ばれる無次元数で、入力パラメータη、σ、σ_ｋを用いて以下のように表すことができる。

【0057】

【数11】

・・・（１３）

【0058】

（突起曲率の変換）

【0059】

ここで、非特許文献１記載のＧＳＥモデルにおいて、粗面及び平滑面についての真実接触面積Ａ_ｒは以下（１４）式で表される。

【0060】

【数12】

・・・（１４）

【0061】

また、垂直抗力Ｗは以下（１５）式で表される。

【0062】

【数13】

・・・（１５）

【0063】

〈・〉はＭａｃａｕｌｅｙのブラケットで、以下（１６）式のように定義される。

【0064】

【数14】

・・・（１６）

【0065】

ここで、異なる面性状を持つ粗面同士の接触における真実接触面積Ａ_ｒおよび垂直抗力Ｗの導出を行うと、以下（１７）式、及び（１８）式が得られる。

【0066】

【数15】

・・・（１７）

・・・（１８）

【0067】

ここで、ｄは両粗面の基準面間距離であり、これを合成粗さσ_*で無次元化したものを膜厚比と呼ぶ。φ1が、粗面１の突起頂点分布関数であり、φ2が、粗面２の突起頂点分布関数である。

【0068】

ＧＳＥモデルでは粗面の無次元高さξおよび無次元曲率ｇについての二重積分であるが、本発明の実施の形態では、上記（１７）式、及び（１８）式のように、粗面１、及び粗面２のそれぞれの無次元高さξ_ｉおよび無次元曲率ｇ_ｉについての四重積分である。

【0069】

（１７）式、及び（１８）式に前工程で算出した突起頂点分布および入力パラメータを代入すれば、真実接触面積および垂直抗力を求めることができる。さらに、本発明では（１７）式、及び（１８）式に以下の２つの手法で曲率についての補正を行う。

【0070】

第１に、従来の粗さ突起接触モデルはいずれも突起形状を一定の曲率半径を持つ半球と仮定している。しかし、実際の粗さ突起は真球ではなく、先端から徐々に曲率が変化していくものと考えるのが妥当である。そこで、数値的に生成した粗面における突起形状を調べ、それらを平均化した。その結果を図５に示す。突起は先端曲率円よりも広がっていることがわかる。厳密に突起形状の変化を考慮して計算に入れ込むことは困難であるため、ここでは簡易的に先端の曲率半径に一定倍率の補正を掛けるものとする。

【0071】

第２に、Johnson変換における突起曲率の補正である。Johnson分布は粗さの高さ分布を変換するが、それに伴ってそれぞれの突起の曲率も変化する。例えば、図６に示す通り、変換によって突起が高さ方向に引き伸ばされたような形状になった場合、突起の曲率半径は小さくなる。逆に変換によって高さ方向に縮んだ場合は、突起の曲率半径は大きくなる。突起が引き伸ばされるか、縮められるかは突起の位置する高さによって決まる。突起先端の高さがξ_０のとき、変換後の高さξにおける曲率の変化は上記（１１）式のｆ（ξ）に等しい。しかし、先端以外の突起曲率については同様にその位置の高さによって変化するため、突起全体の曲率変化を計算することは非常に困難である。そこで、本発明の実施の形態では突起頂点よりも一定距離下方での高さを用いて変換を行う。

【0072】

以上の２つの突起曲率についての変換を一つの式で表すと、次のようになる。

【0073】

【数16】

・・・（１９）

【0074】

ここで、ｇは変換前の突起曲率、ｇ′は変換後の突起曲率であり、εは実粗面における突起形状の広がりを補正する係数、ｃはJohnson変換における突起曲率変化の基準距離である。

【0075】

補正した突起曲率を用いて（１７）式、及び（１８）式を修正すると以下（２０）式、及び（２１）式となる。

【0076】

【数17】

・・・（２０）

・・・（２１）

【0077】

これらを計算すれば、任意の面性状を持つ粗面同士の接触についての真実接触面積および垂直抗力が算出できる。

【0078】

しかし、（２０）式、及び（２１）式はこのままでは四重積分となり膨大な計算コストが掛かってしまうため、以下（２２）式、及び（２３）式のように演算子を定義する。

【0079】

【数18】

・・・（２２）

・・・（２３）

【0080】

定義したＦ及びＧの演算子を用いれば、真実接触面積Ａ_ｒ及び垂直抗力Ｗは以下（２４）式、及び（２５）式となる。

【0081】

【数19】

・・・（２４）

・・・（２５）

【0082】

実際の数値計算ではこれらの演算子をξについてあらかじめ計算しておき、それを代入することで二重積分と同等のコストで真実接触面積Ａ_ｒ及び垂直抗力Ｗを算出することができる。

【0083】

＜本発明の実施の形態に係る演算装置の構成＞

【0084】

次に、本発明の実施の形態に係る演算装置の構成について説明する。図７に示すように、本発明の実施の形態に係る演算装置１００は、ＣＰＵと、ＲＡＭと、後述する演算処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この演算装置１００は、機能的には図７に示すように入力部１０と、演算部２０と、出力部５０とを備えている。

【0085】

入力部１０は、二つの粗面の各々の接触時におけるパワースペクトル密度Ｐ_i（k）、ヤング率Ｅ_i、ポアソン比ｖ_i、及び見かけの接触面積Ａ_nを受け付ける。また、入力部１０は、二つの粗面の各々の歪度Ｓ_ski及び尖度Ｓ_kuiを受け付ける。

【0086】

演算部２０は、データベース２８と、設定部３０と、分布算出部３２と、接触算出部３４とを含んで構成されている。

【0087】

データベース２８には、上記表１に示すような歪度、及び尖度に対応するJohnson分布のパラメータγ,δ,λ,ζの値が格納されている。

【0088】

設定部３０は、入力部１０で受け付けた二つの粗面の各々の接触時におけるパワースペクトル密度Ｐ_i（k）、ヤング率Ｅ_i、ポアソン比ｖ_i、見かけの接触面積Ａ_n、歪度Ｓ_ski、及び尖度Ｓ_kuiに基づいて、二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する。二乗平均平方根粗さσ、断面曲率σ_ｋi、突起数密度η_iはパワースペクトル密度Ｐ_i（k）に基づいて上記（２）式で算出して、計算条件とする。また、等価ヤング率Ｅ′を上記（６）式で算出して、計算条件とする。

【0089】

分布算出部３２は、設定部３０で設定した計算条件に基づいて、二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布φ_i(ξ)を算出する。ここでは、変換前の突起頂点分布関数φ_０(ξ)として、上記（１２）式で表されるＧＳＥモデルを用いて、粗さ分布から、変換後の突起頂点分布φ_i(ξ)を算出する。ＧＳＥモデルは、突起頂点分布としての、基準面からの突起頂点の高さが各高さであり、かつ、突起頂点の曲率が各曲率である確率を規定する確率密度関数を表す。また、ＧＳＥモデルは、粗さ分布が正規分布によるときの突起頂点分布である。上記（１１）式に従って、ＧＳＥモデルを、粗さ分布が、尖度及び歪度の組み合わせによって定まるパラメータを用いて表される上記（８）式のJohnson分布によるときの突起頂点分布φ_i(ξ)に変換することにより、突起頂点分布φ_i(ξ)を算出する。また、上記（１９）式により、突起頂点分布の変換の際に、粗さ分布がJohnson分布によるときの突起頂点分布の変換に応じて、突起頂点の各高さに対する突起頂点の曲率を補正するように、突起頂点分布φ_i(ξ)を算出する。

【0090】

接触算出部３４は、上記（２４）式、及び（２５）式に従って、設定部３０で計算条件とされた等価ヤング率Ｅ′及び見かけの接触面積Ａ_nと、二つの粗面の各々について算出された突起頂点分布φ_i(ξ)とに基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積Ａ_ｒ及び垂直抗力Ｗを算出する。ここでは、上記（２２）式により突起頂点分布について予め計算された積分積算を表す演算子Ｆ、及び上記（２３）式により二つの突起頂点分布の組み合わせについて予め計算された積分積算を表す演算子Ｇを用いて、真実接触面積Ａ_ｒ及び垂直抗力Ｗを算出する。

【0091】

＜第１の実施の形態に係る演算装置の作用＞

【0092】

次に、第１の実施の形態に係る演算装置１００の作用について説明する。入力部１０において二つの粗面の各々の接触時におけるパワースペクトル密度Ｐ_i（k）、ヤング率Ｅ_i、ポアソン比ｖ_i、見かけの接触面積Ａ_n、歪度Ｓ_ski、及び尖度Ｓ_kuiを受け付けると、演算装置１００は、図８に示す演算処理ルーチンを実行する。

【0093】

まず、ステップＳ１００では、入力部１０で受け付けた二つの粗面の各々の接触時におけるパワースペクトル密度Ｐ_i（k）、ヤング率Ｅ_i、ポアソン比ｖ_i、見かけの接触面積Ａ_n、歪度Ｓ_ski、及び尖度Ｓ_kuiに基づいて、二つの粗面の各々の尖度及び歪度を含む面性状を計算条件として設定する。

【0094】

次に、ステップＳ１０２では、入力部１０で受け付けた二つの粗面の各々の歪度Ｓ_ski及び尖度Ｓ_kuiに基づいて、データベース２８からJohnson分布のパラメータγ,δ,λ,ζの値を読み出す。

【0095】

ステップＳ１０４では、設定部３０で設定した計算条件と、ステップＳ１０２で読み出したJohnson分布のパラメータγ,δ,λ,ζとに基づいて、上記（８）式に従って、二つの粗面の各々について、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、上記（１１）式に従って、粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布φ_i(ξ)を算出する。

【0096】

ステップＳ１０６では、上記（２２）式により突起頂点分布について計算された積分積算を表す演算子Ｆを算出し、上記（２３）式により二つの突起頂点分布の組み合わせについて計算された積分積算を表す演算子Ｇを算出する。

【0097】

ステップＳ１０８では、上記（２４）式、及び（２５）式に従って、入力部１０で受け付けたヤング率Ｅ_i及び見かけの接触面積Ａ_nと、二つの粗面の各々について算出された突起頂点分布φ_i(ξ)と、上記ステップＳ１０６で算出された演算子Ｆ、演算子Ｇとに基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積Ａ_ｒ及び垂直抗力Ｗを算出し、出力部５０に出力して処理を終了する。

【0098】

（実験結果）

【0099】

同じ歪度、及び尖度を持つ粗面同士を接触させたときの真実接触面積について、本発明の実施の形態の手法で算出した結果と、直接数値計算で算出した結果との比較を図９に示す。

【0100】

粗面を数値的に生成する手法として、2D-Weierstrass-Mandelbrot関数（以下、ＷＭ関数）を用いる。ＷＭ関数についての詳細な手法は、例えば非特許文献２の手法を用いることができる。

【0101】

［非特許文献２］：Ciavarella, M., Greenwood, J. A., & Paggi, M. (2008). Inclusion of “interaction” in the Greenwood and Williamson contact theory. Wear, 265(5), 729-734.

【0102】

ここで、真実接触面積Ａ'_ｒは変形を無視した場合での粗面同士のオーバーラップしている部分の投影面積を指しており、弾性変形を仮定した場合の真実接触面積Ａ_ｒの２倍に等しい。図９のように本発明の実施の形態の手法で算出した真実接触面積と、直接数値計算で算出した真実接触面積は非常に良い一致を示す。一方、従来技術では同じ歪度、尖度を持つ粗面同士を接触させるような計算を行うことはできない。また、図１０に同じ歪度、尖度を持つ粗面同士を接触させたときの垂直抗力について示す。なお、図９及び図１０におけるＲｅｇは正規分布の場合を表す。

【0103】

以上説明したように、本発明の実施の形態に係る演算装置によれば、尖度及び歪度の組み合わせに対応する、基準面からの高さが各高さである確率を規定する確率密度関数である粗さ分布を取得し、粗さ分布から、基準面からの突起頂点の高さが各高さである確率を規定する確率密度関数である突起頂点分布を算出し、二つの粗面の各々について算出された突起頂点分布に基づいて、二つの粗面の接触時における実際の接触面積である真実接触面積又は垂直抗力を算出することにより、各々指定された歪度及び尖度を有する粗面同士の接触状況を算出することができる。

【0104】

なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

【0105】

例えば、上述した実施の形態では、（１９）式により、突起頂点の各高さに対する突起頂点の曲率を補正する場合を例に説明したが、これに限定されるものではない。突起頂点の曲率を補正しない場合には、上記ｇ′をｇに置き換えて計算して、真実接触面積Ａ_ｒ及び垂直抗力Ｗを算出すればよい。

【0106】

例えば、上述した実施の形態では、真実接触面積Ａ_ｒ及び垂直抗力Ｗを算出する場合を例に説明したが、これに限定されるものではなく、真実接触面積Ａ_ｒ又は垂直抗力Ｗのいずれかを算出するようにしてもよい。

【符号の説明】

【0107】

１０ 入力部
２０ 演算部
２８ データベース
３０ 設定部
３２ 分布算出部
３４ 接触算出部
５０ 出力部
１００ 演算装置

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】