TECHNICAL PAPERS
Identification of flutter parameters for a wing model
Carlos De Marqui JuniorI; Daniela C. RebolhoII; Eduardo M. BeloIII; Flávio D. MarquesIV
Iengineering School of Sao Carlos; University of Sao Paulo; Laboratory of Aeroelasticity; Flight Dynamics and Control;Av. All Rights Reserved. Trabalhador Sancarlense 400; 13566- 590 Sao Carlos, SP. Brazil; [email protected]
[email protected]
[email protected]
[email protected]
ABSTRACT
風洞での剛体翼のフラッター試験用に柔軟なマウントシステムが開発された。 この実験システムで得られた2自由度のフラッターは、構造的な曲げ振動モードとねじり振動モードの組み合わせとして記述することができる。 この実験装置を用いて、後縁フラップをアクチュエータとしたフラッタ抑制のためのアクティブ制御方式を検証することができます。 制御方式を開発する前に、システムの動的および空力的な特性を調査する必要があります。 実験的にモード解析を行い、モード形状と周波数を決定する。 次に、風洞試験を行い、フラッター現象の特徴を把握し、臨界フラッター速度と臨界フラッター周波数を決定します。 また、臨界速度以下の速度範囲について周波数応答関数を求め、速度上昇に伴うピッチモード、プランジモードおよびカップリング傾向の変化を示しています。 これらの試験で得られた時間領域のピッチおよびプランジデータは、拡張固有値実現アルゴリズムが速度増加に伴うフラッターパラメータを特定する能力を評価するために使用されています。 同定処理の結果は、フラッターに関与するモードの周波数と減衰の発展という形で示される。
Keywords: 同定、フラッター、EERA、空力弾性
はじめに
空力弾性現象は航空構造物における弾性、慣性、空力負荷の相互作用から生じる。 弾性体が気流にさらされると、構造変形が追加の空力力を誘発し、これらの力がさらなる構造変形を生み、それがまたより大きな空力力を誘発する。 この相互作用はフラッターのような空力的な不安定性を引き起こす可能性がある（例えば Försching (1979) を参照）。 第二次世界大戦後、飛行速度の向上と構造的な改良により、空力弾性の問題はより重要なものとなった。 空力弾性の歴史的変遷については、Ashley（1970）、Collar（1959）、Garrick and Reed（1981）、Garrick（1976）などに記述されているが、
フラッターは空力弾性の最も代表的なトピックの一つである。 フラッターは空力的な負荷によって構造モードが同時に結合し、励起される複雑な現象である。 より正式には、航空機の部品が臨界速度より速い速度で自己保持された振動的挙動を示す状態をフラッターと呼ぶ（Wright、1991）。 一般に、フラッターは翼や尾翼のような大きな空力荷重を受ける昇降面で発生する。
飛行フラッター試験（Kehoe, 1995）は、航空機の認証において非常に重要な部分である。 このような時間とコストのかかる試験において、航空機が望ましい条件下でフラッターフリーであることを示すために、飛行範囲を安全に拡大する必要がある。
これらの飛行試験の主なタスクは、空力弾性パラメータ（第2段階）を推定することによって可能になる、次の試験速度での安定性を確信を持って予測することである。 測定された試験データからフラッター発生に関連する速度、あるいはその他の空力弾性不安定性を正確に予測する方法の開発は、これらの試験の安全性を高め、さらにはコストを削減するための重要な方法である（Lind, 2003）。 この目的を達成するために、いくつかの方法が開発されてきた。 一般に、これらの方法はシミュレーションからのデータを使用して開発およびテストされるが、あるアプローチを包絡線拡張に確実に使用できるようになる前に、飛行試験からのデータを含む評価を行う必要がある。 Cooper; Emmett and Wright (1993)が開発した包絡線関数もその一つである。 この関数は、インパルス応答関数がシステムの全体的な安定性に関する情報を含んでいるという仮定に基づくものである。 同様に、離散時間自己回帰移動平均モデル（ARMA）は、Jury Stability基準を用い、システムの全体的な安定性も考慮する（Torii and Matsuzaki, 2001）。 もう一つの方法は、Zimmerman-Weissenburgerフラッターマージン法で、減衰追跡の代わりにRouth安定性基準を用いるべきである（Zimmerman and Weissenburger, 1964）。 フラッターメータは，Lind and Brenner (2000)によって開発されたフラッターマージンの予測に用いられるオンラインモデルベースツールである． このツールは実験データと理論モデルを使用してフラッターの発生を予測する。
飛行試験からフラッターパラメータを予測する前述の方法の能力は、Lind（2003）によって評価されている。 この飛行試験は、F-15を空力構造試験翼（ATW）のホストキャリアとして使用して実施された。 この ATW は完全な航空機ではないが，現実的な翼であり，飛行試験中にフラッター速度が達成 されるところまでエンベロープを拡張することができる． 真のフラッター速度がわかっているため、予測されるフラッター速度の評価に使用することができます。 これらの評価から得られた結果は、異なる条件下での各手法の長所と短所を示しています。 例えば、データベースの手法は、低速試験のデータを用いてフラッター速度を正確に予測することはできませんが、対気速度が増加するにつれて良い解に収束していきます。 しかし、モデルベースのフラッターメーターは、低速試験のデータでは保守的ですが、高速試験のデータでは予測は保守的なまま、真のフラッター速度に収束しません。
本論文では、拡張固有値実現アルゴリズム（Extended Eigensystem Realization Algorithm：EERA）を用いて、フラッターパラメータ、すなわち周波数と減衰を同定するアプローチを紹介する。 これらのフラッターパラメータの同定は，風洞試験から得られたデータを解析することで行われる。 風洞試験は、剛体翼モデルと組み合わせたときに2自由度のフラッタを達成するように設計された柔軟な取り付けシステムで実施されています。 この実験システムの風離れ特性や空力弾性特性は、有限要素シミュレーション、実験モード解析、風洞試験により広範囲に決定されています（De Marqui Jr et al.、2004）。 EERA法はEigensystem Realization Algorithm (ERA)を改良したもので、モードを同時に特定できる時間領域のアルゴリズムである(Juang, 1994)。 EERAは入力と出力の時刻歴からブロックハンケル行列を操作してモードパラメータを計算する (Tasker; Bosse and Fisher, 1998)。 これらの部分空間同定法の開発は、多入力多出力振動系のモードパラメータを推定することの難しさに起因しています。 ここ数年，システム同定の分野で亜空間法が注目されているのは，基本的に非反復的で高速だからである (Favoreel et al., 1999)． そのため、収束の問題が発生せず、また、線形代数の安定した技術に基づくだけなので、数値的にもロバストである。 これらの方法は、固有値または特異値分解を用いてデータの実質的なフィルタリングを行い、特に間隔の狭いモードが存在する場合に有効である。 本質的には、データは直交信号空間とヌル部分空間に分離され、そのどちらかをモードパラメータの推定に使用することができる(Tasker; Bosse and Fisher, 1998)。
名称
m = 出力数
n = 自由度
r = 外部加振数
k = サンプル時間
M = 時間窓内のサンプル数
N = 時間窓内のサンプル数
N = 時間窓内のサンプル数
m = 外部加振数
r = 外部加振数
R = 外部加振数
M = 時間窓内のサンプル数
u(k) = 入力ベクトル
x(k) = 状態ベクトル
y(k) = 応答ベクトル
Ad = システム行列
Bd = 入力行列
Cd = 出力行列
Dd = 直接伝達行列
G = ブロックトープリッツ行列8099 行列
I = 恒等行列
R = 左特異ベクトル行列
S = 右特異ベクトル行列
U = ブロックハンケル行列
X = 状態列の行列
Y = ブロックハンケルマトリックス 出力の
0 = ヌル行列
ギリシャ語記号
G = 拡張観測可能行列
å =特異値の行列
下付き文字
sシフト
最初の2n列
上付き文字
–
T transpose
^ orthogonal
pseudoinverse
物理モデル
物理モデルはNACA 0012 airfoil sectionとフレキシブルマウントシステムに関連した剛体長方形の翼である。 柔軟な取り付けシステムは、剛体翼がフラッターに遭遇する、明確に定義された2自由度の力学系を提供します。 図1には、フラッター取り付けシステムの側面図と透視図が示されている。 フラッター搭載システムは、Dansberryら（1993）で開発されたシステムと同様に、4本の円形ロッドと中央の平板支柱で支持された移動プレートからなる。
ロッドと平板は弾性拘束を与え、移動プレートに固定された硬翼モデルはフラッターが発生するとピッチとプランジの2自由度モードで振動することになる。 ロッド、平板、移動板は鋼鉄製で、すべての接続は固定端です。 翼の模型と移動板はアルミニウム製で、後縁フラップはABS樹脂製である。 寸法は、ロッドが直径0.0055 m、ムービングプレートが0.6 ´ 0.3 m、フラットプレートが0.7 ´ 0.1 ´ 0.002 m、翼モデルが0.8´ 0.45 mである。8033>
フラッターマウントシステムのウィンドオフ特性は、平板支柱の寸法、ロッド、移動板と翼モデルの質量に強く影響されます。 平板支柱とロッドの長さと断面を変更することにより、フレキシブルマウント・システムの周波数とモードシェイプが変更されます。 重りを追加して、フレキシブルマウントと翼のモデルの重心をシステムの弾性軸に移動させることで、ピッチモードとプランジモードを切り離すことができます。 システム弾性軸は、平板支柱の垂直中心線と可動板の中心に位置しています。 また、4本のロッドは風洞壁に対して平行なピッチとプランジ変位を保証する。
フレキシブルシステムを設計するために、ソフトウェアAnsysâを使用して有限要素モデルを開発した。 2種類の要素が使用された。 また、梁4とシェル63をそれぞれ棒材と平板支柱に使用しました。 また、フレキシブルマウントのロッドと平板ストラットの基部には、片持ち梁境界条件を採用しました。 FEM から得られた実験システムの寸法と動的特性は、このシステムの空力弾性挙動を利用可能な風洞に合わせることができるまで修正された。 このシステムの空力弾性挙動は、De Marqui Jr, Belo and Marques (2005)に記載されている数学モデルでシミュレーションされた。
実験装置の設計と構築後、風洞フラッター試験前に固有周波数とモードを確認するために実験モード解析が行われた。 この試験では、25Hz以下の周波数を調査し、翼の制御面をロックさせた。 測定点は、システムに弾性的な制約を与える平板ストラットに配置されています。 実験データからモードシェイプと周波数を同定するために，Tsunaki (1999) によって修正された Eigensystem Realization Algorithm (ERA) が採用されている． 最も重要な固有振動数を表1に示す。 1.
表1より、第1曲げおよび第1ねじりモードがよく定義されており、さらに第3モードがそれらよりも高い周波数で定義されていることがわかる。 理論的には、この条件により風洞試験中の2自由度系が保証され、高次モードは風洞試験中に大きく励起されることはない(Dansberry et al., 1993)。 フレキシブルマウントの設計手順の詳細と、より多くの結果は、De Marqui Jr et al. 明らかに、フラッター解析では、これらの特性と空気力学的特性の相互作用を考慮しなければならない。 空力的な力とモーメント、この研究の場合は揚力とピッチモーメントは、古典的な曲げねじりフラッターに関与するモードを刺激することになります。 その結果、構造物の弾性特性と、機械的摩擦がないと仮定した場合の空力的減衰の原因となり、後流渦によって引き起こされる上昇流によって生じる空力的復元荷重が反応し、気流にエネルギーを消散させることになります。
実験システムは、図1に見られるように、マウントシステムと関連した翼に、2つの歪みゲージと3つの加速度計が装備されています。 1台の加速度計(Kistler KBeam 8303A10M4)は平板支柱の中心線に設置され、プランジ加速度を測定しています。 他の2つの加速度ピックアップ(Kistler KBeam 8304B10)は、移動プレートに設置されています。 これらの加速度計で測定した信号からピッチ加速度を計算する。
ひずみゲージは、有限要素解析から求めた最大ひずみ位置の平板支柱の中心線に設置されている。 8033>
後縁フラップの駆動には、可動板下面（図1参照）に設置したブラシレス電気モータ（Thompson BLD2315B10200）を使用しています。 フラップはモータとシャフトで接続されている。 電気モーターにはエンコーダが搭載されており、フラップの実際の角度位置を測定するために使用される。 PIDコントローラーは、モーターによる後縁フラップの位置の正しい制御を保証するために調整された。
Extended Eigensystem Realization Algorithm – EERA
n自由度を持つ任意の線形時不変力学系は、以下の離散時間状態空間方程式でモデル化することができる。
ここで x(k) はk番目のサンプル瞬間における2n次元の状態ベクトル、 u(k) はr次元の入力ベクトル、rは外部励振の数、 y(k) はm次元応答ベクトル、mはシステムの出力または応答の数、Adは2n´2nシステム行列、Bdは2n´r入力行列、Cdはm´2n出力行列、Ddはm´r直接透過行列である。
EERAを用いた同定手順は、入出力の時刻歴からシステム行列Adを決定することである。 フラッターに関連する特徴、すなわち周波数と減衰はシステムマトリックスAdを用いて推定することができる。 EERA法を用いたシステム行列Adの同定は、Tasker; Bosse and Fischer（1998）によって示された理論に基づいて、以下の手順で記述される。
入力（U）および出力（Y）のブロックハンケル行列は、入力および出力時間（Overschee and De Moor, 1996）から直接得ることができる（8033>
ここで、MおよびNは、同定プロセス中に使用する時間窓内のサンプルの数である。
出力のブロックハンケル行列はVerhaegen and Dewilde (1992),
Gは拡張観測性行列、Xは状態列の行列、Gはマルコフパラメータまたはインパルス応答のブロックトープリッツ行列、つまり…で説明されているように表されることが確認される。
定義により、直交行列は、(Van Overschee and De Moor, 1996),
式（3）に式（3）の右項と左項を後乗算して書くことができる。 (5)にそれぞれ乗じ、直交性の定義を用いると、以下の式が得られる。
特異値分解を適用すると：
ここでR（mM ´ mM）は左特異ベクタ行列、S（N ´ N）は右特異ベクタ行列であり、対応する特異値行列である。 これらの行列の列は正規直交である。
の擬似逆行列は式(1)から求めることができる。(7)から得られる：
一方、
この時点で、出力または応答のブロックハンケル行列のシフトした形式を導入できる：
この新しい行列の寸法は、同定処理で使用する出力時間履歴ベクトルの長さに関連しています（時間窓内のサンプルの数）。 ただし、この窓は時間的に1ステップ以上進まなければならない。
式と同様な方法で、この新しい行列を作る。 (3)と同様に、以下のようになる：
ここで、GsとGsはそれぞれ拡張観測性行列とマルコフパラメータのブロックトープリッツ行列のシフト版：
式と同じ導出に従っている：GsとGsは、それぞれ、マルコフパラメータをシフトした行列である：(1)式に使用されている。 (6)と同様にして、
ここで、この式の右辺の項は、元の観測性行列とシフトしたものを比較して、すなわち、Gs＝GAdと簡単に求めることができる。
式（13）の行列YsU^は、
式（14）に式（7）と式（8）を代入すると、
この段階で、必要な特異値の数を決定する基準が規定され得る。 この数は、同定作業に伴う困難さに応じて変更することができる。 この数は同定されたモデルの次元を確立し、それは同定問題の間に修正されなければならない。 特異値の数が2nとして決定されることを考慮すると、特異値行列は次のように表すことができる：
ここで、R2nはRの最初の2n列を含み、S2nはSの最初の2n列を含んで、行列は
と簡便に書くことができる。
行列R2nとS2nは次の関係を満たす：
式（1）の関係を使用することにより、R2nはRの最初の2n列を含み、S2nはSの最初の2n列を含む。 (17)を特異値分解問題に用いると、次のようになる：
そして、S=S-1 (Watkins, 1991)とすると、次のようになる：
それを考慮して、式（19）と式（20）を式に代入する：
そして式に代入する：(20)式に代入する：
そして式に代入する。 (15)
ここで
式（18）から、
式（24）は、式と比較することができ、式は
となり、式は
から
となり、式は
となる。 (13)と比較すると、システム行列は次のように評価できる：
システム行列Adは、システムの最小実現値である。 この行列の次元は2nであり、それはまた同定されたシステムの次元を決定する。 この実現はモード座標の状態方程式に変換でき、固有値を計算することにより固有振動数と減衰を求めることができる。 上記の式は、入力項の存在によってのみERA式と異なっている。 応答がインパルス入力による場合は、ERAで観測された式と同じになります（Juang, 1994）。
フラッター実験検証
フラップのリアルタイム制御の開発とデータ取得には、dSPACE® DS 1103プロセッサボードを使用しました。 このボードは400MHzのPower PC 604eプロセッサ、16 A/Dおよび8 D/Aチャンネルを備えたI/Oインターフェース、およびインクリメンタルエンコーダインターフェースを備えています（DSPACE®, 2001）。 加速度計、ストレインゲージブリッジ、フラップ位置の信号を同時に取得することができます。 データ取得と信号処理のための計算コードは、Matlab/Simulink®で開発されています。 Simulink®のコードは、Matlab®でReal-Time Workshop®コンパイラを使ってコンパイルされ、Cコードになります。 このCコードがdSPACE®ボードにダウンロードされ、信号処理と入出力制御が実行されます。
図2は、データ収集システムの簡略化されたスキームを示しています。 計算システムのゲインは、測定された信号を必要な物理単位（加速度計ではmVからm/s2またはrad/s2、ひずみゲージではmVからmまたはrad）に変換するために使用されます。 後縁フラップの駆動に使用される電気モーターのエンコーダは1000ラインです。 したがって，後縁位置の測定において0.36度の分解能を達成することができる． 実験では、1秒間に1000サンプルの取得レートを採用しています。
最初の実験では、臨界フラッター速度の検証を行いました。 風洞速度を徐々に上げ、dSPACE®システムを用いてピッチ信号とプランジ信号を測定します。 風洞速度は、試験室に設置したベッツマノメーター、気圧計、温度センサーに接続した静止ピトー管で計測した圧力から得られます。 フラッターは臨界流速である25m/sで観測され、振動挙動が計測されます。 図3は実験中に測定されたピッチとプランジの信号をそれぞれ示す。
フラッター現象の特徴の一つは、現象に関与するモード、すなわちこの場合ピッチとプランジが結合していることである。 この条件は図4で検証され、図3で示された時間領域の信号がその周波数内容で示されている。
このテストは臨界速度でのみシステムの挙動を示すものであった。 しかし、いくつかの力学的特性は風洞流速の増加とともに変化する。 これらの変化を確認するために、他の試験を実施した。 基本的に、周波数応答関数は、速度が増加するにつれて最初の曲げモードとねじりモードの発展を示すいくつかの速度で取得されます。 これらの試験で考慮される入力信号は後縁の位置で、出力信号は翼の後縁で測定された加速度です。
A B&K デュアルチャンネル デジタルスペクトルアナライザー 2032型は周波数応答を得るために使用されます。 これらの応答は風洞のオフ状態から臨界速度にできるだけ近い速度まで得られたものである。 信号入力は、dSPACE®システムで生成されたホワイトノイズで、後縁フラップに送られます。 この信号と加速度は、スペクトルアナライザで処理されます。 図5では、風洞速度の増加に伴うモードの変化を確認することができます。 速度ゼロで得られた周波数応答は、最初の曲げモードとねじりモードに関するピークを示し、予想通り、EMAで得られたものと同じ固有振動数を示しています。 臨界速度付近で測定された最後の周波数応答では、フラッターに関与するモード間のカップリングの傾向を確認することができます。 この結合は約1.6Hzの周波数で発生する傾向があり、図4の結果を裏付けている。
中間速度で得られた周波数応答では、ピッチとプランジの周波数の変化が観察される。 またピッチモードとプランジモードのピークは速度ゼロでの周波数応答のピークほどシャープでないことがわかる。 この事実は、減衰増加に対する流体構造間相互作用の影響と見ることができる。 この傾向は、減衰が消滅してフラッターが発生すると予想される臨界速度付近まで予想される。
同定結果
拡張固有値実現アルゴリズム（EERA）を用いて、フラッターに関わるモードに対する風洞増加速度による周波数および減衰値の変化を定量化した。 EERAを用いて風速変化に伴う減衰の変化を調べることで、フラッターの発生時期を予測することができる。 同定に使用するデータは、本論文で紹介した周波数応答関数を得るために実施した空力試験で取得したものです。 周波数領域試験と同時に、入力信号（後縁のフラップ動作）と歪みゲージで計測した信号（ピッチとプランジの変位）を、dSPACE®収集システムを使用して時間領域で取り込みました。 図6～図8は、ある風洞試験で計測された入力信号と出力信号の例です。 図6では、フラップのたわみ量（度）が描かれています。 これは、空力弾性系への加振として働くように、フラップ角に対してランダムに生成された（一様分布の）信号を表しています。 フラップ動作（図6参照）に対するプランジ応答とピッチ応答をそれぞれ図7と図8に示す。
入出力時間領域データを取得した後に識別処理を実施した。 入出力のブロックハンケル行列の次元（MとN=2M）と考慮すべき特異値の数（2n）は、各流速に対して行われた同定ごとに変更された。 この変動は、モードが結合しつつある高い風速で取得されたデータを用いたパラメータの同定が困難であることによって説明できる
同定プロセスで得られた最終結果は、図9で観察することができる。 ピッチ、プランジ周波数と減衰係数の風速に対する変化を示している。 実験結果（前項参照）と同様に、対気速度25m/s付近でフラッタが予測されることがわかる。 各試験において、ピッチ運動とプランジ運動の周波数と減衰係数は、様々な同定パラメータに対する平均値で求められ、異なる同定システム状態行列を導きました。 図9では、点の雲が同定パラメータの変動に関係し、曲線が周波数と減衰の平均値を表している。 周波数の計算では、EERA法が様々な同定パラメータに対して良好な予測を提供できていることがわかります。 それにもかかわらず、減衰係数の同定では、対気速度ごとの値はより分散していました。 また、ピッチモードの減衰係数は、プランジモードに比べ、ばらつきが少ないようです。 この理由はまだ解明されておらず、EERAを用いたフラッター予測の継続的な検討の対象とする必要があります。 これらの結果は、周波数に関するものよりも劣るかもしれませんが、平均減衰値は、古典的な2次元フラッターの物理と一致するカーブを示しています。
結論
風洞実験による空力弾性試験は、フラッターパラメータの同定に使用されました。 風洞実験によりフラッターの特性を調べ、時間領域と周波数領域で現象を観察することができた。 時間領域では、フラッターの自己保持振動挙動が示された。 周波数領域では、風洞の速度上昇に伴うモードの進化が観察されました。 臨界速度では、カップリング傾向が明確に示された。 ピッチとプランジの減衰の変化は、これらの試験において定性的な方法で得られた。
速度の増加に伴うピッチとプランジのモードの変化を定量化するために、同定法を適用した。 拡張固有値実現アルゴリズムが、フラッター特性評価のために行われた試験で得られた入出力データを時間領域で使用して採用された。 この方法は、フラッターパラメータの同定に使用され、速度に関する性能と、プロセス中に起こりうる数値問題を検証するために使用されました。 EERAの使用は、この同定法で得られた結果と、これまでの風洞試験で得られた結果との整合性を考慮すると、適切であったと言える。 特にプランジモードの減衰係数の同定には困難があった。
この研究で示された同定プロセスがオフラインであることを考慮しても、これまでに得られた結果は、風洞試験中のフラッターパラメータのオンライン同定が検討可能であることを示唆している。 8033>
謝辞
著者らはCAPESとFAPESP（サンパウロ州立研究支援財団）から契約番号1999/04980-0と2000/00390-3で財政支援を受けたことを感謝している。 1970, “Aeroelasticity”, Applied Mechanics Reviews, pp.119-129.
Collar, A.R., 1959, “Aeroelasticity Retrospect and Prospect”, The Journal of the Royal Aeronautical Society, Vol.63, No.577, pp.1-15, 1959.航空弾性に関する研究。
Cooper, J.E. and Crowther, W.J., 1999, “Flutter Speed Prediction During Flight Testing Using Neural Networks”, CEAS/AIAA/ICASE/NASA Langley International Forum on Aeroelasticity and Structural Dynamics, pp.255-264。
Cooper, J.E., Emmett, P.R. and Wright, J.R., 1993, “Envelope Function: また，”Flutter Data “については，”A Tool for Analysing Flutter Data”, Journal of Aircraft, Vol.30, No.5, pp.785-790.
Dansberry, B.E., Durham, M.H., Bennett, R.M., Turnock, D.L., Silva, E.A. and Rivera Jr, J.A., 1993, “Physical Properties of the Benchmark Models Program Supercritical Wing”, NASA TM-4457.
De Marqui Jr., C., Belo, E.M. and Marques, F.D., 2005, “A flutter suppression active controller”, Proc I.Mech.E Part G – Journal of Aerospace Engineering, Vol. 219.
De Marqui Jr., C., Belo, E.M., Tsunaki, R.H., Rebolho, D.C. and Marques, F.D., 2004, “Design and Tests of an Experimental Flutter Mount Systems”, Proceedings of the XXII IMAC, Dearborn, MI.
DS1103 PPC Controller Board, 2001, Hardware Reference, www.dspace.de.
Favoreel, W., Huffel, S.V., De Moor, B, Sima, V. and Verhaegen, M., 1999, “Comparative study between three subspace identification algorithms”, Proceedings of the European Control Conference, Karlsruhe-Germany, 31st August-3rd September, 6p.
Försching, H., 1979, “Aeroelastic Problems in Aircraft Desing”, von Karman Institute for Fluid Dynamics, Lecture series 08: A Survey of Aeroelastic Problems.
Garrick, I.E. and Reed, W.H., 1981, “Historical Development of Aircraft Flutter”, Journal of Aircraft, Vol.18, No.11, pp.897-912.など。
Garrick, I.E., 1976, “Aeroelasticity frontiers and beyond”, 13th Von Karman Lecture, Journal of Aircraft, Vol.13, No.9, pp.641-657.
Juang, J.N., 1994, “Applied System Identification”, Prentice Hall PTR, New Jersey, USA.
Kehoe、M.W.、1995、”A Historical Overview of Flight Flutter Testing”、NASA TM-4720.
Ko, J., Kurdila, A.J. and Strganac, T.J., 1997, “Adaptive Feedback Linearization for the Control of a Typical Wing Section with Structural Nonlinearity”, ASME International Mechanical Engineering Congress and Exposition, Dallas, Texas.
Lind, R. and Brenner, M., 2000, “Flutterometer: このような場合、「このような場合、どのようにすればよいのか？
Lind, R., 2003, “Flight-Test Evaluation of Flutter Prediction Methods”, Journal of Aircraft, Vol.40, No.5, pp.964-970.
Mukhopadhyay, V., 1995, “Flutter suppression control law design and testing for the active flexible wing”, Journal of Aircraft, Vol.32, No.1, pp 45-51.
Tasker, F., Bosse, A. and Fisher, S., 1998, “Real-time modal parameters estimation using subspace methods.” (部分空間法を用いた実時間モードパラメータ推定): 理論”. 機械システムと信号処理, Vol.12, No.6, pp.797-808.
鳥居秀樹・松崎洋一郎，2001，”離散時間システムのフラッターマージン評価”，航空機，Vol.38，No.1，pp.42-47.．
Tsunaki, R. H., 1999, “Identificação Automatizada de Modelos Dinâmicos no Espaço de Estados”, Ph.D. Thesis, University of Sao Paulo.
Van Overschee, P. and De Moor, B., 1996, “Subspace identification for linear systems: theory, implementation, applications”, Kluwer Academic Publishers, Boston, United States.より引用。
Verhaegen, M. and Dewilde, P., 1992, “Subspace model identification part 1. 出力誤差状態空間モデル同定アルゴリズムクラス”, 国際制御学会誌, Vol.56, pp.1187-1210.
Waszak, M.R., 1998, “Modeling the benchmark active control technology wind-tunnel model for active control design applications”, NASA TP-1998-206270.
Watkins, D.S., 1991, “Fundamental of Matrix Computations”, New York, Wiley,USA.
Wright, J.R., 1991, “Introduction to Flutter of Winged Aircraft”, von Karman Institute for Fluid Dynamics, Lecture series 01: Elementary Flutter Analysis.日本学術振興会特別研究員奨励賞受賞。
Zimmerman, N.H. and Weissnburger, J.T., 1964, “Prediction of Flutter Onset Speed Based on Flight Testing at Subcritical Speeds”, Journal of Aircraft, Vol.1, No.4, pp.190-202。