Гибкий самолет крыла

В этом примере мы изучаем данные, собранные от маленького гибкого самолета Летающее Крыло, созданного в Необитаемых Воздушных Лабораториях Транспортного средства, Миннесотском университете [1]. Геометрию самолета показывают ниже.

Крыло самолета может подвергнуться большим деформациям при загрузке. Гибкие частоты режима ниже, чем те в общем самолете с более твердыми крыльями. Этот гибкий проект уменьшает материальные затраты, гибкость увеличений и область значений рейса самолета. Однако, если не управляемый, гибкие режимы могут привести к катастрофической аэроупругой нестабильности (дрожание). Разработка эффективных законов о надзоре для подавления этой нестабильности требует точного определения различных режимов изгиба крыла.

Экспериментальный Setup

Цель эксперимента состоит в том, чтобы собрать ответ вибрации самолета в различных местоположениях в ответ на внешнее возбуждение. Самолет приостановлен от деревянной системы координат с помощью одной пружины в ее центре тяжести. Пружина достаточно гибка так, чтобы собственная частота пружинно-массового колебания не вмешивалась в основные частоты самолета. Входная сила прикладывается с помощью Модели 20 Анхолц-Дики электродинамический шейкер около центра самолета.

Двадцать акселерометров PCB-353B16 помещаются вдоль размаха крыла, чтобы собрать ответ вибрации как показано на следующем рисунке.

Входная команда шейкера задана как постоянный амплитудный вход щебета формы $\mathit{A}\sin \left(\omega \left(\mathit{t}\right)\mathit{t}\right)$. Частота щебета варьируется линейно со временем, то есть, $\omega \left(\mathit{t}\right)={\mathit{c}}_0+{\mathit{c}}_1\mathit{t}$. Частотный диапазон, покрытый входным сигналом, составляет 3-35 Гц. Данные собраны двумя акселерометрами (начальные и конечные акселерометры ребра в одном x-местоположении) за один раз. Следовательно 10 экспериментов проводятся, чтобы собрать все 20 ответов акселерометра. Акселерометр и измерения преобразователя силы все производятся на уровне 2 000 Гц.

Подготовка данных

Данные представлены 10 наборами сигналов two-output/one-input, каждый содержащий 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. Каждое поле является структурой с полями 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™. Входной сигнал обработан как bandlimited.

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

Объедините объекты набора данных в один объект данных мультиэксперимента.

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 данных позволяет сжатие больших наборов данных в меньшее количество выборок, а также более легко настройте поведение оценки к соответствующим частотным диапазонам. FRFs инкапсулируются idfrd объекты.

Оцените two-output/one-input функцию частотной характеристики (FRF) для каждого эксперимента данных. Не используйте работу с окнами. Используйте 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

Конкатенация всего FRFs в сингл 20-output/one-input 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 Гц).

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 Гц, где самые критические гибкие режимы изгиба самолета лежат. Следовательно уменьшайте 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 обеспечивает быструю неитеративную оценку. Структура модели в пространстве состояний сконфигурирована можно следующим образом:

Мы настраиваем опции оценки для n4sid использование n4sidOptions.

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

"Подходящие" числа показывают подгонку процента между данными (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 команда. Мы настраиваем опции оценки для ssest использование ssestOptions команда. На этот раз ковариация параметра также оценивается.

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 соответствует FRFs вполне хорошо в области на 7-20 Гц. Границы неопределенности трудны вокруг большинства резонирующих местоположений. Мы оценили модель 24-го порядка, что означает, что существует самое большее 12 колебательных режимов в системе sys2. Используйте modalfit команда, чтобы выбрать собственные частоты в Гц для этих режимов.

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 предложите две частоты очень близко к 7,8 Гц и три близко к 9,4 Гц. Контроль частотных характеристик около этих частот показывает, что местоположение peaks переключает немного через выходные параметры. Эти несоответствия могут быть удалены лучшим управлением процессом эксперимента, выполнив прямую идентификацию временного интервала с входной полосой пропускания, ограниченной узким диапазоном, сосредоточенным на этих частотах, или соответствуя одному колебательному режиму к частотной характеристике вокруг этих частот. Эти альтернативы не исследуются в этом примере.

Модальный расчет параметра

Мы можем теперь использовать модель sys2 извлекать модальные параметры. Контроль FRFs указывает приблизительно на 10 модальных частот, примерно вокруг частот [5 7 10 15 17 23 27 30] Гц. Мы можем сделать эту оценку более конкретной при помощи 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];

Мы используем значения в Freq вектор как руководство для выбора самых доминирующих режимов из модели sys2. Это сделано с помощью modalfit команда.

[fn,dr,ms] = modalfit(sys2,f,length(Freq),'PhysFreq',Freq);

fn собственные частоты (в Гц), dr соответствующие коэффициенты демпфирования и ms нормированные остаточные значения как вектор-столбцы (один столбец для каждой собственной частоты). В процессе экстракции этих модальных параметров, только конюшня, под - используются ослабленные полюса модели. 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 где каждый столбец соответствует форме для одной частоты. Анимируйте режим путем наложения амплитуд формы режима на координаты датчика и варьирования амплитуд в собственной частоте режима. Анимация показывает изгиб без затухания. Как пример, рассмотрите режим на уровне 15,3 Гц.

AnimateOneMode(3, fn, ms, AccelPos); 

Заключения

Этот пример показывает параметрический основанный на идентификации подход модели к модальной оценке параметра. Используя модель в пространстве состояний 24-го порядка, 8 устойчивых колебательных режимов в области частоты были извлечены 6-32 Гц. Модальная информация была объединена со знанием положений акселерометра, чтобы визуализировать различные режимы изгиба. Некоторые из этих режимов показывают в рисунках ниже.

Ссылки

[1] Гупта, Abhineet, Питер Дж. Сейлер и Брайн П. Дановский. "Оснуйте тесты вибрации на гибком самолете Летающее Крыло". AIAA атмосферная конференция бортмехаников, форум AIAA SciTech. (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