(57)【特許請求の範囲】
【請求項１】
機械構造系の、力と変位との間に履歴を有する非線形復元力特性を解析するに際し、
前記履歴を有する非線形復元力特性を等価的な線形特性で表して近似的にモード解析を行うことを特徴とする、機械構造系の履歴を有する非線形復元力特性の解析方法。
【請求項２】
前記等価的な線形特性で非線形復元力特性を表すに際し、
先ず、実系における入力と出力との関係から実験的または模擬的に基本非線形復元力特性を求め、
次いで、前記基本非線形復元力特性に対する正弦波入力と、この正弦波入力が前記基本非線形復元力特性を通過して歪んだ出力の基本波振動数成分との関係から、伝達特性、等価剛性および等価粘性減衰を数値解析的に求め、
前記伝達特性、等価剛性および等価粘性減衰を多自由度運動方程式と連立させて、機械構造系全体の振動数応答特性を求めることを特徴とする、請求項１記載の機械構造系の履歴を有する非線形復元力特性の解析方法。
【請求項３】
前記機械構造系の非線形運動方程式に対してニュートン−ラプソン法で解を求めるに際し、
非線形要素の等価剛性のデータを求めておき、この非線形要素の等価剛性のデータを区間的多項式によって補間することを特徴とする、請求項２記載の機械構造系の履歴を有する非線形復元力特性の解析方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
この発明は、機械構造系の、力と変位の間に履歴を有する非線形復元力特性を解析するための、近似的なモード解析法（近似的モード解析方法）に関するものであり、特に、振動試験等により実験的に、あるいはＦＥＭ（有限要素法）解析等により模擬的に得られた、離散的な荷重ひずみ特性のデータを用いて履歴特性を再構築し、任意の入力に対して、前記得られたデータに沿った出力を与え、機械構造系の振動特性の改善を図る解析方法に関するものである。
【０００２】
【従来の技術】
モード解析法は、機械構造系の振動解析の手法として技術が確立され、ソフトウェアとしても市販されて、多用されている。しかし、モード解析法自体は線形系を仮定した理論の上に成り立っているのに対し、実際の機械構造系は、ボルト接合部、軸受け等の滑り面部、回転部、図19に示す如きエンジンマウントのゴムパッド、ロボットの関節等において、複数の非線形特性を含むのが通常である。
【０００３】
履歴特性を有さない非線形特性が含まれる場合については、等価線形化法によって非線形特性を表し、近似的にモード解析を実現する方法が提示され、工作機械の剛性設計と加工精度向上、配管系の支持特性評価と耐震設計法への反映、エンジンマウントでの支持における最適設計等に寄与してきた。
【０００４】
しかし、ゴムパッドに典型的に代表される復元力特性は、履歴を有する非線形の特性であり、この特性を含む機械構造系の振動数応答特性を設計段階で推定して、その推定値と、実験で求められる振動数応答特性とを対比した場合の誤差を少なくすることは、従来は困難であった。
【０００５】
設計段階で、要素の構造や、要素同士を結合する特性等を考慮して振動特性を求める際には、有限要素法を用いて解析を行うのが一般的である。しかし、振動特性に関しては時刻歴による解析が主体であり、振動数応答特性を求める場合には振動数毎に時刻歴応答を繰り返すことが必要となるため、線形系の応答を求める場合でも膨大な計算時間を要する。また運動方程式から振動数応答特性を求めることも可能である。しかし、非線形特性を含む場合には、その特性を表し、機械構造系全体の振動特性を求める一般的、簡潔な方法が提示されているとは言い難い。
【０００６】
機械構造系の振動実験を行えば、これらの特性を含んだ振動数応答特性が求められる。モード解析法の中では、実験モード解析と呼ばれるものであり、この手法では、実験で振動数応答特性を求め、固有振動数、減衰定数、質量、ばね、減衰係数等を推定して、系の特性改善を図っている。
【０００７】
しかし、非線形特性が含まれる場合には、線形の特性から振動数応答特性が偏って求められる。従って、振動数の偏りや、固有振動数近傍の特性の歪み等により、固有振動数の推定が難しく、減衰定数は過大に評価されることが一般的であり、この結果、特性の改善に誤った対処が行われたり、適切な対処ができなかったりしていた。
【０００８】
履歴を有しない非線形復元力特性に対しては、等価線形化法を多自由度系運動方程式に適用することで多元の非線形連立方程式が構成でき、これを解くことによって、その解を求めること、すなわち複数非線形特性を含む多自由度系の振動数応答特性を求めることが可能とされてきた。また、有限要素法と結合して非線形復元力特性で要素の構造の振動特性を結合するビルディングブロック法（ＢＢ法）を用いることによって、多自由度非線形振動系のモード解析が近似的に、従来のモード解析法の延長上でできるとされていた。
【０００９】
しかしながら上記方法は、自動車のエンジンを車体に搭載する際に用いられるゴムパッドに典型的に見られるような、履歴特性のある非線形復元力特性を含む系には対処できていなかった。また、軸受け等の非線形特性も履歴非線形特性で表される場合が少なくなく、これにも対処できていなかった。
【００１０】
機械系の振動特性解析において、剛性部分に履歴特性がある場合にこれを表示する方法としては、従来、バイリニアモデル、トリリニアモデル等の、履歴特性を区分的な直線の組み合わせによって表すモデルや、Ｒａｍｂｒｅｇ−Ｏｓｇｏｏｄ型などの、数式で表す履歴特性モデルが用いられている。しかし、実機の非線形特性をこれらのモデルの特性によって模擬することを考える場合、以下の如き問題点がある。
【００１１】
バイリニアモデル、トリリニアモデルは、直線の組み合わせによって履歴特性を表す方法である。これらは表示を行う対象となる履歴特性を区間ごとに、バイリニアモデルにおいては2つ、トリリニアモデルでは3つの傾きの直線とみなし、この直線を区間的に繋ぎ合わせることによって履歴特性の形状を模擬するものである。この方法には、表示を行う対象となる履歴特性が複雑な形状のために２〜３本の直線の区間的な組み合わせでは表現しきれない場合があるという問題がある。
【００１２】
さらに、表示を行う対象となる履歴特性を、傾きが２つ〜３つの直線とみなしてしまうため、区間的に傾きが変化する。このため、表示する対象となる履歴特性が滑らかな形状をしている場合でも、傾きが区間的に変化することによって、任意の入力に対するモデルからの出力は滑らかではなくなる。
【００１３】
また、バイリニアモデルによって履歴特性の表示を行う場合、表示の対象となる履歴特性に対し区間的に直線回帰を行う際に、傾きが変化する位置の取り方によって、同様の履歴特性の表示をしたにもかかわらず得られる結果が変わってしまうという問題がある。
【００１４】
一方、Ｒａｍｂｒｅｇ−Ｏｓｇｏｏｄ型で表される履歴特性モデルは、多項式によって履歴特性を表す方法である。これは、骨格曲線によって最大、最小の座標を決め、その座標から下降曲線と上昇曲線とを結ぶことによって履歴特性の形状を模擬するものである。これによって履歴特性を表す場合には、表示を行う対象となる履歴特性に対して最小二乗法等を用いて、履歴特性モデルと合致させる必要がある。しかし振動試験やＦＥＭ解析等から得られる履歴特性が複雑な形状の場合には、曲線適合させた履歴特性モデルと細部にわたって特性を合わせるようパラメータを設定することは難しい。
【００１５】
また、Ｒａｍｂｒｅｇ−Ｏｓｇｏｏｄ型の履歴特性モデルで表される系について、等価線形化法によって近似的に振動特性解析を行う場合には、非線形要素の近似として任意の入力に対する出力のフーリエ級数を求める必要があるが、このモデルを表示するパラメータによっては、解析的にフーリエ級数を求めることが困難となる場合がある。
【００１６】
履歴特性の有無による等価剛性および等価減衰の違いは、履歴特性がない場合は入出力の関係で位相の遅れ、進みがないのに対し、履歴特性が有る場合は位相関係が含まれており、運動方程式と連立するときにこれを考慮する必要があることである。
【００１７】
履歴非線形特性を記述関数（自動制御の分野での等価線形化法）で表すことについては、電気系の分野では、電流と磁場の特性を考慮した場合に、実際の特性に基づき、これを数値的に記述関数で表すことが可能とされていた。しかし、機械構造系に見られるように多自由度系に適用し、モード解析としてその振動数応答特性を表すことについては、位相特性を当初から持つ多元連立非線形方程式の解法も、数学的関心から対象とされてきた領域を別として特には示されてきていなかった。
【００１８】
【発明が解決しようとする課題】
過去の履歴に依存する変位ひずみ特性を持つ非線形系の振動特性解析を行う場合には、先に示した履歴特性のモデルに対して直線回帰や最小二乗法等の当てはめを行うことが必要であった。このため、当てはめる対象となる振動試験、ＦＥＭ解析等から得られる履歴特性のデータに対して、解析に用いる履歴特性モデルの形状を細部にわたって合わせることが難しく、振動特性解析を行う際にこの点が問題となっていた。
【００１９】
【課題を解決するための手段およびその作用・効果】
この発明は、かかる課題を有利に解決することを目的とするものであり、この発明の解析方法は、履歴型非線形復元力特性を有する機械構造系の振動特性を振動数応答特性で表すことに関し、これを実系の非線形特性に基づいて数値的に等価剛性、等価減衰で表し、多自由度系の運動方程式と連立させて求め、近似的なモード解析を可能とする。これによって、機械構造系振動特性の評価、改善に資する。これにより、設計段階で系の特性を予測する一方、実系の振動実験で求められる振動特性に対して、振動系パラメータのより正確、合理的な推定を可能とするものである。
【００２０】
またこの発明は、振動試験等から得られた履歴特性のデータを何らかのモデルに当てはめる方法による履歴特性の表示方法を用いずに、その得られた履歴特性のデータを基にして、任意の入力に対してこの得られた履歴特性に沿った出力をする手法として、電気電子分野において磁気特性の記録に用いるプライザッハモデルを用いて、機械系の荷重ひずみに関する履歴特性モデルを作り出すものである。またこの履歴特性表示法を用いて、試験特性を取り込んだ機械振動特性解析を行う。
【００２１】
このプライザッハモデルによって、振動試験等から得られる履歴特性モデルを作り出す場合、このモデルは多くのリレー素子から成り立っているため、この履歴特性モデルを構成するマイナーループを測定する必要があるが、振動試験によってこの履歴特性のマイナーループを測定するようなことは難しいと考えられる。また、この発明によって得られる履歴特性モデルを振動特性解析に用いる際の課題として、履歴特性が数値的に表されるため等価線形化法を用いる場合、フーリエ変換によって求められる等価剛性を解析的に求めることが出来ないという点がある。そこでこの発明では、これを数値的に求め、解析に用いる。
【００２２】
またこの発明では、系の非線形方程式をＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法（ニュートン−ラプソン法）で解を求める際に、非線形要素の等価剛性のデータを求めておき、これをスプライン関数等の区間的多項式によって補間する。区間的多項式で表されるスプライン関数は微分することができる。これによって、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法を適用する際のヤコビ行列を作る。これを行うことによって反復ごとに得られる解を用いて履歴特性に対する出力をフーリエ級数展開し、等価剛性、等価減衰を求める方法を効率化することができる。
【００２３】
この発明によれば、振動試験、ＦＥＭ解析等から得られる履歴特性のデータを取り込んだ振動特性解析を行うことができ、より実験と密接な解析が可能となる。なお、この発明の履歴特性解析法は、履歴の存在しない非線形特性に対しても、振動試験やＦＥＭ解析等から得られるデータから同様の方法によって解析することが可能である。
【００２４】
これらにより、工作機械等で、要素構造相互間の結合部に存在する非線形特性、軸受け、滑り面部等に介在する非線形特性、エンジン懸架に用いられる支持部品の非線形特性、ロボット関節部に生じる非線形特性等を運動方程式中に取り込んで、振動数応答特性を明らかにすることが可能となり、これらを考慮した特性の改善、対処を合理的に進めることが実現する。履歴のない非線形特性を考慮した近似モード解析に比して、履歴特性を現実の特性に基づいて、数値的に等価剛性、等価減衰として表して考慮することにより、より正確な推定を迅速に行うことが可能となる。
【００２５】
従って、自動車産業等、国際間、企業グループ間で競争力が問われる時代にあって、これまでの課題を解決して、迅速な生産の立ち上げ、問題点発生時の迅速な対応を可能とすることができる。
【００２６】
【発明の実施の形態】
以下に、この発明の実施の形態を実施例として、図面に基づき詳細に説明する。ここに、図１は、この発明の解析方法の一実施例の手順を示すフローチャート、図２は、その実施例の解析方法の実施に用いる解析システムの構成を例示する説明図である。この実施例の解析方法を概説すると、ここでは先ず、図２(a) に示すように、加振試験装置の加振器１で、解析対象とするゴムパッド等の試験片２に対し加振実験をしたり、その試験片に対しＦＥＭ解析等の模擬的解析（コンピュータによるシミュレーション）を行ったりして、変位入力Ｘに対する荷重出力Ｆの、履歴を有する復元力特性のデータを得て、そのデータをオシロスコープ３を介してパーソナル・コンピュータ（ＰＣ）４に入力する。そして、図１に示すように、ＰＣ４によるデータ処理によって、プライザッハモデルを作成し、次いで等価線形化法により等価剛性および等価減衰係数を求め、次いで他自由度系の運動方程式に、その等価剛性および等価減衰係数を代入し、非線形連立方程式を解くことで、例えば図２(b) および(c) に示す如き周波数応答特性を求める。
【００２７】
具体的には、この実施例の解析方法では、履歴特性モデルを構築する際、振動試験、ＦＥＭ解析等から得られる荷重とひずみに関する履歴特性のデータを、入力変位が増加する場合と、入力変位が減少する場合とに分けて処理する。すなわち、変位Ｘを増加時の変位xiと減少時の変位xdとに分け、この2種類の変数xi, xdによって荷重Ｆとの関係を表し、変位の増加、減少の変化に対する荷重の変化ΔＦを求めて行くことによって離散値として履歴特性モデルを作り出す。
【００２８】
図３に示すように、変位と荷重の関係について、変位が増加する場合に荷重が連続な関数
Ｆ＝j(ｘ)                                             (1)
で表される場合、変位ｘに対しての荷重Ｆの増加は次のように表される。
【数１】
0004571322.tif 000002
これは変位と荷重の関係の傾きを表すものであり、同様に変位が減少する場合の荷重と変位の関係が連続な関数
Ｆ＝ｋ(ｘ)                                            (3)
で表される場合、変位ｘに対しての荷重Ｆの減少は次のように表される。
【数２】
0004571322.tif 000003
【００２９】
ここで、変位ｘを増加時の変位と減少時の変位とで独立した変数として、それぞれ増加時の変位を xi 、減少時の変位をxd とおく。これらxi , xd の変化に対するＦの変化は、次のように表す事が出来る。
【数３】
0004571322.tif 000004
これを xd 対して積分すると、式(2)の変位増加時の変位と荷重の傾きとなる。同様に xi に対して積分すると式(4)の変位減少時の変位と荷重の傾きとなる。これを変位の増加、減少の変化の範囲に対して積分して行くことによって、任意の変位に対しての荷重を出力することになる。このη(xi，xd)をプライザッハ分布関数と呼ぶ。
【００３０】
このプライザッハ分布関数を離散値で考えると次のようになる。変位の値xi, xdと荷重の値ｆの三つを与えることによって表現される図４(a)に示すような特性を考える。これをF=ｆ(xi,xd)と表すこととする。
▲１▼  変位がX＜xiにある時、F=−ｆ／２であり、
▲２▼  さらに変位を増加させると、X=xiの点でＦ＝ｆ／２となり、この状態を維持する。
▲３▼  変位が減少した時、X＞xdにある時、Ｆ＝ｆ／２であり、
▲４▼  さらに変位を減少させると、X=xdの点でＦ＝−ｆ／２となり、この状態を維持する。
以上の特性を与えたものを履歴特性の最小要素と考え、図４(b)に示すように履歴特性を、この角型で表される小要素が集合したものと考える。
この履歴特性は増加時と減少時の経路を逆方向に動く事がないと考え、特性の変位xの最大・最小値をXs, -Xsとすると、上記η(xi,xd)の定義域は
【数４】
0004571322.tif 000005
となる。プライザッハ分布関数を xi − xd 平面で見た時、この分布関数の定義域は図４(c)のような三角形の領域となる。三次元空間で見たプライザッハ分布関数を図５に示す。
【００３１】
（変位増加時の荷重の変化）
分布関数の説明をxi − xd 平面で行う。変位をXよりX+ΔX (ΔX＞0 )に増加させた場合、増加時の変位 xiがXとX+ΔXの間にある時の分布関数η(xi, xd)の全てを足し合わせたものが荷重の増分ΔＦとなり
【数５】
0004571322.tif 000006
と表わされる。
【００３２】
図６に、振幅が Xs の履歴特性について、変位の増加時における分布関数の積分領域の変化を示す。図６(a)の▲１▼では荷重Fは最小値Fminであり
【数６】
0004571322.tif 000007
と表される。このFmin と、変位増加部分のη(xi, xd)を積分したものΔＦとを足したものが荷重Ｆとなる。
【数７】
0004571322.tif 000008
積分される斜線部分は、図６(b)で▲１▼の-Xsの位置から▲２▼の０まで増加して行き、図６(c)においてXs▲３▼まで増加し、Ｆの最大値Fmaxとなる。
【００３３】
（変位減少時の荷重の変化）
変位をXよりX−ΔXに減少させた場合、減少時の変位xdがXとX−ΔXの間にある時の分布関数η(xi, xd)の全てを足し合わせたものが荷重の減少ΔＦとなり、
【数８】
0004571322.tif 000009
と表わされる。
【００３４】
図７に、変位の減少時における分布関数の積分領域の変化を示す。図７(a)の▲４▼では荷重Ｆは最大値Fmaxとなり、
【数９】
0004571322.tif 000010
と表される。このFmaxと、変位減少部分のη(xi, xd)を積分したものΔＦとを足したものが荷重Ｆとなる。
【数１０】
0004571322.tif 000011
斜線部分は図７(b)で▲４▼のXsの位置から▲５▼の０まで減少して行き図７(c)において-Xs▲６▼まで減少し、Ｆの最小値Fminとなる。
以上がプライザッハモデルによる履歴特性表示の概念である。
【００３５】
これを離散値で表すと以下のようになる。離散値で表された分布関数の分割数はM個、入力の最大値となる振幅XはN個目の分割に存在すると考える。変位がここからΔX増加したと考えた時、ΔXの存在する分割の位置をNXとすると変位の増加に対する荷重の変化は、図８(a)に示すように考えられ、次のように表される。
【数１１】
0004571322.tif 000012
【数１２】
0004571322.tif 000013
【００３６】
同様に、変位がここからΔX減少したと考えた時、ΔXの存在する分割の位置をNXとすると変位の減少に対する荷重の変化は、図８(b)に示すように考えられ、次のように表される。
【数１３】
0004571322.tif 000014
【数１４】
0004571322.tif 000015
【００３７】
かかる離散値による履歴特性表示の処理のフローチャートを図９にまとめて示す。
【００３８】
この履歴特性の表示するための処理について、図10に示すフローチャートに基づいて詳細に説明する。
（手順1）
増加時の変位xi、減少時の変位xdに関する平面xi−xd平面を考える。振動試験またはＦＥＭ解析等で、対象とするゴムパッド等の試験片の履歴特性のループについて、図11に示すように、等間隔ΔXでループの振幅がX1,X2, ・・・,Xnと増加するデータを取る。この履歴特性データの変位の値-Xn,-Xn-1,・・・,-X1, 0,X1,X2, ・・・,Xnでの荷重の値を求める。変位の値-Xn,-Xn-1,・・・,-X1, 0,X1,X2,・・・,Xnでの荷重のデータがない場合はデータを補間して求める。
【００３９】
（手順2）
この履歴特性のデータを、変位増加時の荷重データと、変位減少時の荷重のデータとに分ける。振幅Xk(k=1,2,…,n)での変位の増加と変位増加時の荷重を、各振幅の最小値-Xkの時の荷重をF-k,2kとし最大値Xkの時の荷重をFk,2kとした時、2k+1個のデータを、
xi(k) =[-Xk,-Xk-1, ・・・,-X1,X0,X1,X2,・・・,Xk ]  (X0=0)  (17)
Fi(k) =[F-k,2k,Fk1,Fk2, ・・・,Fk,2k]      (18)
と表し、また変位の減少に対する変位減少時の荷重は、
xd(k) = [Xk , Xk-1 ,・・・, -X1, X0, -X1, -X2,・・・, -Xk ]  (X0=0)     (19)
Fd(k) =[ Fk,2k,F-k1, F-k2, ・・・, F-k,2k+1 ]      (20)
とする。
【００４０】
（手順3）
変位の増加に対する荷重の変化ΔFiを求める。ΔFiは、
ΔFi(k) = [Fik1 , Fik2 , ・・・, Fik,2k ]
= [Fk1 - F-k,2k, Fk2 - Fk1, Fk2 - Fk1,・・・, Fk,2k - Fk,2k-1 ]  (21)
と表される。
また変位の減少に対する荷重の変化ΔFdを求める。ΔFdは、
ΔFd(k) = [Fdk1 , Fdk2 ,・・・, Fdk,2k ]
= [F-k1 - F-k,2k, F-k2 - F-k1, F-k2 - F-k1,・・・, F-k,2k - F-k,2k-1 ]    (22)
と表される。
【００４１】
（手順4）
この変位の変化に対する荷重の変化量を用いて、k＝1からｋ＝nまで順にΔFi,ΔFdを次のように処理することによって、η(xi,xd)が求められる。
k=1の時η(xi,xd)は3個となり、次のように表される。
η(1)(xi,xd) = [η(X0,X1),η(-X1,X1),η(-X1, X0)]
= [Fi12,Fd11-Fi12,Fd12] = [Fi12, Fd11-η(X0,X1), Fd12]    (23)
k=2の時η(xi,xd)は7個となり、k=1の時に求めたFi12, Fi11 , Fd12 , Fd11を用いて、次のように表される。
η(2)(xi,xd) =[η(X1,X2),η( X0,X2),η(-X1,X2),η(-X2, X2),η(-X2,-X1),η(-X2,X0) ,η(-X2, X1)]
=[Fi24, Fi23-Fi12, Fi22-Fi11, Fi21-{η(X1,X2)+η(X0,X2)+η(-X1,X2),Fd22-Fd11, Fd23-Fd12, Fd24}      (24)
【００４２】
同様に、k=nの時η(xi,xd)は4n-1個となり、k=nの時に求めたFin-1,2n-2, Fin-1,2n-3 ,・・・, Fin-1,1 , Fdn-1,1, ・・・, Fdn-1,2n-3, Fdn-1,2n-2,を用いて次のように表される。
【数１５】
0004571322.tif 000016
【００４３】
（手順5）
これらを次のように表すことによって、履歴特性モデルを表すプライザッハ分布関数表が求まる。
【数１６】
0004571322.tif 000017
【００４４】
図12，13は、変位の振幅が X1からX2へ増加して行く場合の処理を示し、この説明を以下に示す。
[1] 変位が-X1から０まで増加したときの荷重の増加ΔFをFi11とし、-X1 での荷重をF-12，０での荷重をF11とすると、Fi11は図12(a)の斜線部▲１▼＋▲３▼の部分となる。
Fi11=F11 - F-12    (28)
[2] さらに X1までの増加に対する荷重の増加をFi12とし、X1 での荷重をF12とした時、 ▲１▼＋▲２▼+▲３▼となり、図12(b) に示すように、▲２▼の値が求まる。
Fi12=F12-F11   (29)
η(X0, X1) = Fi12   (30)
[3] 次に、変位がX1 から０まで減少した時の荷重の減少分をFd11とし、０での荷重をF-11とすると、図12(c)の▲１▼+▲２▼となり、η(-X1, X0)の値が求まる。
Fd11=-(F-11 - F12)   (31)
η(-X1, X0) = Fd12    (32)
すでに求まっている▲２▼η(X0,X1)より▲１▼η(-X1,X1)の値が求まる。
Fd12=−(F-12 - F-11)  (33)
η(-X1,X1) = Fd11 - Fi12  (34)
【００４５】
[4] 同様に、次の1周期分の振幅X2 (X2>X1) に対して以下の処理を行う。
-X2から-X1 まで増加に対する荷重の変化Fi21は、図13(d) に示すように、-X2 での荷重をF-24 、-X1 での荷重をF21とすると▲４▼〜▲７▼を足した値となる。
Fi21=F21 - F-24  (35)
[5] -X1 から０まで増加に対する荷重の変化Fi22は、０ での荷重をF22 とすると、 ▲１▼，▲３▼，▲４▼〜▲８▼を足した値となる。▲１▼，▲３▼，▲４▼〜▲７▼を足した値は求まるので、これから図13(e) に示すように、γη(-X2,-X1)が求まる。
Fi22=F22 - F21  (36)
η(-X2,-X1) = Fd22 - Fi11  (37)
[6] ０からX1 まで増加に対する荷重の変化Fi23は、X1 での荷重をF23 とすると、▲１▼〜▲９▼を足した値となり、▲１▼〜▲８▼までを足した値は求まるので、これから図13(f) に示すように、▲９▼の値η(-X2,X0)が求まる。
Fi23=F23 - F22  (38)
η(-X2, X0)  = Fi23 - Fi12  (39)
[7] X1 からX2 まで増加に対する荷重の変化Fi23は、X2 での荷重をF24 とすると、▲１▼〜丸10を足した値となり、▲１▼〜▲９▼までを足した値は求まるので、これから図13(g) に示すように、丸10の値η(-X2, X1)が求まる。
【外１】
0004571322.tif 000018
Fi24=F24 - F23     (40)
η(-X2, X1) = Fi24  (41)
【００４６】
[8] X2 からX1 まで減少に対する荷重の変化Fd21は、X1 での荷重をF-21 とすると▲７▼〜丸10を足した値となり、▲８▼〜丸10までを足した値は求まるので、これから図13(h) に示すように、▲７▼の値η(-X2, X2)が求まる。
Fd21=-(F-21 - F24)  (42)
η(-X2, X2) = Fd21 - (η8+η9+η10)     (43)
[9] X1 から０まで減少に対する荷重の変化Fd22は、０での荷重をF-23 とすると▲１▼▲２▼▲６▼▲７▼〜丸10を足した値となり、▲１▼▲２▼▲７▼〜丸10までを足した値は求まるので、これから図13(i) に示すように、▲６▼の値η(-X1,X2)が求まる。
Fd22=-(F-22 - F-21)  (44)
η(-X1,X2) = Fd22 - Fd11  (45)
[10] ０から-X1まで減少に対する荷重の変化Fd23は、-X1 での荷重を-F23 とすると▲１▼〜▲３▼，▲５▼〜丸10を足した値となり、▲１▼〜▲３▼、▲６▼〜丸10までを足した値は求まるので、これから図13(j) に示すように、▲５▼の値η( X0,X2)が求まる。
Fd23=-(F-23 - F-22)  (46)
η( X0,X2) = Fd23 - Fd12  (47)
[11] −X1 から-X2 まで減少に対する荷重の変化Fd24は、-X2 での荷重を-F24 とすると▲１▼〜丸10を足した値となり、▲１▼〜▲３▼、▲５▼〜丸10までを足した値は求まるので、これから図13(k) に示すように、▲４▼の値が求まる。
Fd24=-(F-24 - F-23)  (48)
η(X1, X2) = Fd24     (49)
【００４７】
この作業を振幅X1 ,X2 ,X3 , ・・・( X1 ＜X2 ＜X3 ＜・・・)と増加させていき、所望の変位振幅まで繰り返していくことによって、分布関数表を作り出すことができる。
なお、この発明の方法で振動試験から得られた履歴特性のデータより履歴特性モデルをつくりだす場合、プライザッハ分布関数表をスプライン関数等を用いて2次元補間することによってより詳細な形の履歴特性モデルが作ることができる。
【００４８】
手順1の方法を用い、(35)，(36)，(37)によって表されるＲａｍｂｅｒｇ−Ｏｓｇｏｏｄ型の履歴特性に変位xを0.05刻みで‐1から1まで与え、α= 4,β= 1,γ= 2とし、これから得られる荷重Ｆより分布関数表の作成を行った。この結果を図14に示す。また、その分布関数から得られる履歴特性を図15に示す。
【数１７】
0004571322.tif 000019
【００４９】
以下では、上記プライザッハモデルから得られた履歴特性を持った系の等価線形化法による周波数応答特性解析について述べる。ここでは、図１６に示す手順で、多自由度非線形連立方程式の解を求める。
（手順6）
剛性に履歴を含む非線形系の運動方程式は次のように表される。
【数１８】
0004571322.tif 000020
ここで履歴特性はｆ(x, Bx)で表されると考えて、この運動方程式の解を近似的に
【数１９】
0004571322.tif 000021
と置く。履歴特性に対する出力ｆ(x, Bx)は、φに関しての2πの周期関数となる。
式(54)で表される入力に対するプライザッハモデルの出力を(13)(14)(15)(16)によって求める（図４参照）。
【００５０】
（手順7）
これをＦＦＴによってフーリエ展開し、調和成分のみを残すと次のように表される。
【数２０】
0004571322.tif 000022
これを用いて運動方程式を等価な線形系に置き換えると次のように表される。
【数２１】
0004571322.tif 000023
ここで等価剛性と等価減衰係数は、
【数２２】
0004571322.tif 000024
と表される。
【００５１】
式(59)を求め、これから式(58)を解くことによって近似解を得ることができる。
式(54)の入力に対する、プライザッハモデルで表される履歴特性から得られる出力(55)は離散値となるため、このフーリエ係数(56)(57)は数値的に求める。これによって式(59)が得られる。
本実施例では、式(59)のXの取りうる範囲でのke(X)，ce(X)を予め求める。Xの範囲、ke(X), ce(X)は、
【数２３】
0004571322.tif 000025
で表される。収束計算のたびに求まる式(54)に対して式(55)を求める必要がなくなるため計算効率が改善される。
【００５２】
（手順8）
式(58)から、この調和振動の振幅と位相は次の式によって求まる。
【数２４】
0004571322.tif 000026
この手法を多自由度系に適応する。[Ｍ]を質量行列、[Ｃe(Ｘ)]を等価減衰係数Ｃe(Ｘ)によって構成される行列、[Ｋe(Ｘ)]を等価剛性係数Ｋe(Ｘ)から構成される行列とおくと、運動方程式は、
【数２５】
0004571322.tif 000027
となり、XとFは次のようなベクトルとなる。
【数２６】
0004571322.tif 000028
【００５３】
（手順9）
この連立方程式にＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法を適応する場合、[Ｃe(Ｘ)]と[Ｋe(Ｘ)]を構成するＣe(Ｘ)とＫe(Ｘ)は変位振幅Ｘの関数となるため、このヤコビ行列を求めるには、等価減衰係数ｃe(Ｘ)と等価剛性係数ke(X)との導関数を求める必要がある。
そこで本実施例では、得られた[Ｃe(Ｘ)]，[Ｋe(Ｘ)]に対し区分的な多項式によるスプライン関数によって補間を行ない、このスプライン関数を微分する事によって導関数を求める。
N個の接点x = X0 , X1,・・・, XN-1 ,を持つm次のスプライン関数はαi を定数係数、p(x) をm次の多項式とすると(68)で表される。
【数２７】
0004571322.tif 000029
区分的な多項式は微分可能となるため、m次のスプライン関数の微分は次式によって求められる（図17参照）。
【数２８】
0004571322.tif 000030
この手法を用いて[Ｃe(Ｘ)]と[Ｋe(Ｘ)]の要素の微分が求められ、運動方程式ｆ(VX)のヤコビ行列ｆ′(VX)が作られる。
【数２９】
0004571322.tif 000031
【００５４】
（手順10）
式(65)に対してＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ法を適用すると、
【数３０】
0004571322.tif 000032
k回目のXkに対するk+1回目の近似解Xk+1は
【数３１】
0004571322.tif 000033
と表される。この式(72)を繰り返すことによって解を求めていく。かかる計算を式(65)のωごとに行うことによって、図18に示す如き周波数応答特性が求まる。
【図面の簡単な説明】
【図１】  この発明の解析方法の一実施例の手順を示すフローチャートである。
【図２】  上記実施例の解析方法の実施に用いる解析システムの構成を例示する説明図である。
【図３】  履歴を有する非線形復元力特性を例示する特性図である。
【図４】  (a)は、離散値で表される場合の履歴特性の最小要素を示す説明図、(b)は、その最小要素の集合によって表される履歴特性、(c)は、xd−xi平面での分布関数の存在領域を示す説明図である。
【図５】  離散値で表されるプライザッハ分布関数を三次元空間で見た状態で示す説明図である。
【図６】  変位増加時の荷重の変化を示す説明図である。
【図７】  変位減少時の荷重の変化を示す説明図である。
【図８】  離散値での分布関数表の足し合わせの領域を示す説明図である。
【図９】  離散値による履歴特性表示のための処理を示すフローチャートである。
【図１０】  振動試験データから分布関数表を作成するアルゴリズムを示すフローチャートである。
【図１１】  分布関数作成時におけるデータ取得点を示す説明図である。
【図１２】  振動試験データの処理による履歴特性モデルの作成方法（振幅X1）を示す説明図である。
【図１３】  振動試験データの処理による履歴特性モデルの作成方法（振幅X2）を示す説明図である。
【図１４】  数式モデルデータから作成したプライザッハ分布関数を三次元空間で見た状態で示す説明図である。
【図１５】  数式モデルによる履歴特性データとプライザッハモデルによる履歴特性とを示す説明図である。
【図１６】  多自由度非線形連立方程式の解を求めるための処理を示すフローチャートである。
【図１７】  等価剛性係数ke(X)とその微分値ke′(X)を例示する特性線図である。
【図１８】  履歴を有する一自由度系の周波数応答特性を例示する特性線図である。
【図１９】  非線形復元力特性（KNL）を有する機械構造系の一例としてのエンジンマウントを持つエンジンシェイクモデルを示す説明図である。
【符号の説明】
１  加振器
２  試験片
３  オシロスコープ
４  パーソナル・コンピュータ
【図１】
0004571322.tif 000034
【図２】
0004571322.tif 000035
【図３】
0004571322.tif 000036
【図４】
0004571322.tif 000037
【図５】
0004571322.tif 000038
【図６】
0004571322.tif 000039
【図７】
0004571322.tif 000040
【図８】
0004571322.tif 000041
【図９】
0004571322.tif 000042
【図１０】
0004571322.tif 000043
【図１１】
0004571322.tif 000044
【図１２】
0004571322.tif 000045
【図１３】
0004571322.tif 000046
【図１４】
0004571322.tif 000047
【図１５】
0004571322.tif 000048
【図１６】
0004571322.tif 000049
【図１７】
0004571322.tif 000050
【図１８】
0004571322.tif 000051
【図１９】
0004571322.tif 000052