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

Этот пример показывает, как оценить функции частотной характеристики (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. Можно оценить качество идентифицированных моделей с помощью команд resid и compare. Если вы подтвердили качество модели, можно использовать их для извлечения модальных параметров. Это показывают кратко с помощью средства оценки пространства состояний.

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.