このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
柔軟全翼機のモード解析
この例では、柔軟翼機の曲げモードの計算を示します。翼の振動応答が、翼長の複数の点で収集されます。このデータを使用して、システムの動的モデルを同定します。モーダル パラメーターが、同定されたモデルから抽出されます。モーダル パラメーター データをセンサー位置情報と組み合わせて、翼のさまざまな曲げモードを可視化します。この例では Signal Processing Toolbox™ が必要です。
柔軟翼機
この例では、ミネソタ大学無人航空機研究室で作製された小型の柔軟全翼機から収集したデータについて検討します [1]。航空機のジオメトリを以下に示します。
航空機の翼は、荷重により大きく変形する可能性があります。柔軟なモード周波数は、剛体の翼をもつ一般的な航空機より低くなります。この柔軟な設計により、材料コストが削減され、航空機の動きはより機敏になり、飛行距離も長くなります。ただし、柔軟なモードを制御しないと、空力弾性が非常に不安定になる (揺れが生じる) 可能性があります。こうした不安定さを抑制するための効果的な制御則を設計するには、翼のさまざまな曲げモードを正確に特定する必要があります。
実験の設定
実験の目的は、外部の励起に応じてさまざまな場所で航空機の振動応答を収集することです。航空機は、重心にある 1 個のバネを使用して木枠から吊るしてあります。このバネには十分な柔軟性があり、バネ-質量振動の固有振動数が航空機の基本周波数に干渉することはありません。入力の力は、航空機の中心近くの Unholtz-Dickie Model 20 動電型シェーカー経由で加えられます。
次の図に示すように、20 個の PCB-353B16 加速度計を翼長に沿って配置し、振動応答を収集します。
シェーカーの入力コマンドは、 という形式の定振幅チャープ入力として指定します。チャープ周波数は時間とともに線形に変化し、 となります。入力信号でカバーされる周波数範囲は 3 ~ 35 Hz です。データは 2 個の加速度計 (1 つの x 位置の前方と後方のエッジにある加速度計) で同時に収集されます。つまり、10 回の実験を行って、20 個全部の加速度計の応答が収集されます。加速度計と力変換器の測定値はすべて、2000 Hz でサンプリングされます。
データの準備
データは 2 出力/1 入力の信号の 10 セットで表現され、それぞれに 600K のサンプルが含まれます。データは MathWorks サポート ファイル サイトで入手できます。免責事項を参照してください。データ ファイルをダウンロードし、データを MATLAB® ワークスペースに読み込みます。
url = 'https://www.mathworks.com/supportfiles/ident/flexible-wing-GVT-data/FlexibleWingData.mat'; websave('FlexibleWingData.mat',url); load FlexibleWingData.mat MeasuredData
変数 MeasuredData
は、フィールド Exp1
、Exp2
、...、Exp10
をもつ構造体です。各フィールドは、2 つの応答および対応する入力の力データを含むフィールド y
および u をもつ構造体です。最初の実験のデータをプロットします。
fs = 2000; % data sampling frequency Ts = 1/fs; % sample time y = MeasuredData.Exp1.y; % output data (2 columns, one for each accelerometer) u = MeasuredData.Exp1.u; % input force data t = (0:length(u)-1)' * Ts; figure subplot(211) plot(t,y) ylabel('Outputs (m/s^2)') legend('Leading edge','Trailing edge') subplot(212) plot(t,u) ylabel('Input') xlabel('Time (seconds)')
モデルの同定用にデータを準備するために、データは iddata
オブジェクトにパッケージ化されます。iddata
オブジェクトは、System Identification Toolbox™ で時間領域データをパッケージ化するための標準的な手段です。入力信号は帯域幅を制限して処理されます。
ExpNames = fieldnames(MeasuredData); Data = cell(1, 10); for k = 1:10 y = MeasuredData.(ExpNames{k}).y; u = MeasuredData.(ExpNames{k}).u; Data{k} = iddata(y, u, Ts, 'InterSample', 'bl'); end
データセット オブジェクトを 1 つの複数実験データ オブジェクトに結合します。
Data = merge(Data{:})
Data = Time domain data set containing 10 experiments. Experiment Samples Sample Time 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 Outputs Unit (if specified) y1 y2 Inputs Unit (if specified) u1
モデルの同定
周波数応答が実際の航空機と可能な限り一致する動的モデルを同定します。動的モデルは、システムの入力と出力の数学的関係を微分または差分方程式としてカプセル化します。動的モデルの例は、伝達関数と状態空間モデルです。System Identification Toolbox では、動的モデルは idtf
(伝達関数の場合)、idpoly
(AR、ARMA モデルの場合)、idss
(状態空間モデルの場合) などのオブジェクトによってカプセル化されます。動的モデルは、時間領域または周波数領域の測定データに対して tfest
や ssest
のような推定コマンドを実行することによって作成できます。
この例ではまず、etfe
コマンドを使用した経験的な伝達関数の推定により、時間領域の測定データを周波数応答データに変換します。次に、推定された FRF を使用して、航空機の振動応答の状態空間モデルを同定します。動的モデルの同定に時間領域データを直接使用することも可能です。ただし、FRF 形式のデータを使用すると、大きいデータセットを少数のサンプルに圧縮できると同時に、関連する周波数範囲に合わせて推定動作を調整することも簡単になります。FRF は idfrd
オブジェクトによってカプセル化されます。
データ実験ごとに、2 出力/1 入力の周波数応答関数 (FRF) を 1 つ推定します。ウィンドウ処理は使用しません。応答の計算に 24,000 個の周波数点を使用します。
G = cell(1, 10); N = 24000; for k = 1:10 % Convert time-domain data into a FRF using ETFE command Data_k = getexp(Data, k); G{k} = etfe(Data_k, [], N); % G{k} is an @idfrd object end
すべての FRF を単一の 20 出力/1 入力の FRF に連結します。
G = cat(1, G{:}); % concatenate outputs of all estimated FRFs G.OutputName = 'y'; % name outputs 'y1', 'y2', ..., 'y20' G.InterSample = 'bl';
推定された周波数応答の感覚をつかむために、出力 1 と 15 (任意に選択) の振幅をプロットします。対象の周波数範囲 (4 ~ 35 Hz) を拡大します。
opt = bodeoptions; % plot options opt.FreqUnits = 'Hz'; % show frequency in Hz opt.PhaseVisible = 'off'; % do not show phase OutputNum = [1 15]; % pick outputs 1 and 15 for plotting clf bodeplot(G(OutputNum, :), opt) % plot frequency response xlim([4 35]) grid on
FRF は少なくとも 9 つの共振周波数を示しています。解析のため、ここでは、航空機の最も重大な柔軟な曲げモードが存在する 6 ~ 35 Hz の周波数スパンに注目します。そのため、FRF をこの周波数領域まで減らします。
f = G.Frequency/2/pi; % extract frequency vector in Hz (G stores frequency in rad/s) Gs = fselect(G, f>6 & f<=32) % "fselect" selects the FRF in the requested range (6.5 - 35 Hz)
Gs = IDFRD model. Contains Frequency Response Data for 20 output(s) and 1 input(s). Response data is available at 624 frequency points, ranging from 37.96 rad/s to 201.1 rad/s. Sample time: 0.0005 seconds Output channels: '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)' Input channels: 'u1' Status: Estimated model
したがって、Gs
には、20 箇所の測定位置の周波数応答測定値が含まれます。次に、Gs
に適合する状態空間モデルを同定します。部分空間推定器 n4sid
は、非反復の推定を素早く行います。状態空間モデルの構造は、次のように構成します。
24 次の連続時間モデルを推定します。この次数は、さまざまな次数を試して FRF に対するモデルの適合結果をチェックすることによって見つかりました。
モデルには直達項が含まれています (D 行列はゼロ以外)。これは、解析を 35 Hz に制限しているものの、翼の帯域幅のほうがはるかに大きいからです (応答は 35 Hz で有意)。
計算時間を短縮するために、パラメーター共分散の計算を抑制します。
FRF の振幅は、周波数範囲内で大きく変化します。低い振幅をより高い振幅と同じように重要視するために、応答の平方根と逆に変化するカスタムの重み付けを適用します。チャネルが 20 個あるため、重みベクトルはそれぞれから取得された重みの平均として取得されます (20 チャネルの FRF は形状が類似しているため、これで問題ありません)。
n4sidOptions
を使用して、n4sid
の推定オプションを設定します。
FRF = squeeze(Gs.ResponseData); Weighting = mean(1./sqrt(abs(FRF))).'; n4Opt = n4sidOptions('EstimateCovariance',false,... 'WeightingFilter',Weighting,... 'OutputWeight',eye(20)); sys1 = n4sid(Gs,24,'Ts',0,'Feedthrough',true,n4Opt); Fit = sys1.Report.Fit.FitPercent'
Fit = 1×20
54.2689 55.3738 84.0188 85.7017 81.9259 80.2293 55.7448 40.6773 58.6429 76.2997 83.5370 84.6555 85.9240 84.8780 76.8937 81.0399 74.8321 79.4069 65.0803 76.3328
"Fit" の数値は、データ (Gs
) とモデル (sys1
) の周波数応答の間の適合の割合を、正規化平方根平均二乗誤差 (NRMSE) 適合度の尺度を使用して示しています。最低と最良の適合が次にプロットされます。
[~,Imin] = min(Fit); [~,Imax] = max(Fit); clf bodeplot(Gs([Imin, Imax],:), sys1([Imin, Imax],:), opt); xlim([6 32]) title('Worst (top) and best (bottom) fits between data and model') grid on legend('Data', 'Model')
モデル sys1
で達成された適合は、モデルのパラメーターを非線形最小二乗法で反復計算することによって、さらに改善できます。これは、ssest
コマンドを使用して達成できます。ssestOptions
コマンドを使用して、ssest
の推定オプションを設定します。今回は、パラメーター共分散も推定します。
ssOpt = ssestOptions('EstimateCovariance',true,... 'WeightingFilter',n4Opt.WeightingFilter,... 'SearchMethod','lm',... % use Levenberg-Marquardt search method 'Display','on',... 'OutputWeight',eye(20)); sys2 = ssest(Gs, sys1, ssOpt); % estimate state-space model (this takes several minutes) Fit = sys2.Report.Fit.FitPercent'
Fit = 1×20
88.8715 88.6587 88.9388 89.8630 87.8597 88.1235 79.7016 82.5930 74.2458 82.6593 90.2663 88.7739 89.9033 89.0571 88.8040 88.2251 86.9307 87.5250 83.2249 82.8676
前回と同様に、最悪と最良の適合をプロットします。また、モデルの周波数応答の 1-sd 信頼領域も可視化します。
[~, Imin] = min(Fit); [~, Imax] = max(Fit); clf h = bodeplot(Gs([Imin, Imax],:), sys2([Imin, Imax],:), opt); showConfidence(h, 1) xlim([6.5 35]) title('Worst (top) and best (bottom) fits between data and refined model') grid on legend('Data', 'Model (sys2)')
調整済みの状態空間モデル sys2
は、7 ~ 20 Hz の領域で FRF に非常によく適合します。共振が最大の位置の周辺で、不確かさ範囲が狭くなっています。ここでは 24 次のモデルを推定しましたが、これは、システム sys2
内に最大で 12 の振動モードがあるということです。modalfit
コマンドを使用して、これらのモードの固有振動数を Hz 単位で取得します。
f = Gs.Frequency/2/pi; % data frequencies (Hz) fn = modalfit(sys2, f, 12); % natural frequencies (Hz) disp(fn)
7.7720 7.7954 9.3176 9.4052 9.4919 15.3462 19.3281 23.0219 27.4035 28.6873 31.7036 64.2624
fn
の値は、2 つの周波数が 7.8 Hz に非常に近く、3 つが 9.4 Hz に近いことを示しています。これらの周波数に近い周波数応答を調べると、出力ごとにピークの位置が少しずれることがわかります。実験プロセスの制御を強化し、これらの周波数を中心とする狭い範囲に制限した入力帯域幅で時間領域同定を直接実行するか、これらの周波数に近い周波数応答に対して単一の振動モードを適合させることによって、こうしたばらつきを除去できる場合があります。こうした代替方法については、この例では扱っていません。
モーダル パラメーターの計算
モデル sys2
を使用してモーダル パラメーターを抽出できるようになりました。FRF を調べると、おおむね [5 7 10 15 17 23 27 30] Hz の周波数に近い 10 個前後のモーダル周波数が確認できます。この評価をより具体的にするために、modalsd
コマンドを使用して、基となるモデルの次数が変更されたときにモーダル パラメーターの安定性をチェックすることができます。この操作には長時間かかります。そのため、結果のプロットはイメージとして直接挿入されます。図を再現するには、以下のコメント ブロック内のコードを実行します。
FRF = permute(Gs.ResponseData,[3 1 2]); f = Gs.Frequency/2/pi; %{ figure pf = modalsd(FRF,f,fs,'MaxModes',20,'FitMethod','lsrf'); %}
プロットと pf
の値を調べると、実際の固有振動数の調整されたリストが提示されます。
Freq = [7.8 9.4 15.3 19.3 23.0 27.3 29.2 31.7];
ここでは、モデル sys2
から最も主要なモードを選択するための指針として、Freq
ベクトルの値を使用します。これを実行するには modalfit
コマンドを使用します。
[fn,dr,ms] = modalfit(sys2,f,length(Freq),'PhysFreq',Freq);
fn
は固有振動数 (Hz 単位)、dr
は対応する減衰係数、ms
は列ベクトル (固有振動数ごとに 1 列ずつ) として正規化された残差です。これらのモーダル パラメーターの抽出プロセスでは、安定していて減衰率の低いモデルの極のみが使用されます。ms
列には、正の虚数部をもつ極のみに対応するデータが含まれます。
モード形状の可視化
さまざまな曲げモードを可視化するには、上記で抽出したモーダル パラメーターが必要です。さらに、測定点の位置に関する情報も必要です。これらの位置 (x-y 値) は、行列 AccePos
の加速度計ごとに記録されます。
AccelPos = [... % see figure 2 16.63 18.48; % nearest right of center 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; % farthest right -16.63 18.48; % nearest left of center -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]; % farthest left
モード形状が含まれている行列 ms
では、各列が 1 つの周波数の形状に対応しています。モード形状の振幅をセンサー座標の上に重ね、モードの固有振動数における振幅を変えることにより、モードをアニメーション化します。アニメーションは、減衰のない曲げを示します。例として、15.3 Hz におけるモードについて考えます。
AnimateOneMode(3, fn, ms, AccelPos);
まとめ
この例は、モーダル パラメーター推定に対するパラメトリック モデル同定ベースのアプローチを示しています。24 次の状態空間モデルを使用して、6 ~ 32 Hz の周波数範囲で 8 個の安定した振動モードを抽出しました。モーダル情報を加速度計の位置の情報と結び付けて、さまざまな曲げモードを可視化しました。これらのモードの一部を以下の図に示します。
参考文献
[1] Gupta, Abhineet, Peter J. Seiler, and Brian P. Danowsky. "Ground Vibration Tests on a Flexible Flying Wing Aircraft." AIAA Atmospheric Flight Mechanics Conference, AIAA SciTech Forum. (AIAA 2016-1753).
function AnimateOneMode(ModeNum, fn, ModeShapes, AccelPos) % Animate one mode. % ModeNum: Index of the mode. % Reorder the sensor locations for plotting so that a continuous, % non-intersecting curve is traced around the body of the aircraft. PlotOrder = [19:-2:11, 1:2:9, 10:-2:2, 12:2:20, 19]; Fwd = PlotOrder(1:10); Aft = PlotOrder(20:-1:11); x = AccelPos(PlotOrder,1); y = AccelPos(PlotOrder,2); xAft = AccelPos(Aft,1); yAft = AccelPos(Aft,2); xFwd = AccelPos(Fwd,1); yFwd = AccelPos(Fwd,2); wn = fn(ModeNum)*2*pi; % Mode frequency in rad/sec T = 1/fn(ModeNum); % Period of modal oscillation Np = 2.25; % Number of periods to simulate tmax = Np*T; % Simulate Np periods Nt = 100; % Number of time steps for animation t = linspace(0,tmax,Nt); ThisMode = ModeShapes(:,ModeNum)/max(abs(ModeShapes(:,ModeNum))); % normalized for plotting z0 = ThisMode(PlotOrder); zFwd = ThisMode(Fwd); zAft = ThisMode(Aft); clf col1 = [1 1 1]*0.5; xx = reshape([[xAft, xFwd]'; NaN(2,10)],[2 20]); yy = reshape([[yAft, yFwd]'; NaN(2,10)],[2 20]); plot3(x,y,0*z0,'-', x,y,0*z0,'.', xx(:), yy(:), zeros(40,1),'--',... 'Color',col1,'LineWidth',1,'MarkerSize',10,'MarkerEdgeColor',col1); xlabel('x') ylabel('y') zlabel('Amplitude') ht = max(abs(z0)); axis([-100 100 10 40 -ht ht]) view(5,55) title(sprintf('Mode %d. FN = %s Hz', ModeNum, num2str(fn(ModeNum),4))); % Animate by updating z-coordinates. hold on col2 = [0.87 0.5 0]; h1 = plot3(x,y,0*z0,'-', x,y,0*z0,'.', xx(:), yy(:), zeros(40,1),'--',... 'Color',col2,'LineWidth',1,'MarkerSize',10,'MarkerEdgeColor',col2); h2 = fill3(x,y,0*z0,col2,'FaceAlpha',0.6); hold off for k = 1:Nt Rot1 = cos(wn*t(k)); Rot2 = sin(wn*t(k)); z_t = real(z0)*Rot1 - imag(z0)*Rot2; zAft_t = real(zAft)*Rot1 - imag(zAft)*Rot2; zFwd_t = real(zFwd)*Rot1 - imag(zFwd)*Rot2; zz = reshape([[zAft_t, zFwd_t]'; NaN(2,10)],[2 20]); set(h1(1),'ZData',z_t) set(h1(2),'ZData',z_t) set(h1(3),'ZData',zz(:)) h2.Vertices(:,3) = z_t; pause(0.1) end end