(57)【要約】【課題】多種多様な歯車装置を試験対象とすることができ、容易な計測によって歯車のかみ合い状態を定量的に推定できる歯車装置のかみ合い状態推定方法を提供する。【解決手段】この歯車装置のかみ合い状態推定方法は、歯車装置の振動を計測するステップと、計測された振動の情報から歯車装置の周波数応答関数を取得するＦＲＦ取得ステップと、歯車装置の負荷条件下の静的伝達誤差を計算するＬＳＴＥ計算ステップとを有する。そして、ＦＲＦ取得ステップは、かみ合い成分の複数の振動応答の特性曲線を抽出する抽出ステップと、歯車装置の模擬周波数応答関数を求めるシミュレーションステップと、実数値ＧＡによりモデル化された歯車装置のパラメータの最適化を行う最適化ステップと、最適化により得られた模擬周波数応答関数と複数の振動応答の特性曲線とを合成して歯車装置の周波数応答関数とするＦＲＦ推定ステップとを有する。【選択図】図１３
2018025462.tif 000002
【特許請求の範囲】
【請求項１】
歯車装置の振動を計測する計測ステップと、
前記計測された振動の情報から前記歯車装置の周波数応答関数を取得するＦＲＦ取得ステップと、
前記取得された周波数応答関数と前記計測された振動の情報とから前記歯車装置の負荷条件下の静的伝達誤差を計算するＬＳＴＥ計算ステップと、
を有し、
前記ＦＲＦ取得ステップは、
前記計測された振動の情報からかみ合い成分の複数の振動応答の特性曲線を抽出する抽出ステップと、
前記歯車装置をモデル化し、部分構造合成法により前記モデル化された前記歯車装置の数値シミュレーションを行って模擬周波数応答関数を求めるシミュレーションステップと、
前記模擬周波数応答関数の関数形状が、前記かみ合い成分の複数の振動応答の特性曲線に合うように、実数値ＧＡにより前記モデル化された歯車装置のパラメータの最適化を行う最適化ステップと、
前記最適化されたパラメータが適用された前記数値シミュレーションによって得られた模擬周波数応答関数と重なるように、前記かみ合い成分の複数の振動応答の特性曲線を、各絶対値を変更し且つ連続的に接続して前記周波数応答関数とするＦＲＦ推定ステップとを有する、
歯車装置のかみ合い状態推定方法。
【請求項２】
負荷条件下の歯面接触解析によって、歯車のかみ合い状態を特定する複数のパラメータを変化させながら求められた負荷条件下の静的伝達誤差と、前記ＬＳＴＥ計算ステップで計算された負荷条件下の静的伝達誤差とを比較して、前記歯車装置のかみ合い状態を推定するステップをさらに含む、
請求項１記載の歯車装置のかみ合い状態推定方法。
【請求項３】
前記ＬＳＴＥ計算ステップは、前記振動の情報のうち、かみ合い周波数の非整数次の成分を除外した振動情報を用いて負荷条件下の静的伝達誤差を計算する、
請求項１又は請求項２に歯車装置のかみ合い状態推定方法。
【請求項４】
コンピュータに、
入力された歯車装置の振動の情報から前記歯車装置の周波数応答関数を取得するＦＲＦ取得手順、
前記取得された周波数応答関数と前記入力された振動の情報とから前記歯車装置の負荷条件下の静的伝達誤差を計算するＬＳＴＥ計算手順と、
を実行させるためのプログラムであって、
前記ＦＲＦ取得手順は、
前記入力された振動の情報からかみ合い成分の複数の振動応答の特性曲線を抽出する抽出手順と、
前記歯車装置をモデル化し、部分構造合成法により前記モデル化された前記歯車装置の数値シミュレーションを行って模擬周波数応答関数を求めるシミュレーション手順と、
前記模擬周波数応答関数の関数形状が、前記かみ合い成分の複数の振動応答の特性曲線に合うように、実数値ＧＡにより前記モデル化された歯車装置のパラメータの最適化を行う最適化手順と、
前記最適化されたパラメータが適用された前記数値シミュレーションによって得られた模擬周波数応答関数と重なるように、前記かみ合い成分の複数の振動応答の特性曲線を、各絶対値を変更し且つ連続的に接続して前記周波数応答関数とするＦＲＦ推定手順とを有するプログラム。
【請求項５】
コンピュータに、
負荷条件下の歯面接触解析によって、歯車のかみ合い状態を特定する複数のパラメータを変化させながら求められた負荷条件下の静的伝達誤差と、前記ＬＳＴＥ計算手順で計算された負荷条件下の静的伝達誤差とを比較して、前記歯車装置のかみ合い状態を推定する手順、
をさらに実行させるための請求項４記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、歯車装置のかみ合い状態推定方法およびプログラムに関する。
【背景技術】
【０００２】
歯車装置は、歯車対がかみ合って回転運動を伝達する装置である。歯車装置には振動騒音という課題がある。
【０００３】
歯車装置の振動騒音は、歯対がかみ合うときの励振が原因となって発生する。励振の要素には、（１）歯車対のかみ合い剛性の変動による係数励振と、（２）理論歯面からの実歯面形状の偏差による強制変位とがある。そして、これらの２つの励振が総合された量として、負荷条件下の静的伝達誤差（loaded static transmission error、以下「ＬＳＴＥ」と記す。）がある。ＬＳＴＥは、歯車のかみ合い起振力を表す指標として歯車装置の設計及び評価に広く用いられている。
【０００４】
ＬＳＴＥには、歯車の諸元及び伝達荷重の他に、歯車自体の精度、歯車周辺の各部品の精度、歯車および各部品の組付精度、ならびに、負荷による歪みといった各因子の偏差が重なり合って、互いにかみ合う一対の歯車の伝達誤差としてμｍ単位で表れる。従来、例えば歯車装置の設計時において、上記各因子の歯車への影響が考慮されることはあった。しかしながら、歯車装置が実際に運転されているときの歯車の伝達誤差を、μｍ単位で予測することは困難であった。
【０００５】
また、本願発明に関連する技術として、特許文献１には、歯車の歯面形状を検査する歯車検査装置が開示されている。この歯車検査装置は、基準歯車に検査歯車をかみ合わせ基準歯車に負荷を加える。そして、歯車検査装置は、検査歯車を回転したときに生じる振動加速度を計測する。さらに、特許文献１の技術では、基準歯車と検査歯車について周波数応答関数が別途計算される。そして、計測された振動加速度と計算された周波数応答関数とから、検査歯車の歯面形状が求められる。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００７−００３５０７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
特許文献１の歯車検査装置は、検査歯車の個々の歯面形状を求めるのに有効である。しかしながら、特許文献１の技術は、予め歯面形状を詳細に計測した基準歯車が必要である。さらに、特許文献１の技術は、１つの歯車を単体で検査する必要があり、ハウジングの内部で一組又は複数組の歯車がかみ合って構成される歯車装置をそのまま検査することはできない。
【０００８】
また、特許文献１の技術では、計測された振動加速度のうち、歯車のかみ合い周波数の非整数次成分から周波数応答関数の形状が求められる（特許文献１の段落００４８を参照）。そして、この周波数応答関数により、かみ合い成分の特徴を除外した検査が行われる（特許文献１の図７、８のグラフの欠損部分を参照、横軸３０、６０、９０の部分がかみ合い成分に相当）。ここで、かみ合い成分とは、歯対がかみ合う期間を一周期としたかみ合い周波数とその高調波とを含む周波数成分を意味する。一方、歯車装置の振動騒音には歯車のかみ合い成分の特徴が大きく影響する。従って、特許文献１の技術には、歯車装置の振動騒音に関する解析を行う場合に、かみ合い成分の特徴が除外されてしまうという欠点がある。
【０００９】
なお、特許文献１には、検査歯車の振動加速度からかみ合い成分を抽出し、周波数応答関数の形状を求めることが説明されている（特許文献１の段落００８９を参照）。しかし、この方法では周波数応答関数を定量的に求めることはできず、それ故、この周波数応答関数を用いた歯車の励振力の解析は定性的なものに留まっている（特許文献１の図２５−図２８およびその説明を参照）。
【００１０】
ところで、歯車装置の運転時の特性について詳細な解析を行うには、歯車装置のＬＳＴＥが求められるとよい。しかしながら、ＬＳＴＥの計測は、専用の試験機を必要とする。このため、試験機を、様々なサイズ又は様々な軸配置を有する産業用の歯車装置に広く適用してＬＳＴＥを計測することは困難である。
【００１１】
また、たとえＬＳＴＥが求められたとしても、それをもとにして、歯車自体および歯車周辺の各部品（軸、軸受けなど）の精度、歯車および各部品の組付精度、並びに、負荷による歪みといった様々な要素による影響を推定できる手法は明らかにされていない。このため、ＬＳＴＥが求められただけでは、歯車装置の振動騒音に影響を与える歯車の特徴量など、歯車の詳細な解析を行うことはできない。
【００１２】
本発明は、多種多様な歯車装置を試験対象とすることができ、容易な計測によって歯車のかみ合い状態を定量的に推定できる歯車装置のかみ合い状態推定方法およびそのプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１３】
本発明に係る歯車装置のかみ合い状態推定方法は、
歯車装置の振動を計測する計測ステップと、
前記計測された振動の情報から前記歯車装置の周波数応答関数を取得するＦＲＦ取得ステップと、
前記取得された周波数応答関数と前記計測された振動の情報とから前記歯車装置の負荷条件下の静的伝達誤差を計算するＬＳＴＥ計算ステップと、
を有し、
前記ＦＲＦ取得ステップは、
前記計測された振動の情報からかみ合い成分の複数の振動応答の特性曲線を抽出する抽出ステップと、
前記歯車装置をモデル化し、部分構造合成法により前記モデル化された前記歯車装置の数値シミュレーションを行って模擬周波数応答関数を求めるシミュレーションステップと、
前記模擬周波数応答関数の関数形状が、前記かみ合い成分の複数の振動応答の特性曲線に合うように、実数値ＧＡにより前記モデル化された歯車装置のパラメータの最適化を行う最適化ステップと、
前記最適化されたパラメータが適用された前記数値シミュレーションによって得られた模擬周波数応答関数と重なるように、前記かみ合い成分の複数の振動応答の特性曲線を、各絶対値を変更し且つ連続的に接続して前記周波数応答関数とするＦＲＦ推定ステップとを有する方法とした。
【００１４】
本発明の係るプログラムは、
コンピュータに、
入力された歯車装置の振動の情報から前記歯車装置の周波数応答関数を取得するＦＲＦ取得手順と、
前記取得された周波数応答関数と前記入力された振動の情報とから前記歯車装置の負荷条件下の静的伝達誤差を計算するＬＳＴＥ計算手順と、
を実行させるためのプログラムであって、
前記ＦＲＦ取得手順は、
前記入力された振動の情報からかみ合い成分の複数の振動応答の特性曲線を抽出する抽出手順と、
前記歯車装置をモデル化し、部分構造合成法により前記モデル化された前記歯車装置の数値シミュレーションを行って模擬周波数応答関数を求めるシミュレーション手順と、
前記模擬周波数応答関数の関数形状が、前記かみ合い成分の複数の振動応答の特性曲線に合うように、実数値ＧＡにより前記モデル化された歯車装置のパラメータの最適化を行う最適化手順と、
前記最適化されたパラメータが適用された前記数値シミュレーションによって得られた模擬周波数応答関数と重なるように、前記かみ合い成分の複数の振動応答の特性曲線を、各絶対値を変更し且つ連続的に接続して前記周波数応答関数とするＦＲＦ推定手順とを有する。
【発明の効果】
【００１５】
本発明によれば、多種多様な歯車装置を試験対象とすることができ、容易な計測によって歯車のかみ合い状態を定量的に推定することができる。
【図面の簡単な説明】
【００１６】
【図１】本実施形態のＬＳＴＥの推定の仕組みを示す説明図である。
【図２】本実施形態の歯車対のかみ合い状態の推定の仕組みを示す説明図である。
【図３】試験対象の歯車装置と本実施形態に係るかみ合い状態推定装置の一例を示す一部破断の側面図である。
【図４】かみ合い状態推定装置の計測系を示す図である。
【図５】計測ステップで得られる計測結果の一例であり、図５（ａ）は軸のねじり方向の振動加速度の計測結果を、図５（ｂ）は入力軸の回転検知の結果を、それぞれ示す。
【図６】振動加速度の周波数スペクトルを示すグラフであり、図６（ａ）は、歯車の回転速度を一定にしたときの振動加速度の周波数スペクトルを、図６（ｂ）は、回転速度を連続的に変化させたとききの振動加速度の周波数スペクトルを、それぞれ示す。
【図７】同期平均の仕組みを示す説明図であり、図７（ａ）は計測結果の信号波形を、図７（ｂ）は同期平均を、図７（ｃ）は同期平均のフーリエ変換を、それぞれ示す。
【図８】計測された振動加速度の周波数スペクトルの一例を示すグラフであり、図８（ａ）は同期平均前の周波数スペクトルを、図８（ｂ）は入力軸（高速軸）の同期平均の周波数スペクトルを、図８（ｃ）は中間軸の同期平均の周波数スペクトルを、それぞれ示す。
【図９】振動加速度のかみ合い成分を示すグラフであり、図９（ａ）は或る回転速度の振動加速度に含まれるかみ合い成分を、図９（ｂ）は振動加速度の測定結果から抽出されたかみ合い成分の振動応答の特性曲線を、それぞれ示す。
【図１０】ＦＲＦの数値シミュレーションに適用される歯車装置のモデル化の一例を示す説明図である。
【図１１】歯車対のモデル化の一例を示す説明図である。
【図１２】数値シミュレーションにより計算されたＦＲＦの一例を示すグラフである。
【図１３】最適化されたＦＲＦの数値シミュレーションの結果を示すグラフである。
【図１４】ＦＲＦ推定ステップの処理を示す説明図である。
【図１５】推定されたＬＳＴＥを示すグラフである。
【図１６】理論歯面からの実歯面形状の偏差の表現を示す図であり、図１６（ａ）−（ｃ）は第１要素から第３要素を、図１６（ｄ）は第１要素から第３要素の総合を示す。
【図１７】ＬＴＣＡにより計算されたＬＳＴＥと振動加速度の計測結果から推定されたＬＳＴＥとを比較するグラフである。
【図１８】理論歯面からの実歯面形状の偏差の推定結果を示すグラフである。
【図１９】歯面の接触面圧および歯当たりの推定結果を示すグラフである。
【図２０】歯当たりの観察結果を示す図である。
【図２１】計算装置のＣＰＵ３６が実行するかみ合い状態推定処理の手順を示すフローチャートである。
【図２２】図２１のＦＲＦ取得処理の詳細を示すフローチャートである。
【発明を実施するための形態】
【００１７】
以下、本発明の実施形態について図面を参照して詳細に説明する。
【００１８】
本発明の実施形態に係る歯車装置のかみ合い状態推定方法は、歯車装置の実際の運転時における一対の歯車のかみ合い状態を定量的に推定および解析するものである。このかみ合い状態推定方法は、産業用の減速機など多種多様な歯車装置を試験対象とし、容易な計測により、歯車のかみ合い状態を推定および解析できる。ここで、かみ合い状態とは、歯車装置の実際の運転時における、歯車の理論歯面からの実歯面形状の偏差、歯当たり、歯面の面圧分布などを意味する。
【００１９】
先ず、概要を説明する。実施形態のかみ合い状態推定方法は、歯車装置の振動計測の結果からＬＳＴＥを推定する第１段階と、推定したＬＳＴＥからかみ合い状態を推定する第２段階とに分けられる。以下では、歯車装置として減速機を試験する場合について説明する。
【００２０】
＜第１段階（ＬＳＴＥの推定方法）の概要＞
図１は、本実施形態のＬＳＴＥの推定の仕組みを示す。図中、ΔＴＥは減速機のＬＳＴＥ、Ｈは減速機の周波数応答関数（frequency response function、以下「ＦＲＦ」と記す）、Ｘは減速機の振動を示す。図１の上段に示されるように、減速機の振動Ｘは、主にＬＳＴＥとＦＲＦとによって発生する。従って、本実施形態のＬＳＴＥの推定方法は、図１の下段に示されるように、計測した振動をＦＲＦで除算することでＬＳＴＥを推定する。
【００２１】
推定に用いるＦＲＦは、歯車装置を振動モデルにより表わし、数値シミュレーションを行って取得される。このＦＲＦは、振動の計測結果と適合するように、数値シミュレーションのパラメータが最適化されて求められる。パラメータの最適化を行う際、数値シミュレーションにより得られたＦＲＦと、振動の計測結果から抽出された振動応答の曲線との一致の度合いが目的関数にされる。しかし、この場合、目的関数は多峰性を持つ。このため、本実施形態では、実数値ＧＡを用いてパラメータの最適化が達成される。
【００２２】
本実施形態のＦＲＦの推定方法（ｃ）と、その他の手法（ａ）、（ｂ）との、利点および欠点を表１に示す。表１のように、本実施形態のＦＲＦの推定方法（ｃ）は、ＦＲＦを定量的に且つ高い精度で推定できる。
【表１】
2018025462.tif 000003
【００２３】
＜第２段階（かみ合い状態の推定方法）の概要＞
図２は、本実施形態の歯車対のかみ合い状態の推定の仕組みを示す。第２段階のかみ合い状態の推定方法では、第一段階で推定されたＬＳＴＥと、負荷条件下の歯面接触解析（loaded tooth contact analysis、以下「ＬＴＣＡ」と記す）とが利用される。
【００２４】
ＬＴＣＡでは、歯車のかみ合う歯をモデル化し、歯車の歯面形状の理論値に偏差を加えてＬＳＴＥを計算することができる。
【００２５】
本実施形態のかみ合い状態の推定方法は、図２に示すように、実歯面形状の偏差を調整しながらＬＴＣＡを繰り返し計算し、且つ、第１段階で推定されたＬＳＴＥとＬＴＣＡで計算されたＬＳＴＥとを比較する。そして、両方のＬＳＴＥが最も適合するように歯車の実歯面形状が見つけられる。両方のＬＳＴＥが最も適合する実歯面形状を探すため、本実施形態では、第一段階と同じく実数値ＧＡを用いて、モデル化された歯車装置のパラメータの最適化が行われる。
【００２６】
このように求められた実歯面形状は、歯車自体および歯車周辺の各部品の精度、歯車および各部品の組付精度、並びに、負荷による歪みといった各要素による影響が加わったときの歯車の歯面形状を表している。また、このときのＬＴＣＡの結果から、歯車装置の運転時における歯車の歯当たり、歯面の面圧分布等のかみ合い状態がわかる。
【００２７】
続いて、上記第１段階および第２段階の各処理の詳細を説明する。第１段階の処理は、計測ステップと、同期平均ステップと、ＦＲＦ取得ステップと、ＬＳＴＥ計算ステップとを含む。さらに、ＦＲＦ取得ステップは、振動応答のかみ合い成分の取得ステップと、ＦＲＦの数値シミューションステップと、数値シミュレーションの最適化ステップと、ＦＲＦ推定ステップとを含む。
【００２８】
＜試験対象の歯車装置と本実施形態に係るかみ合い状態推定装置の構成＞
図３は、試験対象の歯車装置と本実施形態に係るかみ合い状態推定装置の一例を示す。図４は、かみ合い状態推定装置の計測系を示す図である。
【００２９】
試験対象の歯車装置は、歯車（はすば歯車）Ｇ１、Ｇ２を有する２段減速の減速機１０である。本実施形態に係るかみ合い状態推定装置は、減速機１０の入力軸１１に駆動力を加えるモータ２１と、減速機１０の出力軸１２に負荷を加える抵抗器２２と、計測系の構成を有する。計測系の構成は、図４に示すように、入力軸１１の１回転を検知する回転検知器（例えば光電式検出器）３１と、出力軸１２のねじり方向の加速度を計測する加速度センサ３２と、加速度センサ３２の出力を外部に取り出すスリップリング３３とを備える。また、計測系の構成は、回転検知器３１の出力および加速度センサ３２の出力を記録するデータロガー３５と、記録されたデータを解析する計算装置３６とを備える。
【００３０】
計算装置３６は、コンピュータであり、ＣＰＵ（Central Processing Unit）３６ａと、本実施の形態に係るかみ合い状態推定プログラム４１を記録した記憶装置３６ｂと、作業用のメモリ３６ｃとを備える。また、計算装置３６は、データロガー３５から計測された振動情報を入力するインタフェース３６ｄと、利用者から各種パラメータ等の入力を行う入力装置３６ｅと、利用者へ画像出力或いは印刷出力を行う出力装置３６ｆとを備える。
【００３１】
＜計測ステップ＞
計測ステップでは、加速度センサ３２より出力軸１２のねじり振動の振動加速度が計測され、回転検知器３１により入力軸１１の回転が計測される。振動加速度を計測する部位は、歯車Ｇ１、Ｇ２を収容して回転軸を支持する歯車箱、又は他の部位を採用してもよいが、ねじり振動を計測することで、歯車の励振による振動を高い精度で計測できる。計測時の運転条件としては、出力軸１２のトルクを一定とし、入力軸１１の回転速度を広く変化するという条件が採用される。
【００３２】
図５（ａ）、（ｂ）は、計測ステップで得られる計測結果の一部を示す。図５（ａ）のグラフには入力軸１１の回転速度を一定としたときのねじり振動の振動加速度が表れている。図５（ｂ）のグラフには入力軸１１の１回転ごとに１パルスの回転トリガ信号が現れている。入力軸１１の回転速度を広く変化する過程、回転速度が各値に固定された状態で、図５（ａ）、（ｂ）のような計測結果が得られる。
【００３３】
図６（ａ）、（ｂ）は、計測結果から得られる振動加速度の周波数スペクトルを示すグラフである。図６（ａ）は、入力軸１１の回転速度が一定のときに得られた振動加速度の周波数スペクトルである。図６（ｂ）は、縦軸に入力軸１１の回転速度、横軸に周波数、第３軸に振動加速度をとった三次元グラフであり、何れか１つの回転速度に固定したときの信号を抽出したものが図６（ａ）のグラフ線を表す。入力軸１１の回転速度を広く変化する過程で得られた全ての計測結果を集めると、図６（ｂ）のようなグラフになる。
【００３４】
＜同期平均ステップ＞
図７（ａ）−（ｃ）は、同期平均処理の仕組みを示す。
【００３５】
同期ステップは、計測された振動加速度の時間波形の同期平均を行う。時間波形とは時間を横軸にとった波形を意味する。同期平均は、時間領域での平均化処理であり、着目する歯車の回転周期で切り分けられた複数の時間波形を、始端を合せて平均する処理である。
【００３６】
例えば、入力軸に固定される歯車に着目した場合、計測されたトリガ信号の周期が歯車の回転周期を表す。この場合、図７（ａ）に示すように、トリガ信号の周期で複数の振動加速度の時間波形が切り分けられ、図７（ｂ）に示すように、複数の時間波形が始端を合せて平均化される。これにより、着目する歯車の振動加速度の信号が得られる。この信号からは、着目する歯車に非同期の振動成分が消去されている。また、図７（ｃ）に示すように、時間波形をフーリエ変換すれば着目する歯車の振動加速度の周波数スペクトルが得られる。
【００３７】
同様に、出力軸に固定される歯車、或いは、中間軸に固定される歯車に着目する場合には、トリガ信号と歯車の歯数とから、これらの歯車の回転周期を求めることができる。従って、求められた回転周期で切り分けられた振動加速度の時間波形が平均化されることで、着目する各歯車の振動加速度の信号が得られる。
【００３８】
図８（ａ）は同期平均前の振動加速度の周波数スペクトルの一例を示す。図８（ｂ）、（ｃ）は入力軸（高速軸）の同期平均の周波数スペクトルと中間軸の同期平均の周波数スペクトルとをそれぞれ示す。
【００３９】
図８（ｂ）、（ｃ）に示すように、同期平均によって、振動加速度の時間波形から着目する歯車の信号を分離して取り出すことができる。グラフ上に矢印で示したように、同期平均の結果には、個々の歯車の振動を表す信号が抽出されている。
【００４０】
以下では、入力軸（高速軸）の同期平均を行った振動加速度の信号を用いて、入力軸に固定された歯車を分析する場合について説明する。なお、別の歯車を分析する場合には、対象の歯車に対応する同期平均された振動加速度の信号を用いて、同様の処理を行えばよい。
【００４１】
＜振動応答のかみ合い成分の取得ステップ＞
図９（ａ）、（ｂ）は、振動加速度のかみ合い成分を示すグラフである。図９（ａ）は、或る回転速度のときの振動加速度に含まれるかみ合い成分を示す。図９（ｂ）は、歯車のかみ合い成分の振動応答の特性波形を示す。
【００４２】
図９（ａ）に示すように、或る回転速度のときの振動加速度の周波数スペクトルが得られると、かみ合い周波数ｆｚおよびその高調波成分２ｆｚ、３ｆｚの加速度値が求められる。かみ合い周波数とは、或る一対の歯のかみ合いの開始から次の一対の歯のかみ合いの開始までを１周期とした周波数である。かみ合い周波数ｆｚとその高調波成分２ｆｚ、３ｆｚ・・・と含めてかみ合い成分と呼ぶ。
【００４３】
歯車の回転速度が変わると、かみ合い成分の周波数も連続的に変化する。よって、回転速度を連続的に変化させて計測された複数の振動加速度の信号波形の各々から、かみ合い成分の加速度値を抽出すると、図９（ｂ）に示すように、各かみ合い成分に対応する振動応答の特性曲線が求められる。
【００４４】
なお、図６（ｂ）に示したように、振動加速度の計測時には入力軸の回転速度に上限値と下限値とが生じる。従って、図９（ｂ）のように振動応答の特性曲線が得られる周波数範囲は、かみ合い成分ごとに異なる上限値と下限値とが生じる。例えば、かみ合い周波数ｆｚの特性曲線は約５０Ｈｚ−５３０Ｈｚの範囲で得られ、高調波成分２ｆｚの特性曲線は約１００Ｈｚ−１０６０Ｈｚの範囲で得られ、高調波成分３ｆｚの特性曲線は約１５０Ｈｚ−１６００Ｈｚの範囲で得られる。これら３つの特性曲線は、単独では、横軸全ての周波数範囲をカバーできないが、３つの特性曲線を合せると、停止状態に近い低い周波数範囲を除いて、横軸全ての周波数範囲をカバーできる。
【００４５】
＜ＦＲＦの数値シミュレーション＞
図１０は、ＦＲＦの数値シミュレーションに適用される歯車装置のモデル化の一例を示す。
【００４６】
本実施形態では、ＦＲＦを推定するために、先ず、数値シミュレーションを行ってＦＲＦ（本発明に係る模擬周波数応答関数に相当）を計算する。数値シミュレーションでは、歯車装置を線形時不変系で表し、歯車対のかみ合い部にＬＳＴＥによる強制変位が作用するモデルが採用される。計算方法は部分構造合成法が用いられる。
【００４７】
歯車装置が定数係数系のモデルで表され、歯車対のかみ合い部に調和励振力による強制振動が作用すると見なすと、運動方程式は式（１）で表すことができる。
【００４８】
【数１】
2018025462.tif 000004
ここで、Ｍは質量行列、Ｃは減衰行列、Ｋは剛性行列、ωは角周波数を表す。
【００４９】
さらに、変位ベクトルＸと起振力ベクトルＥについて複素数を用いて式（２）のように表し、式（１）に代入すれば式（３）が得られる。式（３）を解くことでＦＲＦが求められる。
【００５０】
質量行列Ｍと剛性行列Ｋを作成する方法は、部分構造合成法の一つである特性行列合成法（文献［１］長松昭男、大熊政明、"部分構造合成法"、培風館、1991年、pp.76-78, 155-156を参照）を用いる。減速機１０を含む振動系を全系として、これを回転軸系（入力軸１１、出力軸１２など）と歯車箱１３という二つの分系に分け、歯車箱１３の軸受支持部を二つの分系の結合部とする。分系それぞれについて特性行列を作成すれば、有限要素法における個々の有限要素の特性行列を統合するときの手順と同様にして、全系の特性行列を作成することができる。
【００５１】
回転軸系（入力軸１１、出力軸１２、中間軸など）は、一つの節点につき並進と回転の６自由度を考慮したばね質量モデルとする。軸受の剛性は各軸方向のばね、軸の剛性はねじり、伸縮、曲げに対応するばねに置き換える。
【００５２】
図１１は、歯車対のモデル化の一例を示す。歯車対としては、各歯車の中心とピッチ点に配置される二つの節点ｎ０、ｎ１が剛体の要素ｅ１で接続され、各歯車本体の質量が中心の節点ｎ１に与えられたモデルを採用できる。また、ピッチ点には駆動歯車Ｇ１と被動歯車Ｇ２それぞれに対応する二つの節点ｎ０１、ｎ０２（図１１（ｂ）を参照）があり、これらはかみ合う歯のばねで接続されたモデルとする。このようなモデルでは、歯のばねは、作用平面Ｓ１上の歯の直角方向のかみ合い剛性と、作用平面Ｓ１の垂直方向を軸とする回転の剛性をもつ。また、一対の歯車Ｇ１、Ｇ２のＬＳＴＥは、歯のばねの両端に作用する符号が反対で大きさが等しい二つの調和励振力Ｆ１、Ｆ２により表現される。
【００５３】
歯車箱１３は有限要素法を用いてモデル化する。ただし有限要素モデルの特性行列はそのままでは大き過ぎるので、特性行列を再構成して、自由度を回転軸系との結合部である軸受支持部のみに減少させる。再構成の方法は、固有値解析を行って固有モードを求め、結合部の小自由度特性行列を構成する方法（文献［１］を参照）を用いる。
【００５４】
減衰行列Ｃは、歯車装置１０に働く減衰を比例ヒステリシス減衰（文献［２］長松昭男、"モード解析"、培風館、1985年、pp.86-88を参照）と考えて、式（４）のように表す。
【００５５】
【数２】
2018025462.tif 000005
ここで、Ｄはヒステリシス減衰行列、ａとｂは比例定数を表す。
【００５６】
図１２は、数値シミュレーションにより計算されたＦＲＦの一例を示す。上述の数値シミュレーションにより計算されたＦＲＦは、ほぼ停止状態の周波数範囲をカバーする一方、精度が低くなるため高い周波数範囲をカバーできない。また、数値シミュレーションによるＦＲＦの関数形状は、減速機のモデル化のパラメータによって変化し、すぐには振動応答の特性曲線と合致しない。なお、図１２では、かみ合い成分の振動応答の特性曲線を重ねて示している。
【００５７】
＜数値シミュレーションの最適化ステップ＞
詳細は後述するが、後段のＦＲＦの推定ステップにおいては、振動加速度から得られるかみ合い成分の振動応答が、数値シミュレーションにより得られたＦＲＦと重ねられて１つの曲線に合成される。この合成の際、数値シミュレーションによるＦＲＦは振動加速度から得られるかみ合い成分の振動応答の曲線形状と一致することが求められる。従って、数値シミュレーションの正確さは、かみ合い起振力の推定結果に大きく影響する。
【００５８】
そこで、最適化ステップでは、数値シミュレーションによるＦＲＦが、かみ合い成分の振動応答に最も適合するように、モデル化された減速機１０のパラメータの最適化が行われる。パラメータの最適化は、一対の歯車Ｇ１、Ｇ２のかみ合い剛性、軸受、および、カップリングの並進と回転の６つの自由度の各方向における剛性およびヒステリシス減衰の比例定数の各パラメータを対象に行われる。質量行列は減速機の図面等から十分な精度で作成することができるので最適化の対象にしない。以上により最適化するパラメータの数は３２個である。
【００５９】
最適化問題において、目的関数が多峰性を持つ場合に必要となる大域的探索の効果的な手法として遺伝的アルゴリズム（genetic algorithm：ＧＡ）が知られている。さらに、連続変数を最適化する問題には、遺伝子表現に実数値を用いる実数値ＧＡが変数の連続性を考慮している点で有効とされる（文献［５］小林重信、"実数値GAのフロンティア"、人工知能学会論文誌、2009年、Vol.24, No.1、pp.147-162と、文献［６］金田行雄、笹井理生、古橋武、"計算科学講座 第9巻 超多自由度系の最適化"、共立出版、2013年、pp.1-117を参照）。
【００６０】
実数値ＧＡにおける一つの個体は、最適化するパラメータの数と同数の自由度ｎをもつ実数値ベクトルξで表される。この個体に対応する振動モデル（モデル化した減速機１０）において、最適化するパラメータｎ個の中でのｉ番目の値Ｑｉは、式（５）を用いて求める。
【００６１】
【数３】
2018025462.tif 000006
ここで、Ｐｉはｉ番目のパラメータの基準値、ξｉは個体のベクトルξのｉ番目の要素を表す。個体のベクトルξの初期値は、各部の剛性を対象としたとき−１〜０．５の幅で、ヒステリシス減衰の係数を対象としたとき−５〜−１の幅で、一様乱数により決定する。
【００６２】
個体の集団には分割集団モデルを用いる。この方法は分散遺伝的アルゴリズム（distributed genetic algorithm：ＤＧＡ）と呼ばれ、ＧＡにおける個体の集団を複数の分割集団に分割し、各分割集団内で独立してＧＡを行う。ＤＧＡは、通常のＧＡに比べて計算時間の短縮や高品質な解が得られるとされている（文献［７］廣安知之、三木光範、上浦二郎、"実験計画法を用いた分散遺伝的アルゴリズムのパラメータ推定"、情報処理学会論文誌、数理モデル化と応用、2002年、Vol.43、No.SIG10(TOM7)と、文献［６］を参照）。ここでの１分割集団あたりの個体数は６４、分割集団の数は４である。
【００６３】
ある世代の集団から次の世代の集団を生成する手順には、世代交代モデルＪＧＧ（文献［５］を参照）を用いる。ＪＧＧでは、まず集団からｎｐ個の個体をランダムに非復元抽出して、交叉に用いる親個体群とする。次に親個体群に交叉を繰り返し適用して、ｎｃ個の子個体を生成する。次に子個体群から評価が上位のｎｐ個の個体を選び、集団に加えることで次の世代集団とする。親個体群から子個体群を生成する交叉の方法には、ＲＥＸｓｔａｒを用いる（文献［５］を参照）。ＲＥＸｓｔａｒによる交叉では、親個体群の重心を基準点として、そこからマクロにみて評価が高くなる方向に移動した点が中心となり、一様分布をもとに計算される子個体生成分布に従って子個体がランダムに生成される。ここでの集団から抽出される親個体の数ｎｐは１分割集団あたりｎ＋１である。ここでの子個体の生成数ｎｃは１分割集団あたり６４である。
【００６４】
最適化ステップは、生成された子個体を用いて数値シミュレーションを行い、得られたＦＲＦをかみ合い成分の振動応答と比較することで、個体の評価を行う。評価には、次のφ１とφ２との二つの目的関数を用いることができる。
【数４】
2018025462.tif 000007
【００６５】
式（６）のＥハットは、与えられたＸハットとＨに対してφ１が最小になるように値を決める。φ１はＸハットからのＨのずれを複素平面上の距離を用いて表したものである。もう一方のφ２は、ＸハットとＨそれぞれの周波数による変化を表すＸハット’とＨ’のずれを複素平面上の距離を用いて表したものである。式（６）と式（７）の波括弧内の第２項は、｜Ｈ｜≪｜Ｘハット｜となるＨに対して目的関数の値がほぼ同じになり適切に評価できなくなることを防ぐために追加されている。ＨとＸハットが完全に一致するときφ１とφ２の両方が零になるので、φ１とφ２が小さい個体について高い評価が与えられる。
【００６６】
二つの目的関数を用いる多目的最適化のアルゴリズムにはＮＳＧＡ−ＩＩ（文献［８］Deb, K., Pratap, A., Agarwal, S. and Meyarivan, T., "A fast and elitist multiobjective genetic algorithm: NSGA-II"、IEEE Transactions on Evolutionary Computation, 2002, Vol.6, No.2 (2002), pp.182-197を参照）を利用できる。ＮＳＧＡ−ＩＩはパレート優越関係と多様性維持メカニズムを用いて個体評価を行う方法であり、ＧＡによる多目的最適化の分野で幅広く用いられている（文献［６］を参照）。パレート優越関係に基づく個体評価では、個体群の中で他の解に優越されない非劣解に高い評価を与えて１から順番にランクを割り当てて、このランクを個体評価における第一基準とする。多様性維持メカニズムによる個体評価では、同じランクが与えられた他の個体と比較して目的関数の値が離れている個体ほど高い評価を与えて、これを個体評価の第二基準とする。
【００６７】
このアルゴリズムは、評価値が上位の子個体を集団に加えた後に、各分割集団内で選択した複数個の個体を別の分割集団へ移動させる移住操作を行う。どの分割集団からどの分割集団へ移住するかの選択、および移住する個体の選択はランダムに決定する（文献［７］を参照）。
【００６８】
以上の手順で世代交代を繰り返し行った後、各分割集団からパレート優越関係でランク１の個体を取り出し、取り出された個体ベクトルの集団の重心からの距離が最も小さい個体が、最適化の結果として採用される。ここでの世代交代の回数は５００回である。
【００６９】
図１３は、パラメータの最適化後のＦＲＦの数値シミュレーションの結果を示す。パラメータの最適化を行うことで、数値シミュレーションによるＦＲＦが、振動加速度から求められるかみ合い成分の振動応答の特性曲線によく一致することが分かる。これにより、ＦＲＦの推定の精度を高めることができる。
【００７０】
＜ＦＲＦ推定ステップ＞
図１４（ａ）−（ｄ）は、ＦＲＦ推定ステップの処理を示す説明図である。
ＦＲＦ推定ステップは、かみ合い成分の複数の振動応答の特性曲線を、縦軸方向に変位する調整を行って、数値シミュレーションにより計算されたＦＲＦに重ね、その後、各かみ合い成分の振動応答の特性曲線が連続的に接続されるように合成する。連続的に接続できない箇所については、両方の特性曲線の平均値などを適宜採用すればよい。このような合成により、１本の曲線としてＦＲＦを求めることができる。図１４の例では、かみ合い周波数ｆｚおよび高調波成分２ｆｚ−Ｍｆｚにそれぞれ対応する多数の特性曲線を合成している。
【００７１】
このような複数の特性曲線の合成は、振動応答の振幅の特性曲線（図１４（ａ）、（ｃ）を参照）と、振動応答の位相の特性曲線と（図１４（ｂ）、（ｄ）を参照）で行うことができる。これにより、周波数ごとに振幅と位相とが１つに決まったＦＲＦが推定される。
【００７２】
数値シミュレーションで得たＦＲＦは、精度が低くなるため高い周波数範囲をカバーできないが、振動応答の特性曲線を合成して得られたＦＲＦは、高いい周波数範囲までカバーされている。
【００７３】
＜ＬＳＴＥ計算ステップ＞
続いて、ＬＳＴＥ計算ステップは、推定されたＦＲＦと、計測された振動加速度とから、対象とする歯車のＬＳＴＥの推定値を計算する。ＬＳＴＥの計算は、式（9）による。
【数５】
2018025462.tif 000008
ここで、Ｘｍは着目する歯車のｍ次の回転次数成分の振動応答、ΔＴＥｍは着目する歯車におけるＬＳＴＥのｍ次の回転次数成分、Ｈは推定されたＦＲＦ、ｆｒは着目する歯車の回転周波数、変数にハットマークは変数のフーリエ変換を示す。Ｘｍ、ΔＴＥｍは、ｍが着目する歯車の歯数の整数倍のときにかみ合い成分に関する量となる。
【００７４】
図１５（ａ）−（ｃ）は、推定されたＬＳＴＥを示すグラフを示す。これらのグラフは、歯の直角方向の値に換算されたＬＳＴＥを示している。また、これらのグラフは、９歯の歯車のＬＳＴＥを示している。
【００７５】
式（９）に従えば、Ｚを着目する歯車の歯数として、歯数の整数倍の回転次数ｍ＝Ｚ，２Ｚ，…ＭＺ（Ｍは整数）のΔＴＥｍをそれぞれ計算し、ΔＴＥｍのｍについての総和を計算することで、主に着目する歯車のかみ合い成分の特徴を表わすＬＳＴＥを求めることができる。図１５（ｂ）は、着目する歯車について、かみ合い成分の１次から７次の回転次数ｍを使用して求めたＬＳＴＥを示している。このようなＬＳＴＥには、歯車装置の振動騒音に大きく影響する歯車の特徴が表れている。着目する歯車について、歯数の整数倍の回転次数ｍ＝Ｚ，２Ｚ，…ＭＺ（Ｍは整数）の振動応答のみを使用し、非整数次成分を除外して求められたＬＳＴＥは、図１５（ｂ）に示されるように、着目する歯車の各歯共通の特徴を表す。
【００７６】
また、周波数成分を示すｍとして、かみ合い成分の両側に現れる側帯波の周波数成分を含めてＬＳＴＥを計算することもできる。図１５（ａ）は、着目する歯車について、回転次数ｍとして、かみ合い成分の１〜７次と、側帯波の正負１〜４次を含めてＬＳＴＥを計算した結果を示している。図１５（ａ）に示すように、着目する歯車について、非整数次の成分を含めて計算されたＬＳＴＥは、着目する歯車の個々の歯ごとの特徴も表す。
【００７７】
図１５（ｃ）は、図１５（ａ）のＬＳＴＥを歯車の回転次数を横軸にとった周波数スペクトルとして表したグラフである。グラフの横軸の値は、歯車の回転周波数を乗算することで周波数に換算される。このように、ＬＳＴＥを周波数スペクトルとして表わすことで、歯車の振動騒音の周波数とＬＳＴＥの大きさとの関係を確認することができる。
【００７８】
なお、ＦＲＦは歯車装置の振動の周波数特性である。よって、伝達トルク又は歯車の合成誤差が多少変化してもＦＲＦは変わらない。従って、ＦＲＦが既に求められている場合には、ＦＲＦを推定する手順を省略して、ＬＳＴＥの推定を行ってもよい。
【００７９】
続いて、第２段階の処理であるかみ合い状態の推定ステップについて詳細に説明する。
【００８０】
＜かみ合い状態の推定ステップ＞
かみ合い状態の推定ステップは、これまでの手順によって推定されたＬＳＴＥと、ＬＴＣＡ（負荷条件下の歯面接触解析：loaded tooth contact analysis）によって計算されるＬＳＴＥとを比較する。そして、両方のＬＳＴＥが一致するような実歯面形状が探索される。このときの実歯面形状は、歯車自体および歯車周辺の各部品の精度、歯車および各部品の組付精度、ならびに、負荷による歪みといった各要素による影響が加わったときの歯車の歯面形状を表している。また、このときのＬＴＣＡの結果から、歯当たり、歯面の面圧分布等のかみ合い状態が分かる。
【００８１】
続いて、ＬＴＣＡによる計算の一例を示す。
【００８２】
はすば歯車のかみ合う歯対ｉにおける接触線上の分布荷重は、次の連立積分方程式を解くことによって求められる（文献［３］久保愛三、梅澤清彦、"誤差をもつ円筒歯車の荷重伝達特性に関する研究（第１報，基礎的考察）"、日本機械学會論文集、1977年、Vol.43, No.371、pp.2771-2783と、文献［４］松村茂樹、梅澤清彦、北條春夫、"誤差を考慮したはすば歯車の軽負荷における回転方向振動解析（第２報，回転方向振動シミュレーションの開発）"、日本機械学会論文集Ｃ編、1994年、Vol.60, No.575、pp.2436-2443を参照）。
【００８３】
【数６】
2018025462.tif 000009
ここで、ｘは接触線上の変形観測点の座標、ξは接触線上の荷重作用点の座標、Ｐｉ（ξ）は接触線上での分布荷重、Ｋｂ（ｘ、ξ）は歯の曲げせん断たわみの影響関数、Ｋｃ（ｘ、ξ）は歯の接触による近寄りの影響関数、Δは被動歯車の駆動歯車からの回転角の遅れを作用線上の変位に換算した値、ｅｉ（ｘ）は歯対ｉの接触線上における合成誤差（理論インボリュートヘリコイドを０として歯が凹になる方向が正）、βｂは基礎円筒ねじれ角である。
【００８４】
歯対ｉについて、式（１０）を用いて作用線上の各かみ合い位置で被動歯車の遅れΔと分布荷重Ｐｉ（ξ）の関係を求めれば、歯対ｉに関する伝達荷重Ｗｉは式（１１）によって表され、歯車対の伝達荷重Ｗは式（１２）によって表される。
【００８５】
【数７】
2018025462.tif 000010
【００８６】
式（９）〜（１２）により求めた各かみ合い位置におけるΔとＷの関係をもとに、複数の歯対の同時かみ合いを考慮して、伝達荷重Ｗを一定として各かみ合い位置でのΔを計算すれば、ＬＳＴＥが求められる。
【００８７】
ＬＴＣＡは、このような計算原理に基づいてＬＳＴＥを計算する。
【００８８】
さらに、本実施形態では、ＬＴＣＡにより計算されるＬＳＴＥが、振動加速度の計測によって推定されたＬＳＴＥに最も適合するように、ＬＴＣＡの計算条件に関するパラメータの最適化が行われる。ここでのパラメータの最適化は、ＬＴＣＡに入力する実歯面形状の偏差を対象として行われる。
【００８９】
図１６（ａ）−（ｄ）は、理論歯面からの実歯面形状の偏差の表現を示す。図１６（ａ）の歯幅の位置毎の歯形形状は、歯幅方向の複数の位置に５本の歯形形状を与える。各歯形形状は、１６点の座標における、偏差の値により表され、パラメータの数は８０個である。図１６（ｂ）の共通の歯形形状は作用線方向の座標の２次関数によって表され、１次と２次の項の係数を与えるので、パラメータの数は２個である。図１６（ｃ）の共通の歯すじ形状は歯幅方向の座標の４次関数によって表され、１次から４次までの各項の係数を与えるので、パラメータの数は４個である。図１６（ｄ）の理論歯面からの実歯面形状の偏差は、図１６（ａ）から（ｃ）までを重ね合わせて合成することで得られる。従って、最適化するパラメータの数は８６である。
【００９０】
パラメータの最適化計算は、ＦＲＦ数値シミュレーションと同様に、実数値ＧＡを用いて行われる。実数値ＧＡにおいて、一つの個体は、最適化するパラメータの数と同数の自由度ｎをもつ実数値ベクトルξで表される。この個体に対応する計算モデルにおいて、最適化するパラメータｎ個の中でのｉ番目の値Ｑｉは、式（１３）を用いて求められる。
【数８】
2018025462.tif 000011
ここで、Ｐｉはｉ番目のパラメータの基準値、ξｉは個体のベクトルξのｉ番目の要素を表す。個体のベクトルξの初期値は、−１〜１の幅で、一様乱数により決定する。
【００９１】
個体の集団には分割集団モデルを用いる。この方法は分散遺伝的アルゴリズム（distributed genetic algorithm：ＤＧＡ）と呼ばれ、ＧＡにおける個体の集団を複数の分割集団に分割し、各分割集団内で独立してＧＡを行う。ＤＧＡは、通常のＧＡに比べて計算時間の短縮や高品質な解が得られるとされている（文献［７］、［６］を参照）。ここでの１分割集団あたりの個体数は１７２、分割集団の数は４である。
【００９２】
ある世代の集団から次の世代の集団を生成する手順には、世代交代モデルＪＧＧ（文献［５］を参照）を用いる。ＪＧＧでは、まず集団からｎｐ個の個体をランダムに非復元抽出して、交叉に用いる親個体群とする。次に親個体群に交叉を繰り返し適用して、ｎｃ個の子個体を生成する。次に子個体群から評価が上位のｎｐ個の個体を選び、集団に加えることで次の世代集団とする。親個体群から子個体群を生成する交叉の方法には、ＲＥＸｓｔａｒを用いる（文献［５］を参照）。ＲＥＸｓｔａｒによる交叉では、親個体群の重心を基準点として、そこからマクロにみて評価が高くなる方向に移動した点が中心となり、一様分布をもとに計算される子個体生成分布に従って子個体がランダムに生成される。ここでの集団から抽出される親個体の数ｎｐは１分割集団あたりｎ＋１である。ここでの子個体の生成数ｎｃは１分割集団あたり１７２である。
【００９３】
ここでの最適化の処理では、生成された個体による形状偏差を用いてＬＴＣＡによる計算を行い、得られたＬＳＴＥを振動加速度の計測により推定されたＬＴＣＡと比較することで、個体の評価が行われる。評価には、次のφ１とφ２との目的関数を用いることができる。
【００９４】
【数９】
2018025462.tif 000012
ここで、ｔｉは時間、ΔＡｊはｊ番目のトルク条件におけるＬＳＴＥのＬＴＣＡ計算結果、ΔＶｊはｊ番目のトルク条件におけるＬＳＴＥの振動計測による推定結果である。ΔＡｊとΔＶｊが完全に一致するときφ１とφ２の両方が零になるので、φ１とφ２が小さい個体について高い評価が与えられる。ここでのトルク条件の数は４水準である。
【００９５】
二つの目的関数を用いる多目的最適化のアルゴリズムにはＮＳＧＡ−ＩＩ（文献［８］を参照）を利用できる。ＮＳＧＡ−ＩＩはパレート優越関係と多様性維持メカニズムを用いて個体評価を行う方法であり、ＧＡによる多目的最適化の分野で幅広く用いられている（文献［６］を参照）。パレート優越関係に基づく個体評価では、個体群の中で他の解に優越されない非劣解に高い評価を与えて１から順番にランクを割り当てて、このランクを個体評価における第一基準とする。多様性維持メカニズムによる個体評価では、同じランクが与えられた他の個体と比較して目的関数の値が離れている個体ほど高い評価を与えて、これを個体評価の第二基準とする。
【００９６】
このアルゴリズムは、評価値が上位の子個体を集団に加えた後に、各分割集団内で選択した複数個の個体を別の分割集団へ移動させる移住操作を行う。どの分割集団からどの分割集団へ移住するかの選択、および移住する個体の選択はランダムに決定する（文献［７］を参照）。
【００９７】
以上の手順で世代交代を繰り返し行った後、各分割集団からパレート優越関係でランク１の個体を取り出し、取り出された個体ベクトルの集団の重心からの距離が最も小さい個体が、最適化の結果として採用される。ここでの世代交代の回数は５９０回である。
【００９８】
このような処理によって、実際の運転条件下における歯車の各歯面形状およびかみ合い状態が推定される。
【００９９】
＜歯車対のかみ合い状態の解析例＞
続いて、上述の方法によって推定される歯車対のかみ合い状態を示す定量的なデータを示す。
【０１００】
図１７（ａ）−（ｄ）は、ＬＴＣＡにより計算されたＬＳＴＥと、振動加速度の計測結果から推定されたＬＳＴＥとの比較を示す。このように、複数のトルクが加えられた状態で、両方のＬＳＴＥが一致しているのが分かる。
【０１０１】
図１８は、理論歯面からの実歯面形状の偏差の推定結果を示す。図１９は、接触面圧および歯当たりの推定結果を示す。図２０は、振動計測後に確認のために行った歯当たりの観察結果を示す。
【０１０２】
図１９において、接触面圧が正の値の領域が歯面が接触している歯当たりの領域である。この結果はトルク３２０Ｎｍのときのデータである。図２０の観察結果から、接触領域が歯の中心よりも左側に偏っていることから片当たりが起きていることがわかる。この傾向は図１９の推定結果とよく一致する。
【０１０３】
＜実施形態効果＞
以上のように、本実施形態の歯車装置のかみ合い状態推定方法によれば、専用の試験機等を必要とせず、一般的な振動計測機器を用いた容易な計測の結果から、減速機の実際の運転条件下における歯車のＬＳＴＥを定量的に且つ高い精度で推定することができる。従って、本実施形態によれば、様々なサイズ又は様々な軸配置を有する産業用の歯車装置に広く適用して、試験を行うことができる。
【０１０４】
さらに、本実施形態によれば、推定されたＬＳＴＥに基づいて、歯車のかみ合い状態が定量的に且つ詳細に推定される。そして、推定された歯車のかみ合い状態から、例えば、歯車自体および歯車周辺の各部品の精度、歯車および各部品の組付精度、負荷による歪みといった各要素の実歯面形状への影響を分析することができる。従って、このような分析結果を、設計工程或いは製造工程にフィードバックして、様々な問題点を明らかにしたり問題点の改善を図ることができる。
【０１０５】
また、本実施形態の歯車装置のかみ合い状態推定方法によれば、歯車のかみ合い成分の振動応答の特性曲線を用いてＦＲＦが定量的に推定される。そして、このＦＲＦを用いてＬＳＴＥ、および歯車のかみ合い状態が推定される。従って、推定される歯車のかみ合い状態には、減速機の振動騒音に大きく影響するかみ合い成分の特徴が定量的に表れる。このため、本実施形態の方法は、減速機の振動騒音の改善を図るための調査分析に特に有用である。
【０１０６】
次に、本実施形態の計算装置３６に備わるかみ合い状態推定プログラム４１について説明する。
【０１０７】
図２１は、計算装置３６のＣＰＵ３６ａが実行するかみ合い状態推定処理の手順を示すフローチャートである。図２２は、図２１のＦＲＦ取得処理の詳細を示すフローチャートである。
【０１０８】
図２１のかみ合い状態推定処理は、記憶装置３６ｂに記憶されたかみ合い状態推定プログラム４１に従って、ＣＰＵ３６ａが実行する処理であり、入力装置３６ｅを介した利用者からの要求に基づいて開始される。
【０１０９】
この処理が開始されると、図２１に示すように、ＣＰＵ３６ａは、先ず、データロガー３５から計測された振動情報を入力する処理（ステップＳ１）を実行する。次いで、ＣＰＵ３６ａは、同期平均処理（ステップＳ２）と、ＦＲＦ取得処理（ステップＳ３）と、ＬＳＴＥ計算処理（ステップＳ４）と、かみ合い状態の推定処理（ステップＳ５）とを実行する。
【０１１０】
さらに、図２２に示すように、ＦＲＦ取得処理では、ＣＰＵ３６ａは、振動応答のかみ合い成分の取得処理（ステップＳ１１）と、ＦＲＦの数値シミューション処理（ステップＳ１２）と、数値シミュレーションの最適化処理（ステップＳ１３）と、ＦＲＦ推定処理（ステップＳ１４）とを実行する。
【０１１１】
ここで、同期平均処理、ＬＳＴＥ計算処理、および、かみ合い状態の推定処理とは、前述した同期平均ステップ、ＬＳＴＥ計算ステップ、および、かみ合い状態の推定ステップの各処理を示す。また、振動応答のかみ合い成分の取得処理、ＦＲＦの数値シミューション処理、数値シミュレーションの最適化処理、および、ＦＲＦ推定処理は、前述した振動応答のかみ合い成分の取得ステップ、ＦＲＦの数値シミューション、数値シミュレーションの最適化ステップ、および、ＦＲＦ推定ステップの各処理を示す。
【０１１２】
また、ＦＲＦ取得処理は、本発明に係るＦＲＦ取得ステップならびにＦＲＦ取得手順の一例に相当する。ＬＳＴＥ計算処理は、本発明に係るＬＳＴＥ計算ステップならびにＬＳＴＥ計算手順の一例に相当する。振動応答のかみ合い成分の取得処理は、本発明に係る抽出ステップならびに抽出手順の一例に相当する。ＦＲＦの数値シミューション処理は、本発明に係るシミュレーションステップならびにシミュレーション手順の一例に相当する。数値シミュレーションの最適化処理は、本発明に係る最適化ステップならびに最適化手順の一例に相当する。ＦＲＦ推定処理は、本発明に係るＦＲＦ推定ステップならびにＦＲＦ推定手順の一例に相当する。かみ合い状態の推定処理は、本発明に係るかみ合い状態を推定するステップならびに手順の一例に相当する。
【０１１３】
各ステップＳ２〜Ｓ５の処理において設定変更可能なパラメータは、その都度、利用者から入力装置３６ｅを介して入力を要求するように処理されてもよい。或いは、予めＣＰＵ３６ａがデフォルトのデータを設定し、利用者が変更要求を入力することで、変更できるように処理されてもよい。
【０１１４】
上記のかみ合い状態推定処理プログラム４１は、コンピュータ読み取り可能な記憶媒体として、ＲＯＭやハードディスク等の他、フラッシュメモリ等の不揮発性メモリ、ＣＤ−ＲＯＭ等の可搬型記録媒体に格納されていてもよい。記録媒体は、非一時的な記録媒体であってもよい。また、かみ合い状態推定処理プログラム４１であるプログラムデータを所定の通信回線を介して提供する媒体としては、キャリアウェーブ（搬送波）が適用されてもよい。
【０１１５】
以上、本実施形態について説明した。しかしながら、本発明は、上記の実施形態に限られるものではない。例えば、上記実施形態では、ＦＲＦの数値シミュレーションに適用される減速機のモデル化の一例を示したが、モデル化の方法はこの例に限られるものではない。また、実数値ＧＡにより最適化される減速機の振動モデルのパラメータ、および、実数値ＧＡにより最適化される歯面の偏差を示すパラメータも、上記実施の形態で示したものに制限されるものではない。
【０１１６】
また、上記実施の形態では、ＬＴＣＡを利用してＬＳＴＥの推定値から歯車対のかみ合い状態を推定および解析した例を示したが、かみ合い状態を推定する方法は、これに限られない。例えば、ＬＴＣＡを使わない実験的な方法として、歯車の合成誤差、および、その他の様々な誤差が明らかになっている多数の歯車装置の伝達誤差を専用の試験機を用いて計測し、様々な誤差に起因する伝達誤差の計測結果を集めたデータベースを用いてもよい。この場合、データベースの中から、試験対象の歯車装置と同一形態の歯車装置であって、試験対象の歯車装置のＬＳＴＥの推定値と合致するＬＳＴＥを有する歯車装置のデータを抽出する。これにより、試験対象の歯車装置のかみ合い状態を、データベースから抽出された歯車装置のデータによって推定および解析することができる。
【符号の説明】
【０１１７】
１０  減速機（歯車装置）
３１  回転検知器
３２  加速度センサ
３５  データロガー
３６  計算装置
４１  プログラム
Ｇ１、Ｇ２  歯車
【図１】
2018025462.tif 000013
【図２】
2018025462.tif 000014
【図３】
2018025462.tif 000015
【図４】
2018025462.tif 000016
【図５】
2018025462.tif 000017
【図６】
2018025462.tif 000018
【図７】
2018025462.tif 000019
【図８】
2018025462.tif 000020
【図９】
2018025462.tif 000021
【図１０】
2018025462.tif 000022
【図１１】
2018025462.tif 000023
【図１２】
2018025462.tif 000024
【図１３】
2018025462.tif 000025
【図１４】
2018025462.tif 000026
【図１５】
2018025462.tif 000027
【図１６】
2018025462.tif 000028
【図１７】
2018025462.tif 000029
【図１８】
2018025462.tif 000030
【図１９】
2018025462.tif 000031
【図２０】
2018025462.tif 000032
【図２１】
2018025462.tif 000033
【図２２】
2018025462.tif 000034