(57)【特許請求の範囲】
【請求項１】
外力が作用したときに互いに共形接触を行う２つの接触物体の接触面における面圧を、コンピュータを用いて計算する方法であって、
前記接触物体の接触面の形状、弾性係数、ポアソン比、および作用する外力を入力して記憶手段に記憶させる入力過程と、入力された値を用いて前記面圧を計算する解析過程と、計算した結果を出力する出力過程とでなり、
前記解析過程として、
前記２つの接触物体の互いに接する面の間に、基準曲面を設定してこの基準曲面をセルが縦横に並ぶメッシュに分割し、接触物体間の距離を計算する点である計算点をセル上に定めるメッシュ分割過程と、
各セルに単位分布荷重が作用したときの前記２つの接触物体の表面の前記各計算点に対応する点に生じる変位である影響量を算出する変位影響算出過程と、
各セルにつき面圧と前記単位分布荷重による影響量の積の合計が接触領域内部において接触物体間の、前記単位分布荷重が与えられることにより互いに近づく方向に移動する距離である接近量と、前記単位分布荷重が与えられる前の距離である初期隙間の差に釣り合い、かつ各セルの面圧の合計が外力と釣り合うという釣り合い条件式から、各セルの面圧を算出する接触面圧計算過程と、
この解析過程の途中で、前記メッシュに分割する大きさを修正してその修正結果により、メッシュ分割後の前記各過程を再度行うメッシュ修正過程とを含み、
前記変位影響算出過程および接触面圧計算過程では、計測点間の距離の算出に上記基準曲面上の距離を用い、かつこの基準曲面上の位置の違いによる法線方向の修正を行って面圧を算出する、
ことを特徴とする共形接触下の接触面圧計算方法。
【請求項２】
請求項１において、前記２つの接触物体が転がり軸受の内輪とボールとであって、
前記メッシュ分割過程では、前記基準曲面を円筒面とし、
前記変位影響量算出過程では、前記円筒面からなる基準曲面上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をｘ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面上の法線方向に離れる座標位置をｚ、前記ｚ軸方向の変位をｗとして、
前記外力が法線方向の力である場合は、次式（１）で表される単位集中荷重作用時の変位解析解により前記影響量を計算する、
【数１】
0005026225.tif 000018
θ：前記座標位置ｘをｙ軸回りの半径ＲO の極座標で表す角度
ただし、（１）式により与えられる変位ｗは位置（θ，ｙ）の要素が位置（θ′，ｙ′）に及ぼす変位、
共形接触下の接触面圧計算方法。
【請求項３】
請求項２において、前記接触面圧計算過程では、前記（１）式で計算されるｚ軸方向の変位をｗを全セルについて重ね合わせる式である次式（２）を、前記釣り合い条件式の一つとして用いて各セルの面圧分布を求め、
【数２】
0005026225.tif 000019
かつ、前記釣り合い条件式の他の一つとして次式（３）を用いて前記接近量αを求める、
【数３】
0005026225.tif 000020
共形接触下の接触面圧計算方法。
【請求項４】
請求項１ないし請求項３のいずれか１項において、前記メッシュ修正過程は、前記接触面圧計算過程の後に行う共形接触下の接触面圧計算方法。
【請求項５】
外力が作用したときに互いに共形接触を行う２つの接触物体の接触面における面圧および表面下応力を、コンピュータを用いて計算する方法であって、
請求項１ないし請求項４のいずれか１項に記載の接触面圧の計算方法と、この計算方法で各セルの面圧を計算した後に、表面下応力を計算する過程を含み、
この表面下応力計算過程として、
前記外力が法線方向の力である場合は、次式（４）〜（９）で与えられる単位集中荷重作用時の表面下の応力を用いて、個々のセルの全域にわたり積分して、かつ全セルの影響の和をとって表面下応力を求める、
【数４】
0005026225.tif 000021
共形接触下の接触面圧および表面下応力計算方法。
【請求項６】
請求項５において、前記２つの接触物体が転がり軸受の内輪とボールとであって、
前記メッシュ分割過程では、前記基準曲面を円筒面とし、
前記変位影響量算出過程では、前記円筒面からなる基準曲面上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をθ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面上の法線方向に離れる座標位置をｚとして、
位置（θ，ｙ）の要素ｅが位置（θ’，ｙ’，ｚ）に及ぼす単位面圧負荷時の剪断応力τzyを、次式（１０）で求める、
【数５】
0005026225.tif 000022
共形接触下の接触面圧および表面下応力計算方法。
【請求項７】
外力が作用したときに互いに共形接触を行う２つの接触物体の接触面における面圧を計算する装置であって、
前記接触物体の接触面の形状、弾性係数、ポアソン比、および作用する外力を入力して記憶手段に記憶させる入力処理手段と、入力された値を用いて前記面圧を計算する解析手段と、計算した結果を出力する出力処理手段とでなり、
前記解析手段として、
前記２つの接触物体の互いに接する面の間に、基準曲面を設定してこの基準曲面をセルが縦横に並ぶメッシュに分割し、接触物体間の距離を計算する点である計算点をセル上に定めるメッシュ分割手段と、
各セルに単位分布荷重が作用したときの前記２つの接触物体の表面の前記各計算点に生じる変位である影響量を算出する変位影響算出手段と、
各セルにつき面圧と前記単位分布荷重による影響量の積の合計が接触領域内部において接触物体間の、前記単位分布荷重が与えられることにより互いに近づく方向に移動する距離である接近量と、前記単位分布荷重が与えられる前の距離である初期隙間の差に釣り合い、かつ各セルの面圧の合計が外力と釣り合うという釣り合い条件式から、各セルの面圧を算出する接触面圧計算手段と、
この接触面圧計算手段による前記各セルの面圧の算出よりも後に、前記メッシュ分割手段で前記メッシュに分割する大きさを修正してその修正結果により、前記変位影響算出手段および前記接触面圧計算手段による算出を再度行わせるメッシュ修正手段とを含み、
前記変位影響算出手段および接触面圧計算手段は、計測点間の距離の算出に上記基準曲面上の距離を用い、かつこの基準曲面上の位置の違いによる法線方向の修正を行って面圧を算出する、
ことを特徴とする共形接触下の接触面圧計算装置。
【請求項８】
請求項７において、前記２つの接触物体が転がり軸受の内輪とボールとであって、
前記メッシュ分割手段は、前記基準曲面を円筒面とし、
前記変位影響量算出手段は、前記円筒面からなる基準曲面上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をｘ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面上の法線方向に離れる座標位置をｚ、前記ｚ方向の変位をｗとして、
前記外力が法線方向の力である場合は、次式（１）で表される単位集中荷重作用時の変位解析解により前記影響量を計算する、
【数６】
0005026225.tif 000023
θ：前記座標位置ｘをｙ軸回りの半径ＲO の極座標で表す角度
ただし、（１）式により与えられる変位ｗは位置（θ，ｙ）の要素が位置（θ′，ｙ′）に及ぼす変位、
共形接触下の接触面圧計算装置。
【請求項９】
請求項８において、前記接触面圧計算手段は、前記（１）式で計算されるｚ軸方向の変位ｗを全セルについて重ね合わせる式である次式（２）を、前記釣り合い条件式の一つとして用いて各セルの面圧分布を求め、
【数７】
0005026225.tif 000024
かつ、前記釣り合い条件式の他の一つとして次式（３）を用いて前記接近量αを求める、
【数８】
0005026225.tif 000025
共形接触下の接触面圧計算装置。
【請求項１０】
外力が作用したときに互いに共形接触を行う２つの接触物体の接触面における面圧および表面下応力を計算する装置であって、
請求項７ないし請求項９のいずれか１項に記載の接触面圧の計算装置と、表面下応力計算手段とを含み、
この表面下応力計算手段は、
前記外力が法線方向の力である場合は、次式（４）〜（９）で与えられる単位集中荷重作用時の表面下の応力を用いて、個々のセルの全域にわたり積分して、かつ全セルの影響の和をとって表面下応力を求める、
【数９】
0005026225.tif 000026
共形接触下の接触面圧および表面下応力計算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、接触領域が平坦と見なせない接触物体間の接触面圧および表面下応力の計算を、有限要素法や境界要素法を使用することなく行う計算方法および装置に関する。
【背景技術】
【０００２】
接触面圧の計算問題は下記の二つの分野に分けられる。
一つは非共形接触問題（Counter-formal（またはNon-Conforming）Problem ）であり、接触領域が接触部の物体半径に比べて小さい場合の問題である。
他の一つは、共形接触問題（Conformal （またはConforming）Problem ）であり、接触領域が接触部の物体半径に比べて大きい場合の問題である。
【０００３】
非共形接触問題（Counter-formal Problem）では、接触領域はほぼ平面上にあると見なせるため、比較的解き易い。
特に接触面が二次関数で表される曲面の場合は、ヘルツ（Hertz ）接触理論などが広く用いられており、数表を用いて電卓レベルで計算が可能である。二次関数で表すことができない曲面の場合は、半無限体に集中荷重が作用する場合の解析解（ブジネスク（Boussinesq）の解の重ね合わせにより、数値的に面圧分布を求める手法が知られている。また、これらと類似の解法で各種文献（例えば非特許文献１や非特許文献２）が出ており、この手法をベースとした面圧計算ソフトウェアも市販されている。
【特許文献１】特開昭５９−６９５１９号公報
【非特許文献１】ハーネット（Hartnett, M. J. ）著 ,「転がり要素軸受の接触応力解析」（“The Analysis of Contact Stresses in Rolling Element Bearings”）, ASMEJ. of Lubr., 101, Jan.p105,(1979),
【非特許文献２】栫井著, 「三次元弾性接触問題の境界要素解析」, トライボロジスト, 36,6,(1991), 429-433
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかし、共形接触問題（Conformal Problem ）については、上記のような解析手法は適用できないため、有限要素法や境界要素法を用いて解く以外に計算方法がない。
有限要素法や境界要素法を使用して計算する場合、計算時間が多大になるという課題がある。
【０００５】
このような課題を解消する方法，装置として、本出願人は、次に示す方法，装置を先に提案した（特願２００７−０２６５７７号）。
この提案例は、半無限体に集中荷重が作用したときの変位や表面下応力の解析解を共形接触問題に適用して面圧および表面下応力を計算する手法であり、法線方向荷重作用時はブジネスク（Boussinesq）の解を用いて重ね合わせる。この際に、接触２物体の曲面近傍に基準となる曲面を設置して、その基準曲面をメッシュ分割し、２物体間の距離を計測する点をメッシュ上に設置し、各セルに分布荷重が作用したときの２物体間計測点への変位量の影響を算出する。その位置での接触判定によって接触領域を決め、かつ各セルの面圧の合計が外力と釣合うように２物体間の接近量を決めて、接触面圧を算出する。
【０００６】
上記提案例によると、有限要素法や境界要素法を使用する方法に比べて、比較的簡単に短時間で、かつ精度良く、共形接触における接触面圧の計算、さらに表面下の応力の計算をすることができる。
しかし、上記提案例では、３方向の変位を考慮していたため、まだ計算方法が複雑で、かつ数値積分を行うために計算時間がかかる。そのため、より高速な計算を可能にすることが望まれる。
【０００７】
この発明の目的は、有限要素法や境界要素法を使用することなく、非共形接触問題で行われている数値計算方法を拡張して共形接触問題にも適用できるようにし、かつ上記提案例を改善して、１方向の変位の考慮で計算できて、簡単に短時間で、かつ精度良く、共形接触における接触面圧の計算、さらに表面下の応力の計算をできるようにした方法および装置を提供することを目的とする。
【課題を解決するための手段】
【０００８】
この発明の共形接触下の接触面圧および表面下応力計算方法は、半無限体に集中荷重が作用したときの変位や表面下応力の解析解を共形接触問題に適用して、面圧および表面下応力を計算する手法であって、解析解として、法線方向荷重作用時はブジネスク（Boussinesq）の解を用いる。
【０００９】
具体的には、この発明の共形接触下の接触面圧計算方法は、外力が作用したときに互いに共形接触（接触領域がほぼ平面上にあると見なせない接触）を行う２つの接触物体Ｍ１，Ｍ２の接触面における面圧を、コンピュータを用いて計算する方法であって、
前記接触物体の接触面の形状、弾性係数、ポアソン比、および作用する外力を入力して記憶手段に記憶させる入力過程（Ｓ１）と、入力された値を用いて前記面圧を計算する解析過程（Ｓ２）と、計算した結果を出力する出力過程（Ｓ３）とでなり、
前記解析過程（Ｓ２）として、次の各過程（Ｔ１，Ｔ２，Ｔ４，Ｔ５）を含む。
【００１０】
すなわち、前記２つの接触物体Ｍ１，Ｍ２の互いに接する面の間に接触面形状に近い形状の基準曲面Ｓを設定してこの基準曲面ＳをセルＣが縦横に並ぶメッシュに分割し、接触物体Ｍ１，Ｍ２間の距離を計算する点である計算点Ａをセル上に定めるメッシュ分割過程（Ｔ１）と、
各セルＣに単位分布荷重が作用したときの前記接触物体Ｍ１，Ｍ２の表面の前記各計算点Ａに対応する点に生じる変位である影響量を算出する変位影響算出過程（Ｔ２）と、
各セルＣにつき面圧Ｐと前記単位分布荷重による影響量の積の合計が接触領域内部において接触物体間の接近量と初期隙間の差に釣り合い、かつ各セルＣの面圧Ｐの合計が外力Ｆと釣り合うという釣り合い条件式から、各セルの面圧Ｐを算出する接触面圧計算過程（Ｔ３）とを含む。なお、前記「接近量」は、前記単位分布荷重が与えられることにより互いに近づく方向に移動する距離であり、前記「初期隙間」は、前記単位分布荷重が与えられる前の距離である。
【００１１】
前記変位影響算出過程（Ｔ２）および接触面圧計算過程（Ｔ３）では、計測点間の距離の算出に上記基準曲面上の距離を用いる。また、基準曲面上の位置の違いによる法線方向の修正を行って面圧を算出する。すなわち、影響係数の算出に基準曲面の法線方向の修正を行って面圧を算出する。
【００１２】
なお、上記基準曲面Ｓは、２つの接触物体Ｍ１，Ｍ２の互いに接する面の接触面形状に近い形状であることが好ましい。前記変位影響算出過程（Ｔ２）において、各セルＣに作用する単位分布荷重は、例えば、均一な面圧とするが、１次関数または２次関数に従って変化する面圧を仮定しても良い。接触物体Ｍ１，Ｍ２間の距離は、離散的な計測点Ａで評価するものとするが、計測点Ａは、例えばメッシュの中心にとる。
前記変位影響算出過程（Ｔ２）では、単位集中荷重作用時の変位解析解（法線方向の荷重に対しては例えばブジネスク（Boussinesq）の解を１つの要素全域にわたり積分して、要素１つの面圧分布によって生じる２つの接触物体Ｍ１，Ｍ２の面上の点（上記計測点Ａに対応）の変位量を求める。上記の「要素」は、セルのことである。
【００１３】
この発明方法は、前記２つの物体Ｍ１，Ｍ２が、転がり軸受の内輪とボールとである場合に適用することができる。その場合、前記メッシュ分割過程（Ｔ１）では、前記基準曲面Ｓを円筒面とする。
前記変位影響量算出過程（Ｔ２）では、前記円筒面からなる基準曲面Ｓ上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をｘ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面Ｓ上の法線方向に離れる座標位置をｚ、これら３軸方向の変位を、それぞれｕ，ｖ，ｗとして、
前記外力が法線方向の力である場合は、次式（１）で表される単位集中荷重作用時の変位解析解ｗを、全てのセルについて重ね合わせることで前記影響量を計算する。
【００１４】
【数１０】
0005026225.tif 000002
θ：前記座標位置ｘをｙ軸回りの半径ＲO の極座標で表す角度
ただし、（１）式により与えられる変位ｗは位置（θ，ｙ）の要素が位置（θ′，ｙ′）に及ぼす変位。
【００１５】
すなわち、この発明方法において、前記変位影響量算出過程（Ｔ２）では、単位集中荷重作用時の変位解析解（Boussinesqの解）を１つの要素全域にわたり積分して、要素１つの面圧分布によって生じる２つの接触物体Ｍ１，Ｍ２の表面上の点（上記計測点Ａに対応）の変位量を求める。
位置（θ，ｙ）の要素ｅが位置（θ’，ｙ’）に及ぼす単位面圧負荷時の変位は（１）式で与える。このとき、基準曲面Ｓの法線方向が位置により変わるためCOS(θ−θ’) がかかってくる。
ここで重要なのは、(1) 要素間のθ方向距離として基準曲面上の距離（Ｒ０θ−Ｒ０θ’）を用いる点、(2) 基準曲面Ｓの法線方向の修正を、上記COS(θ−θ’) で行うこと（すなわち、影響係数の算出に基準曲面Ｓの法線方向の修正を行って面圧を算出すること）の２点である。
上記（１）式は、解析的に積分可能であり、数値積分は不要となる。
【００１６】
前記接触面圧計算過程（Ｔ３）は、釣り合い条件式から各セルの面圧Ｐを算出する過程であり、例えば次のように計算する。
すなわち、次の（２）式が釣り合うように各セルの面圧分布を決める。
【００１７】
【数１１】
0005026225.tif 000003
【００１８】
また、次の（３）式が釣り合うように接近量をαを求める。
【００１９】
【数１２】
0005026225.tif 000004
【００２０】
なお、（２）式は、各セルにつき面圧と前記単位分布荷重による影響量の積の合計が接触領域内部において接触物体間の接近量と初期隙間の差に釣り合うという釣り合い条件式であり、また（１）式で計算される変位ｗを全要素について重ね合わせる式である。（３）式は、各セルの面圧の合計が外力と釣り合うという釣り合い条件式である。
【００２１】
この発明方法によると、半無限体に集中荷重が作用したときの変位の解析解を共形接触問題に適用し、また接触面圧計算時の計測点間の距離の算出に上記基準曲面上の距離を用い、かつ影響係数の算出に基準曲面の法線方向の修正を行って面圧を算出するようにしたため、比較的簡単に、したがって短時間で、共形接触における接触面圧の計算を行うことができる。
また、有限要素法を用いて接触面圧を計算した場合と同程度に精度の良い計算結果が得られることが、この発明方法と有限要素法との計算結果を比較することで、確認できた。
この発明は、このように、従来では共形接触には半無限体が適用できないと考えられていた共形接触問題に、ブジネスクの半無限体近似の解を適用し、良好な面圧の解析結果を得ることができる。
なお、この発明方法は、非共形接触の場合にも適用可能であるが、非共形接触の場合はより簡単な計算で面圧を求めることができるため、共形接触の場合に効果的となる。
【００２２】
この発明方法において、表面下応力を計算する場合は、上記のいずれかの方法で各セルＣの面圧Ｐを計算した後に、表面下応力を計算する。
この表面下応力計算過程において、
セルに対する荷重が法線方向の力である場合は、次式（４）〜（９）で与えられる単位集中荷重作用時の表面下の応力を用いて、個々のセルＣの全域にわたり積分して、かつ全セルＣの影響の和をとって表面下応力を求める。
【００２３】
【数１３】
0005026225.tif 000005
【００２４】
セルに対する荷重が法線方向の力である場合は、次のように計算する。前記２つの接触物体Ｍ１，Ｍ２が転がり軸受のボールと内輪とである場合、前記メッシュ分割過程（Ｔ１）では、前記のように前記基準曲面Ｓを円筒面とする。前記変位影響量算出過程（Ｔ２）では、前記円筒面からなる基準曲面Ｓ上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をθ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面上の法線方向に離れる座標位置をｚとして、
位置（θ，ｙ）の要素ｅが位置（θ’，ｙ’，ｚ）に及ぼす単位面圧負荷時の剪断応力τzyを、次式（１０）で求める。
【００２５】
【数１４】
0005026225.tif 000006
【００２６】
このように表面下応力を計算するについても、計算が比較的簡単で短時間で行え、また精度の良い表面下応力の計算が行える。
【００２７】
この発明方法において、前記解析過程の途中で、前記メッシュに分割する大きさを修正してその修正結果により、メッシュ分割後の各過程を再度行うメッシュ修正過程（Ｔ６）を含むようにしても良い。
接触物体の接触面の形状によっては、接触領域の縁がギザギザになるなどして、離散化の誤差が大きくなる場合があるが、解析途中でメッシュを修正することで、離散化の誤差を低減させることができる。
【００２８】
なお、この発明方法は、前記２つの接触物体Ｍ１，Ｍ２が、転がり軸受における内輪または外輪である軌道輪と、ボール等の転動体とである場合に、効果的に適用できる。転がり軸受の軌道輪とボール等の転動体との接触においても、接触領域がほぼ平面上にあると見なせず、共形接触の問題として面圧等を求めることが必要であるため、この発明方法を効果的に適用することができる。
【００２９】
この発明の共形接触下の接触面圧計算装置は、外力が作用したときに互いに共形接触を行う２つの接触物体の接触面における面圧を、計算する装置であって、
前記接触物体Ｍ１，Ｍ２の接触面の形状、弾性係数、ポアソン比、および作用する外力Ｆを入力して記憶手段６に記憶させる入力処理手段４と、入力された値を用いて前記面圧を計算する解析手段７と、計算した結果を出力する出力処理手段５とでなり、
前記解析手段７として、次の各手段８，９，１１，１２を備える。
【００３０】
すなわち、前記２つの接触物体Ｍ１，Ｍ２の互いに接する面の間に、基準曲面Ｓを設定してこの基準曲面ＳをセルＣが縦横に並ぶメッシュに分割し、接触物体Ｍ１，Ｍ２間の距離を計算する点である計算点Ａをセル上に定めるメッシュ分割手段８と、
各セルＣに均一な単位分布荷重が作用したときの前記接触物体Ｍ１，Ｍ２の表面の前記各計算点Ａに対応する点の変位への影響量を算出する変位影響算出手段９と、
各セルＣにつき面圧Ｐと前記単位分布荷重による影響量の積の合計が接触領域内部において接触物体間の接近量と初期隙間の差に釣り合、かつ各セルＣの面圧Ｐの合計が外力と釣り合うという釣り合い条件式から、各セルの面圧Ｐを算出する接触面圧計算手段１０とを含む。
なお、前記「接近量」は、前記単位分布荷重が与えられることにより互いに近づく方向に移動する距離であり、前記「初期隙間」は、前記単位分布荷重が与えられる前の距離である。
前記変位影響算出手段９および接触面圧計算手段１０は、計測点間の距離の算出に上記基準曲面Ｓ上の距離を用い、かつこの基準曲面Ｓ上の位置の違いによる法線方向の修正を行って面圧を算出する。すなわち、影響係数の算出に基準曲面Ｓの法線方向の修正を行って面圧を算出する。
前記メッシュ修正手段１２は、前記接触面圧計算手段１０による前記各セルＣの面圧Ｐの算出よりも後に、前記メッシュ分割手段８で前記メッシュに分割する大きさを修正してその修正結果により、前記変位影響算出手段９および前記接触面圧計算手段１０による算出を再度行わせる。
【００３１】
この発明の表面下応力計算装置において、前記２つの接触物体Ｍ１，Ｍ２が転がり軸受の内輪とボールとである場合に、次の構成としても良い。
前記メッシュ分割手段８は、前記基準曲面Ｓを円筒面とする。
前記変位影響量算出手段９は、前記円筒面からなる基準曲面Ｓ上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をｘ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面上の法線方向に離れる座標位置をｚ、これら３軸方向の変位を、それぞれｕ，ｖ，ｗとして、
前記外力が法線方向の力である場合は、次式（１）で表される単位集中荷重作用時の変位解析解ｗを、全てのセルについて重ね合わせることで前記影響量を計算する。
【００３２】
【数１５】
0005026225.tif 000007
θ：前記座標位置ｘをｙ軸回りの半径ＲO の極座標で表す角度
ただし、（１）式により与えられる変位ｗは位置（θ，ｙ）の要素が位置（θ′，ｙ′）に及ぼす変位。
【００３３】
この場合に、前記接触面圧計算手段１０は、前記（１）式で計算されるｚ軸方向の変位ｗを全セルについて重ね合わせる式である次式（２）を、前記釣り合い条件式の一つとして用いて各セルの面圧分布を求める。
【００３４】
【数１６】
0005026225.tif 000008
【００３５】
また、前記釣り合い条件式の他の一つとして、次式（３）を用いて前記接近量αを求める。
【００３６】
【数１７】
0005026225.tif 000009
【００３７】
この構成の面圧計算装置によると、この発明方法を実施して、共形接触を行う接触面における面圧を、比較的簡単に短時間で計算でき、かつ精度良く計算することができる。
【００３８】
この発明の面圧計算装置において、さらに表面下応力の計算機能を付加する場合は、次の表面下応力計算手段１３を設けることで、表面下応力計算装置を構成する。
前記表面下応力計算手段１３は、セルに対する荷重が法線方向の力である場合は、次式（４）〜（９）で与えられる単位集中荷重作用時の表面下の応力を用いて、個々のセルの全域にわたり積分して、かつ全セルの影響の和をとって表面下応力を求める。
【００３９】
【数１８】
0005026225.tif 000010
【００４０】
この構成の場合、この発明の表面下応力計算方法を実施して、共形接触を行う接触面における表面下応力を、比較的簡単に短時間で計算でき、かつ精度良く計算することができる。
【発明の効果】
【００４１】
この発明の共形接触下の接触面圧計算方法および計算装置によると、半無限体に集中荷重が作用したときの変位の解析解を共形接触問題に適用し、また接触面圧計算時の計測点間の距離の算出に上記基準曲面上の距離を用い、かつ影響係数の算出に基準曲面の法線方向の修正を行って面圧を算出するようにしたため、比較的簡単に、したがって短時間で、共形接触における接触面圧の計算が行え、しかも精度の良い計算結果を得ることができる。特に、解析解として、法線方向荷重作用時はブジネスクの解を重ね合わせるようにした場合は、有限要素法を用いて計算した場合と同程度に精度の良い計算結果を得ることができる。特に、前記解析過程の途中で、前記メッシュに分割する大きさを修正してその修正結果により、メッシュ分割後の各過程を再度行うメッシュ修正過程を含むため、接触物体の接触面の形状によっては、接触領域の縁がギザギザになるなどして、離散化の誤差が大きくなる場合があるが、解析途中でメッシュを修正することで、離散化の誤差を低減させることができる。
この発明の共形接触下の接触面圧，表面下応力の計算方法および計算装置によると、接触面圧に加えて、さらに表面下応力についても、比較的簡単に短時間で、また精度良く、共形接触における表面下応力の計算を行うことができる。特に、メッシュ修正手段を設けたため、解析途中でメッシュを修正することにより、離散化の誤差を低減させることができる。
【発明を実施するための最良の形態】
【００４２】
この発明の一実施形態を図１ないし図８と共に説明する。図１は、この発明方法の各過程の流れ図を、図２はこの発明方法を実施する装置の概念構成のブロック図を、図３は接触面圧計算過程の具体例を、図４は接触モデルを、図５は基準曲面の例をそれぞれ示す。
【００４３】
この共形接触下の接触面圧および表面下応力計算方法は、外力が作用したときに互いに共形接触を行う２つの接触物体Ｍ１，Ｍ２（図４）の接触面における面圧等を、コンピュータを用いて計算する方法である。２つの接触物体Ｍ１，Ｍ２は、例えば、転がり軸受におけるボールと内輪とである。
この計算方法は、前記接触物体Ｍ１，Ｍ２の接触面の形状、弾性係数（ここでは、ヤング率Ｅと、ポアソン比ν）、および作用する外力Ｆを入力して記憶手段に記憶させる入力過程（Ｓ１）と、入力された値を用いて前記面圧を計算する解析過程（Ｓ２）と、計算した結果を出力する出力過程（Ｓ３）とでなる。解析過程（Ｓ２）は、順次行われる次のメッシュ分割過程（Ｔ１）、変位影響算出過程（Ｔ２）、接触面圧計算過程（Ｔ３）、および表面下応力計算過程（Ｔ４）を含む。
【００４４】
計算装置となるコンピュータ１（図２）は、中央処理装置（ＣＰＵ）およびメモリ等の記憶手段（いずれも図示せず）を有し、このコンピュータ１のハードウェア（オペレーションシステムを含む）、およびこのコンピュータ１に実行させる接触面圧・表面下応力計算プログラム（図示せず）により、図２に概念構成で示す各手段が構成される。
入力処理手段４は、入力装置２から、接触物体Ｍ１，Ｍ２の接触面の形状、弾性係数、ポアソン比、および作用する外力Ｆを入力して入力情報記憶手段６に記憶させる手段である。ここで言う「形状」は、大きさを含む概念である。例えば、前記接触物体Ｍ１，Ｍ２の接触面の形状のデータとしては、ボールからなる接触物体Ｍ１の半径Ｒ1 と、内輪からなる接触物体Ｍ２の軌道面の溝半径Ｒ2 軸受のPCD 、接触角とを入力する。入力処理手段４は、上記の形状等の情報の他に、後述の解析手段７でメッシュ分割する枡目の各軸方向の大きさやセル個数（ＮＡ，ＮＢ）等を、入力装置２から入力して入力情報記憶手段６に記憶させる機能を備えるものとしても良い。また、入力処理手段４は、入力装置２から入力された情報に対して、解析手段７で演算するための前処理を行う機能を有するものとしても良い。入力装置２は、キーボード等であっても、また通信手段や記憶素子等であっても良い。入力情報記憶手段６は、コンピュータ１の上記記憶手段における所定の記憶領域である。出力処理手段５は、解析手段７で計算した結果を、液晶ディスプレイ，プリンタ，あるいは通信機器等からなる出力装置３に出力する手段である。
【００４５】
解析手段７は、入力処理手段４で入力されて入力情報記憶手段６に記憶されたデータを用いて前記面圧および表面下応力を計算する手段であり、メッシュ分割手段８、変位影響算出手段９、接触面圧計算手段１０、および表面下応力計算手段１１を備える。解析手段７に備えられた上記各手段８〜１２は、それぞれ図１（Ｂ）の流れ図におけるメッシュ分割過程（Ｔ１）、変位影響算出過程（Ｔ２）、接触面圧計算過程（Ｔ３）、および表面下応力計算過程（Ｔ４）を実行する手段であり、その実行に必要な計算式や設定値等を記憶している。上記各手段８〜１１の具体的な内容は、次に述べる上記各過程（Ｔ１）〜（Ｔ４）の説明の通りである。
解析手段７は、上記の他に、後に説明するメッシュ修正手段１２を有している。
【００４６】
図１（Ａ）の解析過程（Ｓ２）を、同図（Ｂ）および図４〜図６と共に説明する。
この解析仮定（Ｓ２）は、法線方向加重が作用したときのハーネット（Hartnett) の手法を拡張して共形接触問題（Conformal Problem ）を解く手法である。
【００４７】
まず、メッシュ分割過程（Ｔ１）として、図４のように、接触物体Ｍ１，Ｍ２の間に基準曲面Ｓを設定する。この基準曲面Ｓを、図５，図６のようにセルＣが縦横に並ぶメッシュに分割し、接触物体Ｍ１，Ｍ２間の距離を計算する点である計算点ＡをセルＣ上に定める。
また、上記基準曲面Ｓは、２つの接触物体Ｍ１，Ｍ２の接触面形状に近い形状の曲面とする。
接触物体Ｍ１，Ｍ２は、例えば図１１に示す転がり軸受におけるボール２３と内輪２１とである。同図の転がり軸受は、内輪２１と外輪２２の断面円弧面状の軌道面間に複数のボール２３を介在させ、ボール２３を保持する保持器２４を設けたアンギュラ玉軸受、または深溝玉軸受等の玉軸受である。
接触物体Ｍ１，Ｍ２が、上記のように転がり軸受におけるボールと内輪とである場合、基準曲面Ｓは円筒面とみなす。同図では、基準曲面Ｓは、内輪である接触物体Ｍ２の軌道面とボールである接触物体Ｍ１の外球面との間に介在する円筒面であって、その円筒中心軸は内輪である接触物体Ｍ２の軸心と平行としている。
【００４８】
なお、基準曲面Ｓについて、接触面形状の正確な形状は本来わからないが、ボールと内輪の接触では接触面形状は円筒面の一部に近いとみなせる（解析上もこの基準面を採用している）。また、凸球面と凹球面の接触では接触面形状はほぼ球面とみなせる。従来の「接触面が平坦と見なせる」場合は平面である。複雑な曲面同士の接触では上記のような簡単な曲面で表現できない場合が出てくる。基準面形状をどのようにするかは解析精度に影響してくるため、できるだけ接触面形状に近い方が良いが、複雑にすると計算が難しくなるというデメリットがある。従って「基準曲面Ｓ」は解析する人がケースバイケースで適宜決定していく。
【００４９】
計算点Ａは、接触物体Ｍ１，Ｍ２間の距離を離散的な点で評価するものとするためにセルＣの位置を代表させる点である。計算点Ａは、セルＣ内の適宜の位置に設定すれば良いが、ここではセルＣの中心にとる。
メッシュ分割する枡目の大きさ、つまりセルＣの大きさは、適宜設定すればよいが、入力処理手段４からの入力等によって自由に設定可能としても良い。また、基準曲面Ｓの広さは、接触物体Ｍ１，Ｍ２の接触領域を全て含むように、接触領域より広めに適宜設定する。ここでは、セルＣがＸ軸方向にＮＡ個、Ｙ軸方向にＮＢ個並ぶ基準曲面Ｓとする。
【００５０】
変位影響算出過程（Ｔ２）では、メッシュ分割過程（Ｔ１）の後、各セルＣに均一な分布荷重Ｐが作用したときの接触物体Ｍ１，Ｍ２の表面の前記各計算点Ａに対応する点の変位への影響量を算出する。この場合に、ここでは、個々のセルＣには全体に均一な面圧Ｐが作用すると仮定するが、１次関数または２次関数等の所定の関数に従いセルＣ内の位置によって変わる面圧が作用すると仮定しても良い。
【００５１】
このとき、互いに直交する３軸方向の座標位置をｘ，ｙ，ｚ、これら３軸方向の変位を、それぞれｕ，ｖ，ｗとして計算する。なお、ｘ，ｙ軸方向の変位ｕ，ｖは、後に、単位集中荷重作用時の変位解析解（Boussinesqの解）を説明するために記載したが、この発明では、これらｘ，ｙ軸方向の変位ｕ，ｖを用いず、ｚ軸方向の変位ｗのみを用いて計算する。
２つの接触物体Ｍ１，Ｍ２が、上記のように転がり軸受のボールと内輪であって、基準面Ｓを円筒面とする場合、この円筒面からなる基準曲面Ｓ上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をｘ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面上の法線方向に離れる座標位置をｚ、これら３軸方向の変位を、それぞれｕ，ｖ，ｗとする。なお、基準面Ｓ上のｘ座標位置は、角度θにより表現される極座標Ｒo θで表される。
【００５２】
前記変位影響算出過程（Ｔ２）では、単位集中荷重作用時の変位解析解（Boussinesqの解）を１つの要素全域にわたり積分して、要素１つの面圧分布によって生じる２接触物体Ｍ１，Ｍ２の表面上の点（上記計測点Ａに対応）の変位量を求める。
すなわち、位置（θ，ｙ）の要素ｅが位置（θ’，ｙ’）に及ぼす単位面圧負荷時の変位を、次の（１）式で与える。このとき、基準曲面の法線方向が位置により変わるためCOS(θ−θ’) がかかってくる。なお、ここでは、一つのセルＣを一つの要素として計算する。
ここで重要なのは、(1) 要素間のθ方向距離として基準面上の距離（Ｒ０θ−Ｒ０θ’）を用いる点、および(2) 基準面の法線方向の修正を、COS(θ−θ’) の乗算により行うこと（すなわち、影響係数の算出に基準曲面Ｓの法線方向の修正を行って面圧を算出すること）の２点である。
【数１９】
0005026225.tif 000011
【００５３】
なお、前記の式（１）は、単位集中荷重作用時の変位解析解（Boussinesqの解）から、後述のように導かれる。
【００５４】
前記接触面圧計算過程（Ｔ３）では、各セルＣにつき面圧と前記単位分布荷重による影響量の積の合計が接触領域内部において接触物体Ｍ１，Ｍ２間の接近量と初期隙間の差に釣り合い、かつ各セルＣの面圧の合計が外力と釣り合ういう釣り合い条件式から、各セルの面圧を算出する。
【００５５】
接触面圧計算過程（Ｔ３）は、具体的には、次の（２）式が釣り合うように面圧分布を決め、また（３）式が釣り合うように、接近量αを求める。（２）式は、各セルＣijにつき面圧Ｐklと前記単位分布荷重による影響量ｆijklの積の合計が接触領域内部において接触物体Ｍ１，Ｍ２間の接近量と初期隙間の差Ｄijに釣り合うという条件式であり、（１）式で求められる変位ｗを全要素について重ね合わせる式である。（３）式は各セルＣijの面圧Ｐijの合計が外力Ｆと釣り合ういう条件式である。
なお、Ｃijは、ｘ軸方向にｉ番目でｙ軸方向にｊ番目にあるセルＣを示す。Ｑklは、ｘ軸方向にｋ番目でｙ軸方向にｌ番目の点Ａを示す（なお、ｉ：１〜ＮＡ、ｊ：１〜ＮＢ、ｋ：１〜ＮＡ、ｌ：１〜ＮＢ）。
【００５６】
【数２０】
0005026225.tif 000012
【００５７】
【数２１】
0005026225.tif 000013
【００５８】
上記の（２）式および（３）式が釣り合うように各セルＣの面圧値を決めれば、離散化された全体の面圧分布が求まる。
なお、上記のように（２）式を解く過程で、２つの接触物体Ｍ１，Ｍ２間の接近量と初期隙間の差Ｄijが正の値となるセルＣijについては、そのセルＣijの面圧は零としており、この接近量と初期隙間の差Ｄijが零または負の値となるセルＣijの並び範囲が、接触領域である。
（２）式は、解析的に積分可能であり、数値積分は不要となる。
【００５９】
図３に、（２）式、（３）式を、コンピュータにおける接触面圧計算手段１０によって解くアルゴリズムの一例を示す。
まず、α値（：２つの接触物体Ｍ１，Ｍ２の接近量）の初期値を設定する（ステップＨ１）。
そのα値に応じて各セルＣijの接触面圧を設定する（ステップＨ２）。この場合に、ｉ＝１〜ＮＡの範囲、およびｊ＝１〜ＮＢの範囲で全てのセルＣijの接触面圧Ｐijを設定する。その設定した各セルＣijの接触面圧Ｐijで（２）式が成り立つか否かを判定し、成り立たない場合は、セルＣijの接触面圧の設定過程（ステップＨ２）に戻る（Ｈ３）。ステップＨ２では、各セルＣijの接触面圧Ｐijを再度設定し直し、その設定したセルＣijの接触面圧Ｐijで（２）式が成り立つか否かを判定し、成り立たない場合はセルＣijの接触面圧Ｐijを設定過程（Ｈ２）に戻る。このような処理を（２）式が成り立つまで繰替し、（２）式が成り立つと、ステップＨ４に進む。
ステップＨ４では、（３）式が成り立つか否かを判定し、成り立たない場合は、ステップＨ１に戻ってα値を、再度設定し直す。
この後、前記と同様にしてステップＨ２〜Ｈ４，Ｈ１の処理を繰り返し、（３）式が成り立つと、処理を終了する。
このようにして求めた、式（２）、（３）が成り立つときの、ステップＨ２で設定した各セルＣijの接触面圧Ｐijが、接触面圧計算手段１０の出力となる接触面圧の値である。  なお、セルＣijの接触面圧の最初の値は何でもよいが、できるだけ実際の面圧に近い値を設定した方が早く収束する。
【００６０】
この接触面圧計算方法によると、このように、半無限体に集中荷重が作用したときの変位の解析解を共形接触問題に適用し、また接触面圧計算時の計測点間の距離の算出に上記基準曲面Ｃ上の距離を用い、かつ影響係数の算出に基準曲面Ｃの法線方向の修正を行って面圧を算出するようにしたため、比較的簡単に、したがって短時間で、共形接触における接触面圧の計算を行うことができる。
また、有限要素法を用いて接触面圧を計算した場合と同程度に精度の良い計算結果が得られることが、この実施形態の方法と有限要素法との計算結果を比較することで、確認できた。
この実施形態は、このように、従来では共形接触には半無限体が適用できないと考えられていた共形接触問題に、ブジネスクの半無限体近似の解を適用し、良好な面圧の解析結果を得ることができる。
なお、この発明方法は、非共形接触の場合にも適用可能であるが、非共形接触の場合はより簡単な計算で面圧を求めることができるため、共形接触の場合に効果的となる。
【００６１】
表面下応力計算過程（Ｔ４）では、上記のように求めた面圧分布に対して、次式（４）〜（９）で与えられる単位集中荷重作用時の表面下の応力を用いて、個々のセルＣの全域にわたり積分して、かつ全セルＣの影響の和をとって表面下応力が求まる。
【００６２】
【数２２】
0005026225.tif 000014
【００６３】
セルＣに対する荷重が法線方向の力である場合は、次のように計算する。前記２つの接触物体Ｍ１，Ｍ２が転がり軸受のボールと内輪とである場合、前記メッシュ分割過程（Ｔ１）では、前記のように前記基準曲面Ｓを円筒面とする。前記変位影響量算出過程（Ｔ２）では、前記円筒面からなる基準曲面Ｓ上の、円筒中心軸に垂直な平面で切断した断面の円弧に沿う方向の座標位置をθ、円筒中心軸に沿う方向の座標位置をｙ、基準曲面上の法線方向に離れる座標位置をｚとして、
位置（θ，ｙ）の要素ｅが位置（θ’，ｙ’，ｚ）に及ぼす単位面圧負荷時の剪断応力τzyを、次式（１０）で求める。
【００６４】
【数２３】
0005026225.tif 000015
【００６５】
この場合も、積分は解析的に求まるため、数値積分は不要である。
【００６６】
このように表面下応力を計算するについても、計算が比較的簡単で短時間で行え、また精度の良い表面下応力の計算が行える。
【００６７】
なお、上記実施形態において、接触物体Ｍ１，Ｍ２の形状によっては、接触領域の縁がギザギザにあることがあるが、そのような場合は、解析途中で、例えば接触面圧計算過程（Ｔ３）の後、メッシュ分割の修正、つまりセルＣの各軸ｘ，ｙ方向の大きさを修正してその修正結果により、メッシュ分割後の各過程を再度行わせるメッシュ修正過程（Ｔ５）を含むようにしても良い。
図２のメッシュ修正手段１２は、このようなメッシュ修正過程（Ｔ５）を実行する手段である。
【００６８】
接触物体Ｍ１，Ｍ２が、転がり軸受におけるボールと内輪である場合における、この実施形態、前述の提案例（特願２００７−０２６５７７号）、および有限要素法（ＦＥＭ）の計算結果を、図９に示す。３者とも良く一致しており、どれも使用可能レベルである。一方、従来手法（接触面が平坦という前提の解析）では、図１０のようであり、正しくない結果となる。
【００６９】
前述の提案例の手法では、３方向の変位を考慮していたため、計算時間が長くかかっていたが（ＦＥＭよりは早い）が、この実施形態は、これを１方向とすることができた。また、全順の提案例は数値積分で行っていたが、この実施形態では解析積分に置き換えることができ、計算の簡素化と高速化が達成できた。
計算の所要時間の例を示すと、前述の提案例の手法では１９７．０８秒であったのに対して、実施形態の手法では１１２．４５秒となり、従来手法の５７％であった。
【００７０】
つぎに、上記の式（１）とブジネスク（Boussinesq）の解との関係を説明する。
【００７１】
作用荷重が法線方向である場合の単位集中荷重作用時の変位解析解であるブジネスク（Boussinesq）の解は、次式（１１）〜（１３）によって現される。
なお、図７はブジネスクの解の座標系を示す。同図のような半無限体平面の原点Ｏに単位荷重が作用するとき、任意の点であるＱ点のｘ，ｙ，ｚ方向の変位を、それぞれｕ，ｖ，ｗとする。
【００７２】
【数２４】
0005026225.tif 000016
【００７３】
この式から、次のようにして、上記の式（１）が導かれる。
【数２５】
0005026225.tif 000017
【００７４】
これを要素eが点Qに及ぼす変位として相対距離を用いて積分したものが(1)式である。本手法では(1)式のｘには基準面上の距離を用いるので、ｘの代わりにR0θを用いる。なお、積分に当たっては以下の変数変換
ｘ＝Ｒ0θ  ⇒  ｄｘ＝Ｒ0ｄθ  
を用いて（１）式が得られる。
【図面の簡単な説明】
【００７５】
【図１】（Ａ）はこの発明の一実施形態に係る共形接触下の接触面圧および表面下応力計算方法を示す流れ図、（Ｂ）はその解析過程の流れ図である。
【図２】この発明の一実施形態に係る共形接触下の接触面圧および表面下応力計算装置の概念構成を示すブロック図である。
【図３】接触面圧計算過程の具体例の流れ図である。
【図４】接触モデルの一例の説明図である。
【図５】基準曲面およびメッシュ分割の形態例を示す説明図である。
【図６】基準曲面の展開図である。
【図７】ブジネスクの解の座標系の説明図である。
【図８】接触物体間の初期隙間の説明図である。
【図９】接触物体の重荷重時の面圧分布を、この実施形態の方法、提案例の方法、およびＦＥＭ解析で求めた結果を示すグラフである。
【図１０】同接触物体の重荷重時の面圧分布を、市販の軽荷重時用のソフトウェアで計算した結果の説明図である。
【図１１】同面圧を求める対象となる転がり軸受の部分断面図である。
【符号の説明】
【００７６】
１…コンピュータ
４…入力処理手段
５…出力処理手段
６…入力情報記憶手段
７…解析手段
８…メッシュ分割手段
９…変位影響算出手段
１０…接触面圧計算手段
１１…表面下応力計算手段
Ｍ１，Ｍ２…接触物体
Ｓ…基準曲面
Ｃ…セル
【図１】
0005026225.tif 000027
【図２】
0005026225.tif 000028
【図３】
0005026225.tif 000029
【図４】
0005026225.tif 000030
【図５】
0005026225.tif 000031
【図６】
0005026225.tif 000032
【図７】
0005026225.tif 000033
【図８】
0005026225.tif 000034
【図９】
0005026225.tif 000035
【図１０】
0005026225.tif 000036
【図１１】
0005026225.tif 000037