Коррекция матрицы демпфирования с использованием экспериментальных значений коэффициентов модального демпфирования Текст научной статьи по специальности «Физика»

УДК 534.1

Коррекция матрицы демпфирования с использованием экспериментальных значений коэффициентов модального демпфирования

© С.Н. Дмитриев1, Р.К. Хамидуллин2

1МГТУ им. Н.Э. Баумана, Москва, 105005, Россия 2ОАО «ВПК «НПО машиностроения», г. Реутов Московской обл., 143966, Россия

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

Ключевые слова: динамика, конечный элемент, демпфирование, модальный анализ.

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

Теоретическая часть. Колебания линейной системы с п степенями свободы [1, 2] описываются матричным дифференциальным уравнением следующего вида:

Ми + £и + К и = Р, (1)

где М — матрица массы; Б — матрица демпфирования; К — матрица жесткости; и — вектор узловых перемещений; Р — вектор внешних

сил. В расчетах обычно предполагается [2], что демпфирование является классическим, т. е. пропорциональным матрице жесткости К (случай внутреннего демпфирования [1]):

Б = аК, (2)

или матрице массы М (случай внешнего демпфирования [2]):

Б = РМ (3)

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

-Шф + Кф = 0, (4)

где X — собственное число; ф — собственный вектор.

Собственные числа \ и векторы ф(г) для разных номеров тонов колебаний г собираются соответственно в диагональную матрицу собственных чисел Л (5) и матрицу форм Ф (6), где они располагаются в порядке возрастания собственных чисел:

Л = Х2, ..., Хп); (5)

Ф = [ф(1), ф(2), ..., ф(и)]. (6)

Найденные собственные числа и векторы удовлетворяют соотношениям ортогональности с матрицами жесткости и массы (К- и М-орто-гональности) [1]:

ФТКФ = Л; (7)

ФТМФ = Е, (8)

где Е — единичная матрица.

Переход к координатам форм осуществляется согласно следующему соотношению:

и = Фя, (9)

где я = [я1, я2, ..., яИ]Т — вектор координат форм. После подстановки (9) и использования соотношений ортогональности (7) и (8), уравнение (1) приводится к виду

Я + Щ + Ля = 0. (10)

Здесь диагональная матрица коэффициентов модального демпфирования

Я = ФТБФ (11)

для случая (2), когда матрица демпфирования пропорциональна матрице жесткости:

Я = аЛ, (12)

а в случае (3), когда матрица демпфирования пропорциональна матрице масс,

Я = (13)

Известно [3], что гипотеза классического демпфирования противоречит опытным данным. Величина конструкционного демпфирования может зависеть не только от частоты, но и от амплитуды и формы колебаний. Несмотря на высказанные выше соображения, по экспериментальным данным обычно строится динамическая модель в виде системы несвязанных между собой уравнений (гипотеза Базиля [3]), а демпфирование предполагается пропорциональным модальной скорости. Матрица Я в этой модели является диагональной. Как указано в [3], гипотеза Базиля о диагональности матрицы демпфирования (в модальных координатах) в общем случае предполагает, что эта матрица представляет собой линейную комбинацию матриц М и К. Однако данное утверждение не вполне точно. Очевидно, что в общем случае Я не может быть представлена и виде суммы (12) и (13), так как п различных коэффициентов демпфирования нельзя задать, варьируя лишь два параметра: а и р.

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

Если бы опытным путем были найдены все п коэффициентов модального демпфирования, то восстановление матрицы Б можно было бы осуществить, проводя обратное преобразование уравнения (10) к виду (1). Подставим в (10) выражения для вектора координат форм (14) и для матрицы собственных чисел (15):

q = ФтМи (14)

и

Л = ФтКФ. (15)

В результате подстановки получим

ФТМ Ii + R ФМи + ФТК ФФТМи = Фт P.

(16)

Выражение (16) умножим слева на транспонированную обращенную матрицу форм Ф-т = МФ, учтем что ФТМ = Ф-1. Тогда

Mü + Diu + Ки = P,

где матрица демпфирования

D = М ФЯФТМ (17)

В частных случаях классического демпфирования матрица R, найденная по D (17), принимает вид (12) или (13). Покажем, что для случая (12), когда матрица демпфирования пропорциональна матрице жесткости, имеет место тождество

D = М ФRФTM = аК. (18)

Умножим (17) слева на ФТ, а справа — на Ф. В силу соотношений M- и K-ортогональности получим R = аЛ. Следовательно, выполняется соотношение (2), а значит, и (18). Аналогичным образом можно показать, что случай, когда матрица демпфирования пропорциональна матрице масс, также охватывается предложенной общей записью (17).

Теперь рассмотрим ситуацию, когда часть коэффициентов матрицы модального демпфирования найдена опытным путем.

Разделим матрицу коэффициентов модального демпфирования на две части — измеренные коэффициенты демпфирования R1 и коэффициенты демпфирования R2, которые предполагается задавать нормативно (допустим, пропорционально матрице жесткости). Формы колебаний также разделим на две части Ф1 и Ф2 в соответствии с измеренными и расчетными коэффициентами демпфирования. Представим R в виде суммы (19):

R = R1 + R2 =

" Rj 0" " 0 0 "

+

_ 0 0 _ _ 0 R2 _

(19)

Подставим (19) в формулу (18) для матрицы D:

D = М ФRФ1М = М Ф^Ф1 М + М ФТ?2ФТ М =

= М [Ф! Ф2 ]

" Rj 0" "ФТ"

_ 0 0 _ _ФТ _

М + М [Ф1 Ф2 ]

"0 0 " "ФТ"

_0 R2 _ _Ф2_

М =

= М Ф1 R^T М + М Ф2 R^T М = D1 + D2.

(20)

Первое слагаемое выражения (20) мы определяем на основе найденных опытным путем коэффициентов модального демпфирования R. При определении второго слагаемого предположим, что соответствующие ему коэффициенты модального демпфирования имеют вид (12) или (13), т. е. для высших тонов колебаний демпфирование является классическим, например, в случае (12):

D2 = aM Ф2 Л2 ФТ2М. (21)

Сложность заключается в том, что высшие тона колебаний для системы с достаточно большим числом степеней свободы не определяются ни расчетным путем, ни экспериментально. Иными словами, Л2 и Ф2 неизвестны. Найти D2 можно, подставив в выражение (20) формулу (18), тем самым получив выражение для скорректированной матрицы демпфирования:

D = aK = M Ф^ФТМ + M Ф2Я2 Ф]Ы = M Ф^Ф^М + D2. (22)

Откуда вторая составляющая матрицы демпфирования

D2 =aK - aMФ1Л1Ф?М. (23)

Собственные числа и векторы для низших тонов колебаний являются известными. Таким образом, скорректированная по результатам опыта матрица D вычисляется по формуле

D = aK + МФ1 (R1 - аЛ1) Ф^М. (24)

Отметим еще раз, что в этой формуле a — априорно заданный нормативный коэффициент демпфирования, а Ri — диагональная матрица модальных коэффициентов демпфирования, найденных опытным путем для части тонов колебаний.

Расчетная часть. Проверка расчетных формул, полученных в теоретической части статьи, была проведена на конечно-элементной модели балки. В экспериментальной части статьи для балки с такими геометрическими параметрами приведены результаты частотных испытаний. Балка имела поперечное сечение в форме швеллера № 10 и была вывешена через два рым-болта М10 на резиновых жгутах. Геометрические характеристики швеллера № 10 были приняты согласно ГОСТ 8240-89. Погонная плотность pF = 8,59 кг/м, момент инерции сечения при изгибе в плоскости меньшей жесткости Jx = 174 • 10-8 м4. Другие характеристики швеллера в расчете не использовались. Длина швеллера L = 1,53 м. Масса рым-болта М10, согласно ГОСТ 4751-73, — 0,012 кг. Жесткость резиновых амортизаторов определялась перед проведением испытаний путем сравнения их длины в свободном и нагру-

женном состояниях балки, разница составила 0,2 м. По этим данным жесткость к = 320 Н/м. Подвеска прикреплялась на расстоянии 0,285 м, а усилия штоков были приложены на расстояниях по 0,1 м от каждого края. Масса мембраны и штока силовозбудителя, согласно технической документации, составляет 0,36 кг.

Расчет проводился методом конечного элемента. Число узловых точек — 201, матрица жесткости соответствовала классической теории изгиба Бернулли—Эйлера, матрица масс — диагональная. Для поверки достаточности числа разбиений модель была протестирована на шар-нирно опертой по концам балке. Первые шесть собственных частот совпали с аналитическим решением с точностью до четырех значащих цифр. Рым-болты, штоки и подвеска учитывались как сосредоточенные массы и жесткости. Подвеска является достаточно податливой, такие условия закрепления имитируют свободную балку. В расчете первые две частоты получаются на два порядка меньше, чем третья, поэтому они рассматриваются как нулевые. Значения ненулевых собственных частот приведены в табл. 1.

Таблица 1

Частота колебаний балочной модели (расчет методом конечного элемента)

Показатель Номер тона

3 4 5 6

Собственная частота, Гц 102,43 288,21 558,00 944,26

Испытания также проводились в частотном диапазоне, позволяющем определить третью и последующие частоты.

Процедура формирования матрицы демпфирования (24), оформленная в виде программного модуля m_demp на языке фортран, проверялась путем обратного преобразования по формуле (11) с последующим сравнением матриц Я, заложенных в расчет и полученных после обратного преобразования. Относительная погрешность составила приблизительно 10-17. Разница в откликах на внезапно приложенную в центре балки силу в 2 Н, рассчитанных по конечно-элементной модели (1) и с разложением по собственным формам (10) с последующим переходом к перемещениям (9), была замечена только в шести-семи-значащих разрядах. В расчет были заложены коэффициенты демпфирования, найденные в экспериментальной части статьи.

Дополнительно в расчетной части была апробирована фортран-программа о1кНк, разработанная авторами для расчетного моделирования вибрационных испытаний при многоточечном возбуждении. Эта программа использует математический аппарат теории передаточных

функций [1, 2, 4], которые вычисляются по частотам, формам колебаний и коэффициентам демпфирования конечно-элементной модели (1). Зависимости для передаточных функций с учетом демпфирования приведены в [4]. Программа реализована в среде Compaq Visual Fortran. Version 6.6. Расчет проводился на синусоидальное возбуждение от двух симметрично расположенных силовозбудителей. В результате были определены действительная и мнимая части откликов в окрестности резонансных частот ненулевых тонов с первого по четвертый (два низших тона считались нулевыми и в опыте не определялись). В процессе расчета было обнаружено также, что для возбуждения первого и третьего ненулевых тонов возбудители должны совершать движение в фазе, а второго и четвертого — в противофазе. В дальнейшем это подтвердилось при модальных испытаниях. В качестве примера на рис. 1 приводятся расчетные графики мнимой (а) и действительной (б) частей амплитудно-частотной характеристики (АЧХ) в окрестности первой резонансной частоты для точки 3 (на расстоянии 1,5 см от левого конца балки), на которых значения даны в частотном диапазоне 102,0.. .102,8 Гц. По оси ординат указаны значения ускорений в точке 1. Расчетные АЧХ для этого и последующих тонов колебаний качественно совпадают и по значениям относительно близки к экспериментальным, приводимым ниже, хотя

м/с2

4,0 Е + 02 3,5Е + 02 3,0 Е + 02 2,5 Е + 02 2,0 Е + 02 1,5Е + 02 1,0Е + 02 0.5Е+01 0,0 Е + 00

а

2,0 Е + 02 1,5£+02 \,0Е + 02 0,5 Е + 02 0,0£ + 02 -0,5 Е + 02 -1,0£+02 -\,5Е + 02 -1,0£+02

102,0 Гц „ 102,8 Гц

о

Рис. 1. Расчетная амплитудно-частотная характеристика в окрестности первого ненулевого тона: а — ее мнимая часть; б — действительная

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

Экспериментальная часть. Модальные испытания проходили в виброзале ОАО «ВПК «НПО машиностроения». Возбуждение вибраций осуществлялось двумя электродинамическими силовозбудителями 20Ш20С № 645 и № 656 системы «Ргс^ега» с максимальной силой возбуждения по 20 кгс. Испытания проводились в диапазоне частот 50.. .1000 Гц. Амплитудное значение сил составляло 2 Н. Замеры параметров колебаний осуществлялись с помощью 16 датчиков типа 1СР модели 3097А2. Обработку результатов, построение полученных форм колебаний и АЧХ проводили с помощью аппаратуры LMS Scadas III. Внешний вид экспериментальной установки показан на рис. 2.

Рис. 2. Экспериментальная установка

При обработке результатов испытаний и определении параметров демпфирования объекта принимались следующие допущения:

• не учитывалась амплитудозависимость характеристик;

• измерительная ось каждого из вибродатчиков совпадает с направлением возбуждения вибрации;

• частота дискретизации достаточна для точного определения параметров колебаний;

• влиянием шумовых составляющих сигналов можно пренебречь.

В работе применялось два метода определения логарифмических

декрементов: резонансной кривой и свободных затухающих колебаний.

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

Ар

5 = п-

рез

где ррез — резонансная частота; Ар — ширина резонансного пика, на которой мнимая составляющая АЧХ уменьшается в 2 раза.

В качестве примера на рис. 3 показана АЧХ в районе первого из-гибного тона швеллера, необходимая для определения логарифмического декремента.

-15

104,0 104,1 104,2 104,4 104,4 /104,5 104,6 104,7 104,8 Гц

104,48

Рис. 3. Экспериментальная АЧХ в окрестности первого изгибного тона

швеллера

15,91 104,38 // 31,82 \ 104,57

III — 104,48 (-Г I I

—/

Метод определения характеристик демпфирования по затухающим колебаниям заключается в получении осциллограмм свободных колебаний исследуемой системы, по темпу убывания которых определяется величина декремента колебаний 5 [5]. На практике оценка логарифмического декремента осуществляется с помощью виброграммы свободных колебаний исследуемого тела по формуле

5 = - 1п

п

1 Г и Л

V иг+п у

где и, и + п — амплитуды колебаний в начале и в конце интервала, состоящего из п колебаний.

Результаты проведенных экспериментальных исследований (значения амплитуд, собственных частот и декрементов колебаний для первых четырех полученных изгибных тонов колебаний) приведены в табл. 2

Таблица 2

Экспериментальные значения собственных частот и коэффициентов демпфирования

Показатель Номер тона

3 4 5 6

Собственная частота (по АЧХ), Гц 104,47 287,52 549,25 866,06

Собственная частота (по фигурам Лиссажу), Гц 104,44 287,60 549,21 866,26

Максимальная перегрузка (по АЧХ) g 31,8 14,4 21,5 7,9

Декремент колебаний (по АЧХ в точке 1) х 10_3 5,7 9,5 5,2 12,4

Декремент колебаний (по экспоненциальному затуханию в точке 1) х10 6,0 9,5 5,5 12,2

Собственные частоты определялись двумя методами: построением АЧХ и уточнением частоты по фигурам Лиссажу (по сдвигу фазы). По определению, на резонансной частоте сдвиг фазы между возбуждающей силой и вынужденными колебаниями составляет 90° [6]. При построении АЧХ шаг дискретизации принимался равным 0,02 Гц.

При расчете декрементов по затухающим колебаниям в расчет принимались лишь первые 30 периодов колебаний после выключения задающей нагрузки. Расчетные значения частот, приведенные в табл. 1, близки к найденным экспериментально.

Использование опытных данных для коррекции матрицы демпфирования. Для пересчета логарифмического декремента колебаний в коэффициент а перед матрицей жесткости в (12) воспользуемся формулой [1]

5

а = — 2п2/

Значения коэффициентов, использовавшихся для коррекции матрицы демпфирования Б, приведены в табл. 3.

Таблица 3

Значения коэффициентов демпфирования, использованных в расчете

Показатель Номер тона

3 4 5 6

Принимаемое в расчете значение логарифмического декремента колебаний ><10 5,9 9,5 5,4 12,3

Коэффициент а 2,86 • 10_6 1,67 • 10_6 4,98 • 10_7 7,19 • 10_7

Для первых двух тонов (смещения как жесткого целого) коэффициенты демпфирования принимались равными нулю. Коэффициенты из табл. 3 были использованы для третьего—шестого тонов колебаний. Для тонов выше шестого а = 0,144 • 10-5 — среднее значение по результатам четырех опытов. Эти коэффициенты были использованы для получения матрицы демпфирования.

Заключение. Разработана и проверена процедура коррекции матрицы демпфирования на основе опытных данных. При коррекции использованы коэффициенты демпфирования, найденные при модальных испытаниях балочной модели.

Апробирована программа моделирования испытаний при многоточечном возбуждении. Полученные расчетным путем АЧХ близки к экспериментальным.

ЛИТЕРАТУРА

[1] Ильин М.М., Колесников К.С., Саратов Ю.С. Теория колебаний. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2003, 271 с.

[2] Вибрации в технике. Справочник. Челомей В.Н., ред. T. 1. Москва, Машиностроение, 1978, 352 с.

[3] Валяев В.И., Кожевникова Г.Л. Внутреннее и конструкционное демпфирование в механических системах и его влияние на динамические характеристики авиационных конструкций. Обзор. Москва, ЦАГИ, ОНТИ, 1973, № 402, 74 с.

[4] Дмитриев С.Н., Хамидуллин Р.К. Уточненная формула для вычисления коэффициентов передаточной матрицы в задачах статистической динамики.

«Наука и образование: электронное научно-техническое издание». 2013, № 3, doi: 10.7463/0313.0543198

[5] Писаренко Г.С., Матвеев В.В., Яковлев А.П. Методы определения характеристик демпфирования колебаний упругих систем. Киев, Наукова думка, 1976, 88 с.

[6] Иориш Ю.И. Виброметрия. Измерения вибрации и ударов. Общая теория, методы и приборы. Москва, Государственное научно-техническое издательство машиностроительной литературы, 1963, 772 с.

Статья поступила в редакцию 21.05.2013

Ссылку на эту статью просим оформлять следующим образом:

С.Н. Дмитриев, Р.К. Хамидуллин. Коррекция матрицы демпфирования с использованием экспериментальных значений коэффициентов модального демпфирования.

Инженерный журнал: наука и инновации, 2013, вып. 3. URL: http://engjournal.ru/ catalog/machin/rocket/619.html

Дмитриев Сергей Николаевич родился в 1950 г., окончил МВТУ им. Н.Э. Баумана в 1973 г. Канд. техн. наук, доц. кафедры «Аэрокосмические системы» МГТУ им. Н.Э. Баумана. Автор более 70 работ в области динамики конструкций летательных аппаратов. e-mail: dim.sm2@yandex.ru

Хамидуллин Руслан Камилевич родился в 1989 г., окончил МГТУ им. Н.Э. Баумана. в 2012 г. Инженер ОАО «ВПК «НПО машиностроения». Автор двух работ в области исследования вибраций. e-mail: brus555555@yandex.ru