柔軟翼機

この例では,ミネソタ大学無人航空機研究室で作製された小型の柔軟全翼機から収集したデータについて検討します[1]。航空機のジオメトリを以下に示します。

航空機の翼は,荷重により大きく変形する可能性があります。柔軟なモード周波数は,剛体の翼をもつ一般的な航空機より低くなります。この柔軟な設計により,材料コストが削減され,航空機の動きはより機敏になり,飛行距離も長くなります。ただし,柔軟なモードを制御しないと,空力弾性が非常に不安定になる(揺れが生じる)的可能性があります。こうした不安定さを抑制するための効果的な制御則を設計するには,翼のさまざまな曲げモードを正確に特定する必要があります。

実験の設定

実験の目的は,外部の励起に応じてさまざまな場所で航空機の振動応答を収集することです。航空機は,重心にある1個のバネを使用して木枠から吊るしてあります。このバネには十分な柔軟性があり,バネ——質量振動の固有振動数が航空機の基本周波数に干渉することはありません。入力の力は,航空機の中心近くのUnholtz-Dickie模型20動電型シェーカー経由で加えられます。

次の図に示すように,20個のpcb - 353 b16转椅加速度計を翼長に沿って配置し,振動応答を収集します。

シェーカーの入力コマンドは, $\mathit{一个}\mathrm{罪}（\omega （\mathit{t}）\mathit{t}）$という形式の定振幅チャープ入力として指定します。チャープ周波数は時間とともに線形に変化し, $\omega （\mathit{t}）＝{\mathit{c}}_0+{\mathit{c}}_1\mathit{t}$となります。入力信号でカバーされる周波数範囲は3 ~ 35赫兹です。データは2個の加速度計(1つのx位置の前方と後方のエッジにある加速度計)で同時に収集されます。つまり10回の実験を行って,20個全部の加速度計の応答が収集されます。加速度計と力変換器の測定値はすべて,2000 Hzでサンプリングされます。

データの準備

データは2出力/ 1入力の信号の10セットで表現され,それぞれに600 kのサンプルが含まれます。データはMathWorksサポートファイルサイトで入手できます。免責事項を参照してください。データファイルをダウンロードし,データをMATLAB®ワークスペースに読み込みます。

url =“//www.tatmou.com/金宝appsupportfiles/ident/flexible-wing-GVT-data/FlexibleWingData.mat”；websave (“FlexibleWingData.mat”url);负载FlexibleWingData.matMeasuredData

変数MeasuredDataは,フィールドExp1、Exp2、...、Exp10をもつ構造体です。各フィールドは2つの応答および対応する入力の力データを含むフィールドyおよびuをもつ構造体です。最初の実験のデータをプロットします。

fs = 2000;数据采样频率Ts = 1 / f;%样品时间y = MeasuredData.Exp1.y;%输出数据(2列，每个加速度计一列)u = MeasuredData.Exp1.u;输入力数据%t = (0:length(u)-1)' * t;图subplot(211) plot(t,y) ylabel(“输出(m / s ^ 2)”)传说(“前沿”,“后缘”) subplot(212) plot(t,u) ylabel(“输入”)包含(的时间(秒)）

モデルの同定用にデータを準備するために,データはiddataオブジェクトにパッケージ化されます。iddataオブジェクトは,系统辨识工具箱™で時間領域データをパッケージ化するための標準的な手段です。入力信号は帯域幅を制限して処理されます。

ExpNames =字段名(MeasuredData);数据= cell(1, 10);为k = 1:10 y = MeasuredData.(ExpNames{k}).y;k u = MeasuredData。(ExpNames {}) .u;数据{k} = iddata(y, u, t, t)“InterSample”,“提单”）;结束

データセットオブジェクトを1つの複数実験データオブジェクトに結合します。

Data =合并(数据{:})

Data =包含10个实验的时域数据集。样本时间Exp1 600001 0.0005 Exp2 600001 0.0005 Exp3 600001 0.0005 Exp4 600001 0.0005 Exp5 600001 0.0005 Exp6 600001 0.0005 Exp7 600001 0.0005 Exp8 600001 0.0005 Exp9 600001 0.0005 Exp10 600001 0.0005输出单元(如果指定)y1 y2输入单元(如果指定)u1

モデルの同定

周波数応答が実際の航空機と可能な限り一致する動的モデルを同定します。動的モデルは,システムの入力と出力の数学的関係を微分または差分方程式としてカプセル化します。動的モデルの例は,伝達関数と状態空間モデルです。系统辨识工具箱では,動的モデルはidtf(伝達関数の場合),idpoly(AR, ARMAモデルの場合),中的难点(状態空間モデルの場合)などのオブジェクトによってカプセル化されます。動的モデルは,時間領域または周波数領域の測定データに対して特遣部队や党卫军のような推定コマンドを実行することによって作成できます。

この例ではまず,etfeコマンドを使用した経験的な伝達関数の推定により,時間領域の測定データを周波数応答データに変換します。次に,推定され润扬悬索桥たを使用して,航空機の振動応答の状態空間モデルを同定します。動的モデルの同定に時間領域データを直接使用することも可能です。润扬悬索桥ただし,形式のデータを使用すると,大きいデータセットを少数のサンプルに圧縮できると同時に,関連する周波数範囲に合わせて推定動作を調整することも簡単になります。频はidfrdオブジェクトによってカプセル化されます。

データ実験ごとに2出力/ 1入力の周波数応答関数(降维)を1つ推定します。ウィンドウ処理は使用しません。応答の計算に24000個の周波数点を使用します。

G = cell(1,10);N = 24000;为k = 1:10%使用ETFE命令将时域数据转换为频响Data_k = getexp(Data, k); / /数据G{k} = etfe(Data_k， []， N);% G{k}是一个@idfrd对象结束

润扬悬索桥すべてのを単一の20出力/ 1入力润扬悬索桥のに連結します。

G = cat(1, G{:});连接所有估计频响函数的输出G.OutputName =“y”；% name输出'y1'， 'y2'，…,“日元”G.InterSample =“提单”；

推定された周波数応答の感覚をつかむために,出力1と15(任意に選択)の振幅をプロットします。対象の周波数範囲(4 ~ 35赫兹)を拡大します。

选择= bodeoptions;%绘图选项opt.FreqUnits =“赫兹”；%表示频率，单位为Hzopt.PhaseVisible =“关闭”；%不显示相位输出:[1 15];%选择输出1和15绘图clf bodeplot(G(OutputNum，:)， opt)%绘制频率响应网格xlim (35 [4])在

频は少なくとも9つの共振周波数を示しています。解析のため,ここでは,航空機の最も重大な柔軟な曲げモードが存在する6 ~ 35赫兹の周波数スパンに注目します。润扬悬索桥そのため,をこの周波数領域まで減らします。

f = G.Frequency / 2 /π;%提取频率矢量，Hz (G存储频率，rad/s)G = fselect(G, f>6 & f<=32)% "fselect"选择请求范围内的频响(6.5 - 35 Hz)

Gs = IDFRD模型。包含20输出和1输入的频率响应数据。响应数据可在624个频率点获得，范围从37.96 rad/s到201.1 rad/s。样品时间:0.0005秒输出通道:‘y(1)’,‘y(2)’,‘y(3)’,‘y(4)’,‘y(5)’,‘y(6)’,‘y(7)’,‘y(8)’,‘y(9)’,‘y(10)’,‘y(11)’,‘y(12)’,‘y(13)’,‘y(14)’,‘y(15)’,‘y(16)’,‘y(17)’,‘y(18)’,‘y(19)’,‘y(20)输入通道:“u1”状态:估计模型

したがって,Gsには20箇所の測定位置の周波数応答測定値が含まれます。次に,Gsに適合する状態空間モデルを同定します。部分空間推定器n4sidは,非反復の推定を素早く行います。状態空間モデルの構造は,次のように構成します。

n4sidOptionsを使用して,n4sidの推定オプションを設定します。

频=挤压(Gs.ResponseData);权重=平均(1. /√(abs(降维)))。”;n4Opt = n4sidOptions (“EstimateCovariance”假的,.．.“WeightingFilter”权重,.．.“OutputWeight”、眼睛(20));sys1 = n4sid (Gs, 24岁,“t”0,“引线”,真的,n4Opt);适合= sys1.Report.Fit.FitPercent '

适合=1×2057.0200 57.9879 85.0160 86.3815 80.4879 80.4430 58.2216 45.2692 61.5057 76.7612 84.7305 86.2600 86.4266 85.0251 76.9208 82.1191 74.7982 79.6837 67.9078 76.7249

“适合”の数値は,データ(Gs)とモデル(sys1)の周波数応答の間の適合の割合を,正規化平方根平均二乗誤差(NRMSE)適合度の尺度を使用して示しています。最低と最良の適合が次にプロットされます。

[~, Imin] = min(配合);[~, Imax] = max(配合);clf bodeplot (Gs (Imin, Imax,:), sys1 ((Imin, Imax):),选择);xlim(32[6])标题(“最坏(上)和最好(下)符合数据和模型”网格)在传奇(“数据”,“模型”）

モデルsys1で達成された適合は,モデルのパラメーターを非線形最小二乗法で反復計算することによって,さらに改善できます。これは,党卫军コマンドを使用して達成できます。ssestOptionsコマンドを使用して,党卫军の推定オプションを設定します。今回は,パラメーター共分散も推定します。

ssOpt = ssestOptions (“EstimateCovariance”,真的,.．.“WeightingFilter”, n4Opt。WeightingFilter,.．.“SearchMethod”,“lm”,.．.%使用Levenberg-Marquardt搜索方法“显示”,“上”,.．.“OutputWeight”、眼睛(20));sys2 = ssest(Gs, sys1, ssOpt);%估计状态空间模型(这需要几分钟)适合= sys2.Report.Fit.FitPercent '

适合=1×2089.7225 89.5185 89.7260 90.4986 88.5522 88.8727 81.3225 83.5975 75.9215 83.1763 91.1358 89.7310 90.6844 89.8444 89.6685 89.1467 87.8532 88.0385 84.2898 83.3578

前回と同様に,最悪と最良の適合をプロットします。また,モデルの周波数応答の1-sd信頼領域も可視化します。

[~， Imin] = min(Fit);[~， Imax] = max(Fit);clf h = bodeplot (Gs (Imin, Imax,:), sys2 ((Imin, Imax):),选择);showConfidence(h, 1) xlim([6.5 35]) title(“最坏(上)和最好(下)符合数据和细化模型之间的关系”网格)在传奇(“数据”,“模型(sys2)”）

調整済みの状態空間モデルsys2は7 ~ 20 Hzの領域润扬悬索桥でに非常によく適合します。共振が最大の位置の周辺で,不確かさ範囲が狭くなっています。ここでは24次のモデルを推定しましたが,これは,システムsys2内に最大で12の振動モードがあるということです。modalfitコマンドを使用して,これらのモードの固有振動数を赫兹単位で取得します。

f = Gs.Frequency / 2 /π;数据频率% (Hz)Fn = modalfit(sys2, f, 12);固有频率% (Hz)disp (fn)

7.7721 7.7953 9.3147 9.4009 9.4910 15.3463 19.3291 23.0219 27.4164 28.7256 31.7014 63.3034

fnの値は2つの周波数が7.8赫兹に非常に近く,3つが9.4赫兹に近いことを示しています。これらの周波数に近い周波数応答を調べると,出力ごとにピークの位置が少しずれることがわかります。実験プロセスの制御を強化し,これらの周波数を中心とする狭い範囲に制限した入力帯域幅で時間領域同定を直接実行するか,これらの周波数に近い周波数応答に対して単一の振動モードを適合させることによって,こうしたばらつきを除去できる場合があります。こうした代替方法については,この例では扱っていません。

モーダルパラメーターの計算

モデルsys2を使用してモーダルパラメーターを抽出できるようになりました。频を調べると,おおむね[5 7 10 15 17日23日27日30]赫兹の周波数に近い10個前後のモーダル周波数が確認できます。この評価をより具体的にするために,modalsdコマンドを使用して,基となるモデルの次数が変更されたときにモーダルパラメーターの安定性をチェックすることができます。この操作には長時間かかります。そのため,結果のプロットはイメージとして直接挿入されます。図を再現するには,以下のコメントブロック内のコードを実行します。

频=排列(Gs。ResponseData,中[3 1 2]);f = Gs.Frequency / 2 /π;% {数字润扬悬索桥pf = modalsd (f, f, MaxModes, 20日“FitMethod”,“lsrf”);%}

プロットとpfの値を調べると,実際の固有振動数の調整されたリストが提示されます。

频率= [7.8 9.4 15.3 19.3 23.0 27.3 29.2 31.7];

ここでは,モデルsys2から最も主要なモードを選択するための指針として,频率ベクトルの値を使用します。これを実行するにはmodalfitコマンドを使用します。

(fn,博士、女士)= modalfit (sys2, f,长度(频率),“PhysFreq”、频率);

fnは固有振動数(Hz単位),博士は対応する減衰係数,女士は列ベクトル(固有振動数ごとに1列ずつ)として正規化された残差です。これらのモーダルパラメーターの抽出プロセスでは,安定していて減衰率の低いモデルの極のみが使用されます。女士列には,正の虚数部をもつ極のみに対応するデータが含まれます。

モード形状の可視化

さまざまな曲げモードを可視化するには,上記で抽出したモーダルパラメーターが必要です。さらに,測定点の位置に関する情報も必要です。これらの位置(x - y値)は行列AccePosの加速度計ごとに記録されます。

AccelPos = [.．.%见图216.63 - 18.48;离中心最近的右边16.63 - 24.48;27.90 - 22.22;27.90 - 28.22;37.17 - 25.97;37.17 - 31.97;46.44 - 29.71;46.44 - 35.71;55.71 - 33.46;55.71 - 39.46;%最远的权利-16.63 - 18.48;最靠近中间偏左的百分比-16.63 - 24.18;-27.90 - 22.22;-27.90 - 28.22;-37.17 - 25.97;-37.17 - 31.97;-46.44 - 29.71;-46.44 - 35.71;-55.71 - 33.46;-55.71 - 39.46);%最远的离开

モード形状が含まれている行列女士では,各列が1つの周波数の形状に対応しています。モード形状の振幅をセンサー座標の上に重ね,モードの固有振動数における振幅を変えることにより,モードをアニメーション化します。アニメーションは,減衰のない曲げを示します。例として,15.3赫兹におけるモードについて考えます。

AnimateOneMode(3, fn, ms, AccelPos);

まとめ

この例は,モーダルパラメーター推定に対するパラメトリックモデル同定ベースのアプローチを示しています。24次の状態空間モデルを使用して,6 ~ 32个赫兹の周波数範囲で8個の安定した振動モードを抽出しました。モーダル情報を加速度計の位置の情報と結び付けて,さまざまな曲げモードを可視化しました。これらのモードの一部を以下の図に示します。

参考文献

Gupta, Abhineet, Peter J. Seiler和Brian P. Danowsky。柔性飞行翼飞机地面振动试验AIAA大气飞行力学大会，AIAA科技论坛。张仁(2016 - 1753)。

函数AnimateOneMode(ModeNum, fn, ModeShapes, AccelPos)%动画一个模式。% ModeNum:模式索引。重新排序传感器位置，以便绘制一个连续的、%不相交曲线绕机身轨迹。PlotOrder = [19: 2:11 1:2:9 10: 2:2, 12:2:20, 19);Fwd = PlotOrder (1:10);尾= PlotOrder (20: 1:11);x = AccelPos (PlotOrder, 1);y = AccelPos (PlotOrder 2);船尾xAft = AccelPos (1);yAft = AccelPos(船尾,2);xFwd = AccelPos(录象,1);yFwd = AccelPos(前轮驱动,2);wn = fn (ModeNum) * 2 *π;%模式频率，单位为rad/秒T = 1 / fn (ModeNum);%模态振荡周期Np = 2.25;%要模拟的周期数达峰时间= Np * T;%模拟Np周期元= 100;%动画的时间步数t = linspace(0,达峰时间,Nt);ThisMode = ModeShapes (:, ModeNum) / max (abs (ModeShapes (:, ModeNum)));标绘归一化%z0 = ThisMode (PlotOrder);zFwd = ThisMode (Fwd);zAft = ThisMode(尾);CLF col1 = [1 1 1]*0.5;xx =重塑([[xAft, xFwd]';南(2,10)],[2]20);yy =重塑([[yAft, yFwd]';南(2,10)],[2]20);plot3 (x, y, 0 * z0,“- - -”0 * z0, x, y,,“。”， xx(:)， yy(:)， 0 (40,1)，“——”,.．.“颜色”col1,“线宽”, 1“MarkerSize”10“MarkerEdgeColor”, col1);包含(“x”) ylabel (“y”) zlabel (“振幅”) ht = max(abs(z0));轴([-100 100 10 40 -ht ht])视图(5,55)title(sprintf())% d’模式。FN = %s Hz'、ModeNum num2str (fn (ModeNum), 4)));通过更新z坐标来实现动画。持有在Col2 = [0.87 0.5 0];h1 = plot3 (x, y, 0 * z0,“- - -”0 * z0, x, y,,“。”， xx(:)， yy(:)， 0 (40,1)，“——”,.．.“颜色”、col2“线宽”, 1“MarkerSize”10“MarkerEdgeColor”、col2);h2 = fill3 (x, y, 0 * z0 col2,“FaceAlpha”, 0.6);持有从为k = 1:Nt Rot1 = cos(wn*t(k));Rot2 =罪(wn * t (k));z_t = real(z0)*Rot1 - imag(z0)*Rot2;zAft_t = real(zAft)*Rot1 - image (zAft)*Rot2;zFwd_t = real(zFwd)*Rot1 - imag(zFwd)*Rot2;zz =重塑([[zAft_t, zFwd_t]';南(2,10)],[2]20);集(h1 (1),“ZData”z_t)组(h1 (2),“ZData”z_t)组(h1 (3),“ZData”，zz(:) h2. vertex (:，3) = z_t;暂停(0.1)结束结束