WO2022065500A1

WO2022065500A1 - 振動解析方法、プログラム、記憶媒体

Info

Publication number: WO2022065500A1
Application number: PCT/JP2021/035491
Authority: WO
Inventors: 孝広近藤; 大樹住川
Original assignee: 国立大学法人九州大学
Priority date: 2020-09-28
Filing date: 2021-09-28
Publication date: 2022-03-31
Also published as: US20230366774A1; EP4220105A1; JP2022054596A

Abstract

この振動解析方法は、局所的に非線形性を有する大規模系の振動を解析する方法であって、線形状態量に対する方程式に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式から線形状態量の高次モードの影響を補正して消去する工程（１）と、前記実モード方程式から、元の大規模系の解に対する影響の大きなモードである支配的モードを抽出し、影響の小さなモードである副次的モードについては、近似解を利用してその影響を前記非線形状態量に対する方程式に補正項として取り込んで消去し、低次元モデルを導出する工程（２）と、前記低次元モデルを利用して周波数応答を計算する工程（３）と、を備える。

Description

振動解析方法、プログラム、記憶媒体

　本発明は、主に局所的強非線形性を有する大規模系の高性能振動解析手法に関する。より具体的には、新型複素モード解析を利用した低次元化法に関する。本願は、２０２０年９月２８日に、日本に出願された特願２０２０－１６１７２７号に基づき優先権を主張し、その内容をここに援用する。

　機械は、互いに相対運動する多数の要素から構成されており、基礎支持部や要素間の結合部に、軸受や歯車等のような強い非線形性を有する要素を持っている。しかも、設計時の応力解析や振動解析においては、有限要素法などを利用することによって、実機に近い精密な大規模自由度系としてモデル化されることが多い。したがって、ほぼ全ての機械は、局所的に強い非線形性を有する大規模自由度系であるといえる。

　一方、近年、機械に対する高性能化の要求が益々強まっている。その要求を実現するには、設計段階において系内に不可避的に含まれる非線形性の影響を十分に考慮した振動解析を実施する必要がある。しかも、複雑多岐にわたる非線形振動の全貌を明らかにするには、単に特定の系パラメータに対して数値シミュレーションを行うだけでは不十分で、安定判別を含めた周波数特性の全容を解明しなければならない。

　これまで、多自由度非線形系の振動解析には、アルゴリズムに規則性があり幅広い問題に適用可能であることから、線形固有モード展開に基づくガラーキン法（例えば、下記特許文献１参照）がよく用いられて来た。

日本国特開２０１７－２１１８７７号公報

　しかしながら、ガラーキン法には、採用モード数を増やして高精度の解析を行おうとすると計算時間やメモリ量が爆発的に増大することや、非線形性が局所的にしか存在しないモデルであってもモード座標に変換後は非線形性が系全体に分散しモード分離の効果がなくなることなどの重大な欠点がある。

　このような状況を打開するために、少数のモードによる低次元モデルを用いて高精度の解析を行うための研究が精力的になされてきたが、未だ根本的な解決には至ってない。その主な問題点としては、採用するモードを適切に選択しないと、非線形系に特有のモード間の連成作用により、解析の精度（とりわけ安定判別の精度）が著しく悪化することが挙げられる。このように、大規模非線形系に対する実用的な振動解析手法はまだ存在しないといっても過言ではなく、この問題に対処し得る高性能な振動解析手法を開発することが、機械設計に関わる数値計算上の喫緊の課題となっている。

　本発明は、背景で述べた問題の抜本的な解決を目指したものであり、局所的に強い非線形性を有する大規模系を対象に、非線形性の影響を適切に反映した高精度の低次元モデルを合理的に構成する手法を提案している。

　提案手法は、主に次の二つの工程から構成される。
（１）線形状態量に対する方程式［後述の式（１）および式（４）］に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式［後述の式（２）］から線形状態量の高次モードの影響を補正して消去する工程。
（２）低次モードに対する実モード方程式から元の大規模系の解に対する影響の大きなモード（支配的モード）を抽出し、影響の小さなモード（副次的モード）については近似解を利用してその影響を非線形状態量に対する方程式［後述の式（２）］に補正項として取り込んで消去する低次元モデル導出工程、および周波数応答を計算する過程で低次モードの中から支配的モードを（なるべく合理的に）抽出する工程。

　この低次元モデルを利用すれば、これまでは不可能であった大規模非線形系の安定判別を含む振動解析を高速かつ高精度で行うことが可能になる。しかも、提案手法は非常に高い汎用性を有しており、局所的に強い非線形性を有する大規模振動系全般に適用可能である。

本解析手法の概要を示す図である。本解析手法のフローチャートである。低次元化法適用例の解析モデル概略図である。フルモデルにおける１次から５次の主共振領域の周波数応答図である。フルモデルにおける１次から２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝５９、ｎ_ａ＝５）における１次から５次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝５９、ｎ_ａ＝５）における１次から２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝５９、ｎ_ａ＝１５）における１次から５次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝５９、ｎ_ａ＝１５）における１次から２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝５、ｎ_ａ＝５）における１次から５次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝１０、ｎ_ａ＝５）における１次から５次の主共振領域の周波数応答図である。低次モード数ｎ_ｌの１次の固有振動数への影響を示す図である。低次モード数ｎ_ｌの２次の固有振動数への影響を示す図である。低次モード数ｎ_ｌの３次の固有振動数への影響を示す図である。低次モード数ｎ_ｌの４次の固有振動数への影響を示す図である。低次モード数ｎ_ｌの５次の固有振動数への影響を示す図である。フルモデルにおける２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝２０，δ＝５×１０^－４）における２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝２０，δ＝１×１０^－４）における２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝２０，δ＝１×１０^－５）における２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝１５，δ＝７×１０^－４）における１次から２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝１５，δ＝１×１０^－４）における１次から２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝１５，δ＝２×１０^－５）における１次から２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝５９９，ｎ_ａ＝２０）における１次から２次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝３０，δ＝１×１０^－３）における４次から５次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝３０，δ＝５×１０^－４）における４次から５次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝３０，δ＝２×１０^－４）における４次から５次の主共振領域の周波数応答図である。低次元モデル（ｎ_ｌ＝５９９，ｎ_ａ＝２０）における４次から５次の主共振領域の周波数応答図である。本解析手法を実施する装置の一例を示すブロック図である。

　以下、上述した二つの工程における具体的な実施事項について説明する。
　なお、本明細書において、数式以外の本文中では、項や数式の表現が限られる。このため、例えば、本明細書に画像として取り込まれて項（０－１）のように表される項を、本明細書の本文中では、「^ＣΨ^～ _１ ^Ｔ」と表記する。本明細書に画像として取り込まれて項（０－２）のように表される項を、本明細書の本文中では、「ｙ^・」と表記する。なお、「ｙ^・」は、ｙを時間により１階微分した項を意味する。

〔解析対象および基礎式〕
　本発明で提案する手法は、全自由度がＮで局所的に強い非線形性を有する大規模自由度の振動系を解析対象とする。その基礎式は、一般的に次のように表すことができる。

　ｙ_１は、非線形力が作用しない状態量の物理座標をまとめたｎ次元ベクトル、ｙ_２は、非線形力が作用する状態量の物理座標をまとめたｍ次元ベクトルである。以下では、ｙ_１を線形状態量、ｙ_２を非線形状態量と呼び、便宜的に式（１）を線形状態量の方程式、式（２）を非線形状態量の方程式と呼ぶ。ｎ_２は、ｙ_２およびｙ^・ _２に関する強非線形関数であり、ｆ_１，ｆ_２は、周期２π／ωの調和強制外力とする。Ｍは、質量要素に関する係数行列（質量行列）である。Ｃは、減衰要素に関する係数行列（減衰行列）である。Ｋは、ばね要素に関する係数行列（剛性行列）である。提案手法の汎用性を確保するため、Ｍ_１１，Ｃ_１１，Ｋ_１１に対して正定性を仮定しない。ｍ＜＜ｎであるから、非線形性は局所的に存在することが分かる。このように、提案手法は、局所的に強い非線形性を有する大規模振動系全般に適用可能である。

　低次元化法の定式化にあたり、式（１）の係数行列を次のような転置行列に置き換え、外力を零とした系を同時に取り扱う。

　ここに、右上添え字「Ｔ」は転置記号である。式（４）のｙ_２は、式（１）および式（２）を満足する解であるが、ｙ_１は、式（１）および式（２）を満足する解とは一般に異なるので、式（４）ではｙ^＊ _１と表記している。以下では、式（４）を式（１）の双対系と呼ぶ。

＜＜線形状態量の方程式に対する新型複素モード解析の適用＞＞〔１階の微分方程式への変換〕
　新型複素モード解析に基づく低次元化法を定式化するため、式（１）および式（４）を、次のような１階の微分方程式に変換する。なおＩ_ｎは、ｎ次の単位行列である。

〔一般固有値問題〕
　本低次元化法では、まず、式（５）および式（６）でｘ^・ _２＝０、ｘ_２＝０と置くことによって非線形状態量を固定し、ｇ_１＝０とした自由振動系から、次のような一般固有値問題を導出する。

　式（８）および式（９）から求められる固有値λは一致し、一般に複素数になる。また、対応する固有ベクトルＸ_１およびＸ^＊ _１も複素ベクトルになる。式（８）を基準にすれば、Ｘ_１が右複素固有ベクトル、Ｘ^＊ _１が左複素固有ベクトルとなり、式（９）を基準にすれば、Ｘ^＊ _１が右複素固有ベクトル、Ｘ_１が左複素固有ベクトルとなる。振動学的には、固有値λの虚部が固有角振動数に対応し、Ｘ_１およびＸ^＊ _１は左右の複素拘束モードに対応する。

〔低次モードの複素拘束モード〕
　第１段階として、式（８）および式（９）から複素固有値と左右の複素固有ベクトル（複素拘束モード）を求める。ところが、一般固有値問題の次元２ｎが大きいとき、すべての固有値と固有ベクトルを求めるのは困難である（大規模一般固有値問題の解法は、数値計算分野に残る未解決の難問であり、事実上不可能に近い）。ただし、低次から順に（絶対値の小さいものから順に）比較的少数の固有値と対応する固有ベクトルを求めることは、現時点でも可能である。
　一方、線形系・非線形系ともに、モード次数が高くなるにつれて振動特性に及ぼす影響は一般に小さくなることが知られている。

　そこで、本低次元化法では、解の精度に大きな影響を及ぼす可能性のある１次からｎ_ｌ（＜＜ｎ）次までの比較的低次の固有値と固有ベクトルのみを、一般固有値問題から求め、ｎ_ｌ＋１次以上の高次モードの影響に関しては、ｎ_ｌ次までの低次モードを用いて近似する。ｎ_ｌの値については、解析すべき最大振動数の数倍程度の振動数範囲内に存在する固有振動数（固有値の虚部）の最大次数を目安に設定すれば良い。

　この手続きに必要な関係式は次の通りである。ただし、下添え字「ｌ」および「ｈ」はそれぞれ低次モードおよび高次モードに関連する変数または物理量であることを示す。また、左上添え字「Ｃ」が付された物理量は複素数であることを示し、同じく「Ｒ」および「Ｉ」は原則としてそれぞれその実部および虚部を表す。
　まず、式（８）の１次からｎ_ｌ（＜＜ｎ）次までの固有値と対応する左右の固有ベクトルを求める。これらの固有値は、全て実部が負の複素数であるものとして、ｐ（＝１，‥，ｎ_ｌ）次の複素共役な固有値を次のように表す。

　ここに、ω_ｐ，１およびζ_ｐ，１はそれぞれｐ次モードの不減衰固有角振動数および減衰比を表す。さらに、これら２個ｎ_ｌ組（計２ｎ_ｌ個）の固有値をひとまとめにして、次のように行列表示する。

　ここに、diag[ ]は対角行列を表す。
　この複素固有値に対応する式（８）のｎ_ｌ組（２ｎ_ｌ本）の右複素固有ベクトルをひとまとめにした右複素拘束モード行列を、^ＣΦ^＾ _１ｌとする。同じく、ｎ_ｌ組（２ｎ_ｌ本）の左複素固有ベクトルをひとまとめにした左複素拘束モード行列を、^ＣΨ^＾ _１ｌとする。ただし、^ＣΦ^＾ _１ｌと^ＣΨ^＾ _１ｌは、次式を満足するように正規化されているものとする。

　このとき、^ＣΦ^＾ _１ｌ，^ＣΨ^＾ _１ｌはそれぞれ次のように表される。なお^Ｒφ_ｐ，１，^Ｉφ_ｐ，１は、それぞれｐ次モードの右複素固有ベクトルにおける実部、虚部である。また、^Ｒψ_ｐ，１，^Ｉψ_ｐ，１は、それぞれｐ次モードの左複素固有ベクトルにおける実部、虚部である。

〔低次モードの複素拘束モードの実数化〕
　計算の高速化のために、低次モードに対する左右の複素拘束モード行列^ＣΦ^＾ _１ｌおよび^ＣΨ^＾ _１ｌから、次のようにして右実拘束モード行列^ＲΦ^＾ _１ｌ、および左実拘束モード行列^ＲΨ^＾ _１ｌを導出する。

〔低次モードに対する実モード座標の導入〕
　線形状態量の物理座標ｘ_１およびｘ^＊ _１は、次のように低次モードに基づくｘ_１ｌおよびｘ^＊ _１ｌと高次モードに基づくｘ_１ｈおよびｘ^＊ _１ｈの和で表すことができる。なおｙ_１ｌおよびｙ^＊ _１ｌは、それぞれ低次モードにおける線形状態量である。ｙ_１ｈおよびｙ^＊ _１ｈは、それぞれ高次モードにおける線形状態量である。

　このうち、ｘ_１ｌとｘ^＊ _１ｌに対しては、それぞれ^ＲΦ^＾ _１ｌおよび^ＲΨ^＾ _１ｌを介して、対応する実モード座標η_１ｌおよびη^＊ _１ｌを次式により導入する。なおξ_１ｌは、低次モードの物理座標に対する実モード座標である。

〔低次モードに対する実モード方程式の導出〕
　式（１５）および式（１７）を式（５）に代入した上で左から^ＲΨ^＾ _１ｌ ^Ｔを乗じて整理すると、次のような、低次モードに対する実モード方程式が導出される。

　双対系の式（６）に対しても、式（１６）および式（１８）を代入した上で左から^ＲΦ^＾ _１ｌ ^Ｔを乗じて整理すると、次のような実モード方程式が導出される。

　このように、係数行列が正定性を満たさない場合でも、質量行列、減衰行列および剛性行列を同時に対角化した上で、実数型２階微分方程式の形でモード方程式を導出することができる。これが新型複素モード解析法の最大の特長である。

〔低次モードを利用した高次モードの影響の補正〕
　非線形状態量の方程式（２）中の線形状態量ｙ_１を、式（１５）に示すように低次モードに基づくｙ_１ｌと高次モードに基づくｙ_１ｈの和で表し、ｙ_１ｌに関しては、式（１７）で導入した低次モードに対応する実モード座標に変換する。さらに、高次モードに基づくｙ_１ｈの影響を補正して省略すると、次式が求められる。

　ここに、Ｍ^～ _２２ｈ，Ｃ^～ _２２ｈ，Ｋ^～ _２２ｈ，ｐ_２１ｈ，ｑ_２１ｈが高次モードの補正項を表す。これらは、いずれも次のように低次モードから求められる。

＜＜低次元モデルの導出＞＞
〔支配的モードと副次的モードの分離〕
　式（１７）および式（１８）で導入した実モード座標η_１ｌおよびη^＊ _１ｌの中には、元の大規模系の解の精度に及ぼす影響が大きなモード座標（支配的モードと呼ぶ）と影響の小さなモード座標（副次的モードと呼ぶ）が存在する。下添え字「ａ」および「ｂ」により、それぞれ支配的モードおよび副次的モードに関連する変数または物理量であることを示す。例えば、η_１ａは支配的モード座標を、η_１ｂは副次的モード座標を表す。

　式（１４）で求めた低次モードに対する左右の実モード行列^ＲΦ^＾ _１ｌ，^ＲΨ^＾ _１ｌを、ｎ_ａ組（２ｎ_ａ本）の支配的モードからなるモード行列^ＲΦ^＾ _１ａ，^ＲΨ^＾ _１ａと、ｎ_ｂ組（２ｎ_ｂ本）の副次的モードからなるモード行列^ＲΦ^＾ _１ｂ，^ＲΨ^＾ _１ｂとに分離し（ｎ_ｌ＝ｎ_ａ＋ｎ_ｂ）、次のように再整理する。ただし、以下では、下添字「＃」は「ａ」または「ｂ」を表す。

　実モード行列の分離に対応して、式（１９）、式（２０）、式（２１）のモード座標も、次のように支配的モードと副次的モードに分離する。

〔副次的モードの近似解〕
　副次的モードに対しては非線形性の影響が小さいのでこの影響を無視し、ｙ_２，ξ_１ｂ，ξ^＊ _１ｂが基本周期２π／ωで調和振動しているものと仮定すると、ｙ^・・ _２＝－ω^２ｙ_２，ξ^・・ _１ｂ＝－ω^２ξ_１ｂ，ξ^・・＊ _１ｂ＝－ω^２ξ^＊ _１ｂと表される。さらに、調和外力ではｆ^・・ _１＝－ω^２ｆ_１であることを考慮すると、式（２４）および式（２５）より、副次的モードの近似解が次のように求められる。

〔副次的モードの消去と低次元モデルの導出〕
　式（２６）中の副次的モード座標ξ_１ｂ，ξ^・ _１ｂ，ξ^・・ _１ｂに式（２７）の近似解を代入し、副次的モードの影響を補正して省略すると、次式を得る（この導出過程で式（２８）の結果を利用している）。

　ここに、Ｍ^～ _２２ｂ，Ｃ^～ _２２ｂ，Ｋ^～ _２２ｂ，ｐ_２１ｂ，ｑ_２１ｂが副次的モードの補正項を表し、次式で与えられる。

　さらに、式（２４）で下添え字「＃」を「ａ」とした式と式（２９）から、次のように（ｎ_ａ＋ｍ）次元に低次元化された方程式が求められる。

　式（３１）で表されるモデルを低次元モデルと呼ぶ。式（３１）を利用して近似解を計算する。

〔線形状態量の物理座標の近似解〕
　低次元モデルの式（３１）から近似解ξ_１ａ，ξ^・ _１ａ，ｙ_２，ｙ^・ _２が求められたとき、必要に応じて線形自由度の近似解を物理座標に変換する。
　この手続きに必要な関係式は次の通りである。

　ｙ_１ａ，ｙ^・ _１ａは支配的モードに基づく物理座標を、ｙ_１ｂ，ｙ^・ _１ｂは副次的モードに基づく物理座標を、ｙ_１ｈ，ｙ^・ _１ｈは高次モードに基づく物理座標を表す。

〔支配的モード座標の抽出法〕
　非線形系において外力の角振動数ωを徐々に変化させながら複雑多岐にわたる周波数応答の全容を解明しようとする場合には、あるωに対して求められた近似解をω±Δωに対する初期値として逐次近似計算することが多い。そこで、あるωにおいて求められた近似解に含まれる情報を利用して、ω±Δωで抽出すべきξ_１ａを選定する。
　この手続きに必要な関係式は次の通りである。
　いま、あるωに対して、ξ_１ａ，ξ^・ _１ａ，ｙ_２，ｙ^・ _２からなる低次元モデルの近似解が求められたとする。この近似解ｙ_２，ｙ^・ _２から、次式により低次モード座標ξ_１ｌの近似解が計算できる。

　この近似解のモード別の振幅を次式で定義できる。

　ω±Δωの計算においては、このモード振幅｜｜ξ_１，ｐ｜｜が、設定した閾値δよりも大きいモードのみをξ_１ａとして抽出する。抽出された支配的モードの本数をｎ_ａ、残りの副次的モードの本数をｎ_ｂとすると、ｎ_ｌ＝ｎ_ａ＋ｎ_ｂである。場合によっては、ｎ_ａの上限値ｎ_{ａ，ｍａｘ}を定めることもある。
　この手法により、解の特性に対応してなるべく少ない本数の支配的モードを適切に選択しながら周波数応答を効率的に計算することが可能となる。この手法ではモード振幅の閾値δおよび支配的モード本数の上限値ｎ_{ａ，ｍａｘ}の設定法が問題となるが、実用的にはδについては徐々に小さくしながら、ｎ_{ａ，ｍａｘ}については徐々に大きくしながら複数の周波数応答を求め、解析結果が変化しなくなった時点で高精度の解が求められたものとみなせばよい。低次元モデルの計算はかなり高速で行えるので、実用上は複数の計算にかかるコストは大きな問題にはならないと考えられる。なお、解析結果が変化しなくなった時点とは、解析結果が全く変化しない場合だけではなく、解析結果を比較した結果、その差が、所定の基準を下回った場合を含む。
　なお、計算を開始する最初のωにおける支配的モード本数ｎ_ａの設定方法については、例えば、ｎ_ａ＝ｎ_ｌ（低次モードの本数）とすればよい。

＜＜周波数応答の計算手順＞＞
　以下に、周波数応答の計算手順を簡潔にまとめる（図１、図２参照）。
（Ｓｔｅｐ１）前処理
　（Ｓｔｅｐ１－１）
　　線形状態量に対する低次モードの複素固有値と複素固有ベクトルの計算［式（８），式（１０），式（１２），式（１３）］および複素固有ベクトルの実数化［式（１４）］
　（Ｓｔｅｐ１－２）
　　低次モードの物理座標に対する実モード座標の導入［式（１７），式（１８）］および実モード方程式の導出［式（１９），式（２０）］
　（Ｓｔｅｐ１－３）
　　非線形状態量の方程式に対する低次モードを利用した高次モードの影響の補正［式（２１）］
（Ｓｔｅｐ２）支配的モードの抽出［式（３６），式（３７）］（Ｓｔｅｐ３）低次元モデルの導出
　（Ｓｔｅｐ３－１）
　　支配的モードと副次的モードの分離［式（２３），式（２４），式（２５），式（２６）］
　（Ｓｔｅｐ３－２）
　　副次的モードの近似解の計算［式（２７），式（２８）］
　（Ｓｔｅｐ３－３)
　　低次元モデルの導出［式（２９），式（３０），式（３１）］（Ｓｔｅｐ４）低次元モデル［式（３１）］の周期解の計算（Ｓｔｅｐ５）計算終了・継続の判定
　角振動数ωが解析範囲の上限または下限に到達すれば計算を終了する。到達していなければ、ω±Δωをωに置き換え、Ｍ_１１，Ｃ_１１，Ｋ_１１がωの関数である場合にはＳｔｅｐ１に戻り、定数の場合にはＳｔｅｐ２に戻り計算を継続する。

＜＜低次元化法適用例＞＞
〔１．解析モデル〕
　新型複素モード解析を利用した低次元化法の有効性、特に、高次モード消去法および支配的モードの抽出法の有効性を確認するため、直線状梁構造物の面内曲げ振動を対象に本低次元化法を適用する。解析モデルの模式図を図３に、計算に用いたパラメータを表１に示す。

　このモデル１０は、長さ１．５ｍ、外径３０ｍｍ、内径２５ｍｍの一様中空円形断面の梁１の両端を、非線形要素２（粘性減衰器と漸硬型３次ばね）によって支持したものであり、左端から０．５ｍの位置には、断片線形ばね（ガタ）３を配置している。また、外力の振動数が比較的高い領域においても、振幅の大きな振動が発生して非線形性の影響が強く表れるように、遠心力タイプの調和強制外力Ｆを左端から１ｍの位置に作用させている。
　梁１の物性値は、鉄鋼材料を想定している。梁１については３０個の微小要素に等分割し、有限要素法により質量行列、剛性行列を導出した。この場合の自由度は、線形状態量がｎ＝５９、非線形状態量がｍ＝３であり、総自由度がＮ＝６２となる。

　これは大規模系と呼ぶには相応しくない自由度であるが、低次元化法を適用しない場合の解を求めることのできる限界に近い自由度である。新型複素モード解析の長所が現れるように、減衰行列は分布減衰を与えて非比例減衰系としている。表２に本解析モデルにおいて漸硬型ばね定数を０として線形化した系の１次から１０次の不減衰固有振動数を示す。

　図４から図９に本解析モデルの周波数応答図を示す。図４、図６、図８は１次から５次までの主共振領域を、図５、図７、図９は１次および２次の主共振領域の拡大図を示している。図４および図５は低次元化法を適用しない場合（フルモデル）の結果、図６から図９は低次元化法を適用した場合（低次元モデル）の結果である。

　低次元モデルの構成に際しては高次モード消去法および支配的モードの抽出法は適用せず、外力の振動数が拘束モードの固有振動数に近い順に一定本数（ｎ_ａ＝５、１５）を抽出した。低次元モデルの結果では、支配的モード本数（数値は右軸）を、塗り潰した▽印で示したが、図６から図９では一定であることを表している（結果として太い実線のように見えている）。

　周波数応答における応答曲線はすべてガタを配置した点の最大片振幅を描いており、振幅がガタの隙間（０．００２５ｍ）を超えると応答曲線が屈曲している。図中の実線は安定解、破線は不安定解、〇印はサドルノード分岐点、□印はピッチフォーク分岐点、△印はホップ分岐点を表す。

　ｎ_ａ＝５の場合の図６および図７とフルモデルの図４および図５を比較すると概ね良い一致を示しているが、強い非線形性であるガタの影響が現れる領域で応答曲線の形状や振幅のピーク値、安定判別結果が異なっている。図８および図９に示すようにｎ_ａ＝１５まで支配的モード数を増加させるとフルモデルと一致した結果が得られる。

〔２．高次モード消去法の適用〕
　本解析モデルに高次モード消去法を適用して、解析に用いる低次モード数ｎ_ｌの解析結果への影響を調査した。図１０および図１１は支配的モード数をｎ_ａ＝５と固定して低次モード数をｎ_ｌ＝５，１０とした場合の結果である。高次モード消去法を適用しない場合の図６（ｎ_ｌ＝ｎ＝５９）と比較すると、図１０のｎ_ｌ＝５では１次、２次の主共振領域ではほぼ変わらない応答が得られているが、振動数が高くなるにつれてずれが大きくなる。図１１のｎ_ｌ＝１０とすればこの解析範囲においてはほぼ一致した結果が得られていることがわかる。

　このように解析を行う振動数の範囲によって必要な低次モード数は変化し、低次モード数を増加させるとより高い振動数範囲まで高精度に解が求まる。この解析モデルでは５９本中の１０本と少数の低次モードのみで高次モードの影響を高精度に推定できていることから、大規模系の解析に対してこの性質は極めて有効であることが期待される。

　図１２から図１６には低次モードの数の固有振動数への影響を示しており、図１２から図１６は解析振動数範囲の１次から５次の固有振動数である。これらの固有振動数は、高次モード消去法のみを適用して作成した低次元モデル（ｎ_ｌ＝ｎ_ａ）から求めた。点は低次元モデル、破線はフルモデルの固有振動数である。

　図１２から図１６の比較から、低次の固有振動数ほど精度がよく高次になるにつれて低精度である。また、低次モード数の増加に伴って精度が向上することがわかる。このように低次モード数に対して固有振動数の精度と周波数応答解析結果の精度が同様の傾向を示しており、低次モード数ｎ_ｌの値の設定法としては、計算コストの観点から、固有振動数の精度を目安にすればよいと考えられる。

〔３．支配的モードの抽出法の適用〕
　本解析モデルに対して支配的モードの抽出法を適用した。解析範囲は２次の主共振領域に限定した。低次元化法には高次モード消去法も併せて適用しており、低次モード数はｎ_ｌ＝２０とした。図１７はフルモデルの結果、図１８から図２０はそれぞれ閾値をδ＝５×１０^－４，１×１０^－４，１×１０^－５として抽出法を適用した結果であり、支配的モード数が外力の振動数に対して変化していることがわかる。図１８から図２０にかけて閾値を徐々に小さくしていくと、モード数が増加し解析に必要なモードが抽出されることでフルモデルの解に近づいていくことがわかる。図２０のδ＝１×１０^－５ではフルモデルと一致した結果が得られている。

　いずれの閾値においても、解の振幅が小さく非線形性の影響が表れていない領域では支配的モード数は比較的少なく、振幅が大きくなるにつれて支配的モード数が増加しており、解の特性に応じて支配的モードを変化させながら効率的に解析できていることが確認された。

〔４．自由度の大きいモデル〕
　本解析モデルにおいて分割数を３００としたモデルに対して、本手法を適用して有効性の検証を行った。この場合の自由度は、線形状態量がｎ＝５９９、非線形状態量がｍ＝３であり、総自由度がＮ＝６０２となる。この自由度ではフルモデルに対する解析は困難となる。解析は１次から２次の主共振領域と４次から５次の主共振領域にわけて行った。

〔４．１．１次、２次の主共振領域〕
　解析振動数範囲は、固有振動数および１章〔１．解析モデル〕の結果から、外力の振動数を２０Ｈｚから２００Ｈｚまでとした。高次モード消去法を適用する際の低次モード数は、最大振動数の２００Ｈｚの５倍程度の６次の固有振動数を目安にｎ_ｌ＝１５とした。
　図２１から図２３はそれぞれ閾値をδ＝７×１０^－４，１×１０^－４，２×１０^－５と徐々に小さくしながら支配的モードの抽出法を適用した結果である。図２１と図２２を比較すると、１次および２次共振の応答に差異がみられる。図２２と図２３を比較すると、応答の変化は小さくなったが、２次のピーク付近の解がわずかに異なっていることがわかる。

　図２３の結果が高精度な結果であることを確認するため、高次モード消去法および支配的モードの抽出法は適用せずに１章と同様にｎ_ａを徐々に増加させながら低次元モデルの解析を行い、応答に変化がなくなった時点での結果を本解析モデルの正解とみなした。その結果を図２４に示しており、低次モード数がｎ_ｌ＝５９９、支配的モード数がｎ_ａ＝２０である。この結果は、図２３と非常に良い一致を示していることがわかる。

〔４．２．４次、５次の主共振領域〕
　解析振動数範囲は４５０Ｈｚから１０００Ｈｚまでとして、低次モード数は最大振動数の１０００Ｈｚの３倍程度の１０次の固有振動数を目安にｎ_ｌ＝３０とした。前節と同様に閾値を徐々に小さくしながら支配的モードの抽出法を適用した。図２５から図２７はそれぞれ閾値をδ＝１×１０^－３，５×１０^－４，２×１０^－４とした結果である。
　前節と同様に図２５から図２７を比較すると低次元モデルの応答の変化が小さくなっていくことがわかる。また、図２７と、図２８のｎ_ｌ＝５９９、ｎ_ａ＝２０の結果との比較から、この振動数範囲においても本低次元化法によってフルモデルとほとんど変わらない結果を効率的に計算できている。ただし、解析を行う振動数範囲によって低次モード数や最適な閾値は異なる。
　以上のように、本解析モデルに対する具体的な数値計算を通して、本発明の有効性が確認された。

　なお上記解析手法は、図２９に示すような解析装置１００によって実施されてもよい。
　解析装置１００は、パーソナルコンピュータ、サーバー、又は産業用コンピュータ等の装置によって実現される。解析装置１００は、処理部１１０と、記憶部１２０と、入力部１３０と、出力部１４０と、を含む。処理部１１０は、線形状態量に対する方程式に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式から線形状態量の高次モードの影響を補正して消去する。処理部１１０は、実モード方程式から、元の大規模系の解に対する影響の大きなモードである支配的モードを抽出し、影響の小さなモードである副次的モードについては、近似解を利用してその影響を非線形状態量に対する方程式に補正項として取り込んで消去し、低次元モデルを導出する。処理部１１０は、低次元モデルを利用して周波数応答を計算する。処理部１１０は、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）などのハードウェアプロセッサが記憶部１２０に格納されたプログラム（ソフトウェア）を実行することにより実現される。また、これらの機能部のうち一部または全部は、ＬＳＩ（Ｌａｒｇｅ　Ｓｃａｌｅ　Ｉｎｔｅｇｒａｔｉｏｎ）やＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、ＧＰＵ（Ｇｒａｐｈｉｃｓ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）などのハードウェア（回路部；ｃｉｒｃｕｉｔｒｙを含む）によって実現されてもよいし、ソフトウェアとハードウェアの協働によって実現されてもよい。プログラムは、予めＨＤＤ（Ｈａｒｄ　Ｄｉｓｋ　Ｄｒｉｖｅ）やフラッシュメモリなどの記憶装置（非一過性の記憶媒体を備える記憶装置）に格納されていてもよいし、ＤＶＤやＣＤ－ＲＯＭなどの着脱可能な記憶媒体（非一過性の記憶媒体）に格納されており、記憶媒体がドライブ装置に装着されることでインストールされてもよい。入力部１３０は、例えば、キーボードやマウス、タッチパネル等によって実現される。出力部１４０は、例えば、ディスプレイやプリンタ、タッチパネル等によって実現される。解析手法の実現に際し、閾値など、設定必要な情報に関しては、例えば、予め記憶部１２０に記憶されていてもよく、入力部１３０から研究者が入力してもよい。解析結果は、出力部１４０に出力してもよい。

Claims

　局所的に非線形性を有する大規模系の振動を解析する方法であって、
　線形状態量に対する方程式に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式から線形状態量の高次モードの影響を補正して消去する工程（１）と、
　前記実モード方程式から、元の大規模系の解に対する影響の大きなモードである支配的モードを抽出し、影響の小さなモードである副次的モードについては、近似解を利用してその影響を前記非線形状態量に対する方程式に補正項として取り込んで消去し、低次元モデルを導出する工程（２）と、
　前記低次元モデルを利用して周波数応答を計算する工程（３）と、を備える振動解析方法。
　角振動数ωについて、角振動数ωが解析範囲の上限または下限に到達するまで、角振動数ω±Δωをωに置き換え、前記工程（１）から前記工程（３）を繰り返す、請求項１に記載の振動解析方法。
　角振動数ω±Δωをωに置き換え前記工程（１）から前記工程（３）を繰り返すときに、角振動数ω±Δωにおける前記支配的モードと前記副次的モードとを、角振動数ωについての前記低次元モデルに基づいて分離する、請求項２に記載の振動解析方法。
　角振動数ωにおける前記低次元モデルから求められるモード振幅と、所定の閾値の関係に基づいて、角振動数ω±Δωにおける前記支配的モードと前記副次的モードとを分離する、請求項３に記載の振動解析方法。
　前記新型複素モード解析では、線形状態量に対する低次モードの複素固有値と複素固有ベクトルの計算および複素固有ベクトルの実数化をし、低次モードの物理座標に対する実モード座標の導入および実モード方程式の導出をする、請求項１から４のいずれか１項に記載の振動解析方法。
　局所的に非線形性を有する大規模系の振動を解析する装置としてのコンピュータを、
　線形状態量に対する方程式に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式から線形状態量の高次モードの影響を補正して消去する第１処理部と、
　前記実モード方程式から、元の大規模系の解に対する影響の大きなモードである支配的モードを抽出し、影響の小さなモードである副次的モードについては、近似解を利用してその影響を前記非線形状態量に対する方程式に補正項として取り込んで消去し、低次元モデルを導出する第２処理部と、
　前記低次元モデルを利用して周波数応答を計算する第３処理部と、して機能させるプログラム。
　請求項６に記載のプログラムが記憶された記憶媒体。