(57)【要約】【課題】任意の周波数特性をもつばね・減衰要素と，任意の非線形性をもつばね・減衰要素が共存するシステムの過渡的および定常的な振動・運動挙動の解析を可能にする。【解決手段】任意の周波数特性を持つ、入力変位Ｘ（ω）と出力荷重Ｆ（ω）の伝達関数Ｈ（ω）をラプラス変換してｓ平面上の関数Ｈ（ｓ）に変換する。次に、Ｈ（ｓ）をｚ変換によりｚ平面上の関数Ｈ（ｚ）に変換する。そして、このＨ（ｚ）を逆ｚ変換し、入力と出力の離散時間データを示す漸化式を得る。この漸化式を利用した、運動方程式を時刻応答解析する。【選択図】　　　　図３
2004101304.tif 000002
【特許請求の範囲】
【請求項１】
システムの過渡的および定常的な挙動を解析する時刻歴応答解析プログラムであって、
周波数特性を有するシステムの構成要素の振動および減衰特性を周波数領域の伝達関数として表現するステップと、
得られた周波数伝達関数をラプラス変換のｓ平面上の伝達関数に変換するステップと、
得られたｓ平面上の伝達関数をｚ変換のｚ平面上の伝達関数に変換するステップと、
得られたｚ平面上の伝達関数に逆ｚ変換を施して、離散時間領域における入力変位と出力の関係式を導くステップと、
得られた関係式と、他の非線形要素の動特性を考慮して、システムの挙動を示す微分方程式を時間領域で数値積分するステップと、
をコンピュータに実行させ、システムの挙動の時刻歴応答を解析する時刻歴応答解析プログラム。
【請求項２】
周波数特性をもつばね・減衰要素と、非線形性をもつばね・減衰要素が共存する動的システムの過渡的および定常的な振動挙動を解析する時刻歴応答解析プログラムであって、
動的システムの構成要素において，そのばね・減衰特性（動剛性）が周波数依存性を持つ場合に，その動剛性を周波数領域の伝達関数（出力荷重と入力変位の比）として表現するサブプログラムと、
得られた周波数伝達関数をラプラス変換のｓ平面上の伝達関数に変換するサブプログラムと、
得られたｓ平面上の伝達関数をｚ変換のｚ平面上の伝達関数に変換するサブプログラムと、
得られたｚ平面上の伝達関数に逆ｚ変換を施して、離散時間領域における入力変位と出力荷重の関係式を導くサブプログラムと、
得られた関係式と他の非線形要素の動特性を考慮して、動的システムの挙動の方程式を時間領域で数値積分するサブプログラムと、
を有し、
これらサブプログラムを組み合わせてコンピュータに実行させることによって、動的システムの振動挙動を解析する時刻歴応答解析プログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
機械システムなどの過渡的および定常的な挙動を解析する時刻歴応答解析プログラムに関する。
【０００２】
【従来の技術】
従来より、機械システムや、電子回路システムにおいて、その動作の最適化を図るなどの目的で、過渡的および定常的な振動現象などを解析することが必要であり、多くの解析法が提案されている。
【０００３】
（１）機械システム（特許文献１）や電子回路（特許文献２）の時刻歴応答解析法
特許文献１には、常微分方程式で示される油圧・制御・車両系の機械システムの過渡振動を数値積分法で解析することが示されている。また、特許文献２には、電子回路の過渡振動現象を高速・高精度に解析する数値積分法が示されている。
【０００４】
（２）線形な機械システムの周波数応答解析法（特許文献３、４）
特許文献３、４には、車両や一般構造物を対象として、線形要素から構成される機械システムの周波数応答を解析することが示されている。また、実験モード解析の結果を使って、機械システムに作用する外力や、構造体の質量行列・剛性行列を同定する方法が示されている。
【０００５】
（３）市販の機構・振動解析ソフト
また、機械システムの振動・運動の時刻歴応答解析を行う汎用ソフトが市販されており、制御要素として２次遅れ系の伝達関数要素を備えているものが知られている。市販の機構・振動解析ソフトとしては、ＬＭＳ　Ｉｎｔｅｒｎａｔｉｏｎａｌ　社の機構解析プログラム「ＤＡＤＳ（Ｄｙｎａｍｉｃ　Ａｎａｌｙｓｉｓ　ａｎｄ　Ｄｅｓｉｇｎ　Ｓｙｓｔｅｍ）　Ｒｉｖｉｓｉｏｎ　８．０」等がある。
【０００６】
また、物体変位・振動解析方法を利用した汎用のエンジン振動解析ソフトについては、特許文献５に記載がある。
【０００７】
【特許文献１】
特開２００２−８３２５２号公報
【特許文献２】
特開２００１−３５７０９３号公報
【特許文献３】
特開２００２−７３７０３号公報
【特許文献４】
特開２００１−３５０７４１号公報
【特許文献５】
特開２０００−３０５９２２号公報
【０００８】
【発明が解決しようとする課題】
しかし、上述の従来の手法では、次のような問題がある。
【０００９】
（１）機械システム（特許文献１）や電子回路（特許文献２）の時刻歴応答解析法
非線形なばね・減衰要素を含む機械システムの過渡的な振動・運動解析は可能であるが、ゴム部品など周波数依存性の大きなばね・減衰要素を含む機械システムの解析は不可能である。
【００１０】
（２）線形な機械システムの周波数応答解析法（特許文献３、４）
ゴム部品など周波数依存性の大きなばね・減衰要素を含む機械システムの周波数応答解析は可能であるが、過渡的な振動・運動の解析や，非線形なばね・減衰要素を含む機械システムの振動解析は不可能である。
【００１１】
（３）市販の機構・振動解析ソフト
２次遅れ系の伝達関数要素を利用すれば、限られた周波数特性のばね・減衰要素は解析可能である。しかし、次数が限定されているため、任意の周波数特性をもつばね・減衰要素に対応できない。
【００１２】
本発明の目的は、任意の周波数特性をもつばね・減衰要素と，任意の非線形性をもつばね・減衰要素をが共存するシステムの過渡的および定常的な振動挙動の解析を可能にすることにある。
【００１３】
【課題を解決するための手段】
本発明は、システムの過渡的および定常的な挙動を解析する時刻歴応答解析プログラムであって、周波数特性を有するシステムの構成要素の振動および減衰特性を周波数領域の伝達関数として表現するステップと、得られた周波数伝達関数をラプラス変換のｓ平面上の伝達関数に変換するステップと、得られたｓ平面上の伝達関数をｚ変換のｚ平面上の伝達関数に変換するステップと、得られたｚ平面上の伝達関数に逆ｚ変換を施して、離散時間領域における入力変位と出力の関係式を導くステップと、得られた関係式と、他の非線形要素の動特性を考慮して、システムの挙動を示す微分方程式を時間領域で数値積分するステップと、をコンピュータに実行させ、システムの挙動の時刻歴応答を解析することを特徴とする。
【００１４】
このように、本発明によれば、周波数領域で表現された伝達関数から、ラプラス変換とｚ変換の理論を利用して離散時間領域における入力と出力の関係式を導く。これによって、システムの挙動を示す微分方程式の時刻歴応答解析を行う数値積分の中で、伝達関数要素に発生する出力の計算が可能となる。従って、任意の周波数依存性を有する要素を含むシステムの周波数応答解析が可能であるとともに、過渡的な振動の解析や、非線形なばね・減衰要素を含むシステムの振動解析が可能となる。
【００１５】
また、本発明は、周波数特性をもつばね・減衰要素と、非線形性をもつばね・減衰要素が共存する動的システムの過渡的および定常的な振動挙動を解析する時刻歴応答解析プログラムであって、動的システムの構成要素において，そのばね・減衰特性（動剛性）が周波数依存性を持つ場合に，その動剛性を周波数領域の伝達関数（出力荷重と入力変位の比）として表現するサブプログラムと、得られた周波数伝達関数をラプラス変換のｓ平面上の伝達関数に変換するサブプログラムと、得られたｓ平面上の伝達関数をｚ変換のｚ平面上の伝達関数に変換するサブプログラムと、得られたｚ平面上の伝達関数に逆ｚ変換を施して、離散時間領域における入力変位と出力荷重の関係式を導くサブプログラムと、得られた関係式と他の非線形要素の動特性を考慮して、動的システムの挙動の方程式を時間領域で数値積分するサブプログラムと、を有し、これらサブプログラムを組み合わせてコンピュータに実行させることによって、動的システムの振動挙動を解析することを特徴とする。
【００１６】
このように、本発明によれば、周波数領域で表現された伝達関数から、ラプラス変換とｚ変換の理論を利用して離散時間領域における入力変位と出力荷重の関係式を導く。これによって、運動方程式の時刻歴応答解析を行う数値積分の中で、伝達関数要素に発生する出力荷重の計算が可能となる。従って、ゴム部品など周波数依存性の大きなばね・減衰要素を含む機械システムの周波数応答解析が可能であるとともに、過渡的な振動・運動の解析や、非線形なばね・減衰要素を含む機械システムの振動解析が可能である。
【００１７】
【発明の実施の形態】
以下、本発明の実施形態について、説明する。
【００１８】
本発明は、コンピュータプログラムに関するものであり、通常の汎用コンピュータがプログラムの実行に利用される。プログラムは各種記録媒体の他、通信などで提供することも可能であり、提供形態が限定されるものではない。
【００１９】
まず、図１に、解析対象の機械システムの構成図を示す。構造体１０は、単独あるいは複数の部品で構成される弾性体および剛体要素である。外力ｆｅは、時刻ｔの関数ｆｅ（ｔ）として構造体に作用する。非線形要素１２は、軸受油膜やばね・減衰要素など、構造体の変位ｘと速度ｄｘ／ｄｔに依存して出力荷重（反力）ｆＮを構造体１０に返す要素である。その出力荷重ｆＮは次式のように、構造体１０の変位ｘと速度ｄｘ／ｄｔ（図においては、時間微分をドットで表している）に関する非線形関数として表現される。
【数１】
2004101304.tif 000003
【００２０】
伝達関数要素１４は、本発明で新たに解析を可能にする要素である。伝達関数要素１４は、ゴムの支持部品のように、ばね特性が周波数ωの伝達関数Ｈ（ω）で表現されるばね要素（動剛性要素）で、構造体の変位に依存して、出力荷重（反力）ｆを構造体１０に返す。
【００２１】
このシステムでは、構造体１０に印加される力は、ｆｅ、ｆＮ、ｆであり、構造体１０の変位がｘである。そこで、このシステムのダイナミクスを表す運動方程式は、次のように表される。
【数２】
2004101304.tif 000004
【００２２】
ここで、ｘは構造体の変位ベクトル、Ｍは構造体の質量行列、Ｃは減衰行列、Ｋは剛性行列を表す。
【００２３】
この運動方程式を時刻歴応答解析する場合に、周波数領域で表現された伝達関数要素１４の出力荷重（反力）ｆをどのように計算するかが、ポイントとなる。
【００２４】
「伝達関数要素の定義」
ゴム部品などの伝達関数Ｈ（ω）は、入力変位ｘと出力荷重ｆの周波数スペクトルＸ（ω），Ｆ（ω）の比として次式のように定義される。
【数３】
2004101304.tif 000005
【００２５】
ここで、周波数スペクトルＸ（ω），Ｆ（ω）および伝達関数Ｈ（ω）は、複素数として表現する。
【００２６】
図２に伝達関数の一例を示す。一般に，ゴム部品などの伝達関数Ｈ（ω）は加振実験やＦＥＭ解析（有限要素解析）で求めることができる。
【００２７】
「伝達関数要素の解析方法」
図３に、伝達関数要素１４の定式化の手順を示す。
【００２８】
（ｉ）伝達関数の変換その１（ラプラス変換のｓ平面上の関数へ変換）
まず、式（３）により、周波数ωの関数として表現された伝達関数Ｈ（ω）を、次式で表現されるラプラス変換のｓ平面上の関数に変換する。
【数４】
2004101304.tif 000006
【００２９】
ここで、ｎは伝達関数の次数、Ａｎ〜Ａ０およびＢｎ〜Ｂ０は定数、Ａｎ≠０である。
【００３０】
特に、ＦＥＭ解析や実験モード解析などにより、予め伝達関数Ｈ（ω）がωの多項式関数として求められている場合は、次の変換式で容易に式（４）の伝達関数を得ることができる。
【数５】
2004101304.tif 000007
【００３１】
また、加振実験により、伝達関数Ｈ（ω）が離散データとして与えられている場合には、最小二乗法などの同定手段により式（４）の各係数を求める。
【００３２】
（ｉｉ）伝達関数の変換その２（ｚ平面上の関数へ変換）
式（４）で表現される伝達関数を、時刻歴応答解析の数値積分で扱えるようにするため、離散時間が表現できるｚ平面上の伝達関数に変換する。ここで、数値積分の時間ステップをｄｔとする。
【００３３】
式（６）に示すｓ→ｚの変換公式（Ｔｕｓｔｉｎ変換）を式（４）に代入することにより、式（７）に示すｚ平面上の伝達関数を得ることができる。
【数６】
2004101304.tif 000008
【数７】
2004101304.tif 000009
【００３４】
なお、このｚ平面への変換は、Ｔｕｓｔｉｎ変換に限定されることなく、Ｚｅｒｏ次ホールド変換などを利用することもできる。このＺｅｒｏ次ホールド変換の場合、式（６）に代えて、ｓ＝（ｚ−１）／ｄｔが採用される。
【００３５】
（ｉｉｉ）時刻歴応答解析における伝達関数の漸化式
式（７）の伝達関数から、数値積分で扱える離散時間の漸化式を導く。入力変位ｘのｚ変換をＸ（ｚ）、出力荷重ｆのｚ変換をＦ（ｚ）とすると、入・出力の関係は、式（７）の伝達関数Ｈ（ｚ）を用いて、次のように表される。
【数８】
2004101304.tif 000010
【００３６】
式（７）を式（８）に代入することにより、式（９）が得られる。
【数９】
2004101304.tif 000011
【００３７】
式（９）に逆ｚ変換を適用すると、式（１０）が得られる。
【数１０】
2004101304.tif 000012
【００３８】
ここで、ｘｉ−ｋ，ｆｉ−ｋ（ｋ＝０，１，２，．．．，ｎ）は、入力と出力の離散時間データを表す。
【００３９】
式（１０）は、入・出力の過去ｎ個分の履歴データ、時刻ｉの入力ｘｉ、および伝達関数の係数から、時刻ｉの出力ｆｉが計算できることを示す。よって、式（１０）は、時刻歴応答解析の数値積分において、伝達関数要素の入力変位ｘから出力荷重ｆを計算する漸化式として用いることができる。
【００４０】
このように、周波数領域の伝達関数をラプラス変換、ｚ変換、逆ｚ変換することで、入力変位ｘから出力荷重ｆを計算する漸化式が得られる。
【００４１】
「運動方程式の数値積分法」
式（１０）の漸化式を利用して、運動方程式（２）を時刻歴応答解析する。式（２）の右辺にある非線形要素の力ｆＮと伝達関数要素の力ｆは、時刻ｔによって陽に決定される力ではなく、構造体１０の変位と速度に依存する未知の力となる。数値積分の発散防止のため、ニューマークβ法を利用した陰解法の数値積分法により、時刻歴応答解析を行う。
【００４２】
数値積分法の流れ図を図４に示す。数値積分の時間刻みをｄｔとする。時刻ｔｉ−１までの加速度、速度、変位が既知で、時刻ｔｉの各状態量を求める場合を想定する。
【００４３】
まず、時刻を初期設定（ｔｉ＝ｔ０）する（Ｓ１１）。次に、加速度を推定する（Ｓ１２）。例えば、図４に示すように、前の時刻の加速度（ｄ／ｄｔ）２ｘｉ−１と同じ値とする。
【００４４】
この加速度の推定値を用いて、次のニューマークβ法の式から次の時刻の変位ｘｉと速度ｄｘ／ｄｔを計算する（Ｓ１３）。
【数１１】
2004101304.tif 000013
【数１２】
2004101304.tif 000014
【００４５】
ここで、βとγ定数である。例えば、β＝１／６，γ＝１／２に設定する。
【００４６】
次に、式（１１）、（１２）の変位と速度を使って、式（１）の非線形要素１２の出力荷重（反力）ｆＮｉを計算する（Ｓ１４）。
【００４７】
また、変位と過去の状態両履歴から、式（１０）により伝達関数要素１４の出力荷重ｆｉを計算する（Ｓ１５）。
【００４８】
続いて、以上のように求めた速度ｄｘ／ｄｔ、加速度（ｄ／ｄｔ）２ｘ、各要素の出力荷重ｆＮｉ、ｆｉを式（２）の運動方程式に代入することにより、式（１３）のように、新たな加速度の推定値（ｄ／ｄｔ）２ｘｉ＊を、式（１３）で求める（Ｓ１６）。
【数１３】
2004101304.tif 000015
【００４９】
この加速度の推定値（ｄ／ｄｔ）２ｘｉ＊が前の推定値（ｄ／ｄｔ）２ｘｉ（今回は、この推定値は、前の時刻と同一としている）の差が所定の閾値Ｔｏｌ以下かを判定する（Ｓ１７）。すなわち、加速度の推定値（ｄ／ｄｔ）２ｘｉ＊が前の推定値（ｄ／ｄｔ）２ｘｉとほとんど変わらないかを判定する。このＳ１７の判定で、ＮＯ、すなわち有意差があれば、（ｄ／ｄｔ）２ｘｉ＊＝（ｄ／ｄｔ）２ｘｉとする（Ｓ１８）。そして、Ｓ１３に戻り、式（１１）からの上記計算を加速度の推定値（ｄ／ｄｔ）２ｘｉ＊が収束するまで繰り返す。
【００５０】
Ｓ１７の判定で、ＹＥＳであれば、加速度が収束したと判定して、時刻ｔｉの各状態量（加速度（ｄ／ｄｔ）２ｘｉ、速度ｄｘｉ／ｄｔ、変位ｘｉ）を決定し（Ｓ１９）、時刻をｔｉ＝ｔｉ＋ｄｔとして、次のステップに進める（Ｓ２０）。
【００５１】
そして、時刻ｔｉが終了時刻（ｅｎｄＴｉｍｅ）に至ったかを判定し（Ｓ２１）、ＮＯであればＳ１２に戻り、時刻ｔｉ＝ｔｉ＋ｄｔでの処理を行う。一方、Ｓ２１の判定で、ＹＥＳの場合に処理を終了する。
【００５２】
このようにして、各時刻における加速度、速度、変位を順次決定することができる。本実施形態は、次のような特徴がある。
【００５３】
（１）周波数領域で表現された伝達関数から、ラプラス変換とｚ変換の理論を利用して、離散時間領域における入力変位と出力荷重の関係式を導くことにより、運動方程式の時刻歴応答解析を行う数値積分の中で、伝達関数要素に発生する出力荷重の計算を可能にした。
【００５４】
（２）（１）の伝達関数要素の入力変位・出力荷重の関係式と、ニューマークβ法の陰解法を融合させることにより、他の非線形ばね・減衰要素の出力荷重も考慮して、運動方程式の数値積分を安定・高精度に実行できる。
【００５５】
（３）（１）（２）の具体例としては，ゴムマウントの周波数特性を考慮したエンジン・車両系の振動伝達解析、クラッチや流体継ぎ手の周波数特性を考慮したドライブトーレーン系の振動伝達解析などがある。
【００５６】
さらに、共振回路などの電気回路における振動解析にも本実施形態の手法を適用することができる。
【００５７】
「解析法の適用事例」
（ｉ）１自由度系の基本モデル
本解析法の妥当性を検証するため、図５に示す１自由度系の過渡振動を計算し、その結果を厳密解と比較した。このモデルは、１つのばねと、１つのダッシュポットを並列接続し、これに質点を接続した構成である。条件は、Ｍ＝０．１ｋｇ（重力＝０．１×９．８ｋｇｆ）、Ｃ＝０．０００２ｋｇｆｓ／ｍｍ、Ｋ＝１０ｋｇｆ／ｍｍである。また、式（４）の伝達関数の定義において、ｎ＝１，Ａ１＝０．０，Ａ０＝１．０，Ｂ１＝０．０００２，Ｂ０＝１０．０としている。
また、外力は質点に作用する重力のみであり、時刻０における初期値はｘ＝ｄｘ／ｄｔ＝０である。
【００５８】
図５の計算結果において、厳密解を点線で記載したが、実線で示した本解析結果と有意の差がないため、図５において点線は見えていない。このように、過渡振動の計算結果か、本実施形態の計算により、正しい過渡振動の計算が行えることが分かる。
【００５９】
（ｉｉ）ゴム支持を模擬した簡易モデル
実際のゴム部品への適用性を確認するため、図６に示すように、質点をゴムで支持した場合の過渡振動を計算した。本モデルでは、２次の伝達関数でゴムの伝達特性（動剛性）を表現している。外力は質点に作用する重力のみである。また、時刻０における初期値はｘ＝ｄｘ／ｄｔ＝０である。また、Ｍ＝１ｋｇ（重力＝１×９．８ｋｇｆ）、２次伝達関数の各係数は、Ａ０＝１．０，Ａ１＝１．２０１×ｅ−５，Ａ２＝４．００８×ｅ−８，Ｂ０＝５８．７２，Ｂ１＝４．６７７ｅ−３，Ｂ２＝５．１６４ｅ−７である。
【００６０】
この計算によって、図７に示すように、ゴム部品の過渡振動を適切に計算することができる。この計算結果は、実際のゴム部品の挙動を正しく模擬していることが確認されている。
【００６１】
（ｉｉｉ）エンジンの実働振動モデル
実際の機械システムへの適用性を確認するため、図８に示すエンジンシステムの実働振動を計算した。この例は、４気筒エンジンについての振動試験である。
エンジンのクランク軸には４つのピストンが接続され、シリンダ内の爆発によって、クランク軸が回転される。クランク軸の回転は、減速機を介し出力され、モータによってトルクが吸収される。エンジンブロックおよび減速機ケースがマウントゴムによって接地されている。
【００６２】
本モデルでは、図９に示すように、複数のマウントの中のエンジンブロックを接地するゴムマウントの一つ（上下方向の動剛性）に上記事例（ｉｉ）の伝達関数要素を適用し、エンジンの火入れ運転状態のマウント振動を計算した。なお、残りのゴムマウントは、端あるばね・減衰要素とした。
【００６３】
図１０に計算結果による変位ｘの時間変化（時刻歴波形）を示す。マウント振動の時刻歴波形は、エンジン１６回転分の過渡挙動を示しており、ゴムの伝達特性（減衰特性）によって振動が次第に周期定常状態に落ち着く様子を模擬している。なお、この振動計算は、特開２０００−３０５９２２号公報（特許文献５）に記載の物体変位・振動解析方法を利用した汎用のエンジン振動解析ソフトに、本発明の解析法に基づく新解析要素を実装することにより行った。
【００６４】
【発明の効果】
以上説明したように、本発明によれば、周波数領域で表現された伝達関数から、ラプラス変換とｚ変換の理論を利用して離散時間領域における入力と出力の関係式を導く。これによって、システムの挙動を示す微分方程式の時刻歴応答解析を行う数値積分の中で、伝達関数要素に発生する出力の計算が可能となる。従って、任意の周波数依存性を有する要素を含むシステムの周波数応答解析が可能であるとともに、過渡的な振動の解析や、非線形なばね・減衰要素を含むシステムの振動解析が可能となる。特に、機械システムにおける変位、振動について、精度の高い推定が可能になる。
【図面の簡単な説明】
【図１】伝達関数要素を含む機械システムを示す図である。
【図２】伝達関数の例を示す図である。
【図３】伝達関数要素の解析の流れを示す図である。
【図４】時刻歴応答解析の流れ図である。
【図５】１自由度系の基本モデルの構成および過渡振動の計算結果を示す図である。
【図６】ゴム支持を模擬した簡易モデルの構成および過渡振動の計算結果を示す図である。
【図７】ゴム支持簡易モデルの計算結果により変位の時間変化を示す図である。
【図８】エンジンの実働振動モデルの構成を示す図である。
【図９】マウントの取り付けを説明する図である。
【図１０】エンジン実働振動モデルの解析結果における変位の時間変化を示す図である。
【符号の説明】
１０　構造体、１２　非線形要素、１４　伝達関数要素。
【図１】
2004101304.tif 000016
【図２】
2004101304.tif 000017
【図３】
2004101304.tif 000018
【図４】
2004101304.tif 000019
【図５】
2004101304.tif 000020
【図６】
2004101304.tif 000021
【図７】
2004101304.tif 000022
【図８】
2004101304.tif 000023
【図９】
2004101304.tif 000024
【図１０】
2004101304.tif 000025