ホーム > 特許ランキング > 国立大学法人九州大学
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2024-06-28
(45)【発行日】2024-07-08
(54)【発明の名称】振動解析方法、プログラム、記憶媒体
(51)【国際特許分類】
G01M 7/02 20060101AFI20240701BHJP
G01H 17/00 20060101ALI20240701BHJP
【ＦＩ】
G01M7/02 H
G01H17/00 Z
【請求項の数】
7
(21)【出願番号】P 2020161727
(22)【出願日】2020-09-28
(65)【公開番号】P2022054596
(43)【公開日】2022-04-07
【審査請求日】2023-08-07
(73)【特許権者】
【識別番号】504145342
【氏名又は名称】国立大学法人九州大学
(74)【代理人】
【識別番号】100106909
【弁理士】
【氏名又は名称】棚井 澄雄
(74)【代理人】
【識別番号】100149548
【弁理士】
【氏名又は名称】松沼 泰史
(74)【代理人】
【識別番号】100141139
【弁理士】
【氏名又は名称】及川 周
(74)【代理人】
【識別番号】100188558
【弁理士】
【氏名又は名称】飯田 雅人
(72)【発明者】
【氏名】近藤 孝広
(72)【発明者】
【氏名】住川 大樹
【審査官】森口 正治
(56)【参考文献】
【文献】特開昭６３－２３１２３３（ＪＰ，Ａ）
【文献】特開２００５－２４９６８７（ＪＰ，Ａ）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０１Ｍ ７／０２－７／０６
Ｇ０１Ｈ １／００－１７／００
(57)【特許請求の範囲】
【請求項1】
局所的に非線形性を有する大規模系の振動を解析する方法であって、線形状態量に対する方程式に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式から線形状態量の高次モードの影響を補正して消去する工程（１）と、
前記実モード方程式から、元の大規模系の解に対する影響の大きなモードである支配的モードを抽出し、影響の小さなモードである副次的モードについては、近似解を利用してその影響を前記非線形状態量に対する方程式に補正項として取り込んで消去し、低次元モデルを導出する工程（２）と、
前記低次元モデルを利用して周波数応答を計算する工程（３）と、を備える振動解析方法。
【請求項2】
角振動数ωについて、角振動数ωが解析範囲の上限または下限に到達するまで、角振動数ω±Δωをωに置き換え、前記工程（１）から前記工程（３）を繰り返す、請求項１に記載の振動解析方法。
【請求項3】
角振動数ω±Δωをωに置き換え前記工程（１）から前記工程（３）を繰り返すときに、角振動数ω±Δωにおける前記支配的モードと前記副次的モードとを、角振動数ωについての前記低次元モデルに基づいて分離する、請求項２に記載の振動解析方法。
【請求項4】
角振動数ωにおける前記低次元モデルから求められるモード振幅と、所定の閾値の関係に基づいて、角振動数ω±Δωにおける前記支配的モードと前記副次的モードとを分離する、請求項３に記載の振動解析方法。
【請求項5】
前記新型複素モード解析では、線形状態量に対する低次モードの複素固有値と複素固有ベクトルの計算および複素固有ベクトルの実数化をし、低次モードの物理座標に対する実モード座標の導入および実モード方程式の導出をする、請求項１から４のいずれか１項に記載の振動解析方法。
【請求項6】
局所的に非線形性を有する大規模系の振動を解析する装置としてのコンピュータを、線形状態量に対する方程式に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式から線形状態量の高次モードの影響を補正して消去する第１処理部と、
前記実モード方程式から、元の大規模系の解に対する影響の大きなモードである支配的モードを抽出し、影響の小さなモードである副次的モードについては、近似解を利用してその影響を前記非線形状態量に対する方程式に補正項として取り込んで消去し、低次元モデルを導出する第２処理部と、
前記低次元モデルを利用して周波数応答を計算する第３処理部と、して機能させるプログラム。
【請求項7】
請求項６に記載のプログラムが記憶された記憶媒体。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、主に局所的強非線形性を有する大規模系の高性能振動解析手法に関する。より具体的には、新型複素モード解析を利用した低次元化法に関する。
【背景技術】
【0002】
機械は、互いに相対運動する多数の要素から構成されており、基礎支持部や要素間の結合部に、軸受や歯車等のような強い非線形性を有する要素を持っている。しかも、設計時の応力解析や振動解析においては、有限要素法などを利用することによって、実機に近い精密な大規模自由度系としてモデル化されることが多い。したがって、ほぼ全ての機械は、局所的に強い非線形性を有する大規模自由度系であるといえる。
【0003】
一方、近年、機械に対する高性能化の要求が益々強まっている。その要求を実現するには、設計段階において系内に不可避的に含まれる非線形性の影響を十分に考慮した振動解析を実施する必要がある。しかも、複雑多岐にわたる非線形振動の全貌を明らかにするには、単に特定の系パラメータに対して数値シミュレーションを行うだけでは不十分で、安定判別を含めた周波数特性の全容を解明しなければならない。
【0004】
これまで、多自由度非線形系の振動解析には、アルゴリズムに規則性があり幅広い問題に適用可能であることから、線形固有モード展開に基づくガラーキン法（例えば、下記特許文献１参照）がよく用いられて来た。
【先行技術文献】
【特許文献】
【0005】
【文献】特開２０１７－２１１８７７号公報
【発明の概要】
【発明が解決しようとする課題】
【0006】
しかしながら、ガラーキン法には、採用モード数を増やして高精度の解析を行おうとすると計算時間やメモリ量が爆発的に増大することや、非線形性が局所的にしか存在しないモデルであってもモード座標に変換後は非線形性が系全体に分散しモード分離の効果がなくなることなどの重大な欠点がある。
【0007】
このような状況を打開するために、少数のモードによる低次元モデルを用いて高精度の解析を行うための研究が精力的になされてきたが、未だ根本的な解決には至ってない。その主な問題点としては、採用するモードを適切に選択しないと、非線形系に特有のモード間の連成作用により、解析の精度（とりわけ安定判別の精度）が著しく悪化することが挙げられる。このように、大規模非線形系に対する実用的な振動解析手法はまだ存在しないといっても過言ではなく、この問題に対処し得る高性能な振動解析手法を開発することが、機械設計に関わる数値計算上の喫緊の課題となっている。
【0008】
本発明は、背景で述べた問題の抜本的な解決を目指したものであり、局所的に強い非線形性を有する大規模系を対象に、非線形性の影響を適切に反映した高精度の低次元モデルを合理的に構成する手法を提案している。
【課題を解決するための手段】
【0009】
提案手法は、主に次の二つの工程から構成される。
（１）線形状態量に対する方程式［後述の式（１）および式（４）］に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式［後述の式（２）］から線形状態量の高次モードの影響を補正して消去する工程。
（２）低次モードに対する実モード方程式から元の大規模系の解に対する影響の大きなモード（支配的モード）を抽出し、影響の小さなモード（副次的モード）については近似解を利用してその影響を非線形状態量に対する方程式［後述の式（２）］に補正項として取り込んで消去する低次元モデル導出工程、および周波数応答を計算する過程で低次モードの中から支配的モードを（なるべく合理的に）抽出する工程。
（１）線形状態量に対する方程式［後述の式（１）および式（４）］に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式［後述の式（２）］から線形状態量の高次モードの影響を補正して消去する工程。
（２）低次モードに対する実モード方程式から元の大規模系の解に対する影響の大きなモード（支配的モード）を抽出し、影響の小さなモード（副次的モード）については近似解を利用してその影響を非線形状態量に対する方程式［後述の式（２）］に補正項として取り込んで消去する低次元モデル導出工程、および周波数応答を計算する過程で低次モードの中から支配的モードを（なるべく合理的に）抽出する工程。
【発明の効果】
【0010】
この低次元モデルを利用すれば、これまでは不可能であった大規模非線形系の安定判別を含む振動解析を高速かつ高精度で行うことが可能になる。しかも、提案手法は非常に高い汎用性を有しており、局所的に強い非線形性を有する大規模振動系全般に適用可能である。
【図面の簡単な説明】
【0011】
【図1】本解析手法の概要を示す図である。
【図2】本解析手法のフローチャートである。
【図3】低次元化法適用例の解析モデル概略図である。
【図4】フルモデルにおける１次から５次の主共振領域の周波数応答図である。
【図5】フルモデルにおける１次から２次の主共振領域の周波数応答図である。
【図6】低次元モデル（ｎｌ＝５９、ｎａ＝５）における１次から５次の主共振領域の周波数応答図である。
【図7】低次元モデル（ｎｌ＝５９、ｎａ＝５）における１次から２次の主共振領域の周波数応答図である。
【図8】低次元モデル（ｎｌ＝５９、ｎａ＝１５）における１次から５次の主共振領域の周波数応答図である。
【図9】低次元モデル（ｎｌ＝５９、ｎａ＝１５）における１次から２次の主共振領域の周波数応答図である。
【図10】低次元モデル（ｎｌ＝５、ｎａ＝５）における１次から５次の主共振領域の周波数応答図である。
【図11】低次元モデル（ｎｌ＝１０、ｎａ＝５）における１次から５次の主共振領域の周波数応答図である。
【図12】低次モード数ｎｌの１次の固有振動数への影響を示す図である。
【図13】低次モード数ｎｌの２次の固有振動数への影響を示す図である。
【図14】低次モード数ｎｌの３次の固有振動数への影響を示す図である。
【図15】低次モード数ｎｌの４次の固有振動数への影響を示す図である。
【図16】低次モード数ｎｌの５次の固有振動数への影響を示す図である。
【図17】フルモデルにおける２次の主共振領域の周波数応答図である。
【図18】低次元モデル（ｎｌ＝２０，δ＝５×１０－４）における２次の主共振領域の周波数応答図である。
【図19】低次元モデル（ｎｌ＝２０，δ＝１×１０－４）における２次の主共振領域の周波数応答図である。
【図20】低次元モデル（ｎｌ＝２０，δ＝１×１０－５）における２次の主共振領域の周波数応答図である。
【図21】低次元モデル（ｎｌ＝１５，δ＝７×１０－４）における１次から２次の主共振領域の周波数応答図である。
【図22】低次元モデル（ｎｌ＝１５，δ＝１×１０－４）における１次から２次の主共振領域の周波数応答図である。
【図23】低次元モデル（ｎｌ＝１５，δ＝２×１０－５）における１次から２次の主共振領域の周波数応答図である。
【図24】低次元モデル（ｎｌ＝５９９，ｎａ＝２０）における１次から２次の主共振領域の周波数応答図である。
【図25】低次元モデル（ｎｌ＝３０，δ＝１×１０－３）における４次から５次の主共振領域の周波数応答図である。
【図26】低次元モデル（ｎｌ＝３０，δ＝５×１０－４）における４次から５次の主共振領域の周波数応答図である。
【図27】低次元モデル（ｎｌ＝３０，δ＝２×１０－４）における４次から５次の主共振領域の周波数応答図である。
【図28】低次元モデル（ｎｌ＝５９９，ｎａ＝２０）における４次から５次の主共振領域の周波数応答図である。
【図29】本解析手法を実施する装置の一例を示すブロック図である。
【発明を実施するための形態】
【0012】
以下、上述した二つの工程における具体的な実施事項について説明する。
なお、本明細書において、数式以外の本文中では、項や数式の表現が限られる。このため、例えば、本明細書に画像として取り込まれて項（０－１）のように表される項を、本明細書の本文中では、「ＣΨ～ １ Ｔ」と表記する。本明細書に画像として取り込まれて項（０－２）のように表される項を、本明細書の本文中では、「ｙ・」と表記する。なお、「ｙ・」は、ｙを時間により１階微分した項を意味する。
なお、本明細書において、数式以外の本文中では、項や数式の表現が限られる。このため、例えば、本明細書に画像として取り込まれて項（０－１）のように表される項を、本明細書の本文中では、「ＣΨ～ １ Ｔ」と表記する。本明細書に画像として取り込まれて項（０－２）のように表される項を、本明細書の本文中では、「ｙ・」と表記する。なお、「ｙ・」は、ｙを時間により１階微分した項を意味する。
【0013】
【数1】
【0014】
〔解析対象および基礎式〕
本発明で提案する手法は、全自由度がＮで局所的に強い非線形性を有する大規模自由度の振動系を解析対象とする。その基礎式は、一般的に次のように表すことができる。
本発明で提案する手法は、全自由度がＮで局所的に強い非線形性を有する大規模自由度の振動系を解析対象とする。その基礎式は、一般的に次のように表すことができる。
【0015】
【数2】
【0016】
ｙ１は、非線形力が作用しない状態量の物理座標をまとめたｎ次元ベクトル、ｙ２は、非線形力が作用する状態量の物理座標をまとめたｍ次元ベクトルである。以下では、ｙ１を線形状態量、ｙ２を非線形状態量と呼び、便宜的に式（１）を線形状態量の方程式、式（２）を非線形状態量の方程式と呼ぶ。ｎ２は、ｙ２およびｙ・
２に関する強非線形関数であり、ｆ１，ｆ２は、周期２π／ωの調和強制外力とする。Ｍは、質量要素に関する係数行列（質量行列）である。Ｃは、減衰要素に関する係数行列（減衰行列）である。Ｋは、ばね要素に関する係数行列（剛性行列）である。提案手法の汎用性を確保するため、Ｍ１１，Ｃ１１，Ｋ１１に対して正定性を仮定しない。ｍ＜＜ｎであるから、非線形性は局所的に存在することが分かる。このように、提案手法は、局所的に強い非線形性を有する大規模振動系全般に適用可能である。
【0017】
低次元化法の定式化にあたり、式（１）の係数行列を次のような転置行列に置き換え、外力を零とした系を同時に取り扱う。
【0018】
【数3】
【0019】
ここに、右上添え字「Ｔ」は転置記号である。式（４）のｙ２は、式（１）および式（２）を満足する解であるが、ｙ１は、式（１）および式（２）を満足する解とは一般に異なるので、式（４）ではｙ＊
１と表記している。以下では、式（４）を式（１）の双対系と呼ぶ。
【0020】
＜＜線形状態量の方程式に対する新型複素モード解析の適用＞＞
〔１階の微分方程式への変換〕
新型複素モード解析に基づく低次元化法を定式化するため、式（１）および式（４）を、次のような１階の微分方程式に変換する。なおＩｎは、ｎ次の単位行列である。
〔１階の微分方程式への変換〕
新型複素モード解析に基づく低次元化法を定式化するため、式（１）および式（４）を、次のような１階の微分方程式に変換する。なおＩｎは、ｎ次の単位行列である。
【0021】
【数4】
【0022】
〔一般固有値問題〕
本低次元化法では、まず、式（５）および式（６）でｘ・ ２＝０、ｘ２＝０と置くことによって非線形状態量を固定し、ｇ１＝０とした自由振動系から、次のような一般固有値問題を導出する。
本低次元化法では、まず、式（５）および式（６）でｘ・ ２＝０、ｘ２＝０と置くことによって非線形状態量を固定し、ｇ１＝０とした自由振動系から、次のような一般固有値問題を導出する。
【0023】
【数5】
【0024】
式（８）および式（９）から求められる固有値λは一致し、一般に複素数になる。また、対応する固有ベクトルＸ１およびＸ＊
１も複素ベクトルになる。式（８）を基準にすれば、Ｘ１が右複素固有ベクトル、Ｘ＊
１が左複素固有ベクトルとなり、式（９）を基準にすれば、Ｘ＊
１が右複素固有ベクトル、Ｘ１が左複素固有ベクトルとなる。振動学的には、固有値λの虚部が固有角振動数に対応し、Ｘ１およびＸ＊
１は左右の複素拘束モードに対応する。
【0025】
〔低次モードの複素拘束モード〕
第１段階として、式（８）および式（９）から複素固有値と左右の複素固有ベクトル（複素拘束モード）を求める。ところが、一般固有値問題の次元２ｎが大きいとき、すべての固有値と固有ベクトルを求めるのは困難である（大規模一般固有値問題の解法は、数値計算分野に残る未解決の難問であり、事実上不可能に近い）。ただし、低次から順に（絶対値の小さいものから順に）比較的少数の固有値と対応する固有ベクトルを求めることは、現時点でも可能である。
一方、線形系・非線形系ともに、モード次数が高くなるにつれて振動特性に及ぼす影響は一般に小さくなることが知られている。
第１段階として、式（８）および式（９）から複素固有値と左右の複素固有ベクトル（複素拘束モード）を求める。ところが、一般固有値問題の次元２ｎが大きいとき、すべての固有値と固有ベクトルを求めるのは困難である（大規模一般固有値問題の解法は、数値計算分野に残る未解決の難問であり、事実上不可能に近い）。ただし、低次から順に（絶対値の小さいものから順に）比較的少数の固有値と対応する固有ベクトルを求めることは、現時点でも可能である。
一方、線形系・非線形系ともに、モード次数が高くなるにつれて振動特性に及ぼす影響は一般に小さくなることが知られている。
【0026】
そこで、本低次元化法では、解の精度に大きな影響を及ぼす可能性のある１次からｎｌ（＜＜ｎ）次までの比較的低次の固有値と固有ベクトルのみを、一般固有値問題から求め、ｎｌ＋１次以上の高次モードの影響に関しては、ｎｌ次までの低次モードを用いて近似する。ｎｌの値については、解析すべき最大振動数の数倍程度の振動数範囲内に存在する固有振動数（固有値の虚部）の最大次数を目安に設定すれば良い。
【0027】
この手続きに必要な関係式は次の通りである。ただし、下添え字「ｌ」および「ｈ」はそれぞれ低次モードおよび高次モードに関連する変数または物理量であることを示す。また、左上添え字「Ｃ」が付された物理量は複素数であることを示し、同じく「Ｒ」および「Ｉ」は原則としてそれぞれその実部および虚部を表す。
まず、式（８）の１次からｎｌ（＜＜ｎ）次までの固有値と対応する左右の固有ベクトルを求める。これらの固有値は、全て実部が負の複素数であるものとして、ｐ（＝１，‥，ｎｌ）次の複素共役な固有値を次のように表す。
まず、式（８）の１次からｎｌ（＜＜ｎ）次までの固有値と対応する左右の固有ベクトルを求める。これらの固有値は、全て実部が負の複素数であるものとして、ｐ（＝１，‥，ｎｌ）次の複素共役な固有値を次のように表す。
【0028】
【数6】
【0029】
ここに、ωｐ，１およびζｐ，１はそれぞれｐ次モードの不減衰固有角振動数および減衰比を表す。さらに、これら２個ｎｌ組（計２ｎｌ個）の固有値をひとまとめにして、次のように行列表示する。
【0030】
【数7】
【0031】
ここに、diag[ ]は対角行列を表す。
この複素固有値に対応する式（８）のｎｌ組（２ｎｌ本）の右複素固有ベクトルをひとまとめにした右複素拘束モード行列を、ＣΦ＾ １ｌとする。同じく、ｎｌ組（２ｎｌ本）の左複素固有ベクトルをひとまとめにした左複素拘束モード行列を、ＣΨ＾ １ｌとする。ただし、ＣΦ＾ １ｌとＣΨ＾ １ｌは、次式を満足するように正規化されているものとする。
この複素固有値に対応する式（８）のｎｌ組（２ｎｌ本）の右複素固有ベクトルをひとまとめにした右複素拘束モード行列を、ＣΦ＾ １ｌとする。同じく、ｎｌ組（２ｎｌ本）の左複素固有ベクトルをひとまとめにした左複素拘束モード行列を、ＣΨ＾ １ｌとする。ただし、ＣΦ＾ １ｌとＣΨ＾ １ｌは、次式を満足するように正規化されているものとする。
【0032】
【数8】
【0033】
このとき、ＣΦ＾
１ｌ，ＣΨ＾
１ｌはそれぞれ次のように表される。なおＲφｐ，１，Ｉφｐ，１は、それぞれｐ次モードの右複素固有ベクトルにおける実部、虚部である。また、Ｒψｐ，１，Ｉψｐ，１は、それぞれｐ次モードの左複素固有ベクトルにおける実部、虚部である。
【0034】
【数9】
【0035】
〔低次モードの複素拘束モードの実数化〕
計算の高速化のために、低次モードに対する左右の複素拘束モード行列ＣΦ＾ １ｌおよびＣΨ＾ １ｌから、次のようにして右実拘束モード行列ＲΦ＾ １ｌ、および左実拘束モード行列ＲΨ＾ １ｌを導出する。
計算の高速化のために、低次モードに対する左右の複素拘束モード行列ＣΦ＾ １ｌおよびＣΨ＾ １ｌから、次のようにして右実拘束モード行列ＲΦ＾ １ｌ、および左実拘束モード行列ＲΨ＾ １ｌを導出する。
【0036】
【数10】
【0037】
〔低次モードに対する実モード座標の導入〕
線形状態量の物理座標ｘ１およびｘ＊ １は、次のように低次モードに基づくｘ１ｌおよびｘ＊ １ｌと高次モードに基づくｘ１ｈおよびｘ＊ １ｈの和で表すことができる。なおｙ１ｌおよびｙ＊ １ｌは、それぞれ低次モードにおける線形状態量である。ｙ１ｈおよびｙ＊ １ｈは、それぞれ高次モードにおける線形状態量である。
線形状態量の物理座標ｘ１およびｘ＊ １は、次のように低次モードに基づくｘ１ｌおよびｘ＊ １ｌと高次モードに基づくｘ１ｈおよびｘ＊ １ｈの和で表すことができる。なおｙ１ｌおよびｙ＊ １ｌは、それぞれ低次モードにおける線形状態量である。ｙ１ｈおよびｙ＊ １ｈは、それぞれ高次モードにおける線形状態量である。
【0038】
【数11】
【0039】
このうち、ｘ１ｌとｘ＊
１ｌに対しては、それぞれＲΦ＾
１ｌおよびＲΨ＾
１ｌを介して、対応する実モード座標η１ｌおよびη＊
１ｌを次式により導入する。なおξ１ｌは、低次モードの物理座標に対する実モード座標である。
【0040】
【数12】
【0041】
〔低次モードに対する実モード方程式の導出〕
式（１５）および式（１７）を式（５）に代入した上で左からＲΨ＾ １ｌ Ｔを乗じて整理すると、次のような、低次モードに対する実モード方程式が導出される。
式（１５）および式（１７）を式（５）に代入した上で左からＲΨ＾ １ｌ Ｔを乗じて整理すると、次のような、低次モードに対する実モード方程式が導出される。
【0042】
【数13】
【0043】
双対系の式（６）に対しても、式（１６）および式（１８）を代入した上で左からＲΦ＾
１ｌ
Ｔを乗じて整理すると、次のような実モード方程式が導出される。
【0044】
【数14】
【0045】
このように、係数行列が正定性を満たさない場合でも、質量行列、減衰行列および剛性行列を同時に対角化した上で、実数型２階微分方程式の形でモード方程式を導出することができる。これが新型複素モード解析法の最大の特長である。
【0046】
〔低次モードを利用した高次モードの影響の補正〕
非線形状態量の方程式（２）中の線形状態量ｙ１を、式（１５）に示すように低次モードに基づくｙ１ｌと高次モードに基づくｙ１ｈの和で表し、ｙ１ｌに関しては、式（１７）で導入した低次モードに対応する実モード座標に変換する。さらに、高次モードに基づくｙ１ｈの影響を補正して省略すると、次式が求められる。
非線形状態量の方程式（２）中の線形状態量ｙ１を、式（１５）に示すように低次モードに基づくｙ１ｌと高次モードに基づくｙ１ｈの和で表し、ｙ１ｌに関しては、式（１７）で導入した低次モードに対応する実モード座標に変換する。さらに、高次モードに基づくｙ１ｈの影響を補正して省略すると、次式が求められる。
【0047】
【数15】
【0048】
ここに、Ｍ～
２２ｈ，Ｃ～
２２ｈ，Ｋ～
２２ｈ，ｐ２１ｈ，ｑ２１ｈが高次モードの補正項を表す。これらは、いずれも次のように低次モードから求められる。
【0049】
【数16】
【0050】
＜＜低次元モデルの導出＞＞
〔支配的モードと副次的モードの分離〕
式（１７）および式（１８）で導入した実モード座標η１ｌおよびη＊ １ｌの中には、元の大規模系の解の精度に及ぼす影響が大きなモード座標（支配的モードと呼ぶ）と影響の小さなモード座標（副次的モードと呼ぶ）が存在する。下添え字「ａ」および「ｂ」により、それぞれ支配的モードおよび副次的モードに関連する変数または物理量であることを示す。例えば、η１ａは支配的モード座標を、η１ｂは副次的モード座標を表す。
〔支配的モードと副次的モードの分離〕
式（１７）および式（１８）で導入した実モード座標η１ｌおよびη＊ １ｌの中には、元の大規模系の解の精度に及ぼす影響が大きなモード座標（支配的モードと呼ぶ）と影響の小さなモード座標（副次的モードと呼ぶ）が存在する。下添え字「ａ」および「ｂ」により、それぞれ支配的モードおよび副次的モードに関連する変数または物理量であることを示す。例えば、η１ａは支配的モード座標を、η１ｂは副次的モード座標を表す。
【0051】
式（１４）で求めた低次モードに対する左右の実モード行列ＲΦ＾
１ｌ，ＲΨ＾
１ｌを、ｎａ組（２ｎａ本）の支配的モードからなるモード行列ＲΦ＾
１ａ，ＲΨ＾
１ａと、ｎｂ組（２ｎｂ本）の副次的モードからなるモード行列ＲΦ＾
１ｂ，ＲΨ＾
１ｂとに分離し（ｎｌ＝ｎａ＋ｎｂ）、次のように再整理する。ただし、以下では、下添字「＃」は「ａ」または「ｂ」を表す。
【0052】
【数17】
【0053】
実モード行列の分離に対応して、式（１９）、式（２０）、式（２１）のモード座標も、次のように支配的モードと副次的モードに分離する。
【0054】
【数18】
【0055】
〔副次的モードの近似解〕
副次的モードに対しては非線形性の影響が小さいのでこの影響を無視し、ｙ２，ξ１ｂ，ξ＊ １ｂが基本周期２π／ωで調和振動しているものと仮定すると、ｙ・・ ２＝－ω２ｙ２，ξ・・ １ｂ＝－ω２ξ１ｂ，ξ・・＊ １ｂ＝－ω２ξ＊ １ｂと表される。さらに、調和外力ではｆ・・ １＝－ω２ｆ１であることを考慮すると、式（２４）および式（２５）より、副次的モードの近似解が次のように求められる。
副次的モードに対しては非線形性の影響が小さいのでこの影響を無視し、ｙ２，ξ１ｂ，ξ＊ １ｂが基本周期２π／ωで調和振動しているものと仮定すると、ｙ・・ ２＝－ω２ｙ２，ξ・・ １ｂ＝－ω２ξ１ｂ，ξ・・＊ １ｂ＝－ω２ξ＊ １ｂと表される。さらに、調和外力ではｆ・・ １＝－ω２ｆ１であることを考慮すると、式（２４）および式（２５）より、副次的モードの近似解が次のように求められる。
【0056】
【数19】
【0057】
〔副次的モードの消去と低次元モデルの導出〕
式（２６）中の副次的モード座標ξ１ｂ，ξ・ １ｂ，ξ・・ １ｂに式（２７）の近似解を代入し、副次的モードの影響を補正して省略すると、次式を得る（この導出過程で式（２８）の結果を利用している）。
式（２６）中の副次的モード座標ξ１ｂ，ξ・ １ｂ，ξ・・ １ｂに式（２７）の近似解を代入し、副次的モードの影響を補正して省略すると、次式を得る（この導出過程で式（２８）の結果を利用している）。
【0058】
【数20】
【0059】
ここに、Ｍ～
２２ｂ，Ｃ～
２２ｂ，Ｋ～
２２ｂ，ｐ２１ｂ，ｑ２１ｂが副次的モードの補正項を表し、次式で与えられる。
【0060】
【数21】
【0061】
さらに、式（２４）で下添え字「＃」を「ａ」とした式と式（２９）から、次のように（ｎａ＋ｍ）次元に低次元化された方程式が求められる。
【0062】
【数22】
【0063】
式（３１）で表されるモデルを低次元モデルと呼ぶ。式（３１）を利用して近似解を計算する。
【0064】
〔線形状態量の物理座標の近似解〕
低次元モデルの式（３１）から近似解ξ１ａ，ξ・ １ａ，ｙ２，ｙ・ ２が求められたとき、必要に応じて線形自由度の近似解を物理座標に変換する。
この手続きに必要な関係式は次の通りである。
低次元モデルの式（３１）から近似解ξ１ａ，ξ・ １ａ，ｙ２，ｙ・ ２が求められたとき、必要に応じて線形自由度の近似解を物理座標に変換する。
この手続きに必要な関係式は次の通りである。
【0065】
【数23】
【0066】
ｙ１ａ，ｙ・
１ａは支配的モードに基づく物理座標を、ｙ１ｂ，ｙ・
１ｂは副次的モードに基づく物理座標を、ｙ１ｈ，ｙ・
１ｈは高次モードに基づく物理座標を表す。
【0067】
〔支配的モード座標の抽出法〕
非線形系において外力の角振動数ωを徐々に変化させながら複雑多岐にわたる周波数応答の全容を解明しようとする場合には、あるωに対して求められた近似解をω±Δωに対する初期値として逐次近似計算することが多い。そこで、あるωにおいて求められた近似解に含まれる情報を利用して、ω±Δωで抽出すべきξ１ａを選定する。
この手続きに必要な関係式は次の通りである。
いま、あるωに対して、ξ１ａ，ξ・ １ａ，ｙ２，ｙ・ ２からなる低次元モデルの近似解が求められたとする。この近似解ｙ２，ｙ・ ２から、次式により低次モード座標ξ１ｌの近似解が計算できる。
非線形系において外力の角振動数ωを徐々に変化させながら複雑多岐にわたる周波数応答の全容を解明しようとする場合には、あるωに対して求められた近似解をω±Δωに対する初期値として逐次近似計算することが多い。そこで、あるωにおいて求められた近似解に含まれる情報を利用して、ω±Δωで抽出すべきξ１ａを選定する。
この手続きに必要な関係式は次の通りである。
いま、あるωに対して、ξ１ａ，ξ・ １ａ，ｙ２，ｙ・ ２からなる低次元モデルの近似解が求められたとする。この近似解ｙ２，ｙ・ ２から、次式により低次モード座標ξ１ｌの近似解が計算できる。
【0068】
【数24】
【0069】
この近似解のモード別の振幅を次式で定義できる。
【0070】
【数25】
【0071】
ω±Δωの計算においては、このモード振幅｜｜ξ１，ｐ｜｜が、設定した閾値δよりも大きいモードのみをξ１ａとして抽出する。抽出された支配的モードの本数をｎａ、残りの副次的モードの本数をｎｂとすると、ｎｌ＝ｎａ＋ｎｂである。場合によっては、ｎａの上限値ｎａ，ｍａｘを定めることもある。
この手法により、解の特性に対応してなるべく少ない本数の支配的モードを適切に選択しながら周波数応答を効率的に計算することが可能となる。この手法ではモード振幅の閾値δおよび支配的モード本数の上限値ｎａ，ｍａｘの設定法が問題となるが、実用的にはδについては徐々に小さくしながら、ｎａ，ｍａｘについては徐々に大きくしながら複数の周波数応答を求め、解析結果が変化しなくなった時点で高精度の解が求められたものとみなせばよい。低次元モデルの計算はかなり高速で行えるので、実用上は複数の計算にかかるコストは大きな問題にはならないと考えられる。なお、解析結果が変化しなくなった時点とは、解析結果が全く変化しない場合だけではなく、解析結果を比較した結果、その差が、所定の基準を下回った場合を含む。
なお、計算を開始する最初のωにおける支配的モード本数ｎａの設定方法については、例えば、ｎａ＝ｎｌ（低次モードの本数）とすればよい。
この手法により、解の特性に対応してなるべく少ない本数の支配的モードを適切に選択しながら周波数応答を効率的に計算することが可能となる。この手法ではモード振幅の閾値δおよび支配的モード本数の上限値ｎａ，ｍａｘの設定法が問題となるが、実用的にはδについては徐々に小さくしながら、ｎａ，ｍａｘについては徐々に大きくしながら複数の周波数応答を求め、解析結果が変化しなくなった時点で高精度の解が求められたものとみなせばよい。低次元モデルの計算はかなり高速で行えるので、実用上は複数の計算にかかるコストは大きな問題にはならないと考えられる。なお、解析結果が変化しなくなった時点とは、解析結果が全く変化しない場合だけではなく、解析結果を比較した結果、その差が、所定の基準を下回った場合を含む。
なお、計算を開始する最初のωにおける支配的モード本数ｎａの設定方法については、例えば、ｎａ＝ｎｌ（低次モードの本数）とすればよい。
【0072】
＜＜周波数応答の計算手順＞＞
以下に、周波数応答の計算手順を簡潔にまとめる（図１、図２参照）。
（Ｓｔｅｐ１）前処理
（Ｓｔｅｐ１－１）
線形状態量に対する低次モードの複素固有値と複素固有ベクトルの計算［式（８），式（１０），式（１２），式（１３）］および複素固有ベクトルの実数化［式（１４）］
（Ｓｔｅｐ１－２）
低次モードの物理座標に対する実モード座標の導入［式（１７），式（１８）］および実モード方程式の導出［式（１９），式（２０）］
（Ｓｔｅｐ１－３）
非線形状態量の方程式に対する低次モードを利用した高次モードの影響の補正［式（２１）］
（Ｓｔｅｐ２）支配的モードの抽出［式（３６），式（３７）］
（Ｓｔｅｐ３）低次元モデルの導出
（Ｓｔｅｐ３－１）
支配的モードと副次的モードの分離［式（２３），式（２４），式（２５），式（２６）］
（Ｓｔｅｐ３－２）
副次的モードの近似解の計算［式（２７），式（２８）］
（Ｓｔｅｐ３－３)
低次元モデルの導出［式（２９），式（３０），式（３１）］
（Ｓｔｅｐ４）低次元モデル［式（３１）］の周期解の計算
（Ｓｔｅｐ５）計算終了・継続の判定
角振動数ωが解析範囲の上限または下限に到達すれば計算を終了する。到達していなければ、ω±Δωをωに置き換え、Ｍ１１，Ｃ１１，Ｋ１１がωの関数である場合にはＳｔｅｐ１に戻り、定数の場合にはＳｔｅｐ２に戻り計算を継続する。
以下に、周波数応答の計算手順を簡潔にまとめる（図１、図２参照）。
（Ｓｔｅｐ１）前処理
（Ｓｔｅｐ１－１）
線形状態量に対する低次モードの複素固有値と複素固有ベクトルの計算［式（８），式（１０），式（１２），式（１３）］および複素固有ベクトルの実数化［式（１４）］
（Ｓｔｅｐ１－２）
低次モードの物理座標に対する実モード座標の導入［式（１７），式（１８）］および実モード方程式の導出［式（１９），式（２０）］
（Ｓｔｅｐ１－３）
非線形状態量の方程式に対する低次モードを利用した高次モードの影響の補正［式（２１）］
（Ｓｔｅｐ２）支配的モードの抽出［式（３６），式（３７）］
（Ｓｔｅｐ３）低次元モデルの導出
（Ｓｔｅｐ３－１）
支配的モードと副次的モードの分離［式（２３），式（２４），式（２５），式（２６）］
（Ｓｔｅｐ３－２）
副次的モードの近似解の計算［式（２７），式（２８）］
（Ｓｔｅｐ３－３)
低次元モデルの導出［式（２９），式（３０），式（３１）］
（Ｓｔｅｐ４）低次元モデル［式（３１）］の周期解の計算
（Ｓｔｅｐ５）計算終了・継続の判定
角振動数ωが解析範囲の上限または下限に到達すれば計算を終了する。到達していなければ、ω±Δωをωに置き換え、Ｍ１１，Ｃ１１，Ｋ１１がωの関数である場合にはＳｔｅｐ１に戻り、定数の場合にはＳｔｅｐ２に戻り計算を継続する。
【0073】
＜＜低次元化法適用例＞＞
〔１．解析モデル〕
新型複素モード解析を利用した低次元化法の有効性、特に、高次モード消去法および支配的モードの抽出法の有効性を確認するため、直線状梁構造物の面内曲げ振動を対象に本低次元化法を適用する。解析モデルの模式図を図３に、計算に用いたパラメータを表１に示す。
〔１．解析モデル〕
新型複素モード解析を利用した低次元化法の有効性、特に、高次モード消去法および支配的モードの抽出法の有効性を確認するため、直線状梁構造物の面内曲げ振動を対象に本低次元化法を適用する。解析モデルの模式図を図３に、計算に用いたパラメータを表１に示す。
【0074】
【表1】
【0075】
このモデル１０は、長さ１．５ｍ、外径３０ｍｍ、内径２５ｍｍの一様中空円形断面の梁１の両端を、非線形要素２（粘性減衰器と漸硬型３次ばね）によって支持したものであり、左端から０．５ｍの位置には、断片線形ばね（ガタ）３を配置している。また、外力の振動数が比較的高い領域においても、振幅の大きな振動が発生して非線形性の影響が強く表れるように、遠心力タイプの調和強制外力Ｆを左端から１ｍの位置に作用させている。
梁１の物性値は、鉄鋼材料を想定している。梁１については３０個の微小要素に等分割し、有限要素法により質量行列、剛性行列を導出した。この場合の自由度は、線形状態量がｎ＝５９、非線形状態量がｍ＝３であり、総自由度がＮ＝６２となる。
梁１の物性値は、鉄鋼材料を想定している。梁１については３０個の微小要素に等分割し、有限要素法により質量行列、剛性行列を導出した。この場合の自由度は、線形状態量がｎ＝５９、非線形状態量がｍ＝３であり、総自由度がＮ＝６２となる。
【0076】
これは大規模系と呼ぶには相応しくない自由度であるが、低次元化法を適用しない場合の解を求めることのできる限界に近い自由度である。新型複素モード解析の長所が現れるように、減衰行列は分布減衰を与えて非比例減衰系としている。表２に本解析モデルにおいて漸硬型ばね定数を０として線形化した系の１次から１０次の不減衰固有振動数を示す。
【0077】
【表2】
【0078】
図４から図９に本解析モデルの周波数応答図を示す。図４、図６、図８は１次から５次までの主共振領域を、図５、図７、図９は１次および２次の主共振領域の拡大図を示している。図４および図５は低次元化法を適用しない場合（フルモデル）の結果、図６から図９は低次元化法を適用した場合（低次元モデル）の結果である。
【0079】
低次元モデルの構成に際しては高次モード消去法および支配的モードの抽出法は適用せず、外力の振動数が拘束モードの固有振動数に近い順に一定本数（ｎａ＝５、１５）を抽出した。低次元モデルの結果では、支配的モード本数（数値は右軸）を、塗り潰した▽印で示したが、図６から図９では一定であることを表している（結果として太い実線のように見えている）。
【0080】
周波数応答における応答曲線はすべてガタを配置した点の最大片振幅を描いており、振幅がガタの隙間（０．００２５ｍ）を超えると応答曲線が屈曲している。図中の実線は安定解、破線は不安定解、〇印はサドルノード分岐点、□印はピッチフォーク分岐点、△印はホップ分岐点を表す。
【0081】
ｎａ＝５の場合の図６および図７とフルモデルの図４および図５を比較すると概ね良い一致を示しているが、強い非線形性であるガタの影響が現れる領域で応答曲線の形状や振幅のピーク値、安定判別結果が異なっている。図８および図９に示すようにｎａ＝１５まで支配的モード数を増加させるとフルモデルと一致した結果が得られる。
【0082】
〔２．高次モード消去法の適用〕
本解析モデルに高次モード消去法を適用して、解析に用いる低次モード数ｎｌの解析結果への影響を調査した。図１０および図１１は支配的モード数をｎａ＝５と固定して低次モード数をｎｌ＝５，１０とした場合の結果である。高次モード消去法を適用しない場合の図６（ｎｌ＝ｎ＝５９）と比較すると、図１０のｎｌ＝５では１次、２次の主共振領域ではほぼ変わらない応答が得られているが、振動数が高くなるにつれてずれが大きくなる。図１１のｎｌ＝１０とすればこの解析範囲においてはほぼ一致した結果が得られていることがわかる。
本解析モデルに高次モード消去法を適用して、解析に用いる低次モード数ｎｌの解析結果への影響を調査した。図１０および図１１は支配的モード数をｎａ＝５と固定して低次モード数をｎｌ＝５，１０とした場合の結果である。高次モード消去法を適用しない場合の図６（ｎｌ＝ｎ＝５９）と比較すると、図１０のｎｌ＝５では１次、２次の主共振領域ではほぼ変わらない応答が得られているが、振動数が高くなるにつれてずれが大きくなる。図１１のｎｌ＝１０とすればこの解析範囲においてはほぼ一致した結果が得られていることがわかる。
【0083】
このように解析を行う振動数の範囲によって必要な低次モード数は変化し、低次モード数を増加させるとより高い振動数範囲まで高精度に解が求まる。この解析モデルでは５９本中の１０本と少数の低次モードのみで高次モードの影響を高精度に推定できていることから、大規模系の解析に対してこの性質は極めて有効であることが期待される。
【0084】
図１２から図１６には低次モードの数の固有振動数への影響を示しており、図１２から図１６は解析振動数範囲の１次から５次の固有振動数である。これらの固有振動数は、高次モード消去法のみを適用して作成した低次元モデル（ｎｌ＝ｎａ）から求めた。点は低次元モデル、破線はフルモデルの固有振動数である。
【0085】
図１２から図１６の比較から、低次の固有振動数ほど精度がよく高次になるにつれて低精度である。また、低次モード数の増加に伴って精度が向上することがわかる。このように低次モード数に対して固有振動数の精度と周波数応答解析結果の精度が同様の傾向を示しており、低次モード数ｎｌの値の設定法としては、計算コストの観点から、固有振動数の精度を目安にすればよいと考えられる。
【0086】
〔３．支配的モードの抽出法の適用〕
本解析モデルに対して支配的モードの抽出法を適用した。解析範囲は２次の主共振領域に限定した。低次元化法には高次モード消去法も併せて適用しており、低次モード数はｎｌ＝２０とした。図１７はフルモデルの結果、図１８から図２０はそれぞれ閾値をδ＝５×１０－４，１×１０－４，１×１０－５として抽出法を適用した結果であり、支配的モード数が外力の振動数に対して変化していることがわかる。図１８から図２０にかけて閾値を徐々に小さくしていくと、モード数が増加し解析に必要なモードが抽出されることでフルモデルの解に近づいていくことがわかる。図２０のδ＝１×１０－５ではフルモデルと一致した結果が得られている。
本解析モデルに対して支配的モードの抽出法を適用した。解析範囲は２次の主共振領域に限定した。低次元化法には高次モード消去法も併せて適用しており、低次モード数はｎｌ＝２０とした。図１７はフルモデルの結果、図１８から図２０はそれぞれ閾値をδ＝５×１０－４，１×１０－４，１×１０－５として抽出法を適用した結果であり、支配的モード数が外力の振動数に対して変化していることがわかる。図１８から図２０にかけて閾値を徐々に小さくしていくと、モード数が増加し解析に必要なモードが抽出されることでフルモデルの解に近づいていくことがわかる。図２０のδ＝１×１０－５ではフルモデルと一致した結果が得られている。
【0087】
いずれの閾値においても、解の振幅が小さく非線形性の影響が表れていない領域では支配的モード数は比較的少なく、振幅が大きくなるにつれて支配的モード数が増加しており、解の特性に応じて支配的モードを変化させながら効率的に解析できていることが確認された。
【0088】
〔４．自由度の大きいモデル〕
本解析モデルにおいて分割数を３００としたモデルに対して、本手法を適用して有効性の検証を行った。この場合の自由度は、線形状態量がｎ＝５９９、非線形状態量がｍ＝３であり、総自由度がＮ＝６０２となる。この自由度ではフルモデルに対する解析は困難となる。解析は１次から２次の主共振領域と４次から５次の主共振領域にわけて行った。
本解析モデルにおいて分割数を３００としたモデルに対して、本手法を適用して有効性の検証を行った。この場合の自由度は、線形状態量がｎ＝５９９、非線形状態量がｍ＝３であり、総自由度がＮ＝６０２となる。この自由度ではフルモデルに対する解析は困難となる。解析は１次から２次の主共振領域と４次から５次の主共振領域にわけて行った。
【0089】
〔４．１．１次、２次の主共振領域〕
解析振動数範囲は、固有振動数および１章〔１．解析モデル〕の結果から、外力の振動数を２０Ｈｚから２００Ｈｚまでとした。高次モード消去法を適用する際の低次モード数は、最大振動数の２００Ｈｚの５倍程度の６次の固有振動数を目安にｎｌ＝１５とした。図２１から図２３はそれぞれ閾値をδ＝７×１０－４，１×１０－４，２×１０－５と徐々に小さくしながら支配的モードの抽出法を適用した結果である。図２１と図２２を比較すると、１次および２次共振の応答に差異がみられる。図２２と図２３を比較すると、応答の変化は小さくなったが、２次のピーク付近の解がわずかに異なっていることがわかる。
解析振動数範囲は、固有振動数および１章〔１．解析モデル〕の結果から、外力の振動数を２０Ｈｚから２００Ｈｚまでとした。高次モード消去法を適用する際の低次モード数は、最大振動数の２００Ｈｚの５倍程度の６次の固有振動数を目安にｎｌ＝１５とした。図２１から図２３はそれぞれ閾値をδ＝７×１０－４，１×１０－４，２×１０－５と徐々に小さくしながら支配的モードの抽出法を適用した結果である。図２１と図２２を比較すると、１次および２次共振の応答に差異がみられる。図２２と図２３を比較すると、応答の変化は小さくなったが、２次のピーク付近の解がわずかに異なっていることがわかる。
【0090】
図２３の結果が高精度な結果であることを確認するため、高次モード消去法および支配的モードの抽出法は適用せずに１章と同様にｎａを徐々に増加させながら低次元モデルの解析を行い、応答に変化がなくなった時点での結果を本解析モデルの正解とみなした。その結果を図２４に示しており、低次モード数がｎｌ＝５９９、支配的モード数がｎａ＝２０である。この結果は、図２３と非常に良い一致を示していることがわかる。
【0091】
〔４．２．４次、５次の主共振領域〕
解析振動数範囲は４５０Ｈｚから１０００Ｈｚまでとして、低次モード数は最大振動数の１０００Ｈｚの３倍程度の１０次の固有振動数を目安にｎｌ＝３０とした。前節と同様に閾値を徐々に小さくしながら支配的モードの抽出法を適用した。図２５から図２７はそれぞれ閾値をδ＝１×１０－３，５×１０－４，２×１０－４とした結果である。
前節と同様に図２５から図２７を比較すると低次元モデルの応答の変化が小さくなっていくことがわかる。また、図２７と、図２８のｎｌ＝５９９、ｎａ＝２０の結果との比較から、この振動数範囲においても本低次元化法によってフルモデルとほとんど変わらない結果を効率的に計算できている。ただし、解析を行う振動数範囲によって低次モード数や最適な閾値は異なる。
以上のように、本解析モデルに対する具体的な数値計算を通して、本発明の有効性が確認された。
解析振動数範囲は４５０Ｈｚから１０００Ｈｚまでとして、低次モード数は最大振動数の１０００Ｈｚの３倍程度の１０次の固有振動数を目安にｎｌ＝３０とした。前節と同様に閾値を徐々に小さくしながら支配的モードの抽出法を適用した。図２５から図２７はそれぞれ閾値をδ＝１×１０－３，５×１０－４，２×１０－４とした結果である。
前節と同様に図２５から図２７を比較すると低次元モデルの応答の変化が小さくなっていくことがわかる。また、図２７と、図２８のｎｌ＝５９９、ｎａ＝２０の結果との比較から、この振動数範囲においても本低次元化法によってフルモデルとほとんど変わらない結果を効率的に計算できている。ただし、解析を行う振動数範囲によって低次モード数や最適な閾値は異なる。
以上のように、本解析モデルに対する具体的な数値計算を通して、本発明の有効性が確認された。
【0092】
なお上記解析手法は、図２９に示すような解析装置１００によって実施されてもよい。解析装置１００は、パーソナルコンピュータ、サーバー、又は産業用コンピュータ等の装置によって実現される。解析装置１００は、処理部１１０と、記憶部１２０と、入力部１３０と、出力部１４０と、を含む。処理部１１０は、線形状態量に対する方程式に新型複素モード解析を適用して低次モードに対する実モード方程式に変換するとともに、非線形状態量に対する方程式から線形状態量の高次モードの影響を補正して消去する。処理部１１０は、実モード方程式から、元の大規模系の解に対する影響の大きなモードである支配的モードを抽出し、影響の小さなモードである副次的モードについては、近似解を利用してその影響を非線形状態量に対する方程式に補正項として取り込んで消去し、低次元モデルを導出する。処理部１１０は、低次元モデルを利用して周波数応答を計算する。処理部１１０は、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）などのハードウェアプロセッサが記憶部１２０に格納されたプログラム（ソフトウェア）を実行することにより実現される。また、これらの機能部のうち一部または全部は、ＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ）やＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、ＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）などのハードウェア（回路部；ｃｉｒｃｕｉｔｒｙを含む）によって実現されてもよいし、ソフトウェアとハードウェアの協働によって実現されてもよい。プログラムは、予めＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）やフラッシュメモリなどの記憶装置（非一過性の記憶媒体を備える記憶装置）に格納されていてもよいし、ＤＶＤやＣＤ－ＲＯＭなどの着脱可能な記憶媒体（非一過性の記憶媒体）に格納されており、記憶媒体がドライブ装置に装着されることでインストールされてもよい。入力部１３０は、例えば、キーボードやマウス、タッチパネル等によって実現される。出力部１４０は、例えば、ディスプレイやプリンタ、タッチパネル等によって実現される。解析手法の実現に際し、閾値など、設定必要な情報に関しては、例えば、予め記憶部１２０に記憶されていてもよく、入力部１３０から研究者が入力してもよい。解析結果は、出力部１４０に出力してもよい。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】
【図26】
【図27】
【図28】
【図29】