ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

可撓全翼機のモード解析

この例では、可撓翼機の曲げモードの計算を示します。翼の振動応答はその翼長に沿って複数か所で収集されます。データはシステムの動的モデルの識別に使用されます。モーダル パラメーターは識別されたモデルから抽出されます。モーダル パラメーターのデータは、翼のさまざまな曲げモードを可視化するためにセンサーの位置情報と組み合わされます。この例には、Signal Processing Toolbox™ が必要です。

可撓翼機

この例では、ミネソタ大学の Uninhabited Aerial Vehicle Laboratories で組み立てられた小型の可撓全翼機から収集したデータを調査します [1]。航空機の形状は以下のとおりです。

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

実験の設定

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

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

加振器の入力コマンドは、Asin(ω(t)t) という形式の定振幅チャープ入力として指定します。チャープ周波数は時間と共に線形に変化します。つまり、ω(t)=c0+c1t となります。入力信号でカバーされる周波数範囲は 3 ~ 35 Hz です。データは 2 個の加速度計 (1 つの x 位置の前方と後方のエッジにある加速度計) で同時に収集されます。つまり、10 回の実験を行って、20 個全部の加速度計の応答が収集されます。加速度計と力変換器の測定値はすべて、2000 Hz でサンプリングされます。

データの準備

データは 2 出力/1 入力の 10 セットの信号で表現され、それぞれに 600K のサンプルが含まれます。MAT ファイル FlexibleWingData.mat からデータを読み込みます。

load FlexibleWingData.mat MeasuredData

変数 MeasuredData は、フィールド Exp1Exp2、...、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 つずつ推定します。ウィンドウ処理は使用しません。応答の計算に 24000 個の周波数点を使用します。

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 1 つにまとめます。

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 は、非反復の推定をすばやく行います。状態空間モデルの構造は、次のように構成します。

  1. 18 次の連続時間モデルを推定します。この次数は、さまざまな次数を試して FRF に対するモデルの適合結果をチェックすることによって見つかりました。

  2. モデルには直達項が含まれています (D 行列はゼロ以外です)。これは、解析を 35 Hz に制限しているものの、翼の帯域幅の方がはるかに大きいからです (応答は 35 Hz で有意)。

  3. 計算時間を短縮するために、パラメーター共分散の計算を抑制します。

  4. FRF の振幅は、周波数範囲内で大きく変化します。低い振幅をより高い振幅と同じように重要視するために、11 番目の応答の平方根と逆に変化するカスタムの重み付けを適用します。11 番目の出力の選択はどちらかというと任意になりますが、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

   57.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

"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

   89.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],:), 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.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 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