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

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

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

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

Цель эксперимента состоит в том, чтобы собрать виброотдачу самолета в различных местах в ответ на внешнее возбуждение. Самолет подвешивается к деревянной системе координат с помощью одной пружины в центре тяжести. Пружина достаточно гибкая, чтобы естественная частота пружинно-массовых колебаний не мешала основным частотам самолета. Вход сила прикладывается через электродинамический вибросит Unholtz-Dickie Модели 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 реакций акселерометра. Измерения акселерометра и датчика силы производятся при частоте дискретизации 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 объекты. The 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

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

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 объекты.

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

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

Цифры «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)

    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 для извлечения модальных параметров. Проверка 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 нормированные невязки как векторы-столбцы (по одному столбцу на каждую естественную частоту). В процессе извлечения этих модальных параметров используются только стабильные, недостаточно демпфированные полюса модели. The 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] Gupta, Abhineet, Peter J. Seiler, and Brian P. Danowsky. «Наземные вибрационные испытания на гибком летающем самолете». Конференция 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