Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

柔軟全翼機のモード解析

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

柔軟翼機

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

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

実験の設定

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

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

シェーカーの入力コマンドは、A sin(ω(t)t) という形式の定振幅チャープ入力として指定します。チャープ周波数は時間とともに線形に変化し、ω(t)=c0+c1t となります。入力信号でカバーされる周波数範囲は 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 は、フィールド 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 (状態空間モデルの場合) などのオブジェクトによってカプセル化されます。動的モデルは、時間領域または周波数領域の測定データに対して tfestssest のような推定コマンドを実行することによって作成できます。

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

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

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

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

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