Определение частот и декрементов собственных колебаний конструкции по переходным процессам Текст научной статьи по специальности «Математика»

_______УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Том XIX 198 8

№ 1

УДК 629.7.015.4:533.6.013.43

ОПРЕДЕЛЕНИЕ ЧАСТОТ И ДЕКРЕМЕНТОВ СОБСТВЕННЫХ КОЛЕБАНИЙ КОНСТРУКЦИИ ПО ПЕРЕХОДНЫМ ПРОЦЕССАМ

П. Г. Карклэ

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

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

т т

х (О — (Ai е ^ + Ке lt ) = ^ cos (“/ * + ¥/)' 0)

i-i l

где /и — количество частот собственных колебаний, которые удается возбудить и зарегистрировать; Ai = А- + jAi — комплексные амплитуды, р = — 1; А,- = щ + + у»; — собственные значения: А* = щ — jai\ щ — коэффициенты затухания; щ — циклические частоты; В/ = 2 У A'f + А-2 — амплитуды колебаний; <рг—начальные фазы колебаний; tg щ = А'^/А^ ; t — время; х (t) — перемещение, скорость или ускорение в некоторой точке конструкции, углы поворотов рулей или любая координата, описывающая колебания.

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

хк = х(Ш), k = 0, 1, N— 1. (2)

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

является определение собственных значений А.г = а*+ или частот колебаний

(|>£ ГГ,

в герцах /г = —— и декрементов 6г

2" Л

2. Для анализа последовательностей более удобна разностная форма соотношения (1). Рассмотрим случай яг=1, опустив на время индекс I у величин А и Л,:

хк = л: (Ш) = Аехш + А*ех*Ш. (3)

Запишем это же соотношение для моментов времени (6 + 1) Д^ и (к+2) Д*:

*А+1 = еШ + А* еХ*^> (4)

хк+2 = Аехше2ш + А* е'*ш е2к*А‘. (5)

Уравнение (3) вместе с (4) можно разрешить относительно АеХки и А*ех*ш-.

А* ех*ш = (хк+1 - хк еш)1(ех*и - еш),

Аехш = (хкех'и — хк+1)1 (ех*и — еш).

После подстановки в (5) получим

*£-1-2 — Ргхкл 1 " Ро Хк1 (6)

где коэффициенты Р1 и ро равны:

р1 = е^ + еш< Ро = _е»иеш

Полином с этими коэффициентами

Р(1) = 12 — Р21— Ро (7)

имеет корни - .

V __ „ _ * __

71 — ^ » Ъ — 71 — е >

из которых можно получить искомые а и 0).

В общем случае при т> 1 последовательность (2) будет удовлетворять соотношению

хк+2т== Р2т-\хк+2т-г Р2т—2хк+2т-2^ — +Рохк* (8)

что легко показать, рассмотрев 2т+1 уравнений и исключив из первых 2 т коэффи-диенты . Корни полинома

Р(1) = 12т-Р2т-1'12т~1-Р2т-212т~2- - -Ро (9)

позволят определить неизвестные

3. Так называемый «метод прямой идентификации» ([1, 2]) основан именно на

этих формулах. Для измеренных значений хк набирается столько уравнений (8), чтобы

полученную систему можно было разрешить относительно коэффициентов Ро, Ри ■ ■., р2 т — 1- Обычно используют метод наименьших квадратов, набрав максимально возможное количество уравнений в системе. Вычислив корни полинома (9) уг=Г1 + /4',, определяют

к = <ч + М = 1п ц = 1п (о + уяг); (Ю)

а, = 1п (]/г2 + 5?) ; (11)

1 5*

—— агс(§-------------, г = 1, 2, ..., т. (12)

Г1

Простота соотношений (8) казалась привлекательной, но первые же опыты практического применения выявили недостатки подхода. Первый недостаток касается определения числа т. Теоретически его надо выбирать максимально возможным, но чтобы система уравнений была еще не вырождена. Ошибки измерения и вычислений искажают, однако, эту картину, точное значение числу т подобрать оказалось очень трудно. При неправильном выборе т появляются большие погрешности в ро, рг, . •., Р2 т_ 1 ик следовательно, в yi либо за счет плохой обусловленности системы уравнений, либо за счет зависимости результата от номера к, с которого начата обработка. Второй недостаток — высокая чувствительность к случайным помехам (см. [3]).

В работе [4] предложен довольно сложный путь компенсации недостатков метода прямой идентификации. Процедура вычислений состояла из трех частей. В первой части выбиралось число т, но не программой, а человеком, проводящим обработку, на основе анализа матрицы Грамма, формируемой Из сдвинутых по времени последовательностей^, Хк+1........ Х/г+м)' (хк+\> хк+2> — ’ Хк+М +0> (хк+2> хк+3> —• хк+М+2)

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

4. Низкая чувствительность к случайным помехам, свойственная частотным методам исследования, связана с операциями интегрирования, используемыми для выделения компонент гармонического процесса. Похожая процедура иногда применялась для анализа колебаний вида (1). Если вычислить

со

Х(<0)= I* хЦ)е~М1 <и, (13)

о

то при а,ъ<0 выражение для Х(со) будет иметь вид рациональной дроби

*<“) = <?(“)//»(»).

где <2 (со) и Р(со) полиномы, степень которых определяется числом т. В работе [5], например, рассматривается один из таких, методов. Хотя определение Р(ш) и (2(«) связано с некоторыми трудностями, основные ошибки и ограничения появляются при вычислении интеграла (13). Верхний предел заменяется некоторым конечным числом Т, интегрирование заменяется суммированием, а чтобы значения хЦ) при 1—Т не усложняли вида дроби, надо подынтегральное выражение в (14) домножить на функцию вида ■*) (0 = е~*** с известным ц>0, если процесс не успел «затухнуть» сам. Потом это можно учесть, вычитая ц из вычисленных искаженных а,> Метод нагляден и достаточно надежен, широко применялся до разработки программ быстрого преобразования Фурье (БПФ).

Необходимо упомянуть также метод «моментов Пуассона» ([6], [7]). В нем используется интегральное преобразование вида:

Мк [*(*)] ^ = Х°Л = § — £ = 0, I.....

о

ъ р~~^ о

Рк (*о> = ---- •

В принципе можно изготовить прецизионные аналоговые схемы, выполняющие вычисления по указанным формулам, и получить соотношения вида (8), связывающие х°к но реализация алгоритма на цифровой машине пока не представляется целесообразной. Отсутствует также какой-либо опыт практической эксплуатации в отличие от методов 1 и 2. Тем не менее, в работах [6] и [7] показано, что теоретически подход позволяет решать задачу идентификации для сложных процессов, подчиняющихся дифференциальным уравнениям с переменными коэффициентами, с запаздывающим аргументом и т. д. В будущем, возможно, удастся на его основе создать процедуру, пригодную для практического применения.

5. Предлагаемый метод 3 базируется на процедуре БПФ, которая обычно выполняется в соответствии с формулой:

N—1

гі = С ^ хк

2Ш .

--іг>

(14)

*=0

Положим для простоты, что множитель С= 1. В программах это тоже часто делается для сокращения времени счета. Введем обозначение.

2ні

и рассмотрим снова случай т= 1. Для него можем написать:

Л-1

N-1

г1 = А X т? уТк + А*^ т’* *Г*.

й=0 А=0

„ш =

где 71 = е -', 7!

Выполнив суммирование по формуле для геометрической прогрессии, получи...

Б0 + О]

*1 = 4

поскольку 1.

тГ-1

-гГ-1

(15)

'•і—Рі'Ч - Ро

Коэффициенты £)0 и £>1 действительные:

А> = 7* + -4* тГ 71 - АЬ - 71.

£», = Л + А* - - А* 7*'у.

В общем случае при т» 1 имеем сумму слагаемых (15):

=і -(і—РчЧ—Роі

(16)

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

Вычисления будут состоять в следующем. Введенная в ЭВМ последовательность хн, 6=0, 1,...,ЛГ—1 подвергается преобразованию Фурье (14). Новая последовательность комплексных чисел г/ = гг + аппроксимируется формулой (16) методом работы [5], например. Полученные коэффициенты < и р0; используются для формирования ПОЛИНОМОВ (7), корни которых \г И 7; дают значения И 6{ (10) — (12). По коэффициентам Д> г и £>2 ,• можно восстановить Ап

Л =

1

(7г — О

Рі.

+ /■

Ро і Ч~ Р\і о 2 в/

(17)

Заметим, что для оперативных оценок можно пользоваться формулой (15), выделив из всей последовательности 2; несколько значений, ближайших к номеру 1=[ЫAt, где / — частота обрабатываемого тона, интересующего экспериментатора. Чем лучше разделены тона и чем меньше их демпфирование, тем в большей степени допустимо такое упрощение. Задача особенно облегчается, когда информационно-измерительная система оснащена устройством отображения обрабатываемых последовательностей.

6. Особенности методов обработки удобно выявлять по схеме, предложенной в работе [3]. В качестве тестовой используется последовательность

*к = хк +ук,

(18)

где хк — точно соответствует формуле (3), а у к — случайная последовательность, имеющая нулевое математическое ожидание и заданную дисперсию, имитирующая помеху. Варьируется величина е, характеризующая уровень помехи:

е = Бу/Бх,

где Яу и 5* — стандартные отклонения для у к и дг* соответственно.

Задача решается статистически. Для некоторого е при известном хь и 5* генерируется у к с заданным Бу — г 5Х и хк по формуле (18). Выполнится обработка, определяются / и 6, сравниваются с точными, по которым хк формировалась. Для статистической достоверности все это многократно повторяется с новой ун, пока не будут получены оценки смещения декремента Д6, стандартного отклонения и А/ и Б/ для частоты. Эти величины характеризуют в определенной мере методы, позволяют сравнить их друг с другом в одинаковых условиях.

£ -0,3

и,

Рис. 1

В качестве примера последовательностей хк и г; = г1 -|- ]г'1 на рис. 1 приведены два варианта, для случая с = 0,1 и 0,5. Точки соединены условно, для наглядности, к изменяется от 0 до 511.

7. Результаты сравнительного исследования точности определения частот и декрементов различными методами при наличии в сигналах помех приведены на рис. 2—4. Смещение оценки декремента изображено на рис. 2. График показывает, что методам

1 и 2 свойственно систематическое завышение декремента, растущее при повышении уровня шума. Для метода 1 пороговое значение составляет приблизительно 1%. При более интенсивной помехе декремент окажется завышен в два и более раз. Для метода 2 пороговое значение уже около 20%. Метод 3 практически не смещает оценку декремента, точек для него на графике нет. Смещение частоты тоже характерно для методов / и 2, но величина его достаточно мала, им можно пренебречь.

На рис. 3 изображен разброс оценок декремента, по оси ординат отложено стандартное отклонение. Преимущества методов 2 и 3 над методом 1 вполне очевидны. Интересно отметить, что при малых уровнях помехи уменьшение 55 не столь интенсивно, начинают сказываться ошибки вычислений. Операции выполнялись с пятью значащими десятичными цифрами. Разброс оценок частоты приведен на рис. 4. Порядок ошибок здесь существенно меньший, чем в случае декрементов, но характер кривых примерно такой же.

м%

wo-

rn -

/

/

/

/

/

/

/

L /

/

MJL

L

о метод /

* w

/

_L

J

’ W ДО/ 0,1 £=V^

Рис. 2

8. Тем же путем можно сравнить методы и при т> 1, однако новых результатов получить не удается. ’Ошибки оценок при увеличении т растут в основном за счет относительного увеличения уровня помехи для отдельных компонент. Особая ситуация возникает при сближении частот, когда вырождаются все три метода.

Необходимо также добавить, что формула (13) для частоты справедлива, если N

все частоты выше / = подавляются до операции сканирования. В противном

• случае может произойти «перепутывание» частот, характерное для обработки дискретных последовательностей.

ЛИТЕРАТУРА

1. Лебедев В. Г. Алгоритм определения собственных частот и декрементов колебаний по результатам измерений.—В сб. докладов III симпозиума: Колебания упругих конструкций с жидкостью. — М.: ЦНТИ «Волна», 1976.

2. Ibrahim S. R., М i k u 1 с i к F. С. A method for the direct identification of vibration parameters from the free response. — The Shock and Vibration Bulletin, 1977, N 47, Part 4.

3. P a p p a R. S. Some statistical performance characteristics of the ,,ITD” modal identification algorithm. — AIlAA, TP-82-0768, 1982.

4. К a p к л э П. Г. Определение частот и декрементов упругих колебаний конструкции в потоке по ее неустановившимся движениям. — В сб. докладов V симпозиума: Колебания упругих конструкций с жидкостью. — М.: ЦНТИ «Волна», 1984.

5. D a t R., М е u г z е с J. L. Exploitation par lissage mathematique des mesures d’admittance d’un systeme lineaire. — Recherche Aerospatiale, N 1972-4.

6. Saha D. C., Rao G. P. Identification of lumped linear time-varying parameter systems via Poisson moment functionals. — International Journal of Control, 1980, vol. 32, N 4.

7. Saha D. C., Rao G. P. Identification of lumped linear systems in the presence of small unknown delays — the Poisson moment functionals approach. International Journal of Control, 1981, vol. 33, N 5.

Рукопись поступила 8/Х 1986