(57)【要約】  Ｎ個の部材のＭ組の接触面の各節点に対して、予め各部材の部分剛性マトリックスを作成しておき（Ｓ１００，Ｓ１１０）、各節点の接触状態を仮定し（Ｓ１４０）、この仮定の下で、接触パラメータ（分割比，法線方向，摩擦方向）を用いて接触位置マトリックスと接触力マトリックスとを作成する処理（Ｓ１８０）、接触位置マトリックスを用いた方程式と接触力マトリックスと各部材の部分剛性マトリックスを用いた方程式とを連立させて求解する処理（Ｓ１９０）、求解により得られた解において、各節点のうち接触状態の矛盾が生じる節点の一部について矛盾を解消するように接触状態を変更する処理（Ｓ２２０，Ｓ２３０）、を全節点の接触状態に矛盾が生じないようになるまで繰り返す。
2017130811.tif 000002
【特許請求の範囲】
【請求項１】
Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点が非接触の状態の非接触節点であるかすべりが生じているすべり節点であるか固着している固着節点であるかのうちのいずれの接触状態であるかを解析する接触状態解析方法であって、
（１）前記Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておき、
（２）前記Ｎ個の部材の一部に対して力または変位を作用させたときに、
（ａ）前記接触面の一方の要素を構成する節点のうち前記接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか前記非接触節点であるかを仮定すると共に前記コンタクト節点についてはすべり節点であるか固着節点であるかを仮定し、前記コンタクト節点に対しては前記接触面の他方の要素のうち接触しているターゲット要素における前記コンタクト節点の接触位置に関する情報としての分割比と前記コンタクト節点における前記ターゲット要素に対する法線方向とを計算し、前記すべり節点については前記コンタクト節点での前記ターゲット要素における摩擦方向を計算する接触状態仮定ステップと、
（ｂ）前記コンタクト節点の接触状態，分割比，法線方向，摩擦方向を用いて、前記コンタクト節点の変位の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、前記コンタクト節点の力の釣り合いに基づいてコンタクト節点の接触力に関連するベクトルに乗じたときに各節点に接触により生じる力ベクトルとなる接触力マトリックスを作成するマトリックス作成ステップと、
（ｃ）前記接触位置マトリックスを用いた方程式と前記接触力マトリックスと前記線形データから作成される各部材の部分剛性マトリックスを用いた方程式とを連立させて未知数を求解する求解ステップと、
（ｄ）前記求解ステップにおける解において、各節点のうち接触状態の矛盾が生じる節点の少なくとも一部について矛盾を解消するように接触状態を変更する変更ステップと、
を前記求解ステップで求められた全ての節点の接触状態に矛盾が生じないようになるまで前記接触状態仮定ステップと前記マトリックス作成ステップと前記求解ステップと前記変更ステップとを繰り返すことにより、各節点の接触状態を特定する、
ことを特徴とする接触状態解析方法。
【請求項２】
請求項１記載の接触状態解析方法であって、
前記接触位置マトリックスは、前記コンタクト節点が前記固着節点であるときには該コンタクト節点の変位後の位置と前記ターゲット要素における前記コンタクト節点の位置の変位後の位置とが同一であること及び前記コンタクト節点が前記すべり節点であるときには前記コンタクト節点の変位後の位置が前記ターゲット要素の変位後の接触面内であることに基づくマトリックスであり、
前記接触力マトリックスは、前記コンタクト節点の接触力が前記ターゲット要素における前記コンタクト節点の位置の接触力と同じ大きさで反対向きの力であって、前記コンタクト節点が前記すべり節点であるときには前記コンタクト節点の接触力が前記ターゲット要素に対する法線方向の力と前記ターゲット要素の面内方向の摩擦方向の摩擦力との合力であることに基づくマトリックスである、
接触状態解析方法。
【請求項３】
請求項１または２記載の接触状態解析方法であって、
前記接触状態の矛盾は、前記固着節点において接触力の接触面内方向成分と法線方向成分の比が静止摩擦係数より大きいという第１の矛盾、前記すべり節点において仮定した摩擦力方向と解により得られたすべり方向とが逆であるという第２の矛盾、前記コンタクト節点において垂直抗力が負となる第３の矛盾、前記非接触節点において変位が前記接触面の他方側に食い込むという第４の矛盾、のいずれかである、
接触状態解析方法。
【請求項４】
請求項３記載の接触状態解析方法であって、
前記変更ステップにおける接触状態の変更は、前記第１の矛盾が生じている節点に対しては固着節点からすべり節点への変更であり、前記第２の矛盾が生じている節点に対してはすべり節点から固着節点への変更であり、前記第３の矛盾が生じている節点に対してはコンタクト節点から非接触節点への変更であり、前記第４の矛盾が生じている節点に対して非接触節点からコンタクト節点への変更である、
接触状態解析方法。
【請求項５】
請求項１ないし４のうちのいずれか１つの請求項に記載の接触状態解析方法であって、
前記接触状態仮定ステップにより仮定した前記コンタクト節点の分割比，法線方向，摩擦方向と前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向とのいずれかに所定以上の誤差が生じたときには前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向を用いて前記マトリックス作成ステップを実行する、
接触状態解析方法。
【請求項６】
請求項１ないし５のうちのいずれか１つの請求項に記載の接触状態解析方法であって、
前記線形データは、有限要素法により解析されたデータである、
接触状態解析方法。
【請求項７】
請求項１ないし６のうちのいずれか１つの請求項に記載の接触状態解析方法であって、
前記Ｎ個の部材の一部に時間変化する力または変位を作用させる際に、力または変位の作用開始から微小時間経過毎に、作用させた力または変位に対して前記（２）のステップにより特定された各節点の接触状態を記憶する、
接触状態解析方法。
【請求項８】
請求項１ないし７のうちのいずれか１つの請求項に記載の接触状態解析方法を用いて接触状態解析を行なう接触状態解析装置。
【請求項９】
コンピュータを請求項８の接触状態解析装置として機能させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本開示は、接触状態解析方法および接触状態解析装置並びにプログラムに関する。
【背景技術】
【０００２】
従来、この種の接触状態解析方法としては、有限要素法を用いてボルト接合を有する平板構造物の減衰特性を解析するものが提案されている（例えば、非特許文献１参照）。この論文では、有限要素法を接触状態の解析に用いている。まず、線形振動解析を行ない、振動モード形状を基にモード毎の減衰特性の違いを考察する。そして、振動モード形状を強制変位部分とした静的接触解析を行ない、摩擦による散逸エネルギーを算定することによりモード毎の減衰特性を評価している。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】「ボルト接合を有する平板構造物の減衰特性」，平居嵩朗，鞍谷文保，城戸一郎，日本機械学会［Ｎｏ．１２−１２］  Ｄｙｎａｍｉｃｓ  ａｎｄ  Ｄｅｓｉｇｎ  Ｃｏｎｆｅｒｅｎｃｅ  ２０１２  ＵＳＢ論文集［２０１２，９，１８−２１，横浜］
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上述の論文の技術では、大規模な各部材の線形変形問題の有限要素法解析が繰り返し処理に組み込まれているため、解析結果を得るためには長時間を要してしまう。有限要素法解析は線形変形問題の解析に有効であるが、複数部材の求解を同時に行なう場合、計算の規模が大きく処理が多いため、繰り返し処理には不向きである。こうした課題に対して、すべり率を導入し、これを適当に定めることにより、理論値と実験値とを合致させるものも提案されているが、接触面におけるどの部分が固着しており、どの部分がすべっているのかを特定することやすべり量を予測することもできず、本質的な解法と判断することはできない。【０００５】
また、出願人は、２部材の接触面の各部の接触状態を解析する接触状態解析方法については特願２０１４−２２２９８２により提案しているが、３部材以上の各接触面の各部の接触状態の解析については未解決であった。【０００６】
本開示の接触状態解析方法および接触状態解析装置並びにプログラムは、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することを主目的とする。
【課題を解決するための手段】
【０００７】
本開示は、上述の主目的を達成するために以下の手段を採った。【０００８】
本開示の接触状態解析方法は、
Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点が非接触の状態の非接触節点であるかすべりが生じているすべり節点であるか固着している固着節点であるかのうちのいずれの接触状態であるかを解析する接触状態解析方法であって、
（１）前記Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておき、
（２）前記Ｎ個の部材の一部に対して力または変位を作用させたときに、
（ａ）前記接触面の一方の要素を構成する節点のうち前記接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか前記非接触節点であるかを仮定すると共に前記コンタクト節点についてはすべり節点であるか固着節点であるかを仮定し、前記コンタクト節点に対しては前記接触面の他方の要素のうち接触しているターゲット要素における前記コンタクト節点の接触位置に関する情報としての分割比と前記コンタクト節点における前記ターゲット要素に対する法線方向とを計算し、前記すべり節点については前記コンタクト節点での前記ターゲット要素における摩擦方向を計算する接触状態仮定ステップと、
（ｂ）前記コンタクト節点の接触状態，分割比，法線方向，摩擦方向を用いて、前記コンタクト節点の変位の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、前記コンタクト節点の力の釣り合いに基づいてコンタクト節点の接触力に関連するベクトルに乗じたときに各節点に接触により生じる力ベクトルとなる接触力マトリックスを作成するマトリックス作成ステップと、
（ｃ）前記接触位置マトリックスを用いた方程式と前記接触力マトリックスと前記線形データから作成される各部材の部分剛性マトリックスを用いた方程式とを連立させて未知数を求解する求解ステップと、
（ｄ）前記求解ステップにおける解において、各節点のうち接触状態の矛盾が生じる節点の少なくとも一部について矛盾を解消するように接触状態を変更する変更ステップと、
を前記求解ステップで求められた全ての節点の接触状態に矛盾が生じないようになるまで前記接触状態仮定ステップと前記マトリックス作成ステップと前記求解ステップと前記変更ステップとを繰り返すことにより、各節点の接触状態を特定する、
ことを特徴とする。
【図面の簡単な説明】
【０００９】
【図１】本開示の第１実施例の接触状態解析装置２０の構成の概略を示す構成図である。
【図２】接触状態解析プログラムの一例を示すフローチャートである。
【図３】Ｎ個の部材のＭ組の接触面の一例を示す説明図である。
【図４】接触ペアを設定する様子を説明する説明図である。
【図５】第２実施例の接触状態解析装置１２０の構成の概略を示す構成図である。
【図６】第２実施例の接触状態解析プログラムの一例を示すフローチャートである。
【図７】第３実施例の接触状態解析装置２２０の構成の概略を示す構成図である。
【図８】第３実施例の接触状態解析プログラムの一例を示すフローチャートである。
【発明を実施するための形態】
【００１０】
次に、本開示を実施するための形態を実施例を用いて説明する。図１は、本開示の第１実施例の接触状態解析装置２０の構成の概略を示す構成図である。第１実施例の接触状態解析装置２０は、図示するように、一般的な汎用コンピュータ２２にアプリケーションソフトウエアとしての接触状態解析プログラム３０がインストールされたものとして構成されている。接触状態解析プログラム３０は、部分剛性マトリックス構築モジュール３２、接触状態仮定モジュール３４、マトリックス作成モジュール３６、求解モジュール３８、変更モジュール４０、判定モジュール４２、を含んでいる。
【００１１】
部分剛性マトリックス構築モジュール３２は、Ｍ組の接触面により連接しているＮ個の部材の任意の部位に作用する線形変形の範囲内の力（トルクを含む）と変位（ねじり変位を含む）との関係を有限要素法や差分法あるいは境界要素法などを用いて解析し、得られた解を線形データのデータベース４４として記憶し、このデータベース４４から部分剛性マトリックスを作成するモジュールである。ここで、実施例では、「接触面」には、面の一部が接触している面が含まれるだけでなく、面の一部が接触している可能性がある面も含まれる。
【００１２】
接触状態仮定モジュール３４は、Ｍ組の接触面の各要素を構成する節点が接触しているコンタクト節点であるか非接触の非接触節点であるかを仮定すると共にコンタクト節点についてはすべりが生じているすべり節点であるか固着している固着節点であるかを仮定し、コンタクト節点については、コンタクト節点と接触している要素（ターゲット要素）におけるコンタクト節点の接触位置をターゲット要素を構成する節点の内分比により示した分割比や、コンタクト節点におけるターゲット要素に対する法線方向を計算し、更に、すべり節点についてはターゲット要素に対する摩擦方向を計算するモジュールである。本明細書では、分割比と法線方向と摩擦方向を接触パラメータと称する。
【００１３】
マトリックス作成モジュール３６は、コンタクト節点の接触状態（すべり節点であるか固着節点であるかの状態）と接触パラメータ（分割比，法線方向，摩擦方向）を用いて、コンタクト節点の変位の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、コンタクト節点の力の釣り合いに基づいてコンタクト節点の接触力情報ベクトル（固着時の接触力の３方向成分とすべり時の垂直抗力とからなるベクトル）に乗じたときにコンタクト節点に接触により生じる力ベクトルとなる接触力マトリックスを作成するモジュールである。具体的には、接触位置マトリックスは、コンタクト節点が固着節点であるときにはコンタクト節点の変位後の位置とターゲット要素におけるコンタクト節点に対応する点の変位後の位置とが同一であること及びコンタクト節点がすべり節点であるときにはコンタクト節点の変位後の位置がターゲット要素の変位後の接触面内であることに基づいて作成され、接触力マトリックスは、コンタクト節点の接触力がターゲット要素におけるコンタクト節点の位置の接触力と同じ大きさで反対向きの力であって、コンタクト節点がすべり節点であるときにはコンタクト節点に生じる接触力がターゲット要素に対する法線方向の力とターゲット要素の面内方向の摩擦方向の摩擦力の合力であることに基づいて作成される。
【００１４】
求解モジュール３８は、マトリックス作成モジュールにより作成した接触位置マトリックスを用いた方程式と接触力マトリックスと部分剛性マトリックス構築モジュール３２により作成された部分剛性マトリックスを用いた方程式とを連立させて求解するモジュールである。
【００１５】
変更モジュール４０は、求解モジュール３８により得られる解において、各節点のうち接触状態の矛盾が生じる節点の少なくとも一部について矛盾を解消するように接触状態を変更するモジュールである。接触状態の矛盾としては、固着節点において接触力の接触面内方向成分と法線方向成分の比が静止摩擦係数より大きいという第１の矛盾や、すべり節点において仮定した摩擦力方向と解により得られたすべり方向とが略逆であるという第２の矛盾、コンタクト節点において垂直抗力が負となる第３の矛盾、非接触節点において変位が前記接触面の他方側に食い込むという第４の矛盾が想定されている。そして、接触状態の変更としては、第１の矛盾が生じている節点に対しては固着節点からすべり節点への変更や、第２の矛盾が生じている節点に対してはすべり節点から固着節点への変更、第３の矛盾が生じている節点に対してはコンタクト節点から非接触節点への変更、第４の矛盾が生じている節点に対して非接触節点からコンタクト節点への変更となる。
【００１６】
判定モジュール４２は、変更モジュール４０により全ての節点の接触状態に矛盾が生じないようになるまで接触状態仮定モジュール３４の処理やマトリックス作成モジュール３６の処理、求解モジュール３８の処理、変更モジュール４０の処理を繰り返すように各モジュールを起動し、求解モジュール３８で求められた全ての節点の接触状態に矛盾が生じなくなったときに、各節点の変位や接触状態、接触パラメータなどを解データとして記憶装置４６に保存するモジュールである。
【００１７】
上記の各モジュールにより構成された接触状態解析プログラムの一例を示すフローチャートを図２に示す。以下に、図２のフローチャートを用いてＮ個の部材のＭ組の接触面における各要素を構成する節点の接触状態を解析する手法を詳細に説明する。接触状態解析プログラムが実行されると、まず、図３に例示するＮ個の部材がＭ組の接触面により接触している解析対象において、Ｍ組の接触面を構成する各部材に対して有限要素法解析を用いて線形データを求解すると共に線形データをデータベース４４に記憶する（ステップＳ１００）。具体的には、解析対象の各部材の任意の部位に線形変形の範囲内となる所定の力（例えば単位力や単位トルク）や所定の変位（例えば、単位変位や単位ねじり変位）を作用させたときの節点にどのような力が作用しどのような変位が生じるかについての関係を所定の力や所定の変位の作用部位を変化させて各々求解し、所定の力や所定の変位、作用部位、各節点に作用する力、各節点の変位などを線形データとしてデータベース４４に記憶する。
【００１８】
こうして線形データをデータベース４４に記憶すると、記憶した線形データを用いて、Ｎ個の部材に対して節点に対する任意の作用力に対する各要素の変位や力、節点に対する任意の作用変位に対する各要素に作用する力や変位を求める部分剛性マトリックスを作成する（ステップＳ１１０）。
【００１９】
次に、Ｍ組の接触面の摩擦係数を入力し（ステップＳ１２０）、時間の初期化（ｔ＝０）を行なう（ステップＳ１３０）。ここで、摩擦係数は、静止摩擦係数と動摩擦係数とに区分して入力してもよいし、速度や荷重の関数としてもよい。また、簡単のために静止摩擦係数と動摩擦係数とを同一として入力してもよい。
【００２０】
続いて、接触ペアの初期設定を行なう（ステップＳ１４０）。図４は、接触ペアを設定する様子を説明する説明図である。図中、白丸は接触面の一方の面で接触について着目している節点Ｎｃを示し、黒丸は接触面の他方の面の要素のうち節点Ｎｃと接触の可能性のある要素（ターゲット要素候補）Ｅｔ１〜Ｅｔ４を構成する節点Ｎｔ１〜Ｎｔ９を示す。接触ペアの設定は、まず、節点Ｎｃに最も近いターゲット要素候補の節点（図４では節点Ｎｔ５）を探索すると共に、探索した節点を構成節点として含む全ての要素を接触ペアとなりうるターゲット要素候補（図４では要素Ｅｔ１〜Ｅｔ４）として選定する。次に、各ターゲット要素候補Ｅｔ１〜Ｅｔ４に対する節点Ｎｃの位置の情報として分割比ｒ１，ｒ２を計算する。分割比ｒ１，ｒ２は、例えば、図４のターゲット要素候補Ｅｔ１では、節点Ｎｔ１と節点Ｎｔ２に対する節点Ｎｃの内分比を分割比ｒ１とし、節点Ｎｔ２と節点Ｎｔ５に対する節点Ｎｃの内分比を分割比ｒ２とするなどして計算することができる。いま、２つの節点の直線上に着目している節点Ｎｃが存在する場合を考える。着目している節点Ｎｃが２つの節点の中央にあるときには分割比は値０となり、着目している節点Ｎｃが一方の節点と一致するときに分割比は値１又は値−１となる。したがって、着目している節点Ｎｃが２つの節点間に存在するときには分割比は絶対値が値１以下となる。このことから、着目している節点Ｎｃのターゲット要素候補に対する分割比ｒ１，ｒ２の絶対値が共に値１以下となる要素を特定し、この要素を節点Ｎｃのターゲット要素Ｎｔとして特定することができる。なお、節点Ｎｃのターゲット要素候補に対する分割比ｒ１，ｒ２の絶対値が共に値１以下とならない場合、節点Ｎｃは非接触節点となる。このような判定を全ての節点に対して行なう。
【００２１】
接触ペアの初期設定を行なうと、全接触ペアの接触状態を固着と仮定する（ステップＳ１５０）。接触ペアの接触状態としては、コンタクト節点がターゲット要素にすべることなく固着している固着状態と、コンタクト節点がターゲット要素に固着しておらずすべりが生じるすべり状態とがある。実施例では、固着状態のコンタクト節点を固着節点と称し、すべり状態のコンタクト節点をすべり節点と称する。ステップＳ１５０では接触ペアの接触状態の初期設定であるから、適当に各接触ペアに対して固着状態とすべり状態とを仮定してもよいが、実施例では、簡易なものとして全接触ペアの接触状態を固着状態と仮定している。
【００２２】
こうして初期設定が終了すると、時間ｔにおける境界条件を入力し（ステップＳ１６０）、時間ｔの境界条件に対する各節点の状態の解をステップＳ１７０〜ステップＳ２３０までの繰り返し処理を行なうことにより求める。時間ｔにおける境界条件としては、前回の境界条件（時間ｔ−Δｔにおける境界条件）に対してΔｔ時間後の時間ｔで各部材の各部位に作用している力（トルクを含む）をどれだけ変化させたか或いは各部材の各部位に生じさせている変位（回転変位を含む）をどれだけ変化させたかの条件である。繰り返し処理は、各コンタクト節点の法線方向や摩擦方向の算出処理（ステップＳ１７０）、接触位置マトリックスや接触力マトリックスの作成処理（ステップＳ１８０）、接触位置マトリックスを用いた方程式と接触力マトリックスと各部材の部分剛性マトリックスを用いた方程式を連立させて求解する処理（ステップＳ１９０）、求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定する処理（ステップＳ２００）、求解後の各節点の接触状態が妥当であるか否かを判定する処理（ステップＳ２２０）によって構成されている。求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当ではないと判定されたときには求解後の接触パラメータを用いてステップＳ１７０〜Ｓ１９０を再度実行する。求解後の各節点の接触状態が一つでも妥当ではないと判定されたときには、接触状態が妥当ではないと判定された節点のうちの少なくとも一部の節点について接触状態を変更して（ステップＳ２３０）、ステップＳ１７０〜ステップＳ１９０を実行する。
【００２３】
繰り返しの各処理を説明する前に、これらの処理の説明に用いる数式等における基本的な変数について以下に示す。
【００２４】
【数１】
2017130811.tif 000003
【００２５】
【数２】
2017130811.tif 000004
【００２６】
【数３】
2017130811.tif 000005
【００２７】
ステップＳ１７０の法線方向と摩擦方向の算出処理は以下のように行なわれる。法線方向は、コンタクト節点のターゲット要素における法線方向を意味しており、ターゲット要素の各構成節点における法線方向の単位ベクトルを求め、各構成節点の単位ベクトルにコンタクト節点の分割比を用いた形状関数を乗じたものの和として計算することができる。摩擦方向は、コンタクト節点に作用する接触力のターゲット要素の面内成分であり、以下のように計算することができる。なお、摩擦方向は、コンタクト節点がすべり状態のときの摩擦力の計算に必要であるから、コンタクト節点がすべり節点であるときに計算される。コンタクト節点がすべり節点である事象としては、固着節点からすべり節点に変更されたときとすべり節点を維持しているときとを考えることができる。コンタクト節点が固着節点からすべり節点に変更されたときには、ｐ番目の接触ペアにおいて算出される接触力（固着力）を用いると摩擦方向は接触力の面内成分（水平成分）であるから、式（１）として計算することができる。ここで、式（１）では、接触状態として固着状態を考えているから、すべり状態を考えるための垂直成分については除去している。式（１）の右辺の分子は、固着力の水平成分（摩擦力成分）を示しており、右辺の分母は、その絶対値を示している。これにより、摩擦方向は単位ベクトルとして表わされる。一方、コンタクト節点がすべり節点を維持しているときには、接触ペアにおけるコンタクト節点の前回の変位と今回の変位との差とターゲット要素の前回の変位と今回の変位との差とのズレの増分をｄ（ｐ）とすると、ズレの増分ｄ（ｐ）は次式（２）により表わすことができ、ズレの増分ｄ（ｐ）の水平成分の単位ベクトルが摩擦方向となるから、式（３）として計算することができる。
【００２８】
【数４】
2017130811.tif 000006
【００２９】
ステップＳ１８０の接触位置マトリックスの作成は、以下のように行なわれる。コンタクト節点の接触状態が固着状態である場合には、すべりは生じないから変形の前後でのコンタクト節点とターゲット要素におけるコンタクト節点の位置は変化しない。このため、変形後のコンタクト節点の位置ベクトルと変形後のターゲット要素におけるコンタクト節点位置の位置ベクトルは一致する。これを式（４）として表わすことができる。ここで、式（４）の左辺は変形後のコンタクト節点の位置ベクトルであり、右辺は変形後のターゲット要素におけるコンタクト節点位置の位置ベクトルを示している。右辺では形状関数を用いることにより、ターゲット要素におけるコンタクト節点位置の位置ベクトルを表わしているのが解る。コンタクト節点の接触状態がすべり状態である場合には、接触面内での移動は可能であるが接触面外に存在することは許されない。このため、変形後のコンタクト節点の位置ベクトル（Ａ）と変形後のターゲット要素におけるコンタクト節点位置の位置ベクトル（Ｂ）との差のベクトル（Ａ−Ｂ）とコンタクト節点のターゲット要素に対する法線方向の単位ベクトルとの内積は値０となる（式（５））。ベクトル（Ａ），（Ｂ）を式（６），（７）として示す。全ての接触ペアに対して考えると、接触位置マトリックス[Ｍｃｃ］は、式（８）として表わされ、これを変形すると、式（９）となる。
【００３０】
【数５】
2017130811.tif 000007
【００３１】
ステップＳ１８０の接触力マトリックスの作成は、以下のように行なわれる。ｐ番目の接触ペアのコンタクト節点に生じる接触力は、コンタクト節点が固着状態のときには接触力は３方向成分であるから垂直抗力をゼロとして式（１０）のように表わすことができる。一方、コンタクト節点がすべり状態のときには接触力は垂直抗力と摩擦力であるから、接触力の３方向成分をゼロとし、法線方向の単位ベクトルに垂直抗力を乗じたものと摩擦方向の単位ベクトルに摩擦係数と垂直抗力を乗じたものとして式（１１）のように表すことができる。コンタクト節点に生じる接触力は作用反作用の法則により接触ペアをなすターゲット要素にも生じるから、ｐ番目の接触ペアのターゲット要素内のｋ番目の要素内節点番号への接触力の配分はコンタクト節点に生じる接触力に形状関数を乗じたものとなり、式(１２)として表わされる。全ての接触ペアに関する接触力は接触力ベクトル｛ｆｃ｝とまとめることができる。そして、この関係式を拡張することにより、接触による力ベクトルは、全ての接触ペアに関する接触力情報を用いて式（１３）のようにまとめることができ、接触力マトリックス[Ｍｃｆ]を得ることができる。
【００３２】
【数６】
2017130811.tif 000008
【００３３】
ステップＳ１９０の求解する処理は、作成した接触位置マトリックス［Ｍｃｃ］を用いた変位を拘束する方程式と接触力マトリックス［Ｍｃｆ］と各部材の部分剛性マトリックス［Ｋ］を用いた力の釣り合いの方程式とを連立させて求解する処理である。具体的には、外力は内力と接触力の和であり、内力は変位と部分剛性マトリックス［Ｋ］との積であるから、この力の釣り合いは次式（１４）として表わされる。一方、接触における変位の拘束から式（９）が導出されているから、式（１４）と式（９）とを統合して式（１５）を得ることができる。この式（１５）は、未知数と式の数が同じであるから、連立方程式として求解することが可能となる。
【００３４】
【数７】
2017130811.tif 000009
【００３５】
ステップＳ２００の求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定する処理は、前回の解析結果における接触パラメータと今回の解析結果における接触パラメータとが予め定めた誤差の範囲内で一致するか否かにより行なわれる。接触パラメータが妥当でないと判定されたときには今回の解析結果における接触パラメータを用いて再度ステップＳ１７０〜Ｓ１９０を実行する。
【００３６】
ステップＳ２２０の求解後の各節点の接触状態が妥当であるか否かを判定する処理は、以下のように行なわれる。第１に、コンタクト節点が固着節点であるときに接触力の接触面内方向成分と法線方向成分との比が静止摩擦係数より大きいという第１の矛盾が生じているときに妥当ではないと判定する。第２に、コンタクト節点がすべり節点であるときに仮定したすべり方向と求解により求められたすべり方向が概ね反対方向となるという第２の矛盾が生じているときに妥当ではないと判定する。第３に、コンタクト節点の垂直抗力が負の値となる第３の矛盾が生じているときに妥当ではないと判定する。第４に、接触ペアを持たない節点（非接触節点）が接触面の相手側に対して食い込みが生じているという第４の矛盾が生じているときに妥当ではないと判定する。食い込みの判定は、ターゲット要素候補の非接触節点位置の位置から非接触節点へのベクトルと法線方向の単位ベクトルとの内積が予め設定した負の許容値を下回ったときに食い込みと判断することにより行なうことができる。第１から第４のいずれかの矛盾が生じて節点の接触状態が妥当ではないと判定されたときには、妥当ではないと判定された節点うちの少なくとも一部について接触状態を変更し、ステップＳ１７０〜Ｓ２２０の処理を再度行なう。接触状態の変更は、第１の矛盾が生じているときには固着節点からすべり節点に変更することにより行ない、第２の矛盾が生じているときにはすべり節点から固着節点に変更することにより行ない、第３の矛盾が生じているときにはコンタクト節点（接触節点）から非接触節点に変更することにより行ない、第４の矛盾が生じているときには非接触節点からコンタクト節点（接触節点）に変更することにより行なうことができる。矛盾が生じている節点が複数存在する場合、実施例ではその一部についてだけ接触状態を変更するものとした。このようにするのは、解が振動したり発散する可能性があるからである。したがって、矛盾が生じている節点が複数存在する場合には、収束性を高めるために、矛盾が生じている節点のうち矛盾の程度が大きいものから順に一定割合の節点の接触状態を変更したり、第１ないし第４の矛盾の各々に閾値を設けて閾値を超える節点について接触状態を変更したりするのが好ましい。
【００３７】
ステップＳ１６０〜Ｓ２３０の繰り返し処理により解が収束し、ステップＳ２２０で求解後の全ての節点の接触状態に矛盾がなく妥当であると判定すると、各節点に生じる変位や力、接触ペア、接触状態、接触パラメータ（分割比，法線方向，摩擦方向）などを時間ｔの解データとして  記憶装置４６に記憶する（ステップＳ２４０）。そして、時間ｔが予め定めた最終時間ｔ（ｅｎｄ）に至っているか否かを判定し（ステップＳ２５０）、時間ｔが最終時間ｔ（ｅｎｄ）に至っていないときには、微小時間Δｔだけ時間ｔを増加して（ステップＳ２６０）、ステップＳ１６０の微小時間Δｔだけ時間ｔを増加したときの境界条件の入力に戻る。一方、時間ｔが最終時間ｔ（ｅｎｄ）に至っているときには、これで本プログラムを終了する。
【００３８】
以上説明した第１実施例の接触状態解析装置２０では、Ｎ個の部材のＭ組の接触面の一方の面の各節点に対して、非接触節点であるかすべり節点であるか固着節点であるかを仮定し、この仮定の下で、接触パラメータ（分割比，法線方向，摩擦方向）を用いて接触位置マトリックス［Ｍｃｃ］と接触力マトリックス［Ｍｃｆ］とを作成する処理と、接触位置マトリックス［Ｍｃｃ］を用いた方程式と接触力マトリックス［Ｍｃｆ］と各部材の部分剛性マトリックス［Ｋ］を用いた方程式とを連立させて求解する処理と、求解により得られた解において、各節点のうち接触状態の矛盾が生じる節点の一部について矛盾を解消するように接触状態を変更する処理と、を変更する処理で全ての節点の接触状態に矛盾が生じないようになるまで繰り返すことにより、各節点の接触状態を特定することができる。しかも、求解する処理でデータベース化して記憶された線形データにより作成される各部材の部分剛性マトリックスを用いて求解するから、繰り返し処理中に各部材の線形変形問題の解析（大規模線形問題の求解）を行なうことがない。このため、繰り返し処理中に大規模線形問題の解析が行なわれるものに比して解析時間を短縮することができる。また、求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定し、接触パラメータが妥当でないと判定したときには求解後の接触パラメータを用いて再度マトリックスを作成する処理と求解する処理とを実行するから、より適正な解を得ることができる。これらの結果、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる。
【００３９】
また、第１実施例の接触状態解析装置２０では、各節点の接触状態だけでなく、時間ｔに対して解析した各節点に生じる変位や力も記憶するから、各節点に生じる変位や力の時間変化を得ることもできる。
【００４０】
また、第１実施例で説明した図２の接触状態解析プログラムは、コンピュータにインストールすることにより、コンピュータを、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる接触状態解析装置２０として機能させることができる。
【００４１】
第１実施例の接触状態解析装置２０では、「接触面」には、面の一部が接触している面が含まれるだけでなく、面の一部が接触している可能性がある面も含まれるものとしたが、「接触面」には面の一部が接触している面だけが含まれるものとしても構わない。
【００４２】
次に、本開示の第２実施例の接触状態解析装置１２０について説明する。図５は、第２実施例の接触状態解析装置１２０の構成の概略を示す構成図である。第２実施例の接触状態解析装置１２０は、図示するように、一般的な汎用コンピュータ１２２にアプリケーションソフトウエアとしての接触状態解析プログラム１３０がインストールされたものとして構成されている。接触状態解析プログラム１３０は、部分剛性マトリックス構築モジュール１３２、初期状態設定モジュール１３３、接触力ベクトル算出モジュール１３４、求解モジュール１３８、ズレベクトル判定モジュール１３９、パラメータ判定モジュール１４０、接触状態判定モジュール１４２、を含んでいる。
【００４３】
部分剛性マトリックス構築モジュール１３２は、Ｍ組の接触面により連接しているＮ個の部材の任意の部位に作用する線形変形の範囲内の力（トルクを含む）と変位（ねじり変位を含む）との関係を有限要素法や差分法あるいは境界要素法などを用いて解析し、得られた解を線形データのデータベース１４４として記憶し、このデータベース１４４から部分剛性マトリックスを作成するモジュールである。ここで、第２実施例では、「接触面」には、面の一部が接触している面が含まれるだけでなく、面の一部が接触している可能性がある面も含まれる。
【００４４】
初期状態設定モジュール１３３は、Ｎ個の各部材の初期配置を仮定すると共に仮定した初期配置に基づいて各節点の初期のズレベクトルｇ（ｐ）を算出し、さらに、Ｍ組の接触面の各要素を構成する節点が接触しているコンタクト節点であるか非接触の非接触節点であるかを仮定すると共にコンタクト節点についてはすべりが生じているすべり節点であるか固着している固着節点であるかを仮定するモジュールである。
【００４５】
接触力ベクトル算出モジュール１３４は、コンタクト節点については、コンタクト節点と接触している要素（ターゲット要素）におけるコンタクト節点の接触位置をターゲット要素を構成する節点の内分比により示した分割比や、コンタクト節点位置におけるターゲット要素に対する法線方向を計算し、更に、すべり節点についてはコンタクト節点に生じる摩擦力の方向である摩擦方向を計算し、計算した分割比，法線方向，摩擦方向を用いて接触力ベクトルを計算するモジュールである。本明細書では、分割比と法線方向と摩擦方向を接触パラメータと称する。
【００４６】
求解モジュール１３８は、接触力ベクトル算出モジュールにより算出された接触力ベクトルと部分剛性マトリックス構築モジュール１３２により作成された部分剛性マトリックスとを用いて力の釣り合いの方程式を作成し、方程式を連立させて求解するモジュールである。
【００４７】
ズレベクトル判定モジュール１３９は、各節点のズレベクトルｇ（ｐ）を計算した結果、全てのズレベクトルｇ（ｐ）が許容値未満であるか否かを判定し、全てのズレベクトルｇ（ｐ）が許容値未満となっていないときには再度接触力ベクトル算出モジュール１３４と求解モジュール１３８とを実行するモジュールである。
【００４８】
パラメータ判定モジュール１４０は、求解モジュール１３８により得られる解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が予め定めた誤差の範囲内で一致するか否かを判定し、誤差の範囲内で一致しないときには求解後の接触パラメータを用いて再度接触力ベクトル算出モジュール１３４と求解モジュール１３８とズレベクトル判定モジュール１３９とを実行するモジュールである。
【００４９】
接触状態判定モジュール１４２は、求解モジュール１３８により得られる解において、各節点の接触状態に矛盾が生じているか否かを判定し、矛盾が生じているときには矛盾が生じている節点の少なくとも一部について矛盾を解消するように接触状態を変更して接触力ベクトル算出モジュール１３４と求解モジュール１３８とズレベクトル判定モジュール１３９とパラメータ判定モジュール１４０とを実行し、各節点の接触状態に矛盾が生じなくなったときに各節点の変位や接触状態、接触パラメータなどを解データとして記憶装置１４６に保存するモジュールである。接触状態の矛盾としては、固着節点において接触力の接触面内方向成分と法線方向成分の比が静止摩擦係数より大きいという第１の矛盾や、すべり節点において仮定した摩擦力方向と解により得られたすべり方向とが略同方向であるという第２の矛盾、コンタクト節点において垂直抗力が負となる第３の矛盾、非接触節点において変位が前記接触面の他方側に食い込むという第４の矛盾が想定されている。そして、接触状態の変更としては、第１の矛盾が生じている節点に対しては固着節点からすべり節点への変更や、第２の矛盾が生じている節点に対してはすべり節点から固着節点への変更、第３の矛盾が生じている節点に対してはコンタクト節点から非接触節点への変更、第４の矛盾が生じている節点に対して非接触節点からコンタクト節点への変更となる。
【００５０】
上記の各モジュールにより構成された接触状態解析プログラムの一例を示すフローチャートを図６に示す。以下に、図６のフローチャートを用いてＮ個の部材のＭ組の接触面における各要素を構成する節点の接触状態を解析する手法を詳細に説明する。接触状態解析プログラムの各処理を説明する前に、これらの処理の説明に用いる数式等における基本的な変数について以下に示す。
【００５１】
【数８】
2017130811.tif 000010
【００５２】
【数９】
2017130811.tif 000011
【００５３】
接触状態解析プログラムが実行されると、まず、図３に例示するＮ個の部材がＭ組の接触面により接触している解析対象において、Ｍ組の接触面を構成する各部材に対して有限要素法解析を用いて線形データを求解すると共に線形データをデータベース１４４に記憶し、記憶した線形データを用いて部分剛性マトリックスを作成する（ステップＳ３００）。具体的には、まず、解析対象の各部材の任意の部位に線形変形の範囲内となる所定の力（例えば単位力や単位トルク）や所定の変位（例えば、単位変位や単位ねじり変位）を作用させたときの節点にどのような力が作用しどのような変位が生じるかについての関係を所定の力や所定の変位の作用部位を変化させて各々求解し、所定の力や所定の変位、作用部位、各節点に作用する力、各節点の変位などを線形データとしてデータベース１４４に記憶する。そして、記憶した線形データを用いて、Ｎ個の部材に対して節点に対する任意の作用力に対する各要素の変位や力、節点に対する任意の作用変位に対する各要素に作用する力や変位を求める部分剛性マトリックスを作成する。
【００５４】
次に、Ｍ組の接触面の摩擦係数を入力し（ステップＳ３１０）、各部材の配置を設定し（ステップＳ３２０）、時間の初期化（ｔ＝０）を行なう（ステップＳ３３０）。各部材の配置を設定することにより各節点の初期変位は定まる。ここで、摩擦係数は、静止摩擦係数と動摩擦係数とに区分して入力してもよいし、速度や荷重の関数としてもよい。また、簡単のために静止摩擦係数と動摩擦係数とを同一として入力してもよい。各部材の配置としては、（１）各部材が接する配置、（２）各部材がオーバーラップする配置のいずれでもよいが、解の収束性を考慮すると、（２）各部材がオーバーラップする配置とするのが好ましい。なお、オーバーラップが存在しない場合には、後述するように、ズレベクトルｇ（ｐ）が妥当ではないと判断され、各部材の配置が変更される。
【００５５】
続いて、図４を用いて説明したように、接触ペアの初期設定を行なうと共に、初期のズレベクトルｇ（ｐ）を算出する（ステップＳ３４０）。ズレベクトルｇ（ｐ）は、ｐ番目の接触ペアにおけるコンタクト節点とターゲット要素のズレであり、次式（１６）により計算される。
【００５６】
【数１０】
2017130811.tif 000012
【００５７】
初期のズレベクトルｇ（ｐ）を算出すると、ズレベクトルｇ（ｐ）は妥当であるか否かを判定し（ステップＳ３５０）、ズレベクトルｇ（ｐ）が妥当ではないと判定されたときには、変位ベクトルを修正したり各部材の配置を変更し（ステップＳ３６０）、再度ズレベクトルｇ（ｐ）を計算する（ステップＳ３７０）。ズレベクトルｇ（ｐ）が妥当であるか否かの判定は、各接触面においてオーバーラップが存在するか否かにより行なうことができる。接触面においてオーバーラップが全く存在しないときには接触状態であると判別できないため、妥当でないと判定する。この場合、各部材の配置として各部材がオーバーラップする配置に変更したり、変位ベクトルをオーバーラップが生じるように変更し、再度ズレベクトルｇ（ｐ）を計算するのである。なお、ズレベクトルｇ（ｐ）は妥当であるか否かの判定は、ズレベクトルｇ（ｐ）と法線方向の内積の正負によっても判断することができる。
【００５８】
妥当なズレベクトルｇ（ｐ）が設定されると、全接触ペアの接触状態を固着と仮定すると共に全接触ペアの接触力を値０（ゼロ）と仮定する（ステップＳ３８０）。接触ペアの接触状態としては、コンタクト節点がターゲット要素とすべることなく固着している固着状態と、コンタクト節点がターゲット要素に固着しておらずすべりが生じるすべり状態とがある。第２実施例では、固着状態のコンタクト節点を固着節点と称し、すべり状態のコンタクト節点をすべり節点と称する。ステップＳ３８０では接触ペアの接触状態の初期設定であるから、適当に各接触ペアに対して固着状態とすべり状態とを仮定してもよいが、第２実施例では、簡易なものとして全接触ペアの接触状態を固着状態と仮定している。また、接触力についても適当に仮定してもよいが、第２実施例では簡易なものとして値０と仮定している。
【００５９】
こうして初期設定が終了すると、時間ｔにおける境界条件を入力し（ステップＳ３９０）、時間ｔの境界条件に対する各節点の状態の解をステップＳ４００〜ステップＳ４６０までの繰り返し処理を行なうことにより求める。時間ｔにおける境界条件としては、前回の境界条件（時間ｔ−Δｔにおける境界条件）に対してΔｔ時間後の時間ｔで各部材の各部位に作用している力（トルクを含む）をどれだけ変化させたか或いは各部材の各部位に生じさせている変位（回転変位を含む）をどれだけ変化させたかの条件である。繰り返し処理は、各コンタクト節点の法線方向や摩擦方向の算出処理（ステップＳ４００）、接触力ベクトルの算出処理（ステップＳ４１０）、力の釣り合い方程式を求解する処理（ステップＳ４２０）、求解の前後でズレベクトルｇ（ｐ）が妥当であるか否かを判定する処理（ステップＳ４３０）、求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定する処理（ステップＳ４４０）、求解後の各節点の接触状態が妥当であるか否かを判定する処理（ステップＳ４６０）によって構成されている。求解の前後でズレベクトルｇ（ｐ）が妥当ではないと判定されたときには、求解された変位ベクトルを用いて、再度ズレベクトルｇ（ｐ）の算出を行ない、ステップＳ４００〜Ｓ４２０を実行する。求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当ではないと判定されたときには求解後の接触パラメータを用いて（ステップＳ４５０）、ステップＳ４００〜Ｓ４３０を再度実行する。求解後の各節点の接触状態が一つでも妥当ではないと判定されたときには、接触状態が妥当ではないと判定された節点のうちの少なくとも一部の節点について接触状態を変更して（ステップＳ４７０）、ステップＳ４００〜ステップＳ４６０を実行する。
【００６０】
ステップＳ４００の法線方向と摩擦方向の算出処理は以下のように行なわれる。法線方向は、コンタクト節点のターゲット要素における法線方向を意味しており、ターゲット要素の各構成節点における法線方向の単位ベクトルを求め、各構成節点の単位ベクトルにコンタクト節点の分割比を用いた形状関数を乗じたものの和として計算することができる。摩擦方向は、コンタクト節点に作用する接触力のターゲット要素の面内成分であり、以下のように計算することができる。なお、摩擦方向は、コンタクト節点がすべり状態のときの摩擦力の計算に必要であるから、コンタクト節点がすべり節点であるときに計算される。コンタクト節点がすべり節点である事象としては、固着節点からすべり節点に変更されたときとすべり節点を維持しているときとを考えることができる。コンタクト節点が固着節点からすべり節点に変更されたときには、ｐ番目の接触ペアにおいて算出される接触力（固着力）を用いると摩擦方向は接触力の面内成分であるから、次式（１７）として計算することができる。ここで、式（１７）では、すべり状態を考えるための垂直成分については除去している。式（１７）の右辺の分子は、固着力の摩擦力成分を示しており、右辺の分母は、その絶対値を示している。これにより、摩擦方向は単位ベクトルとして表わされる。一方、コンタクト節点がすべり節点を維持しているときには、接触ペアにおけるコンタクト節点の前回の変位と今回の変位との差とターゲット要素の前回の変位と今回の変位との差とのズレの増分をｄ（ｐ）とすると、ズレの増分ｄ（ｐ）は式（１８）により表わすことができ、ズレの増分ｄ（ｐ）の水平成分の単位ベクトルが摩擦方向となるから、式（４）として計算することができる。
【００６１】
【数１１】
2017130811.tif 000013
【００６２】
ステップＳ４１０の接触力ベクトルの算出処理は以下のように行なわれる。まず、ｐ番目の接触ペアのコンタクト節点に生じるコンタクト節点接触力ベクトルｆｃｐ（ｐ）を法線方向に生じる力ｆｎ（ｐ）と２方向の面内方向に生じる力ｆｔ１（ｐ），ｆｔ２（ｐ）とを用いて次式（２０）として表わすと、コンタクト節点が固着状態のときには式（２１）により表わされ、コンタクト節点がすべり状態のときには式（２２）により表わされる。式（２１）中の「Ｋｃ」は接触剛性である。なお、この接触剛性は解析途中で変更しても構わない。そして、ｐ番目の接触ペアのコンタクト節点に生じる接触力は、式（２３）のように表わすことができる。コンタクト節点に生じる接触力は作用反作用の法則により接触ペアをなすターゲット要素にも生じるから、ｐ番目の接触ペアのターゲット要素内のｋ番目の要素内節点番号への接触力の配分はコンタクト節点に生じる接触力に形状関数を乗じたものとなり、式（２４）として表わされる。全ての接触ペアに関する接触力の総和ととることにより接触力ベクトル｛ｆｃ｝は式（２５）のように計算される。
【００６３】
【数１２】
2017130811.tif 000014
【００６４】
ステップＳ４２０の力の釣り合い方程式を求解する処理は以下のように行なわれる。力の釣り合いは次式（２６）として表わされる。式（２６）中、接触力ベクトル｛ｆｃ｝は既知であるため、未知数と式の通が同一となり、連立方程式の求解が可能となる。
【００６５】
【数１３】
2017130811.tif 000015
【００６６】
ステップＳ４３０のズレベクトルｇ（ｐ）が妥当であるか否かを判定する処理は、繰り返し計算により収束するのが前提であるため、その大きさが許容値よりも小さくなる必要がある。したがって、ズレベクトルｇ（ｐ）のいずれかの絶対値が許容値よりも大きいか否かにより、ズレベクトルｇ（ｐ）が妥当であるか否かを判定することができる。ズレベクトルｇ（ｐ）のいずれかの絶対値が許容値よりも大きいときには、再度ステップＳ４００〜Ｓ４２０を実行する。
【００６７】
ステップＳ４４０の求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定する処理は、前回の解析結果における接触パラメータと今回の解析結果における接触パラメータとが予め定めた誤差の範囲内で一致するか否かにより行なわれる。接触パラメータが妥当でないと判定されたときには今回の解析結果における接触パラメータを用いて再度ステップＳ４００〜Ｓ４３０を実行する。
【００６８】
ステップＳ４６０の求解後の各節点の接触状態が妥当であるか否かを判定する処理は、以下のように行なわれる。第１に、コンタクト節点が固着節点であるときに接触力の接触面内方向成分と法線方向成分との比が静止摩擦係数より大きいという第１の矛盾が生じているときに妥当ではないと判定する。第２に、コンタクト節点がすべり節点であるときに仮定したすべり方向と求解により求められたすべり方向が概ね同じ方向となるという第２の矛盾が生じているときに妥当ではないと判定する。第３に、コンタクト節点の垂直抗力が負の値となる第３の矛盾が生じているときに妥当ではないと判定する。第４に、接触ペアを持たない節点（非接触節点）が接触面の相手側に対して食い込みが生じているという第４の矛盾が生じているときに妥当ではないと判定する。食い込みの判定は、ターゲット要素候補への非接触節点の投影位置から非接触節点へのベクトルと法線方向の単位ベクトルとの内積が予め設定した正でない（ゼロ以下）の許容値を下回ったときに食い込みと判断することにより行なうことができる。第１から第４のいずれかの矛盾が生じて節点の接触状態が妥当ではないと判定されたときには、妥当ではないと判定された節点うちの少なくとも一部について接触状態を変更し（ステップＳ４７０）、ステップＳ４００〜Ｓ４４０の処理を再度行なう。接触状態の変更は、第１の矛盾が生じているときには固着節点からすべり節点に変更することにより行ない、第２の矛盾が生じているときにはすべり節点から固着節点に変更することにより行ない、第３の矛盾が生じているときにはコンタクト節点（接触節点）から非接触節点に変更することにより行ない、第４の矛盾が生じているときには非接触節点からコンタクト節点（接触節点）に変更することにより行なうことができる。矛盾が生じている節点が複数存在する場合、第２実施例ではその一部についてだけ接触状態を変更するものとした。このようにするのは、解が振動したり発散する可能性があるからである。したがって、矛盾が生じている節点が複数存在する場合には、収束性を高めるために、矛盾が生じている節点のうち矛盾の程度が大きいものから順に一定割合の節点の接触状態を変更したり、第１ないし第４の矛盾の各々に閾値を設けて閾値を超える節点について接触状態を変更したりするのが好ましい。
【００６９】
ステップＳ４００〜Ｓ４６０の繰り返し処理により解が収束し、ステップＳ４６０で求解後の全ての節点の接触状態に矛盾がなく妥当であると判定すると、各節点に生じる変位や力、接触ペア、接触状態、接触パラメータ（分割比，法線方向，摩擦方向）などを時間ｔの解データとして  記憶装置１４６に記憶する（ステップＳ４８０）。そして、時間ｔが予め定めた最終時間ｔ（ｅｎｄ）に至っているか否かを判定し（ステップＳ４９０）、時間ｔが最終時間ｔ（ｅｎｄ）に至っていないときには、微小時間Δｔだけ時間ｔを増加して（ステップＳ５００）、ステップＳ３９０の微小時間Δｔだけ時間ｔを増加したときの境界条件の入力に戻る。一方、時間ｔが最終時間ｔ（ｅｎｄ）に至っているときには、これで本プログラムを終了する。
【００７０】
以上説明した第２実施例の接触状態解析装置１２０では、まず、Ｎ個の部材がＭ組の接触面により接触している状態で各部材に対して有限要素法解析を用いて線形データを求解すると共に線形データをデータベース１４４に記憶しておく。次に、Ｎ個の部材の配置を仮定して各節点の初期のズレベクトルｇ（ｐ）を計算し、Ｎ個の部材のＭ組の接触面の一方の面の各節点に対して、非接触節点であるかすべり節点であるか固着節点であるかを仮定する。そして、この仮定の下で、接触パラメータ（分割比，法線方向，摩擦方向）を用いて接触力ベクトルを計算する処理と、線形データと接触力ベクトルとを用いて力の釣り合いの方程式を作成すると共に作成した方程式を連立して求解する処理と、を求解においてズレベクトルｇ（ｐ）が許容値未満になるまで繰り返す。そして、求解により得られた解において、各節点のうち接触状態の矛盾が生じなくなるまで矛盾が生じている節点の接触状態を変更して上記の繰り返し処理を繰り返す。これにより、各節点の接触状態を特定する。求解する処理でデータベース化して記憶された線形データを用いて力の釣り合いの方程式を作成して求解するから、繰り返し処理中に各部材の線形変形問題の解析（大規模線形問題の求解）を行なうことがない。このため、繰り返し処理中に大規模線形問題の解析が行なわれるものに比して解析時間を短縮することができる。また、求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定し、接触パラメータが妥当でないと判定したときには求解後の接触パラメータを用いて再度接触力ベクトルを作成する処理と求解する処理とを実行するから、より適正な解を得ることができる。これらの結果、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる。
【００７１】
また、第２実施例の接触状態解析装置１２０では、各節点の接触状態だけでなく、時間ｔに対して解析した各節点に生じる変位や力も記憶するから、各節点に生じる変位や力の時間変化を得ることもできる。
【００７２】
また、第２実施例で説明した図６の接触状態解析プログラムは、コンピュータにインストールすることにより、コンピュータを、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる接触状態解析装置１２０として機能させることができる。
【００７３】
第２実施例の接触状態解析装置１２０では、「接触面」には、面の一部が接触している面が含まれるだけでなく、面の一部が接触している可能性がある面も含まれるものとしたが、「接触面」には面の一部が接触している面だけが含まれるものとしても構わない。
【００７４】
次に、本開示の第３実施例としての接触状態解析装置２２０について説明する。図７は、第３実施例の接触状態解析装置２２０の構成の概略を示す構成図である。第３実施例の接触状態解析装置２２０は、図示するように、一般的な汎用コンピュータ２２２にアプリケーションソフトウエアとしての接触状態解析プログラム２３０がインストールされたものとして構成されている。接触状態解析プログラム２３０は、部分剛性マトリックス構築モジュール２３２、初期状態設定モジュール２３３、コンタクト節点接触力ベクトル算出モジュール２３４、マトリックス作成モジュール、求解モジュール２３８、ズレベクトル判定モジュール２３９、パラメータ判定モジュール２４０、接触状態判定モジュール２４２、を含んでいる。
【００７５】
部分剛性マトリックス構築モジュール２３２は、第２実施例の初期状態設定モジュール１３３と同様に、Ｍ組の接触面により連接しているＮ個の部材の任意の部位に作用する線形変形の範囲内の力（トルクを含む）と変位（ねじり変位を含む）との関係を有限要素法や差分法あるいは境界要素法などを用いて解析し、得られた解を線形データのデータベース１４４として記憶し、このデータベース１４４から部分剛性マトリックスを作成するモジュールである。
【００７６】
初期状態設定モジュール２３３は、Ｎ個の各部材の初期配置を仮定すると共に仮定した初期配置に基づいて各節点の初期のズレベクトルｇ（ｐ）を算出し、さらに、Ｍ組の接触面の各要素を構成する節点が接触しているコンタクト節点であるか非接触の非接触節点であるかを仮定すると共にコンタクト節点についてはすべりが生じているすべり節点であるか固着している固着節点であるかを仮定し、さらに、全接触ペアのコンタクト節点接触力ベクトルのターゲット要素に対する面内方向の成分をゼロと仮定するモジュールである。
【００７７】
コンタクト節点接触力ベクトル算出モジュール２３４は、コンタクト節点については、コンタクト節点と接触している要素（ターゲット要素）におけるコンタクト節点の接触位置をターゲット要素を構成する節点の内分比により示した分割比や、コンタクト節点位置におけるターゲット要素に対する法線方向を計算し、更に、すべり節点についてはコンタクト節点に生じる摩擦力の方向である摩擦方向を計算し、計算した分割比，法線方向，摩擦方向を用いてコンタクト節点接触力ベクトルを計算するモジュールである。
【００７８】
マトリックス作成モジュール２３６は、コンタクト節点の接触状態（すべり節点であるか固着節点であるかの状態）と接触パラメータ（分割比，法線方向，摩擦方向）とを用いて、コンタクト節点の変位の法線方向の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、コンタクト節点の力の釣り合いに基づいてコンタクト節点接触力ベクトル（法線方向成分と面内方向２成分）に乗じたときに接触により生じる力ベクトルとなる接触力マトリックスを作成するモジュールである。具体的には、接触位置マトリックスは、コンタクト節点の変位後の位置がターゲット要素の変位後の接触面内であることに基づいて作成され、接触力マトリックスは、コンタクト節点の接触力がターゲット要素におけるコンタクト節点の位置の接触力と同じ大きさで反対向きの力であって、コンタクト節点がすべり節点であるときにはコンタクト節点に生じる接触力がターゲット要素に対する法線方向の力とターゲット要素の面内方向の摩擦方向の摩擦力の合力であることに基づいて作成される。
【００７９】
求解モジュール２３８は、マトリックス作成モジュール２３６により作成した接触位置マトリックスを用いた方程式と接触力マトリックスと部分剛性マトリックス構築モジュール２３２により作成された部分剛性マトリックスを用いた方程式とを連立させて求解するモジュールである。
【００８０】
ズレベクトル判定モジュール２３９は、各節点のズレベクトルｇ（ｐ）を計算した結果、全てのズレベクトルｇ（ｐ）が許容値未満であるか否かを判定し、全てのズレベクトルｇ（ｐ）が許容値未満となっていないときには再度コンタクト節点接触力ベクトル算出モジュール２３４とマトリックス作成モジュール２３６と求解モジュール２３８とを実行するモジュールである。
【００８１】
パラメータ判定モジュール２４０は、求解モジュール２３８により得られる解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が予め定めた誤差の範囲内で一致するか否かを判定し、誤差の範囲内で一致しないときには求解後の接触パラメータを用いて再度コンタクト節点接触力ベクトル算出モジュール２３４とマトリックス作成モジュール２３６と求解モジュール２３８とズレベクトル判定モジュール２３９とを実行するモジュールである。
【００８２】
接触状態判定モジュール２４２は、求解モジュール２３８により得られる解において、各節点の接触状態に矛盾が生じているか否かを判定し、矛盾が生じているときには矛盾が生じている節点の少なくとも一部について矛盾を解消するように接触状態を変更してコンタクト節点接触力ベクトル算出モジュール２３４とマトリックス作成モジュール２３６と求解モジュール１３８とズレベクトル判定モジュール１３９とパラメータ判定モジュール２４０とを実行し、各節点の接触状態に矛盾が生じなくなったときに各節点の変位や接触状態、接触パラメータなどを解データとして記憶装置２４６に保存するモジュールである。接触状態の矛盾および接触状態の変更については、第２実施例と同一である。
【００８３】
上記の各モジュールにより構成された接触状態解析プログラムの一例を示すフローチャートを図８に示す。以下に、図８のフローチャートを用いてＮ個の部材のＭ組の接触面における各要素を構成する節点の接触状態を解析する手法を詳細に説明する。これらの処理の説明に用いる数式等における基本的な変数は、第２実施例に示したものの他、以下のものを挙げることができる。
【００８４】
接触状態解析プログラムが実行されると、第２実施例と同様に、まず、図３に例示するＮ個の部材がＭ組の接触面により接触している解析対象において、Ｍ組の接触面を構成する各部材に対して有限要素法解析を用いて線形データを求解すると共に線形データをデータベース１４４に記憶し、記憶した線形データを用いて部分剛性マトリックスを作成する（ステップＳ６００）。次に、Ｍ組の接触面の摩擦係数を入力し（ステップＳ６１０）、各部材の配置を設定し（ステップＳ６２０）、時間の初期化（ｔ＝０）を行なう（ステップＳ６３０）。各部材の配置を設定することにより各節点の初期変位は定まる。
【００８５】
続いて、図４を用いて説明したように、接触ペアの初期設定を行なうと共に、上述の式（１６）により初期のズレベクトルｇ（ｐ）を算出する（ステップＳ６４０）。そして、全接触ペアの接触状態を固着と仮定すると共に全接触ペアの接触力を値０（ゼロ）と仮定すると共に、全接触ペアのコンタクト節点接触力ベクトルｆｃｐ（ｐ）における面内方向の２成分をゼロと仮定する（ステップＳ６５０）。
【００８６】
こうして初期設定が終了すると、時間ｔにおける境界条件を入力し（ステップＳ６６０）、時間ｔの境界条件に対する各節点の状態の解をステップＳ６７０〜ステップＳ７４０までの繰り返し処理を行なうことにより求める。時間ｔにおける境界条件については第２実施例と同様である。繰り返し処理は、各コンタクト節点の法線方向や摩擦方向の算出処理（ステップＳ６７０）、面内方向の接触力ベクトルの算出処理（ステップＳ６８０）、接触位置マトリックスと接触力マトリックスを作成する処理（ステップＳ６９０）、接触位置マトリックスと接触力マトリックスを用いた連立方程式を求解する処理（ステップＳ７００）、求解の前後でズレベクトルｇ（ｐ）が妥当であるか否かを判定する処理（ステップＳ７１０）、求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定する処理（ステップＳ７２０）、求解後の各節点の接触状態が妥当であるか否かを判定する処理（ステップＳ７４０）によって構成されている。求解の前後でズレベクトルｇ（ｐ）が妥当ではないと判定されたときには、求解された変位ベクトルを用いて、再度ズレベクトルｇ（ｐ）の算出を行ない、ステップＳ６７０〜Ｓ７１０を実行する。求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当ではないと判定されたときには求解後の接触パラメータを用いて（ステップＳ７３０）、ステップＳ６７０〜Ｓ７２０を再度実行する。求解後の各節点の接触状態が一つでも妥当ではないと判定されたときには、接触状態が妥当ではないと判定された節点のうちの少なくとも一部の節点について接触状態を変更して（ステップＳ７５０）、ステップＳ６７０〜ステップＳ７４０を実行する。
【００８７】
ステップＳ６７０の法線方向と摩擦方向の算出処理は第２実施例のステップＳ４００の法線方向と摩擦方向の算出処理と同一である。
【００８８】
ステップＳ６８０の接触力ベクトルの算出処理は以下のように行なわれる。ｐ番目の接触ペアのコンタクト節点に生じるコンタクト節点接触力ベクトルｆｃｐ（ｐ）を法線方向に生じる力ｆｎ（ｐ）と２方向の面内方向に生じる力ｆｔ１（ｐ），ｆｔ２（ｐ）とを用いて次式（２７）として表わすと、コンタクト節点が固着状態のときには式（２８）により表わされ、コンタクト節点がすべり状態のときには式（２９）により表わされる。式（２８）中の「Ｋｃ」は接触剛性である。なお、この接触剛性は解析途中で変更しても構わない。
【００８９】
【数１４】
2017130811.tif 000016
【００９０】
ステップＳ６９０の接触位置マトリックスの作成は、以下のように行なわれる。コンタクト節点は接触面内での移動は可能であるが接触面外に存在することは許されない。このため、変形後のコンタクト節点の位置ベクトル（Ａ）と変形後のターゲット要素におけるコンタクト節点位置の位置ベクトル（Ｂ）との差のベクトル（Ａ−Ｂ）とコンタクト節点のターゲット要素に対する法線方向の単位ベクトルとの内積は値０となる（式（３０））。ベクトル（Ａ），（Ｂ）を式（３１），（３２）として示す。全ての接触ペアに対して考えると、接触位置マトリックス[Ｍｃｃ］は、式（３３）として表わされ、これを変形すると、式（３４）となる。
【００９１】
【数１５】
2017130811.tif 000017
【００９２】
ステップＳ６９０の接触力マトリックスの作成は、以下のように行なわれる。ｐ番目の接触ペアのコンタクト節点に生じる接触力は式（３５）のように表わすことができる。コンタクト節点に生じる接触力は作用反作用の法則により接触ペアをなすターゲット要素にも生じるから、ｐ番目の接触ペアのターゲット要素内のｋ番目の要素内節点番号への接触力の配分はコンタクト節点に生じる接触力に形状関数を乗じたものとなり、式（３６）として表わされる。全ての接触ペアに関する接触力は接触力ベクトル｛ｆｃｐ｝（式（３７））とまとめることができる。そして、この関係式を拡張することにより、接触による力ベクトルは、全ての接触ペアに関するコンタクト節点接触力ベクトルを用いて式（３８）のようにまとめることができ、接触力マトリックス[Ｍｃｆ]を得ることができる。
【００９３】
【数１６】
2017130811.tif 000018
【００９４】
ステップＳ７００の求解する処理は、作成した接触位置マトリックス［Ｍｃｃ］を用いた変位を拘束する方程式と接触力マトリックス［Ｍｃｆ］と各部材の部分剛性マトリックス［Ｋ］を用いた力の釣り合いの方程式とを連立させて求解する処理である。具体的には、外力は内力と接触力の和であり、内力は変位と部分剛性マトリックス［Ｋ］との積であるから、この力の釣り合いは次式（３９）として表わされる。一方、接触における変位の拘束から式（３４）が導出されているから、式（３９）と式（３４）とを統合して式（４０）を得ることができる。コンタクト節点接触力ベクトルｆｃｐ（ｐ）の面内方向に生じる自由度は既知であるため、この式（４０）は、未知数と式の数が同じとなり、連立方程式として求解することが可能となる。
【００９５】
【数１７】
2017130811.tif 000019
【００９６】
ステップＳ７１０のズレベクトルｇ（ｐ）が妥当であるか否かを判定する処理は、繰り返し計算により収束するのが前提であるため、その大きさが許容値よりも小さくなる必要がある。したがって、ズレベクトルｇ（ｐ）のいずれかの絶対値が許容値よりも大きいか否かにより、ズレベクトルｇ（ｐ）が妥当であるか否かを判定することができる。ズレベクトルｇ（ｐ）のいずれかの絶対値が許容値よりも大きいときには、再度ステップＳ６７０〜Ｓ７００を実行する。
【００９７】
ステップＳ７２０の求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定する処理は、前回の解析結果における接触パラメータと今回の解析結果における接触パラメータとが予め定めた誤差の範囲内で一致するか否かにより行なわれる。接触パラメータが妥当でないと判定されたときには今回の解析結果における接触パラメータを用いて再度ステップＳ６７０〜Ｓ７１０を実行する。
【００９８】
ステップＳ７４０の求解後の各節点の接触状態が妥当であるか否かを判定する処理は、第２実施例のステップＳ４６０の求解後の各節点の接触状態が妥当であるか否かを判定する処理と同様に、第１から第４のいずれかの矛盾が生じているか否かを判定する。第１から第４のいずれかの矛盾が生じて節点の接触状態が妥当ではないと判定されたときには、妥当ではないと判定された節点うちの少なくとも一部について接触状態を変更し（ステップＳ７５０）、ステップＳ６７０〜Ｓ７３０の処理を再度行なう。なお、接触状態の変更については第２実施例と同様である。
【００９９】
ステップＳ６７０〜Ｓ７５０の繰り返し処理により解が収束し、ステップＳ７４０で求解後の全ての節点の接触状態に矛盾がなく妥当であると判定すると、各節点に生じる変位や力、接触ペア、接触状態、接触パラメータ（分割比，法線方向，摩擦方向）などを時間ｔの解データとして  記憶装置２４６に記憶する（ステップＳ７６０）。そして、時間ｔが予め定めた最終時間ｔ（ｅｎｄ）に至っているか否かを判定し（ステップＳ７７０）、時間ｔが最終時間ｔ（ｅｎｄ）に至っていないときには、微小時間Δｔだけ時間ｔを増加して（ステップＳ７８０）、ステップＳ６６０の微小時間Δｔだけ時間ｔを増加したときの境界条件の入力に戻る。一方、時間ｔが最終時間ｔ（ｅｎｄ）に至っているときには、これで本プログラムを終了する。
【０１００】
以上説明した第３実施例の接触状態解析装置２２０では、まず、Ｎ個の部材がＭ組の接触面により接触している状態で各部材に対して有限要素法解析を用いて線形データを求解すると共に線形データをデータベース２４４に記憶しておく。次に、Ｎ個の部材の配置を仮定して各節点の初期のズレベクトルｇ（ｐ）を計算し、Ｎ個の部材のＭ組の接触面の一方の面の各節点に対して、非接触節点であるかすべり節点であるか固着節点であるかを仮定すると共にコンタクト節点接触力ベクトルｆｃｐ（ｐ）の面内方向の成分をゼロと仮定する。そして、この仮定の下で、接触パラメータ（分割比，法線方向，摩擦方向）を用いて接触力ベクトルを計算する処理と、線形データと接触力ベクトルとを用いて接触位置マトリックス［Ｍｃｃ］と接触力マトリックス［Ｍｃｆ］とを作成する処理と、接触位置マトリックス［Ｍｃｃ］を用いた変位の法線方向を拘束する方程式と接触力マトリックス［Ｍｃｆ］と各部材の部分剛性マトリックス［Ｋ］を用いた力の釣り合いの方程式とを連立させて求解する処理と、を求解においてズレベクトルｇ（ｐ）が許容値未満になるまで繰り返す。そして、求解により得られた解において、各節点のうち接触状態の矛盾が生じなくなるまで矛盾が生じている節点の接触状態を変更して上記の繰り返し処理を繰り返す。これにより、各節点の接触状態を特定する。求解する処理でデータベース化して記憶された線形データを用いて力の釣り合いの方程式を作成して求解するから、繰り返し処理中に各部材の線形変形問題の解析（大規模線形問題の求解）を行なうことがない。このため、繰り返し処理中に大規模線形問題の解析が行なわれるものに比して解析時間を短縮することができる。また、求解の前後で接触パラメータ（分割比，法線方向，摩擦方向）が妥当であるか否かを判定し、接触パラメータが妥当でないと判定したときには求解後の接触パラメータを用いて再度接触力ベクトルを作成する処理と求解する処理とを実行するから、より適正な解を得ることができる。これらの結果、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる。
【０１０１】
また、第３実施例の接触状態解析装置２２０では、各節点の接触状態だけでなく、時間ｔに対して解析した各節点に生じる変位や力も記憶するから、各節点に生じる変位や力の時間変化を得ることもできる。
【０１０２】
また、第３実施例で説明した図８の接触状態解析プログラムは、コンピュータにインストールすることにより、コンピュータを、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる接触状態解析装置２２０として機能させることができる。
【０１０３】
第３実施例の接触状態解析装置２２０では、「接触面」には、面の一部が接触している面が含まれるだけでなく、面の一部が接触している可能性がある面も含まれるものとしたが、「接触面」には面の一部が接触している面だけが含まれるものとしても構わない。
【０１０４】
以上に述べたとおり、本開示は以下の項目に記載の接触状態解析方法および接触状態解析装置並びにプログラムを含む。
〔項目１〕
Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点が非接触の状態の非接触節点であるかすべりが生じているすべり節点であるか固着している固着節点であるかのうちのいずれの接触状態であるかを解析する接触状態解析方法であって、
（１）前記Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておき、
（２）前記Ｎ個の部材の一部に対して力または変位を作用させたときに、
（ａ）前記接触面の一方の要素を構成する節点のうち前記接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか前記非接触節点であるかを仮定すると共に前記コンタクト節点についてはすべり節点であるか固着節点であるかを仮定し、前記コンタクト節点に対しては前記接触面の他方の要素のうち接触しているターゲット要素における前記コンタクト節点の接触位置に関する情報としての分割比と前記コンタクト節点における前記ターゲット要素に対する法線方向とを計算し、前記すべり節点については前記コンタクト節点での前記ターゲット要素における摩擦方向を計算する接触状態仮定ステップと、
（ｂ）前記コンタクト節点の接触状態，分割比，法線方向，摩擦方向を用いて、前記コンタクト節点の変位の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、前記コンタクト節点の力の釣り合いに基づいてコンタクト節点の接触力に関連するベクトルに乗じたときに各節点に接触により生じる力ベクトルとなる接触力マトリックスを作成するマトリックス作成ステップと、
（ｃ）前記接触位置マトリックスを用いた方程式と前記接触力マトリックスと前記線形データから作成される各部材の部分剛性マトリックスを用いた方程式とを連立させて未知数を求解する求解ステップと、
（ｄ）前記求解ステップにおける解において、各節点のうち接触状態の矛盾が生じる節点の少なくとも一部について矛盾を解消するように接触状態を変更する変更ステップと、
を前記求解ステップで求められた全ての節点の接触状態に矛盾が生じないようになるまで前記接触状態仮定ステップと前記マトリックス作成ステップと前記求解ステップと前記変更ステップとを繰り返すことにより、各節点の接触状態を特定する、
ことを特徴とする接触状態解析方法。
項目１の接触状態解析方法では、まず、Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておく。即ち、線形データをデータベース化して記憶しておくのである。次に、Ｎ個の部材の一部に対して力または変位を作用させたときに、全ての節点の接触状態に矛盾が生じないようになるまで接触状態仮定ステップとマトリックス作成ステップと求解ステップと変更ステップとを繰り返し、全ての節点の接触状態に矛盾が生じなくなったときの接触状態を各節点の接触状態として特定する。接触状態仮定ステップでは、接触面の一方の要素を構成する節点のうち接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか非接触節点であるかを仮定すると共にコンタクト節点についてはすべり節点であるか固着節点であるかを仮定する。また、コンタクト節点に対しては、接触面の他方の要素のうち接触しているターゲット要素におけるコンタクト節点の接触位置に関する情報としての分割比とコンタクト節点におけるターゲット要素に対する法線方向とを計算する。そして、すべり節点についてはコンタクト節点でのターゲット要素における摩擦方向を計算する。分割比は、ターゲット要素におけるコンタクト節点の接触位置に関する情報であるから、ターゲット要素を含む面内の異なる２方向の内分比として計算することができる。コンタクト節点におけるターゲット要素に対する法線方向は、ターゲット要素を構成する各節点で隣接する２つの要素内節点へのベクトルが構成する平面に対する法線方向の単位ベクトルを求めた後に分割比に応じて各節点での単位法線ベクトルを内挿して計算することができる。すべり節点でのターゲット要素における摩擦方向は、コンタクト節点とターゲット要素の変位の増分の接触面内成分として計算することができる。マトリックス作成ステップでは、コンタクト節点の接触状態，分割比，法線方向，摩擦方向を用いて、コンタクト節点の変位の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、コンタクト節点の力の釣り合いに基づいてコンタクト節点の接触力に関連するベクトルに乗じたときに各節点に接触により生じる力ベクトルとなる接触力マトリックスを作成する。求解ステップでは、接触位置マトリックスを用いた方程式と接触力マトリックスと線形データから作成される各部材の部分剛性マトリックスを用いた方程式とを連立させて未知数を求解する。変更ステップでは、求解ステップにおける解において、各節点のうち接触状態の矛盾が生じる節点の少なくとも一部について矛盾を解消するように接触状態を変更する。変更ステップで全ての節点の接触状態に矛盾が生じないようになるまで接触状態仮定ステップとマトリックス作成ステップと求解ステップと変更ステップとを繰り返すことにより、各節点の接触状態を特定することができる。求解ステップでは、データベース化して記憶された線形データから作成される各部材の部分剛性マトリックスを用いて求解するから、繰り返し処理中には各部材の線形変形問題の解析（大規模線形変形問題の求解）は行なわれない。このため、繰り返し処理中に大規模線形問題の解析が行なわれるものに比して解析時間を短縮することができる。これらの結果、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる。ここで、「接触面」としては、面の一部が接触している面を意味しているものとしてもよいし、面の一部が接触している面と面の一部が接触している可能性がある面との双方を意味しているものとしてもよい。また、線形データは、実測したものでもよいし、有限要素法により解析されたものでもよいし、差分法や境界要素法により解析されたものでもよい。
〔項目２〕
項目１の接触状態解析方法において、前記接触位置マトリックスは、前記コンタクト節点が前記固着節点であるときには該コンタクト節点の変位後の位置と前記ターゲット要素における前記コンタクト節点の位置の変位後の位置とが同一であること及び前記コンタクト節点が前記すべり節点であるときには前記コンタクト節点の変位後の位置が前記ターゲット要素の変位後の接触面内であることに基づくマトリックスであり、前記接触力マトリックスは、前記コンタクト節点の接触力が前記ターゲット要素における前記コンタクト節点の位置の接触力と同じ大きさで反対向きの力であって、前記コンタクト節点が前記すべり節点であるときには前記コンタクト節点の接触力が前記ターゲット要素に対する法線方向の力と前記ターゲット要素の面内方向の摩擦方向の摩擦力との合力であることに基づくマトリックスである。こうすれば、より適正なマトリックスを作成することができる。
〔項目３〕
項目１または２の接触状態解析方法において、前記接触状態の矛盾は、前記固着節点において接触力の接触面内方向成分と法線方向成分の比が静止摩擦係数より大きいという第１の矛盾、前記すべり節点において仮定した摩擦力方向と解により得られたすべり方向とが逆であるという第２の矛盾、前記コンタクト節点において垂直抗力が負となる第３の矛盾、前記非接触節点において変位が前記接触面の他方側に食い込むという第４の矛盾、のいずれかである。
〔項目４〕
項目３の接触状態解析方法において、、前記変更ステップにおける接触状態の変更は、前記第１の矛盾が生じている節点に対しては固着節点からすべり節点への変更であり、前記第２の矛盾が生じている節点に対してはすべり節点から固着節点への変更であり、前記第３の矛盾が生じている節点に対してはコンタクト節点から非接触節点への変更であり、前記第４の矛盾が生じている節点に対して非接触節点からコンタクト節点への変更である。こうすれば、より的確な接触状態の変更を行なうことができる。ここで、第２の矛盾における「逆」は、完全に逆である場合だけでなく、概ね逆である場合も含まれる。
〔項目５〕
項目１〜４のいずれかの接触状態解析方法において、前記接触状態仮定ステップにより仮定した前記コンタクト節点の分割比，法線方向，摩擦方向と前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向とのいずれかに所定以上の誤差が生じたときには前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向を用いて前記マトリックス作成ステップを実行する。こうすれば、より適正に求解することができる。
〔項目６〕
項目１〜５のいずれかの接触状態解析方法において、前記線形データは、有限要素法により解析されたデータである。
〔項目７〕
項目１〜６のいずれかの接触状態解析方法において、前記Ｎ個の部材の一部に時間変化する力または変位を作用させる際に、力または変位の作用開始から微小時間経過毎に、作用させた力または変位に対して前記（２）のステップにより特定された各節点の接触状態を記憶する。こうすれば、時間変化する力または変位に対して微小時間経過する毎の各節点の接触状態を出力することができる。ここで、「（２）のステップ」は、全ての節点の接触状態に矛盾が生じないようになるまで（ａ）接触状態仮定ステップと（ｂ）マトリックス作成ステップと（ｃ）求解ステップと（ｄ）変更ステップとを繰り返すことによって各節点の接触状態を特定するステップである。
〔項目８〕
項目８の接触状態解析装置は、項目１〜７のいずれかの接触状態解析方法を用いて接触状態解析を行う。項目８の接触状態解析装置では、項目１〜７のいずれかの接触状態解析方法を用いるから、項目１〜７のいずれかの接触状態解析方法が奏する効果、即ち、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができるという効果を奏することができる。
〔項目９〕
項目９のプログラムは、コンピュータを項目８の接触状態解析装置として機能させる。したがって、本開示のプログラムでは、項目１〜７のいずれかの接触状態解析方法を用いるから、コンピュータを項目１〜７のいずれかの接触状態解析方法が奏する効果、即ち、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができるという効果を奏する装置として機能させることができる。
〔項目１０〕
Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点が非接触の状態の非接触節点であるかすべりが生じているすべり節点であるか固着している固着節点であるかのうちのいずれの接触状態であるかを解析する接触状態解析方法であって、
（１）前記Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておき、
（２）前記Ｎ個の部材の初期配置を仮定すると共に該仮定した配置に基づいて各節点の初期のズレを算出し、
（３）前記Ｎ個の部材の一部に対して力または変位を作用させたときに、
（ａ）前記接触面の一方の要素を構成する節点のうち前記接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか前記非接触節点であるかを仮定すると共に前記コンタクト節点についてはすべり節点であるか固着節点であるかを仮定する接触状態仮定ステップと、
（ｂ）前記コンタクト節点に対しては前記接触面の他方の要素のうち接触しているターゲット要素における前記コンタクト節点の接触位置に関する情報としての分割比と前記コンタクト節点位置における前記ターゲット要素に対する法線方向とを計算し、前記すべり節点については前記コンタクト節点に生じる摩擦力の方向である摩擦方向を計算するパラメータ計算ステップと、
（ｃ）前記コンタクト節点の接触状態，分割比，法線方向，摩擦方向を用いて、各節点に生じる接触力ベクトルを計算する接触力算出ステップと、
（ｄ）前記線形データと前記接触力ベクトルとを用いて力の釣り合いの方程式を作成すると共に方程式を連立させて未知数を求解する求解ステップと、
（ｅ）前記求解ステップにおける解において節点のズレが許容値未満になるまで前記（ｂ）〜（ｄ）の各ステップを繰り返し実行するズレ起因繰り返しステップと、
（ｆ）前記求解ステップにおける解において、各節点のうち接触状態の矛盾が生じなくなるまで矛盾が生じている節点の接触状態を変更して前記（ｂ）〜（ｅ）の各ステップを繰り返し実行する接触状態起因繰り返しステップと、
を実行することにより、各節点の接触状態を特定する、
ことを特徴とする。
項目１０の接触状態解析方法では、まず、Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておく。即ち、線形データをデータベース化して記憶しておくのである。次に、Ｎ個の部材の初期配置を仮定すると共に、この仮定した配置に基づいて各節点の初期のズレを算出する。続いて、Ｎ個の部材の一部に対して力または変位を作用させたときに、接触面の一方の要素を構成する節点のうち接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか非接触節点であるかを仮定すると共にコンタクト節点についてはすべり節点であるか固着節点であるかを仮定する。次に、各節点のズレが許容値未満になるまで（ｂ）パラメータ計算ステップと（ｃ）接触力算出ステップと（ｄ）求解ステップとを繰り返す。（ｂ）パラメータ計算ステップでは、コンタクト節点に対しては接触面の他方の要素のうち接触しているターゲット要素におけるコンタクト節点の接触位置に関する情報としての分割比とコンタクト節点位置におけるターゲット要素に対する法線方向とを計算し、すべり節点についてはコンタクト節点に生じる摩擦力の方向である摩擦方向を計算する。分割比は、ターゲット要素におけるコンタクト節点の接触位置に関する情報であるから、ターゲット要素を含む面内の異なる２方向の内分比として計算することができる。コンタクト節点位置におけるターゲット要素に対する法線方向は、ターゲット要素を構成する各節点で隣接する２つの要素内節点へのベクトルが構成する平面に対する法線方向の単位ベクトルを求めた後に分割比に応じて各節点での単位法線ベクトルを内挿して計算することができる。すべり節点での摩擦方向は、コンタクト節点とターゲット要素の変位の増分の接触面内成分として計算することができる。（ｃ）接触力算出ステップでは、コンタクト節点の接触状態，分割比，法線方向，摩擦方向を用いて、各節点に生じる接触力ベクトルを計算する。（ｄ）求解ステップでは、記憶しておいた線形データと計算した接触力ベクトルとを用いて力の釣り合いの方程式を作成し、作成した方程式を連立させて未知数を求解する。そして、全ての節点の接触状態に矛盾が生じないようになるまで矛盾が生じている節点の接触状態を変更して（ｂ）パラメータ計算ステップと（ｃ）接触力算出ステップと（ｄ）求解ステップと（ｅ）ズレ起因繰り返しステップとを繰り返し、全ての節点の接触状態に矛盾が生じなくなったときの接触状態を各節点の接触状態として特定する。
項目１０の接触状態解析方法では、求解ステップで、データベース化して記憶された線形データを用いて力の釣り合いの方程式を作成し、作成した方程式を連立させて未知数を求解するから、繰り返し処理中には各部材の線形変形問題の解析（大規模線形変形問題の求解）は行なわれない。このため、繰り返し処理中に大規模線形問題の解析が行なわれるものに比して解析時間を短縮することができる。これらの結果、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる。ここで、「接触面」としては、面の一部が接触している面を意味しているものとしてもよいし、面の一部が接触している面と面の一部が接触している可能性がある面との双方を意味しているものとしてもよい。また、線形データは、実測したものでもよいし、有限要素法により解析されたものでもよいし、差分法や境界要素法により解析されたものでもよい。
〔項目１１〕
項目１０の接触状態解析方法において、前記Ｎ個の部材の初期配置に基づいて前記各節点の初期のズレを算出した結果、前記Ｎ個の部材の初期配置における各接触面においてオーバーラップが存在しないときには、前記Ｎ個の部材の初期配置の変更を行なうと共に前記各節点の初期のズレを算出する。こうすれば、オーバーラップがない場合の収束性の悪化を抑制し、収束性を高めることができる。
〔項目１２〕
項目１０または１１の接触状態解析方法において、前記接触状態の矛盾は、前記固着節点において接触力の接触面内方向成分と法線方向成分の比が静止摩擦係数より大きいという第１の矛盾、前記すべり節点において仮定した摩擦力方向と解により得られたすべり方向とが同方向であるという第２の矛盾、前記コンタクト節点において垂直抗力が負となる第３の矛盾、前記非接触節点において変位が前記接触面の他方側に食い込むという第４の矛盾、のいずれかである。
〔項目１３〕
項目１２の接触状態解析方法において、前記接触状態起因繰り返しステップにおける接触状態の変更は、前記第１の矛盾が生じている節点に対しては固着節点からすべり節点への変更であり、前記第２の矛盾が生じている節点に対してはすべり節点から固着節点への変更であり、前記第３の矛盾が生じている節点に対してはコンタクト節点から非接触節点への変更であり、前記第４の矛盾が生じている節点に対して非接触節点からコンタクト節点への変更である。こうすれば、より的確な接触状態の変更を行なうことができる。ここで、第２の矛盾における「同方向」は、完全に同方向である場合だけでなく、概ね同方向である場合も含まれる。
〔項目１４〕
項目１０〜１３のいずれかの接触状態解析方法において、前記接触状態仮定ステップにより仮定した前記コンタクト節点の分割比，法線方向，摩擦方向と前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向とのいずれかに所定以上の誤差が生じたときには前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向を用いて前記（ｂ）〜（ｆ）の各ステップを繰り返し実行するパラメータ起因繰り返しステップを有する。こうすれば、より適正に求解することができる。
〔項目１５〕
項目１０〜１４のいずれかの接触状態解析方法において、前記線形データは、有限要素法により解析されたデータである。
〔項目１６〕
項目１０〜１５のいずれかの接触状態解析方法において、前記Ｎ個の部材の一部に時間変化する力または変位を作用させる際に、力または変位の作用開始から微小時間経過毎に、作用させた力または変位に対して前記（３）のステップにより特定された各節点の接触状態を記憶する。こうすれば、時間変化する力または変位に対して微小時間経過する毎の各節点の接触状態を出力することができる。
〔項目１７〕
項目１７の接触状態解析装置は、項目１０〜１６のいずれかの接触状態解析方法を用いて接触状態解析を行なう。項目１７の接触状態解析装置では、項目１０〜１６のいずれかの接触状態解析方法を用いるから、項目１０〜１６のいずれかの接触状態解析方法が奏する効果、即ち、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができるという効果を奏することができる。
〔項目１８〕
項目１８のプログラムは、コンピュータを項目１７の接触状態解析装置として機能させる。したがって、項目１８のプログラムでは、項目１０〜１６のいずれかの接触状態解析方法法を用いるから、コンピュータを項目１０〜１６のいずれかの接触状態解析方法が奏する効果、即ち、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができるという効果を奏する装置として機能させることができる。
〔項目１９〕
Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点が非接触の状態の非接触節点であるかすべりが生じているすべり節点であるか固着している固着節点であるかのうちのいずれの接触状態であるかを解析する接触状態解析方法であって、
（１）前記Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておき、
（２）前記Ｎ個の部材の初期配置を仮定すると共に該仮定した配置に基づいて各節点の初期のズレを算出し、
（３）前記Ｎ個の部材の一部に対して力または変位を作用させたときに、
（ａ）前記接触面の一方の要素を構成する節点のうち前記接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか前記非接触節点であるかを仮定すると共に前記コンタクト節点についてはすべり節点であるか固着節点であるかを仮定し、更に各節点の面内方向に生じる力をゼロと仮定する接触状態仮定ステップと、
（ｂ）前記コンタクト節点に対しては前記接触面の他方の要素のうち接触しているターゲット要素における前記コンタクト節点の接触位置に関する情報としての分割比と前記コンタクト節点位置における前記ターゲット要素に対する法線方向とを計算し、前記すべり節点については前記コンタクト節点に生じる摩擦力の方向である摩擦方向を計算するパラメータ計算ステップと、
（ｃ）各節点の面内方向に生じる接触力ベクトルを計算する接触力算出ステップと、
（ｄ）前記コンタクト節点の接触状態，分割比，法線方向，摩擦方向，接触力ベクトルを用いて、前記コンタクト節点の変位の法線方向の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、前記コンタクト節点の力の釣り合いに基づいてコンタクト節点の接触力に関連するベクトルに乗じたときに各節点に接触により生じる力ベクトルとなる接触力マトリックスを作成するマトリックス作成ステップと、
（ｅ）前記接触位置マトリックスを用いた方程式と前記接触力マトリックスと前記線形データから作成される各部材の部分剛性マトリックスを用いた方程式とを連立させて未知数を求解する求解ステップと、
（ｆ）前記求解ステップにおける解において節点のズレが許容値未満になるまで前記（ｂ）〜（ｅ）の各ステップを繰り返し実行するズレ起因繰り返しステップと、
（ｇ）前記求解ステップにおける解において、各節点のうち接触状態の矛盾が生じなくなるまで矛盾が生じている節点の接触状態を変更して前記（ｂ）〜（ｆ）の各ステップを繰り返し実行する接触状態起因繰り返しステップと、
を実行することにより、各節点の接触状態を特定する、
ことを特徴とする。
この項目１９の接触状態解析方法では、まず、Ｎ個の部材の任意の部位に作用する力とその変位との間の関係を線形データとして記憶しておく。即ち、線形データをデータベース化して記憶しておくのである。次に、Ｎ個の部材の初期配置を仮定すると共に、この仮定した配置に基づいて各節点の初期のズレを算出する。続いて、Ｎ個の部材の一部に対して力または変位を作用させたときに、接触面の一方の要素を構成する節点のうち接触面の他方の要素のいずれかに接触している節点であるコンタクト節点であるか非接触節点であるかを仮定すると共にコンタクト節点についてはすべり節点であるか固着節点であるかを仮定し、更に各節点の面内方向に生じる力の自由度をゼロと仮定する。次に、各節点のズレが許容値未満になるまで（ｂ）パラメータ計算ステップと（ｃ）接触力算出ステップと（ｄ）マトリックス作成ステップと（ｅ）求解ステップとを繰り返す。（ｂ）パラメータ計算ステップでは、コンタクト節点に対しては接触面の他方の要素のうち接触しているターゲット要素におけるコンタクト節点の接触位置に関する情報としての分割比とコンタクト節点位置におけるターゲット要素に対する法線方向とを計算し、すべり節点についてはコンタクト節点に生じる摩擦力の方向である摩擦方向を計算する。分割比は、ターゲット要素におけるコンタクト節点の接触位置に関する情報であるから、ターゲット要素を含む面内の異なる２方向の内分比として計算することができる。コンタクト節点位置におけるターゲット要素に対する法線方向は、ターゲット要素を構成する各節点で隣接する２つの要素内節点へのベクトルが構成する平面に対する法線方向の単位ベクトルを求めた後に分割比に応じて各節点での単位法線ベクトルを内挿して計算することができる。すべり節点での摩擦方向は、コンタクト節点とターゲット要素の変位の増分の接触面内成分として計算することができる。（ｃ）接触力算出ステップでは、各節点に生じるコンタクト節点接触力ベクトルにおける面内方向に生じる力の自由度を計算する。（ｄ）マトリックス作成ステップでは、コンタクト節点の接触状態，分割比，法線方向，摩擦方向，コンタクト節点接触力ベクトルを用いて、コンタクト節点の変位の法線方向の拘束に基づいて各節点の変位前の位置ベクトルと各節点の変位ベクトルとの和のベクトルに乗じたときに値０となる接触位置マトリックスを作成すると共に、コンタクト節点の力の釣り合いに基づいてコンタクト節点の接触力に関連するベクトルに乗じたときに各節点に接触により生じる力ベクトルとなる接触力マトリックスを作成する。（ｅ）求解ステップでは、接触位置マトリックスを用いた方程式と接触力マトリックスと線形データから作成される各部材の部分剛性マトリックスを用いた方程式とを連立させて未知数を求解する。そして、全ての節点の接触状態に矛盾が生じないようになるまで矛盾が生じている節点の接触状態を変更して（ｂ）パラメータ計算ステップと（ｃ）接触力算出ステップと（マトリックス作成ステップと（ｅ）求解ステップと（ｆ）ズレ起因繰り返しステップとを繰り返し、全ての節点の接触状態に矛盾が生じなくなったときの接触状態を各節点の接触状態として特定する。
項目１９の接触状態解析方法では、求解ステップで、データベース化して記憶された線形データから作成される各部材の部分剛性マトリックスを用いて求解するから、繰り返し処理中には各部材の線形変形問題の解析（大規模線形変形問題の求解）は行なわれない。このため、繰り返し処理中に大規模線形問題の解析が行なわれるものに比して解析時間を短縮することができる。これらの結果、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができる。ここで、「接触面」としては、面の一部が接触している面を意味しているものとしてもよいし、面の一部が接触している面と面の一部が接触している可能性がある面との双方を意味しているものとしてもよい。また、線形データは、実測したものでもよいし、有限要素法により解析されたものでもよいし、差分法や境界要素法により解析されたものでもよい。
〔項目２０〕
項目１９の接触状態解析方法において、前記接触位置マトリックスは、前記コンタクト節点の変位後の位置が前記ターゲット要素の変位後の接触面内であることに基づくマトリックスであり、前記接触力マトリックスは、前記コンタクト節点の接触力が前記ターゲット要素における前記コンタクト節点の位置の接触力と同じ大きさで反対向きの力であって、前記コンタクト節点が前記すべり節点であるときには前記コンタクト節点の接触力が前記ターゲット要素に対する法線方向の力と前記ターゲット要素の面内方向の摩擦方向の摩擦力との合力であることに基づくマトリックスである。こうすれば、より適正なマトリックスを作成することができる。
〔項目２１〕
項目１９または２０の接触状態解析方法において、前記接触状態の矛盾は、前記固着節点において接触力の接触面内方向成分と法線方向成分の比が静止摩擦係数より大きいという第１の矛盾、前記すべり節点において仮定した摩擦力方向と解により得られたすべり方向とが同方向であるという第２の矛盾、前記コンタクト節点において垂直抗力が負となる第３の矛盾、前記非接触節点において変位が前記接触面の他方側に食い込むという第４の矛盾、のいずれかである。
〔項目２２〕
項目２１の接触状態解析方法において、前記変更ステップにおける接触状態の変更は、前記第１の矛盾が生じている節点に対しては固着節点からすべり節点への変更であり、前記第２の矛盾が生じている節点に対してはすべり節点から固着節点への変更であり、前記第３の矛盾が生じている節点に対してはコンタクト節点から非接触節点への変更であり、前記第４の矛盾が生じている節点に対して非接触節点からコンタクト節点への変更である。こうすれば、より的確な接触状態の変更を行なうことができる。ここで、第２の矛盾における「同方向」は、完全に同方向である場合だけでなく、概ね同方向である場合も含まれる。
〔項目２３〕
項目１９〜２２のいずれかの接触状態解析方法において、前記接触状態仮定ステップにより仮定した前記コンタクト節点の分割比，法線方向，摩擦方向と前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向とのいずれかに所定以上の誤差が生じたときには前記求解ステップにより得られた解による前記コンタクト節点の分割比，法線方向，摩擦方向を用いて前記（ｂ）〜（ｇ）の各ステップを繰り返し実行するパラメータ起因繰り返しステップを有する。こうすれば、より適正に求解することができる。
〔項目２４〕
項目１９〜２３のいずれかの接触状態解析方法において、前記線形データは、有限要素法により解析されたデータである。
〔項目２５〕
項目１９〜２４のいずれかの接触状態解析方法において、前記Ｎ個の部材の一部に時間変化する力または変位を作用させる際に、力または変位の作用開始から微小時間経過毎に、作用させた力または変位に対して前記（３）のステップにより特定された各節点の接触状態を記憶する。こうすれば、時間変化する力または変位に対して微小時間経過する毎の各節点の接触状態を出力することができる。
〔項目２６〕
項目２６の接触状態解析装置は、項目１９〜２５のいずれかの接触状態解析方法を用いて接触状態解析を行なう。項目２６の接触状態解析装置では、項目１９〜２５のいずれかの接触状態解析方法を用いるから、項目１９〜２５のいずれかの接触状態解析方法が奏する効果、即ち、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができるという効果を奏することができる。
〔項目２７〕
項目２７のプログラムは、コンピュータを項目２６の接触状態解析装置として機能させる。したがって、項目２７のプログラムでは、項目１９〜２５のいずれかの接触状態解析方法法を用いるから、コンピュータを項目１９〜２５のいずれかの接触状態解析方法が奏する効果、即ち、Ｎ個の部材のＭ組の接触面の各面における各要素を構成する各節点の接触状態をより的確に且つより迅速に特定することができるという効果を奏することができる。
【０１０５】
以上、本開示を実施するための形態について実施例を用いて説明したが、本開示はこうした実施例に何等限定されるものではなく、本開示の要旨を逸脱しない範囲内において、種々なる形態で実施し得ることは勿論である。
【０１０６】
本出願は、２０１６年１月２５日に出願された日本国特許出願第２０１６−０１１５６９号及び２０１６年３月２５日に出願された日本国特許出願第２０１６−０６１４１８号を優先権主張の基礎としており、引用によりその内容の全てが本明細書に含まれる。
【産業上の利用可能性】
【０１０７】
本開示は、機械製造産業などに利用可能である。
【図１】
2017130811.tif 000020
【図２】
2017130811.tif 000021
【図３】
2017130811.tif 000022
【図４】
2017130811.tif 000023
【図５】
2017130811.tif 000024
【図６】
2017130811.tif 000025
【図７】
2017130811.tif 000026
【図８】
2017130811.tif 000027
【国際調査報告】
2017130811.tif 000028
2017130811.tif 000029
2017130811.tif 000030
2017130811.tif 000031