Модальный анализ летательного аппарата с гибким крылом

В этом примере используются:

В этом примере показано вычисление режимов изгиба летательного аппарата с гибким крылом. Виброустойчивость крыла собирается в нескольких точках вдоль его размаха. Данные используются для идентификации динамической модели системы. Модальные параметры извлекаются из идентифицированной модели. Данные модального параметра объединяются с информацией о положении датчика для визуализации различных режимов изгиба крыла. В этом примере требуется Toolbox™ обработки сигналов.

Летательный аппарат с гибким крылом

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

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

Экспериментальная настройка

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

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

Входная команда шейкера задана как вход чирп постоянной амплитуды формы $A \sin (λ (t)$ t). Частота чирпа изменяется линейно со временем, то $есть (t)_{}_{=}$ c0 + c1t. Диапазон частот, покрываемый входным сигналом, составляет 3-35 Гц. Данные собираются двумя акселерометрами (акселерометрами переднего и заднего краев в одном x-месте) одновременно. Таким образом, было проведено 10 экспериментов для сбора всех 20 откликов акселерометра. Все измерения акселерометра и датчика силы проводятся при частоте 2000 Гц.

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

Данные представлены 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. Каждое поле представляет собой структуру с полями 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 объекты являются стандартным способом упаковки данных временной области в 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

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

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

Идентификация модели

Мы хотим как можно ближе определить динамическую модель, частотная характеристика которой соответствует характеристике фактического самолета. Динамическая модель инкапсулирует математическое соотношение между входами и выходами системы в виде дифференциального или дифференциального уравнения. Примерами динамических моделей являются передаточные функции и модели состояния-пространства. В панели инструментов идентификации системы динамические модели инкапсулируются такими объектами, как idtf (для передаточных функций), idpoly (для моделей AR, ARMA) и idss (для моделей пространства состояний). Динамические модели могут быть созданы путем выполнения команд оценки, таких как tfest и ssest команды на измеренные данные во временной или частотной области.

Для этого примера мы сначала преобразуем измеренные данные временной области в данные частотной характеристики посредством эмпирической оценки передаточной функции с использованием etfe команда. Затем расчетная FRF используется для идентификации модели состояния-пространства вибрационной реакции летательного аппарата. Существует возможность непосредственного использования данных временной области для идентификации динамических моделей. Однако FRF-форма данных позволяет сжимать большие наборы данных в меньшее количество выборок, а также легче корректировать поведение оценки в соответствии с соответствующими частотными диапазонами. FRF инкапсулированы idfrd объекты.

Оценка функции частотной характеристики с двумя выходами/одним входом (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

Объедините все FRF в один 20-выходной/один-входной 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 обеспечивает быструю неидертеративную оценку. Структура модели state-space конфигурируется следующим образом:

Мы оцениваем модель непрерывного времени 18-го порядка. Заказ был найден после некоторых испытаний с различными заказами и проверки результирующего соответствия модели FRF.
Модель содержит член сквозной связи (матрица D ненулевая). Это потому, что мы ограничиваем наш анализ $\leq$ 35 Гц, в то время как ширина полосы крыла значительно больше, чем это (отклик значим при 35 Гц).
Для ускорения вычислений подавляем вычисление ковариации параметров.
Величина FRF значительно изменяется в диапазоне частот. Чтобы гарантировать, что низкие амплитуды получают равное ударение, как и более высокие амплитуды, мы применяем пользовательское взвешивание, которое изменяется обратно как квадратный корень 11-го отклика. Выбор 11-го выхода является несколько произвольным, но работает, поскольку все 20 FRF имеют схожие профили.

Мы настраиваем параметры оценки для 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

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

f = Gs.Frequency/2/pi;     % data frequencies (Hz)
fn = modalfit(sys2, f, 12); % natural frequencies (Hz)
disp(fn)

Значения в fn предполагают две частоты, очень близкие к 7,8 Гц, и три, близкие к 9,4 Гц. Проверка частотных откликов вблизи этих частот показывает, что пики местоположения немного смещаются между выходами. Эти расхождения могут быть устранены посредством лучшего управления экспериментальным процессом, выполнения прямой идентификации во временной области с входной полосой пропускания, ограниченной узким диапазоном, центрированным на этих частотах, или подгонки одиночного колебательного режима к частотной характеристике вокруг этих частот. В этом примере эти альтернативы не рассматриваются.

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

Теперь мы можем использовать модель sys2 для извлечения модальных параметров. Проверка FRF показывает около 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] Гупта, Абхинит, Питер Дж. Зайлер и Брайан П. Дановский. «Наземные вибрационные испытания на летательном аппарате с гибким крылом». Конференция 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

Документация