Естественная частота и демпфирование для моделируемого луча

Возбуждение молота с одним входом/одним выходом

Серия ударов молотка возбуждает 3DOF систему, и датчики регистрируют результирующие смещения. Система пропорционально демпфирована, так что демпфирующая матрица является линейной комбинацией матриц массы и жесткости.

Импортируйте данные для двух наборов измерений, включая сигналы возбуждения, сигналы отклика, временные сигналы и основную истину функции частотной характеристики. Первый набор ответных сигналов, Y1, измеряет перемещение первой массы, а второй, Y2, измеряет вторую массу. Каждый сигнал возбуждения состоит из десяти конкатенированных молотковых влияний, и каждый ответный сигнал содержит соответствующее перемещение. Длительность каждого сигнала влияния составляет 2,53 секунды. Аддитивный шум присутствует в сигналах возбуждения и отклика. Визуализируйте первый канал возбуждения и отклика первого измерения.

[t,fs,X1,X2,Y1,Y2,f0,H0] = helperImportModalData();
X0 = X1(:,1);
Y0 = Y1(:,1);
helperPlotModalAnalysisExample([t' X0 Y0]);

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

winLen = 2.5275*fs; % window length in samples
modalfrf(X0,Y0,fs,winLen,'Sensor','dis')

FRF, рассчитанная с использованием 'H1' по умолчанию устройство оценки содержит три видных peaks в полосах измеренной частоты, соответствующих трем гибким режимам вибрации. Когерентность близка к единице около этих пиков и низка в антирезонансных областях, где отношение сигнал/шум измерения отклика невелико. Когерентность около единицы указывает на высококачественную оценку. The 'H1' оценка оптимальна там, где шум существует только при выходе измерении, тогда как 'H2' оценка оптимальна, когда существует аддитивный шум только на входе [2]. Вычислите 'H1' и 'H2' оценки для этого ФРС.

[FRF1,f1] = modalfrf(X0,Y0,fs,winLen,'Sensor','dis');  % Calculate FRF (H1)
[FRF2,f2] = modalfrf(X0,Y0,fs,winLen,'Sensor','dis','Estimator','H2');

Когда существует значительный шум измерения или возбуждение плохое, параметрические методы могут предложить дополнительные опции для точного извлечения FRF из данных. Метод 'subspace' сначала подбирает модель пространства состояний к данным [3], а затем вычисляет его функцию частотной характеристики. Порядок модели пространства состояний (равный количеству полюсов) и наличие или отсутствие прохода могут быть заданы, чтобы сконфигурировать оценку пространства состояний.

[FRF3,f3] = modalfrf(X0,Y0,fs,winLen,'Sensor','dis','Estimator','subspace','Feedthrough',true);

Здесь FRF3 оценивается путем подбора модели пространства состояний, содержащей передаточный член и оптимальный порядок в области значений 1:10. Сравнение предполагаемых FRF с использованием 'H1', 'H2' и 'subspace' методы к теоретической FRF.

helperPlotModalAnalysisExample(f1,FRF1,f2,FRF2,f3,FRF3,f0,H0);

Оценщики выполняют сравнительно близкий peaks отклика, в то время как 'H2' оценщик переоценивает ответ на антирезонансах. На когерентность не влияет выбор оценщика.

Затем оцените естественную частоту каждого режима, используя алгоритм пикового захвата. Алгоритм пиковой выборки является простой и быстрой процедурой для идентификации peaks в FRF. Это локальный метод, поскольку каждая оценка генерируется из одной функции частотной характеристики. Это также метод с одной степенью свободы (SDOF), поскольку пик для каждого режима рассматривается независимо. В результате для каждого FRF генерируется набор модальных параметров. На основе предыдущего графика задайте частотную область значений от 200 до 1600 Гц, который содержит три пика.

fn = modalfit(FRF1,f1,fs,3,'FitMethod','PP','FreqRange',[200 1600])

fn =

   1.0e+03 *

    0.3727
    0.8525
    1.3707

Естественные частоты составляют приблизительно 373, 853 и 1371 Гц. Постройте график восстановленной FRF и сравните его с измеренными данными с помощью modalfit. FRF восстанавливается с помощью модальных параметров, оцененных из матрицы функции частотной характеристики, FRF1. Функции modalfit снова без выходных аргументов для создания графика, содержащего восстановленный FRF.

modalfit(FRF1,f1,fs,3,'FitMethod','PP','FreqRange',[200 1600])

Восстановленная FRF согласуется с измеренной FRF первых каналов возбуждения и отклика. В следующем разделе рассматриваются два дополнительных местоположения возбуждения.

Вращающееся Возбуждение Молота

Вычислите и постройте графики FRF для характеристик всех трех датчиков с помощью 'H1' по умолчанию оценщик. Задайте тип измерения следующим 'rovinginput' так как у нас есть вращающееся молотковое возбуждение.

modalfrf(X1,Y1,fs,winLen,'Sensor','dis','Measurement','rovinginput')

В предыдущем разделе один набор модальных параметров был вычислен из одного FRF. Теперь оцените модальные параметры с помощью алгоритма комплексной экпоненты методом наименьших квадратов (LSCE). Алгоритмы LSCE и LSRF генерируют один набор модальных параметров путем анализа нескольких ответных сигналов одновременно. Это глобальные методы с несколькими степенями свободы (MDOF), поскольку параметры для всех режимов оцениваются одновременно из нескольких функций частотной характеристики.

Алгоритм LSCE генерирует вычислительные режимы, которые физически не присутствуют в структуре. Используйте диаграмму стабилизации для идентификации физических режимов путем исследования устойчивости полюсов с увеличением количества режимов. Естественные частоты и коэффициенты затухания физических режимов, как правило, остаются на том же месте или являются 'стабильными'. Создайте диаграмму стабилизации и выведите естественные частоты тех полюсов, которые стабильны по частоте.

[FRF,f] = modalfrf(X1,Y1,fs,winLen,'Sensor','dis','Measurement','rovinginput');
fn = modalsd(FRF,f,fs,'MaxModes',20, 'FitMethod', 'lsce'); % Identify physical modes

По умолчанию полюсы классифицируются как стабильные по частоте, если естественная частота для полюсов изменяется менее чем на один процент. Стойкие по частоте полюсы дополнительно классифицируются как стабильные по демпфированию при изменении коэффициента затухания менее чем на пять процентов. Оба критерия могут быть скорректированы под различные значения. Исходя из расположения устойчивых полюсов, выбирайте естественные частоты 373, 852,5 и 1371 Гц. Эти частоты содержатся в выходе fn от modalsd, наряду с естественными частотами других частотно-устойчивых полюсов. Более высокий порядок модели, чем количество физически присутствующих режимов, обычно необходим, чтобы создать хорошие оценки модальных параметров с использованием алгоритма LSCE. В этом случае порядок модели из четырех режимов указывает на три стабильных полюса. Интересующие частоты происходят в первых трех столбцах в 4-й строке fn.

physFreq = fn(4,[1 2 3]);

Оценка естественных частот и демпфирования и графика восстановленных и измеренных FRF. Задайте четыре режима и физические частоты, определенные из схемы устойчивости, 'PhysFreq'. modalfit возвращает модальные параметры только для заданных режимов.

modalfit(FRF,f,fs,4,'PhysFreq',physFreq)

[fn1,dr1] = modalfit(FRF,f,fs,4,'PhysFreq',physFreq)

fn1 =

   1.0e+03 *

    0.3727
    0.8525
    1.3706


dr1 =

    0.0008
    0.0018
    0.0029

Затем вычислите FRF и постройте диаграмму стабилизации для второго набора влияний молотка с датчиком в другом месте. Измените критерий устойчивости на 0,1 процента для частоты и 2,5 процента для демпфирования.

[FRF,f] = modalfrf(X2,Y2,fs,winLen,'Sensor','dis','Measurement','rovinginput');
fn = modalsd(FRF,f,fs,'MaxModes',20,'SCriteria',[0.001 0.025]);

При более строгих критериях большинство полюсов классифицируются как не стабильные по частоте. Полюсы, которые стабильны по частоте и демпфированию, тесно совпадают со средним FRF, что позволяет предположить, что они присутствуют в измеренных данных.

physFreq = fn(4,[1 2 3]);

Извлеките модальные параметры для этого набора измерений и сравните с модальными параметрами для первого набора измерений. Задайте индексы FRF приводной точки, соответствующие месту, где совпадают измерения возбуждения и отклика. Естественные частоты согласуются в пределах части процента, и коэффициенты затухания соглашаются в пределах менее четырех процентов, указывая, что модальные параметры согласованы от измерения до измерения.

[fn2,dr2] = modalfit(FRF,f,fs,4,'PhysFreq',physFreq,'DriveIndex',[1 2])

fn2 =

   1.0e+03 *

    0.3727
    0.8525
    1.3705


dr2 =

    0.0008
    0.0018
    0.0029

Tdiff2 = table((fn1-fn2)./fn1,(dr1-dr2)./dr1,'VariableNames',{'diffFrequency','diffDamping'})

Tdiff2 =

  3x2 table

    diffFrequency    diffDamping
    _____________    ___________

      2.9972e-06      -0.031648 
     -5.9335e-06     -0.0099076 
       1.965e-05      0.0001186 

Параметрический метод для оценки модального параметра может предоставить полезную альтернативу пиковой выборке, и метод LSCE, когда существует шум измерения в FRF или FRF показывает высокую модальную плотность. Подход рациональной функции методом наименьших квадратов (LSRF) подходит для общей передаточной функции знаменателя в FRF с мультивходами, несколькими выходами и, таким образом, получает одну, глобальную оценку модальных параметров [4]. Процедура использования подхода LSRF аналогична процедуре для LSCE. Можно использовать диаграмму стабилизации, чтобы идентифицировать стационарные режимы и извлечь модальные параметры, соответствующие идентифицированным физическим частотам.

[FRF,f] = modalfrf(X1,Y1,fs,winLen,'Sensor','dis','Measurement','rovinginput');
fn = modalsd(FRF,f,fs,'MaxModes',20, 'FitMethod', 'lsrf'); % Identify physical modes using lsfr
physFreq = fn(4,[1 2 3]);
[fn3,dr3] = modalfit(FRF,f,fs,4,'PhysFreq',physFreq,'DriveIndex',[1 2],'FitMethod','lsrf')

fn3 =

  372.6832
  372.9275
  852.4986


dr3 =

    0.0008
    0.0003
    0.0018

Tdiff3 = table((fn1-fn3)./fn1,(dr1-dr3)./dr1,'VariableNames',{'diffFrequency','diffDamping'})

Tdiff3 =

  3x2 table

    diffFrequency    diffDamping
    _____________    ___________

     -7.8599e-06      0.007982  
         0.56255       0.83086  
         0.37799       0.37626  

Окончательное примечание о параметрических методах: метод оценки FRF ('subspace') и метод оценки модального параметра ('lsrf') аналогичны тем, которые используются в System Identification Toolbox для подбора кривой динамических моделей к сигналам временной области или функциям частотной характеристики. Если у вас есть доступный этот тулбокс, можно идентифицировать модели, соответствующие вашим данным, используя такие команды, как tfest и ssest. Оценить качество идентифицированных моделей можно используя compare и resid команды. После проверки качества модели можно использовать их для извлечения модальных параметров. Это кратко показано с использованием оценки пространства состояний.

Ts = 1/fs; % sample time
% Create a data object to be used for model estimation.
EstimationData = iddata(Y0(1:1000), X0(1:1000), 1/fs);
% Create a data object to be used for model validation
ValidationData = iddata(Y0(1001:2000), X0(1001:2000), 1/fs);

Идентифицируйте модель пространства состояний 6-го порядка в непрерывном времени, содержащую исходный термин.

sys = ssest(EstimationData, 6, 'Feedthrough', true)

sys =
  Continuous-time identified state-space model:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
            x1       x2       x3       x4       x5       x6
   x1     4.05    -1765    149.8    -1880   -49.64     -358
   x2     1764  -0.3332     2197   -232.5   -438.3   -128.4
   x3   -152.4    -2198     2.85     4715    255.9    547.5
   x4     1879    228.2    -4713    -15.9    -1216   -28.79
   x5    59.42    440.9   -275.5     1217    35.05    -8508
   x6    363.7    120.2   -545.4   -44.02     8508   -92.45
 
  B = 
            u1
   x1  -0.1513
   x2   -1.911
   x3    4.439
   x4   -3.118
   x5  -0.9416
   x6   -8.039
 
  C = 
               x1          x2          x3          x4          x5          x6
   y1   3.135e-05   2.511e-06   8.634e-06  -1.416e-05   2.218e-06  -6.271e-06
 
  D = 
              u1
   y1  7.564e-09
 
  K = 
               y1
   x1   3.513e+07
   x2  -3.244e+06
   x3  -3.598e+07
   x4  -1.059e+07
   x5   1.724e+08
   x6   7.521e+06
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: yes
   Disturbance component: estimate
   Number of free coefficients: 55
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                    
Estimated using SSEST on time domain data "EstimationData".
Fit to estimation data: 99.3% (prediction focus)           
FPE: 1.235e-16, MSE: 1.189e-16                             

Оцените качество модели, проверив, насколько хорошо она подходит для данных валидации.

clf
compare(ValidationData, sys)  % plot shows good fit

Используйте модель sys для вычисления модальных параметров.

[fn4, dr4] = modalfit(sys, f, 3);

Векторы формы режима ветряного Блейда

Понимание динамического поведения блейдов ветряного двигателя важно для оптимизации эксплуатационной эффективности и прогнозирования отказа блейдов. В этом разделе анализируются данные экспериментального модального анализа для блейда ветряного двигателя и визуализируются формы режимов блейда. Молоток возбуждает блейд турбины в 20 точках, и ссылка акселерометр измеряет отклики в местоположении 18. На основании блейда установлен алюминиевый блок, и блейд возбуждается в закрылковом направлении, перпендикулярном плоской части блейда. FRF собирается для каждого местоположения. Данные FRF любезно предоставлены Лабораторией структурной динамики и акустических систем Массачусетского университета, Лоуэлл. Во-первых, визуализируйте пространственное расположение мест измерения.

Нагрузка и построение расчетов блейд ветряного двигателя для мест 18 и 20. Увеличьте изображение первых нескольких пиков.

[FRF,f,fs] = helperImportModalData();
helperPlotModalAnalysisExample(FRF,f,[18 20]);

Первые два режима появляются как peaks около 37 Гц и 111 Гц. Постройте диаграмму стабилизации для оценки собственных частот. Первые два значения, возвращенные для порядка модели 14, стабильны по частоте и коэффициенту затухания.

fn = modalsd(FRF,f,fs,'MaxModes',20);
physFreq = fn(14,[1 2]);

Затем извлеките формы режима для первых двух режимов с помощью modalfit. Ограничьте подгонку частотой, областью значений от 0 до 250 Гц, исходя из предыдущего графика.

[~,~,ms] = modalfit(FRF,f,fs,14,'PhysFreq',physFreq,'FreqRange',[0 250]);

Формы Mode определяют количественную амплитуду движения для каждого режима структуры в каждом месте. Чтобы оценить вектор формы режима, необходима одна строка или столбец матрицы функции частотной характеристики. На практике это означает, что необходимо либо возбуждение в каждом месте измерения конструкции (в этом случае - гребной молоток), либо что необходимо измерение отклика в каждом месте. Формы режима могут быть визуализированы путем исследования мнимой части FRF. Постройте диаграмму водопада мнимой части матрицы FRF для расположения на одной стороне блейда. Ограничьте область значений максимум 150 Гц для исследования первых двух режимов. Peaks графика представляют формы режима.

measlocs = [3 6 9 11 13 15 17 19 20]; % Measurement locations on blade edge
helperPlotModalAnalysisExample(FRF,f,measlocs,150);

Формы, обозначенные на графике контуром peaks, представляют собой первый и второй изгибающие моменты блейда. Затем постройте график величины векторов формы режима для тех же мест измерения.

helperPlotModalAnalysisExample(ms,measlocs);

В то время как амплитуды масштабируются по-разному (векторы формы режима масштабируются до модала единиц A), контуры формы режима сходятся по форме. Форма первого режима имеет большое перемещение совета и два узла, где амплитуда вибрации равна нулю. Второй режим также имеет большой ход совета и имеет три узла.

Сводные данные

В этом примере анализируются и сравниваются наборы данных моделируемого модального анализа для 3DOF системы, возбуждаемой ровинговым молотком. Он оценил естественную частоту и демпфирование с помощью диаграммы стабилизации и алгоритмов LSCE и LSRF. Модальные параметры были последовательными для двух наборов измерений. В отдельном случае использования формы блейда ветряного двигателя были визуализированы с помощью мнимой части матрицы FRF и векторов формы режима.

Признание

Спасибо доктору Питеру Авитабиле из Лаборатории структурной динамики и акустических систем Массачусетского университета в Лоуэлле за содействие в наборе экспериментальных данных о блейдах ветряного двигателя.

Ссылки

[1] Брандт, Андерс. Анализ шума и вибрации: анализ сигналов и экспериментальные процедуры. Chichester, UK: John Wiley and Sons, 2011.

[2] Vold, Hovard, John Crowley, and G. Thomas Rocklin. Новые способы оценки функций частотной характеристики. Звук и вибрация. Том 18, ноябрь 1984, стр. 34-38.

[3] Питер Ван Оверши и Барт Де Моор. «N4SID: Алгоритмы подпространства для идентификации комбинированных детерминированно-стохастических систем». Автоматика. Том 30, январь 1994 года, стр. 75-93.

[4] Оздемир, А. А., и С. Гумуссой. «Оценка передаточной функции в System Identification Toolbox через векторный подбор кривой». Материалы XX Всемирного конгресса Международной федерации автоматического управления. Тулуза, Франция, июль 2017 года.