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

В этом примере показано, как оценить функцию частотной характеристики (FRFs) и модальные параметры от экспериментальных данных. Первый раздел описывает симулированный эксперимент, который волнует три степени свободы (3DOF), система с последовательностью молотка влияет и записывает получившееся смещение. Функция частотной характеристики, собственные частоты, коэффициенты затухания и векторы формы режима оцениваются для трех режимов структуры. Второй раздел оценивает векторы формы режима от оценок функции частотной характеристики из эксперимента лопатки ветряного двигателя. Турбинная настройка измерения лопатки и получившиеся формы режима визуализируются. Этот пример требует System Identification Toolbox(TM).

Собственная частота и ослабляющий для симулированного луча

Возбуждение Молотка Single-Input/Single-Output

Серия забастовок молотка волнует 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 в измеренном диапазоне частот, соответствуя трем гибким режимам вибрации. Когерентность является близко к одной близости этим peaks, и низко в областях антирезонанса, где отношение сигнал-шум измерения ответа является низким. Когерентность близко к каждый указывает на высококачественную оценку. 'H1' оценка оптимальна, где шум существует только при выходном измерении, тогда как 'H2' средство оценки оптимально, когда существует аддитивный шум только на входе [2]. Вычислите 'H1' и 'H2' оценки для этого FRF.

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

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

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

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

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

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

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

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

   1.0e+03 *

    0.3727
    0.8525
    1.3707

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

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

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

Мобильное возбуждение молотка

Вычислите и постройте FRFs для ответов всех трех датчиков с помощью '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, и 1 371 Гц. Эти частоты содержатся в выходе fn из modalsd, наряду с собственными частотами других устойчивых частотой полюсов. Более высокий порядок модели, чем количество физически существующих режимов обычно необходим, чтобы произвести хорошие модальные оценки параметров с помощью алгоритма LSCE. В этом случае порядок модели четырех режимов указывает на три устойчивых полюса. Частоты интереса возникают в первых трех столбцах в 4-й строке fn.

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

Оцените собственные частоты и затухание и постройте восстановленный и измеренный FRFs. Задайте четыре режима и физические частоты, определенные из схемы устойчивости, '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

Затем вычислите FRFs и постройте диаграмму стабилизации для второго набора ударов молотка с датчиком в другом месте. Измените критерий стабильности в 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 ('подпространство') и модальный метод оценки параметра ('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 любезно обеспечиваются Структурной Динамикой и Системной Лабораторией Акустики в Массачусетском университете, Лоуэлл. Во-первых, визуализируйте пространственное расположение местоположений измерения.

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

[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]);

Формы режима определяют количество амплитуды движения для каждого режима структуры в каждом местоположении. Оценить вектор формы режима, одну строку или столбец матрицы функции частотной характеристики необходимо. На практике это означает, что любой, возбуждение необходимо в каждом местоположении измерения структуры (в этом случае, мобильный молоток), или что измерение ответа необходимо в каждом местоположении. Формы режима могут визуализироваться путем исследования мнимой части 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] Брандт, Андерс. Шум и анализ вибрации: анализ сигнала и экспериментальные процедуры. Чичестер, Великобритания: Джон Вайли и сыновья, 2011.

[2] Vold, Havard, Джон Кроули и Г. Томас Роклин. "Новые Способы Оценить Функции Частотной характеристики". Звук и Вибрация. Издание 18, ноябрь 1984, стр 34-38.

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

[4] Ozdemir, A. A. и С. Гумассой. "Оценка Передаточной функции System Identification Toolbox через Подбор кривой Вектора". Продолжения 20-го Мирового Конгресса Международной федерации Автоматического управления. Тулуза, Франция, июль 2017.