JP3924223B2 - Time history response analysis program - Google Patents

Description

【００２０】
伝達関数要素１４は、本発明で新たに解析を可能にする要素である。伝達関数要素１４は、ゴムの支持部品のように、ばね特性が周波数ωの伝達関数Ｈ（ω）で表現されるばね要素（動剛性要素）で、構造体の変位に依存して、出力荷重（反力）ｆを構造体１０に返す。
【００２１】
このシステムでは、構造体１０に印加される力は、ｆ_e、ｆ_N、ｆであり、構造体１０の変位がｘである。そこで、このシステムのダイナミクスを表す運動方程式は、次のように表される。
【数２】

【００２２】
ここで、ｘは構造体の変位ベクトル、Ｍは構造体の質量行列、Ｃは減衰行列、Ｋは剛性行列を表す。
【００２３】
この運動方程式を時刻歴応答解析する場合に、周波数領域で表現された伝達関数要素１４の出力荷重（反力）ｆをどのように計算するかが、ポイントとなる。
【００２４】
「伝達関数要素の定義」
ゴム部品などの伝達関数Ｈ（ω）は、入力変位ｘと出力荷重ｆの周波数スペクトルＸ（ω），Ｆ（ω）の比として次式のように定義される。
【数３】

【００２５】
ここで、周波数スペクトルＸ（ω），Ｆ（ω）および伝達関数Ｈ（ω）は、複素数として表現する。
【００２６】
図２に伝達関数の一例を示す。一般に，ゴム部品などの伝達関数Ｈ（ω）は加振実験やＦＥＭ解析（有限要素解析）で求めることができる。
【００２７】
「伝達関数要素の解析方法」
図３に、伝達関数要素１４の定式化の手順を示す。
【００２８】
（ｉ）伝達関数の変換その１（ラプラス変換のｓ平面上の関数へ変換）
まず、式（３）により、周波数ωの関数として表現された伝達関数Ｈ（ω）を、次式で表現されるラプラス変換のｓ平面上の関数に変換する。
【数４】

【００２９】
ここで、ｎは伝達関数の次数、Ａ_n〜Ａ₀およびＢ_n〜Ｂ₀は定数、Ａ_n≠０である。
【００３０】
特に、ＦＥＭ解析や実験モード解析などにより、予め伝達関数Ｈ（ω）がωの多項式関数として求められている場合は、次の変換式で容易に式（４）の伝達関数を得ることができる。
【数５】

【００３１】
また、加振実験により、伝達関数Ｈ（ω）が離散データとして与えられている場合には、最小二乗法などの同定手段により式（４）の各係数を求める。
【００３２】
（ｉｉ）伝達関数の変換その２（ｚ平面上の関数へ変換）
式（４）で表現される伝達関数を、時刻歴応答解析の数値積分で扱えるようにするため、離散時間が表現できるｚ平面上の伝達関数に変換する。ここで、数値積分の時間ステップをｄｔとする。
【００３３】
式（６）に示すｓ→ｚの変換公式（Ｔｕｓｔｉｎ変換）を式（４）に代入することにより、式（７）に示すｚ平面上の伝達関数を得ることができる。
【数６】

【数７】

【００３４】
なお、このｚ平面への変換は、Ｔｕｓｔｉｎ変換に限定されることなく、Ｚｅｒｏ次ホールド変換などを利用することもできる。このＺｅｒｏ次ホールド変換の場合、式（６）に代えて、ｓ＝（ｚ−１）／ｄｔが採用される。
【００３５】
（ｉｉｉ）時刻歴応答解析における伝達関数の漸化式
式（７）の伝達関数から、数値積分で扱える離散時間の漸化式を導く。入力変位ｘのｚ変換をＸ（ｚ）、出力荷重ｆのｚ変換をＦ（ｚ）とすると、入・出力の関係は、式（７）の伝達関数Ｈ（ｚ）を用いて、次のように表される。
【数８】

【００３６】
式（７）を式（８）に代入することにより、式（９）が得られる。
【数９】

【００３７】
式（９）に逆ｚ変換を適用すると、式（１０）が得られる。
【数１０】

【００３８】
ここで、ｘ_i-k，ｆ_i-k（ｋ＝０，１，２，．．．，ｎ）は、入力と出力の離散時間データを表す。
【００３９】
式（１０）は、入・出力の過去ｎ個分の履歴データ、時刻ｉの入力ｘ_i、および伝達関数の係数から、時刻ｉの出力ｆ_iが計算できることを示す。よって、式（１０）は、時刻歴応答解析の数値積分において、伝達関数要素の入力変位ｘから出力荷重ｆを計算する漸化式として用いることができる。
【００４０】
このように、周波数領域の伝達関数をラプラス変換、ｚ変換、逆ｚ変換することで、入力変位ｘから出力荷重ｆを計算する漸化式が得られる。
【００４１】
「運動方程式の数値積分法」
式（１０）の漸化式を利用して、運動方程式（２）を時刻歴応答解析する。式（２）の右辺にある非線形要素の力ｆ_Nと伝達関数要素の力ｆは、時刻ｔによって陽に決定される力ではなく、構造体１０の変位と速度に依存する未知の力となる。数値積分の発散防止のため、ニューマークβ法を利用した陰解法の数値積分法により、時刻歴応答解析を行う。
【００４２】
数値積分法の流れ図を図４に示す。数値積分の時間刻みをｄｔとする。時刻ｔ_i-1までの加速度、速度、変位が既知で、時刻ｔ_iの各状態量を求める場合を想定する。
【００４３】
まず、時刻を初期設定（ｔ_i＝ｔ₀）する（Ｓ１１）。次に、加速度を推定する（Ｓ１２）。例えば、図４に示すように、前の時刻の加速度（ｄ／ｄｔ）²ｘ_i-1と同じ値とする。
【００４４】
この加速度の推定値を用いて、次のニューマークβ法の式から次の時刻の変位ｘｉと速度ｄｘ／ｄｔを計算する（Ｓ１３）。
【数１１】

【数１２】

【００４５】
ここで、βとγ定数である。例えば、β＝１／６，γ＝１／２に設定する。
【００４６】
次に、式（１１）、（１２）の変位と速度を使って、式（１）の非線形要素１２の出力荷重（反力）ｆ_Niを計算する（Ｓ１４）。
【００４７】
また、変位と過去の状態両履歴から、式（１０）により伝達関数要素１４の出力荷重ｆ_iを計算する（Ｓ１５）。
【００４８】
続いて、以上のように求めた速度ｄｘ／ｄｔ、加速度（ｄ／ｄｔ）²ｘ、各要素の出力荷重ｆ_Ni、ｆ_iを式（２）の運動方程式に代入することにより、式（１３）のように、新たな加速度の推定値（ｄ／ｄｔ）²ｘ_i*を、式（１３）で求める（Ｓ１６）。
【数１３】

【００４９】
この加速度の推定値（ｄ／ｄｔ）²ｘ_i*が前の推定値（ｄ／ｄｔ）²ｘ_i（今回は、この推定値は、前の時刻と同一としている）の差が所定の閾値Ｔｏｌ以下かを判定する（Ｓ１７）。すなわち、加速度の推定値（ｄ／ｄｔ）²ｘ_i*が前の推定値（ｄ／ｄｔ）²ｘ_iとほとんど変わらないかを判定する。このＳ１７の判定で、ＮＯ、すなわち有意差があれば、（ｄ／ｄｔ）²ｘ_i*＝（ｄ／ｄｔ）²ｘ_iとする（Ｓ１８）。そして、Ｓ１３に戻り、式（１１）からの上記計算を加速度の推定値（ｄ／ｄｔ）²ｘ_i*が収束するまで繰り返す。
【００５０】
Ｓ１７の判定で、ＹＥＳであれば、加速度が収束したと判定して、時刻ｔ_iの各状態量（加速度（ｄ／ｄｔ）²ｘ_i、速度ｄｘ_i／ｄｔ、変位ｘ_i）を決定し（Ｓ１９）、時刻をｔ_i＝ｔ_i＋ｄｔとして、次のステップに進める（Ｓ２０）。
【００５１】
そして、時刻ｔ_iが終了時刻（ｅｎｄＴｉｍｅ）に至ったかを判定し（Ｓ２１）、ＮＯであればＳ１２に戻り、時刻ｔ_i＝ｔ_i＋ｄｔでの処理を行う。一方、Ｓ２１の判定で、ＹＥＳの場合に処理を終了する。
【００５２】
このようにして、各時刻における加速度、速度、変位を順次決定することができる。本実施形態は、次のような特徴がある。
【００５３】
（１）周波数領域で表現された伝達関数から、ラプラス変換とz変換の理論を利用して、離散時間領域における入力変位と出力荷重の関係式を導くことにより、運動方程式の時刻歴応答解析を行う数値積分の中で、伝達関数要素に発生する出力荷重の計算を可能にした。
【００５４】
（２）（１）の伝達関数要素の入力変位・出力荷重の関係式と、ニューマークβ法の陰解法を融合させることにより、他の非線形ばね・減衰要素の出力荷重も考慮して、運動方程式の数値積分を安定・高精度に実行できる。
【００５５】
（３）（１）（２）の具体例としては，ゴムマウントの周波数特性を考慮したエンジン・車両系の振動伝達解析、クラッチや流体継ぎ手の周波数特性を考慮したドライブトーレーン系の振動伝達解析などがある。
【００５６】
さらに、共振回路などの電気回路における振動解析にも本実施形態の手法を適用することができる。
【００５７】
「解析法の適用事例」
（ｉ）１自由度系の基本モデル
本解析法の妥当性を検証するため、図５に示す１自由度系の過渡振動を計算し、その結果を厳密解と比較した。このモデルは、１つのばねと、１つのダッシュポットを並列接続し、これに質点を接続した構成である。条件は、Ｍ＝０．１ｋｇ（重力＝０．１×９．８ｋｇｆ）、Ｃ＝０．０００２ｋｇｆｓ／ｍｍ、Ｋ＝１０ｋｇｆ／ｍｍである。また、式（４）の伝達関数の定義において、ｎ＝１，Ａ１＝０．０，Ａ０＝１．０，Ｂ１＝０．０００２，Ｂ０＝１０．０としている。また、外力は質点に作用する重力のみであり、時刻０における初期値はｘ＝ｄｘ／ｄｔ＝０である。
【００５８】
図５の計算結果において、厳密解を点線で記載したが、実線で示した本解析結果と有意の差がないため、図５において点線は見えていない。このように、過渡振動の計算結果か、本実施形態の計算により、正しい過渡振動の計算が行えることが分かる。
【００５９】
（ｉｉ）ゴム支持を模擬した簡易モデル
実際のゴム部品への適用性を確認するため、図６に示すように、質点をゴムで支持した場合の過渡振動を計算した。本モデルでは、２次の伝達関数でゴムの伝達特性（動剛性）を表現している。外力は質点に作用する重力のみである。また、時刻０における初期値はｘ＝ｄｘ／ｄｔ＝０である。また、Ｍ＝１ｋｇ（重力＝１×９．８ｋｇｆ）、２次伝達関数の各係数は、Ａ₀＝１．０，Ａ₁＝１．２０１×ｅ^-5，Ａ₂＝４．００８×ｅ^-8，Ｂ₀＝５８．７２，Ｂ₁＝４．６７７ｅ^-3，Ｂ₂＝５．１６４ｅ^-7である。
【００６０】
この計算によって、図７に示すように、ゴム部品の過渡振動を適切に計算することができる。この計算結果は、実際のゴム部品の挙動を正しく模擬していることが確認されている。
【００６１】
（ｉｉｉ）エンジンの実働振動モデル
実際の機械システムへの適用性を確認するため、図８に示すエンジンシステムの実働振動を計算した。この例は、４気筒エンジンについての振動試験である。エンジンのクランク軸には４つのピストンが接続され、シリンダ内の爆発によって、クランク軸が回転される。クランク軸の回転は、減速機を介し出力され、モータによってトルクが吸収される。エンジンブロックおよび減速機ケースがマウントゴムによって接地されている。
【００６２】
本モデルでは、図９に示すように、複数のマウントの中のエンジンブロックを接地するゴムマウントの一つ（上下方向の動剛性）に上記事例（ｉｉ）の伝達関数要素を適用し、エンジンの火入れ運転状態のマウント振動を計算した。なお、残りのゴムマウントは、端あるばね・減衰要素とした。
【００６３】
図１０に計算結果による変位ｘの時間変化（時刻歴波形）を示す。マウント振動の時刻歴波形は、エンジン１６回転分の過渡挙動を示しており、ゴムの伝達特性（減衰特性）によって振動が次第に周期定常状態に落ち着く様子を模擬している。なお、この振動計算は、特開２０００−３０５９２２号公報（特許文献５）に記載の物体変位・振動解析方法を利用した汎用のエンジン振動解析ソフトに、本発明の解析法に基づく新解析要素を実装することにより行った。
【００６４】
【発明の効果】
以上説明したように、本発明によれば、周波数領域で表現された伝達関数から、ラプラス変換とz変換の理論を利用して離散時間領域における入力と出力の関係式を導く。これによって、システムの挙動を示す微分方程式の時刻歴応答解析を行う数値積分の中で、伝達関数要素に発生する出力の計算が可能となる。従って、任意の周波数依存性を有する要素を含むシステムの周波数応答解析が可能であるとともに、過渡的な振動の解析や、非線形なばね・減衰要素を含むシステムの振動解析が可能となる。特に、機械システムにおける変位、振動について、精度の高い推定が可能になる。
【図面の簡単な説明】
【図１】 伝達関数要素を含む機械システムを示す図である。
【図２】 伝達関数の例を示す図である。
【図３】 伝達関数要素の解析の流れを示す図である。
【図４】 時刻歴応答解析の流れ図である。
【図５】 １自由度系の基本モデルの構成および過渡振動の計算結果を示す図である。
【図６】 ゴム支持を模擬した簡易モデルの構成および過渡振動の計算結果を示す図である。
【図７】 ゴム支持簡易モデルの計算結果により変位の時間変化を示す図である。
【図８】 エンジンの実働振動モデルの構成を示す図である。
【図９】 マウントの取り付けを説明する図である。
【図１０】 エンジン実働振動モデルの解析結果における変位の時間変化を示す図である。
【符号の説明】
１０ 構造体、１２ 非線形要素、１４ 伝達関数要素。[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a time history response analysis program for analyzing transient and steady behavior of a mechanical system or the like.
[0002]
[Prior art]
Conventionally, it is necessary to analyze transient and steady vibration phenomena for the purpose of optimizing the operation of mechanical systems and electronic circuit systems, and many analysis methods have been proposed. Yes.
[0003]
(1) Time history response analysis method for mechanical systems (Patent Literature 1) and electronic circuits (Patent Literature 2) Patent Literature 1 describes transient vibrations of hydraulic, control, and vehicle mechanical systems represented by ordinary differential equations. It has been shown to analyze with the integration method. Patent Document 2 discloses a numerical integration method for analyzing a transient vibration phenomenon of an electronic circuit with high speed and high accuracy.
[0004]
(2) Frequency response analysis method for linear mechanical system (Patent Documents 3 and 4)
Patent Documents 3 and 4 disclose that the frequency response of a mechanical system composed of linear elements is analyzed for vehicles and general structures. Also, a method for identifying the external force acting on the mechanical system and the mass matrix / rigidity matrix of the structure using the result of the experimental mode analysis is shown.
[0005]
(3) Commercially available mechanism / vibration analysis software and general-purpose software that performs time history response analysis of vibration / motion of mechanical systems are commercially available, and those that have a transfer function element of a second-order lag system as control elements Are known. Commercially available mechanism / vibration analysis software includes a mechanism analysis program “DADS (Dynamic Analysis and Design System) Revision 8.0” by LMS International.
[0006]
Patent Document 5 describes general-purpose engine vibration analysis software using an object displacement / vibration analysis method.
[0007]
[Patent Document 1]
Japanese Patent Laid-Open No. 2002-83252 [Patent Document 2]
JP 2001-357093 A [Patent Document 3]
JP 2002-73703 A [Patent Document 4]
JP 2001-350741 A [Patent Document 5]
JP 2000-305922 A
[Problems to be solved by the invention]
However, the conventional method described above has the following problems.
[0009]
(1) Time history response analysis method for mechanical systems (Patent Literature 1) and electronic circuits (Patent Literature 2) Transient vibration / motion analysis of mechanical systems including nonlinear springs and damping elements is possible, but rubber parts It is impossible to analyze a mechanical system including a spring / damping element having a large frequency dependency.
[0010]
(2) Frequency response analysis method for linear mechanical system (Patent Documents 3 and 4)
Although it is possible to perform frequency response analysis of mechanical systems including springs and damping elements with large frequency dependence such as rubber parts, transient vibration and motion analysis and vibration analysis of mechanical systems including nonlinear springs and damping elements are not possible. Impossible.
[0011]
(3) Commercially available mechanism / vibration analysis software If a transfer function element of a second-order lag system is used, a spring / damping element having a limited frequency characteristic can be analyzed. However, since the order is limited, it cannot cope with a spring / damping element having an arbitrary frequency characteristic.
[0012]
An object of the present invention is to enable analysis of transient and steady vibration behavior of a system in which a spring / damping element having an arbitrary frequency characteristic and a spring / damping element having an arbitrary nonlinearity coexist. is there.
[0013]
[Means for Solving the Problems]
The present invention is a time history response analysis program for analyzing the transient and steady state behavior of a system, and expressing the vibration and damping characteristics of system components having frequency characteristics as a frequency domain transfer function; Transforming the obtained frequency transfer function into a transfer function on the s-plane of Laplace transform, transforming the transfer function on the obtained s-plane into a transfer function on the z-plane of z-transform, and A system in which a transfer function on the z plane is subjected to inverse z transformation to derive a relational expression between input displacement and output in a discrete time domain, the obtained relational expression, and dynamic characteristics of other nonlinear elements are taken into consideration. And a step of numerically integrating a differential equation representing the behavior of the system in the time domain, and executing a time history response of the behavior of the system.
[0014]
Thus, according to the present invention, a relational expression between input and output in the discrete time domain is derived from the transfer function expressed in the frequency domain using the theory of Laplace transform and z transform. This makes it possible to calculate the output generated in the transfer function element in the numerical integration for performing the time history response analysis of the differential equation indicating the behavior of the system. Therefore, it is possible to analyze a frequency response of a system including an element having an arbitrary frequency dependency, and to analyze a transient vibration and a vibration analysis of a system including a nonlinear spring / damping element.
[0015]
Further, the present invention is a time history response analysis program for analyzing transient and steady vibration behavior of a dynamic system in which a spring / damping element having frequency characteristics and a spring / damping element having nonlinearity coexist. Subprogram that expresses dynamic stiffness as a frequency domain transfer function (ratio of output load to input displacement) when the spring / damping characteristics (dynamic stiffness) of a dynamic system component have frequency dependence A subprogram for converting the obtained frequency transfer function into a transfer function on the s-plane of Laplace transform, and a subprogram for converting the obtained transfer function on the s-plane into a transfer function on the z-plane of z-transform A subprogram for deriving the relational expression between the input displacement and the output load in the discrete time domain by applying an inverse z-transform to the obtained transfer function on the z-plane, and the relational expression obtained and other nonlinear elements Considering the characteristics, it has a subprogram that numerically integrates the dynamic system behavior equation in the time domain, and analyzes the vibration behavior of the dynamic system by causing the computer to execute a combination of these subprograms It is characterized by that.
[0016]
As described above, according to the present invention, the relational expression between the input displacement and the output load in the discrete time domain is derived from the transfer function expressed in the frequency domain using the theory of Laplace transform and z transform. This makes it possible to calculate the output load generated in the transfer function element in the numerical integration for performing the time history response analysis of the equation of motion. Therefore, it is possible to analyze the frequency response of mechanical systems including springs and damping elements with large frequency dependence such as rubber parts, as well as transient vibration and motion analysis, and vibrations of mechanical systems including nonlinear springs and damping elements. Analysis is possible.
[0017]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, embodiments of the present invention will be described.
[0018]
The present invention relates to a computer program, and an ordinary general-purpose computer is used to execute the program. The program can be provided by communication in addition to various recording media, and the form of provision is not limited.
[0019]
First, FIG. 1 shows a configuration diagram of a machine system to be analyzed. The structural body 10 is an elastic body or a rigid body element composed of one or a plurality of parts. The external force fe acts on the structure as a function fe (t) at time t. The nonlinear element 12 is an element such as a bearing oil film or a spring / damping element that returns an output load (reaction force) f _N to the structure 10 depending on the displacement x and speed dx / dt of the structure. The output load f _N is expressed as a non-linear function related to the displacement x of the structure 10 and the velocity dx / dt (in the figure, time differentiation is represented by dots) as in the following equation.
[Expression 1]

[0020]
The transfer function element 14 is an element that enables a new analysis in the present invention. The transfer function element 14 is a spring element (dynamic stiffness element) whose spring characteristics are expressed by a transfer function H (ω) having a frequency ω, like a rubber support part, and depending on the displacement of the structure, the output load (Reaction force) Return f to the structure 10.
[0021]
In this system, the forces applied to the structure 10 are f _e , f _N , and f, and the displacement of the structure 10 is x. Therefore, the equation of motion representing the dynamics of this system is expressed as follows.
[Expression 2]

[0022]
Here, x is a displacement vector of the structure, M is a mass matrix of the structure, C is an attenuation matrix, and K is a stiffness matrix.
[0023]
When performing a time history response analysis of this equation of motion, the point is how to calculate the output load (reaction force) f of the transfer function element 14 expressed in the frequency domain.
[0024]
“Definition of Transfer Function Elements”
The transfer function H (ω) of a rubber part or the like is defined as the ratio of the frequency spectrum X (ω), F (ω) of the input displacement x and the output load f as follows:
[Equation 3]

[0025]
Here, the frequency spectrum X (ω), F (ω) and the transfer function H (ω) are expressed as complex numbers.
[0026]
FIG. 2 shows an example of the transfer function. In general, the transfer function H (ω) of a rubber part or the like can be obtained by a vibration experiment or FEM analysis (finite element analysis).
[0027]
"Method of analyzing transfer function elements"
FIG. 3 shows a procedure for formulating the transfer function element 14.
[0028]
(I) Transfer function conversion 1 (Laplace transform converted to function on s plane)
First, the transfer function H (ω) expressed as a function of the frequency ω is converted into a function on the s plane of the Laplace transform expressed by the following equation using Equation (3).
[Expression 4]

[0029]
Here, n is the order of the transfer function, A _{n to} A ₀ and B _{n to} B ₀ are constants, and A _n ≠ 0.
[0030]
In particular, when the transfer function H (ω) is obtained in advance as a polynomial function of ω by FEM analysis or experimental mode analysis, the transfer function of Equation (4) can be easily obtained by the following conversion equation. .
[Equation 5]

[0031]
In addition, when the transfer function H (ω) is given as discrete data by an excitation experiment, each coefficient of Expression (4) is obtained by identification means such as a least square method.
[0032]
(Ii) Transfer function conversion 2 (conversion to function on z plane)
The transfer function expressed by Equation (4) is converted into a transfer function on the z plane that can express discrete time so that it can be handled by numerical integration of time history response analysis. Here, the time step of numerical integration is dt.
[0033]
By substituting the s → z conversion formula (Tustin conversion) shown in Equation (6) into Equation (4), a transfer function on the z plane shown in Equation (7) can be obtained.
[Formula 6]

[Expression 7]

[0034]
Note that the conversion to the z plane is not limited to the Tustin conversion, and Zero-order hold conversion can also be used. In the case of this zero-order hold conversion, s = (z−1) / dt is employed instead of Equation (6).
[0035]
(Iii) Recurrence formula of transfer function in time history response analysis A recurrence formula of discrete time that can be handled by numerical integration is derived from the transfer function of formula (7). Assuming that the z transformation of the input displacement x is X (z) and the z transformation of the output load f is F (z), the relationship between input and output is expressed as follows using the transfer function H (z) of Equation (7): It is expressed as follows.
[Equation 8]

[0036]
By substituting equation (7) into equation (8), equation (9) is obtained.
[Equation 9]

[0037]
Applying the inverse z transform to equation (9) yields equation (10).
[Expression 10]

[0038]
Here, x _ik and f _ik (k = 0, 1, 2,..., N) represent input and output discrete time data.
[0039]
Expression (10) indicates that the output f _{i at} time i can be calculated from the history data for the past n inputs / outputs, the input x _{i at} time _i , and the coefficient of the transfer function. Therefore, Expression (10) can be used as a recurrence expression for calculating the output load f from the input displacement x of the transfer function element in the numerical integration of the time history response analysis.
[0040]
Thus, a recurrence formula for calculating the output load f from the input displacement x is obtained by performing Laplace transform, z transform, and inverse z transform on the transfer function in the frequency domain.
[0041]
"Numerical integration of equations of motion"
Using the recurrence formula of the equation (10), the time history response analysis of the equation of motion (2) is performed. The force f _{N of the} nonlinear element and the force f of the transfer function element on the right side of Expression (2) are not forces determined explicitly by the time t, but are unknown forces that depend on the displacement and speed of the structure 10. . In order to prevent divergence of numerical integration, time history response analysis is performed by the implicit numerical integration method using the Newmark β method.
[0042]
A flowchart of the numerical integration method is shown in FIG. Let dt be the time increment of numerical integration. Assume that the acceleration, speed, and displacement up to time t _i-1 are known, and each state quantity at time t _i is obtained.
[0043]
First, the time is initially set (t _i = t ₀ ) (S11). Next, acceleration is estimated (S12). For example, as shown in FIG. 4, it is set to the same value as the acceleration (d / dt) ² x _i-1 at the previous time.
[0044]
Using the estimated value of acceleration, the displacement xi and velocity dx / dt at the next time are calculated from the next Newmark β method equation (S13).
[Expression 11]

[Expression 12]

[0045]
Where β and γ constants. For example, β = 1/6 and γ = 1/2 are set.
[0046]
Next, the output load (reaction force) f _Ni of the nonlinear element 12 of the equation (1) is calculated using the displacement and the velocity of the equations (11) and (12) (S14).
[0047]
Further, the output load f _i of the transfer function element 14 is calculated from the displacement and the past state history by the equation (10) (S15).
[0048]
Then, above the calculated velocity dx / dt, the acceleration (d / dt) ² x, output load f _Ni of each element, by substituting f _i in the equation of motion of the formula (2), formula (13 ), A new estimated value of acceleration (d / dt) ² x _i * is obtained by equation (13) (S16).
[Formula 13]

[0049]
The difference between the estimated value of acceleration (d / dt) ² x _i * is the previous estimated value (d / dt) ² x _i (this estimated value is the same as that of the previous time). It is determined whether it is equal to or less than Tol (S17). That is, it is determined whether the estimated value (d / dt) ² x _i * of acceleration is almost the same as the previous estimated value (d / dt) ² x _i . If the determination in S17 is NO, that is, if there is a significant difference, (d / dt) ² x _i * = (d / dt) ² x _i is set (S18). Then, returning to S13, the above calculation from the equation (11) is repeated until the estimated acceleration value (d / dt) ² x _i * converges.
[0050]
If YES in S17, it is determined that the acceleration has converged, and each state quantity (acceleration (d / dt) ² x _i , speed dx _i / dt, displacement x _i ) at time t _i is determined. (S19) The time is set to t _i = t _i + dt, and the process proceeds to the next step (S20).
[0051]
Then, it is determined whether the time t _i has reached the end time (endTime) (S21). If NO, the process returns to S12, and the process at time t _i = t _i + dt is performed. On the other hand, if the determination in S21 is YES, the process ends.
[0052]
In this way, the acceleration, speed, and displacement at each time can be sequentially determined. This embodiment has the following features.
[0053]
(1) Using the theory of Laplace transform and z transform from the transfer function expressed in the frequency domain, the time history response analysis of the equation of motion is performed by deriving the relational expression between the input displacement and the output load in the discrete time domain. In the numerical integration to be performed, the output load generated in the transfer function element can be calculated.
[0054]
(2) By combining the input displacement / output load relational expression of the transfer function element in (1) and the implicit method of the Newmark β method, the output load of other nonlinear springs / damping elements is also taken into account. Numerical integration of equations can be performed stably and with high accuracy.
[0055]
(3) As specific examples of (1) and (2), vibration transmission analysis of the engine / vehicle system considering the frequency characteristics of the rubber mount, and vibration transmission analysis of the drive tolane system considering the frequency characteristics of the clutch and fluid joint. and so on.
[0056]
Furthermore, the method of the present embodiment can also be applied to vibration analysis in an electric circuit such as a resonance circuit.
[0057]
"Application examples of analysis methods"
(I) Basic model of one-degree-of-freedom system In order to verify the validity of this analysis method, the transient vibration of the one-degree-of-freedom system shown in FIG. 5 was calculated and the result was compared with the exact solution. In this model, one spring and one dashpot are connected in parallel, and the mass point is connected to this. The conditions are M = 0.1 kg (gravity = 0.1 × 9.8 kgf), C = 0.0002 kgfs / mm, and K = 10 kgf / mm. Further, in the definition of the transfer function of Expression (4), n = 1, A1 = 0.0, A0 = 1.0, B1 = 0.0002, and B0 = 10.0. The external force is only gravity acting on the mass point, and the initial value at time 0 is x = dx / dt = 0.
[0058]
In the calculation result of FIG. 5, the exact solution is indicated by a dotted line, but since there is no significant difference from the analysis result indicated by the solid line, the dotted line is not visible in FIG. Thus, it can be seen that the correct transient vibration can be calculated from the calculation result of the transient vibration or the calculation of the present embodiment.
[0059]
(Ii) Simple model simulating rubber support In order to confirm the applicability to an actual rubber part, as shown in FIG. 6, the transient vibration when the mass point was supported by rubber was calculated. In this model, the transfer characteristic (dynamic rigidity) of rubber is expressed by a secondary transfer function. The external force is only gravity acting on the mass point. The initial value at time 0 is x = dx / dt = 0. Further, M = 1 kg (gravity = 1 × 9.8 kgf), the coefficients of the secondary transfer function are A ₀ = 1.0, A ₁ = 1.201 × e ⁻⁵ , A ₂ = 4.008 × e. ⁻⁸ , B ₀ = 58.72, B ₁ = 4.677e ⁻³ , B ₂ = 5.164e ⁻⁷ .
[0060]
By this calculation, as shown in FIG. 7, the transient vibration of the rubber part can be appropriately calculated. It has been confirmed that this calculation result correctly simulates the behavior of an actual rubber part.
[0061]
(Iii) Actual Vibration Model of Engine In order to confirm applicability to an actual mechanical system, the actual vibration of the engine system shown in FIG. 8 was calculated. This example is a vibration test for a four cylinder engine. Four pistons are connected to the crankshaft of the engine, and the crankshaft is rotated by explosion in the cylinder. The rotation of the crankshaft is output via a reduction gear, and the torque is absorbed by the motor. The engine block and reducer case are grounded by the mount rubber.
[0062]
In this model, as shown in FIG. 9, the transfer function element of the above case (ii) is applied to one of the rubber mounts (vertical dynamic rigidity) for grounding the engine block among the plurality of mounts, The mount vibration in the fired operation state was calculated. The remaining rubber mount was a spring / damping element with ends.
[0063]
FIG. 10 shows the time change (time history waveform) of the displacement x according to the calculation result. The time history waveform of the mount vibration shows a transient behavior for 16 revolutions of the engine, and simulates how the vibration gradually settles into a periodic steady state by the transfer characteristic (damping characteristic) of rubber. This vibration calculation is performed by adding a new analysis element based on the analysis method of the present invention to general-purpose engine vibration analysis software using the object displacement / vibration analysis method described in JP 2000-305922 A (Patent Document 5). Done by implementing.
[0064]
【The invention's effect】
As described above, according to the present invention, a relational expression of input and output in the discrete time domain is derived from the transfer function expressed in the frequency domain using the theory of Laplace transform and z transform. This makes it possible to calculate the output generated in the transfer function element in the numerical integration for performing the time history response analysis of the differential equation indicating the behavior of the system. Therefore, it is possible to analyze a frequency response of a system including an element having an arbitrary frequency dependency, and to analyze a transient vibration and a vibration analysis of a system including a nonlinear spring / damping element. In particular, it is possible to estimate the displacement and vibration in the mechanical system with high accuracy.
[Brief description of the drawings]
FIG. 1 shows a mechanical system including a transfer function element.
FIG. 2 is a diagram illustrating an example of a transfer function.
FIG. 3 is a diagram illustrating a flow of analysis of a transfer function element.
FIG. 4 is a flowchart of time history response analysis.
FIG. 5 is a diagram showing a configuration of a basic model of a one-degree-of-freedom system and a calculation result of transient vibration.
FIG. 6 is a diagram showing a configuration of a simple model simulating rubber support and a calculation result of transient vibration.
FIG. 7 is a diagram showing a change with time of displacement based on a calculation result of a rubber support simple model.
FIG. 8 is a diagram showing a configuration of an actual engine vibration model.
FIG. 9 is a diagram illustrating mounting of a mount.
FIG. 10 is a diagram showing a change with time of displacement in an analysis result of an engine actual vibration model.
[Explanation of symbols]
10 structure, 12 nonlinear element, 14 transfer function element.

Claims (2)

A time history response analysis program for analyzing the transient and steady state behavior of a system,
Expressing the vibration and damping characteristics of the components of the system having frequency characteristics as a transfer function in the frequency domain;
Transforming the obtained frequency transfer function into a transfer function on the s-plane of Laplace transform;
Converting the obtained transfer function on the s-plane to a transfer function on the z-plane of the z-transform,
Applying inverse z-transform to the obtained transfer function on the z-plane to derive a relational expression between input displacement and output in the discrete time domain;
Considering the obtained relational expression and the dynamic characteristics of other nonlinear elements, numerically integrating a differential equation indicating the behavior of the system in the time domain,
Is a time history response analysis program that analyzes the time history response of the behavior of the system. A time history response analysis program that analyzes the transient and steady vibration behavior of a dynamic system in which a spring / damping element with frequency characteristics and a spring / damping element with nonlinearity coexist,
A subprogram that expresses dynamic stiffness as a frequency domain transfer function (ratio of output load to input displacement) when the spring / damping characteristics (dynamic stiffness) of a dynamic system component are frequency dependent. ,
A subprogram for converting the obtained frequency transfer function into a transfer function on the s plane of the Laplace transform;
A subprogram for converting the obtained transfer function on the s-plane into a transfer function on the z-plane of the z-transform,
A subprogram that performs inverse z-transform on the obtained transfer function on the z-plane to derive a relational expression between input displacement and output load in a discrete time domain;
A subprogram that numerically integrates the dynamic system behavior equation in the time domain, taking into account the obtained relational expression and the dynamic characteristics of other nonlinear elements,
Have
A time history response analysis program that analyzes the vibration behavior of a dynamic system by causing a computer to execute a combination of these subprograms.

Images