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

МЕХАНИКА И МАШИНОСТРОЕНИЕ

УДК 519.6

ИСПОЛЬЗОВАНИЕ ВЕЙВЛЕТ-АНАЛИЗА СВОБОДНЫХ КОЛЕБАНИЙ ТРУБОПРОВОДНОЙ СИСТЕМЫ ПРИ ОЦЕНКЕ ХАРАКТЕРИСТИК ДЕМПФИРОВАНИЯ

© 2004 А.Б. Прокофьев, Е.В. Шахматов

Самарский государственный аэрокосмический университет

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

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

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

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

ленным на исследуемом элементе трубопроводной системы. Обработка экспериментальных данных при этом проводится в 3 этапа.

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

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

3. По отфильтрованному сигналу определяется логарифмический декремент затухания.

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

Рис. 1. Амплитудный спектр собственных колебаний трубопровода

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

Вейвлет-преобразование отображает одномерный сигнал в двумерное вейвлетное пространство с координатами "время - масштаб (частота)". Таким образом, исходный сигнал аппроксимируется набором вейвлетов, полученных из "материнского" (или "базового") вейвлета путем сдвига и растяжения:

)=((а,Ь^ )=-1 Г—1, (1)

л/а V а ) к '

где ) - используемая функция

вейвлет-преобразования;

а - коэффициент, определяющий растяжение;

Ь - коэффициент, определяющий сдвиг;

у/0 (^) - материнский вейвлет.

1

Множитель ~Г= в выражении (1) введен

имели постоянную (единичную) норму, т. е.

+ю

^ = 1,

—ю

где " * " указывает на операцию комплексного сопряжения. Кроме того, для вейвлетной

функции ) должно выполняться условие

нулевого среднего:

t d = о

—ю

В основе идеи вейвлет-преобразования лежит представление сигнала х({) в виде

взвешенной суммы базисных функций (к ),

умноженных на коэффициент Ск:

х ( )=£ Ск(к (t)

к

Прямое непрерывное вейвлет-преобра-зование сигнала х^) задается путем вычисления вейвлет-коэффициентов с(а,Ь):

(t —Ь а

С учетом ограниченной области определения сигналов, а также того, что а,Ье R (где R - множество действительных чисел), а Ф 0:

^ — Ь

i /1 _ ь \

C(a,b)= fx(t)^ц -Id

_L Va ^ a )

C (a,b)= f x(t ц

R \ a

\dt

Существует обширный выбор типов вейвлетов. Конкретный выбор анализирующего материнского вейвлета определяется тем, какую информацию необходимо извлечь из сигнала. Наиболее часто в области анализа процессов в гидромеханических системах применяются вейвлеты MHAT (Mexican Hat -мексиканская шляпа) и Морлета [7]. Часто используемый комплексный морлет-вейвлет описывается выражением:

Ц ) =

1

^fb

е2Ш • e

для того, чтобы все вейвлетные функции y(t)

где fB - ширина полосы частот,

fc - центральная частота вейвлета.

a

B

[6] и программно реализовано в расширениях ряда систем компьютерной математики (например, пакет Wavelet Extension для системы Mathcad или пакет Wavelet Toolbox для системы Matlab).

Рассмотрим непрерывное вейвлет-пре-образование с использованием комплексного вейвлета Морлета гармонического сигнала x{t )= cos at.

1

л1¥ва Zo

= a{I, + jI2 ),

+Lj2f ?UI"

где

_ 2jnfcb

e a

cos atdt =

(2) (3)

j = 4rfs 11 = 2

а

cos

2nfc

a

a

e

_(2fc _aa)

+

+ cos

2ж/с

+ a

a

_(2nfc + aa)

Рис. 2. Комплексный вейвлет Морлета, f,=fc=1: а - действительная часть; б - мнимая часть

На рис. 2 представлены графики действительной и мнимой частей комплексного вейвлета Морлета с параметрами

/в = /с = 1 .

Морлет-вейвлет обладает базисом, хорошо локализованным в реальном и фурье-

пространстве, причем с увеличением / растет разрешение в фурье-пространстве, но ухудшается локализация во времени [5]. Морлет-вейвлет представляет собой синусоидальную функцию, моделированную функцией Гаусса. Он обеспечивает наилучшее соотношение между разрешением по времени и разрешением по частоте [7]. Разрешение по времени изменяется в зависимости от частоты таким образом, что большим частотам соответствует большее разрешение по времени. Вейвлет-преобразование дискретных во времени сигналов использует алгоритм Малла

(4)

1 =Jöfs 12 = 2

a

sin

b

2f

_ a

\ a

_ (2fc + aa)2

+

+ sin

2f

+ a

a

(2f _ aa )2

(5)

Соотношения (3), (4) и (5) позволяют получить аналитическое выражение для коэффициентов рассматриваемого вейвлет-пре-образования. Следует отметить, что эти коэффициенты являются комплексными величинами. Таким образом, вейвлетный спектр одномерного сигнала, определяемый коэффициентами С (a,b), представляет собой две

поверхности амплитуды |c(a,b) и фазы

arg С(a,b) коэффициентов вейвлетного преобразования в трехмерном пространстве. В данной работе вейвлетный спектр представляется в виде проекций на плоскость (a, b) с

2

b

b

e

e

b

изолиниями или изоуровнями, позволяющими проследить изменение интенсивности

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

Для получения амплитудного вейвлетно-го спектра гармонической функции с использованием вейвлета Морлета определим модуль комплексного выражения (2):

\С (a,b ] = \A\ ■JtT+lf =—e"f* +0'25aW > ■

2

^ e2*Js™ + e -fam + 2cos2b&.

(6)

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

®прив

а

f

• t = t ■ f

5 прив J i

диск

диск ■

где /диск - частота дискретизации исходного сигнала.

Проекция поверхности амплитуды коэффициентов вейвлетного преобразования гармонического сигнала, построенного с базовым вейвлетом Морлета, представляет собой практически неизменяющуюся с увеличением параметра Ь картину - полосу почти постоянной интенсивности, параллельную оси времени. Сечение поверхности амплитуды коэффициентов плоскостью Ь = СОт1 представлено на рис. 3. График рис. 3 имеет четко выраженный максимум при

2ж/с

а0

а

(7)

Рис. 3. Зависимость амплитуды вейвлет-коэффициентов гармонического сигнала от масштаба преобразования а при Ь=const

Выражение (7) представляет собой приближенное соотношение между масштабами вейвлетного преобразования с применением вейвлета Морлета и частотами гармонического анализа.

Зависимость модуля коэффициентов вейвлетного преобразования от времени строго говоря является периодической функцией, что определяется слагаемым 2 2Ь® в подкоренном выражении соотношения (6).

0,7

Однако при а > ^ ^ в подкоренном выражении и второе слагаемое (е-2жМваа>) и третье (2cos2Ьа^) становятся пренебрежимо малы по сравнению с первым. В этом случае соотношение (6) можно переписать в виде

\c(a,b\ = —e -fB № +0'25a^2 V eMfa(° = 2

a e-fB(f-0,5aaf

(8)

2

С учетом выражения (8) можно уточнить

соотношение (7), взяв производную приравняв ее к 0.

d\C\ _ fB (f -0,5am)1

d\C\ da

и

da

= e

j—Ц + ^■ fBaffc - 0,5a®)\ 144a 2 J

= 0.

(9)

Решение уравнения (9) имеет вид

f +, Ж2/2 +

a

При

fc >>

1 . (10) — ^f соотношение (10) со-

ответствует (7).

Выражение (10) можно переписать в

виде

>„=^■ k(, /

a

(ii)

¡it).

\1, при t > 0, \0, при t < 0;

a - логарифмический декремент затухания.

С(a,b)-

1

ffBa L

hOU i

j2f

t-b1 1 ( t-b

a ¡e fB I a

■ 1(t) ■ Ae ~a cosatdt = B ■ I,

fb

где B =

Ae

/в«2

i 2 j

зуемся результатами, представленными в /8/:

—ax -cx

cos bxdx = ¡л —e 4 \ a

ic + jb )2

erfc

rc + jb^

24 а

+

ic - jb)2

+e

erfc

r c - jb ^

2^ a

При /в = /c = 1 k = 1,025. При /в = 1

и /c = 1,5 k = 1,011.

Если рассматриваемая гармоническая функция имеет амплитуду, отличную от 1 (на-I бе! аб, Х), то для вычисления модуля коэффициентов вейвлет-преобразования необходимо результаты, получаемые по выражениям (6) и (8), умножить на Х.

Рассмотрим теперь непрерывное вейв-лет-преобразование на базе комплексного вейвлета Морлета затухающей периодической функции x(t) ,имеющей нулевые значения при t < 0 :

x(t)= 1(t)Ae~at cos at, где 1(t) - единичная ступенчатая функция Хевисайда:

справедливыми при Rea > 0 .

Выше обозначено: erfc(x) = 1 - er/ (x) - дополнительный интеграл вероятности;

2 x 2

er/ (x) = "^i e - интеграл вероятности. ЫЖ 0

Тогда

I = 0,25 a J f

0,25/в«

er/

0,5 a4fB

e

a

2b .(2—fc a - j-j \ --a

fвa

(1 -

2b .2f .

—2 - j— + ja

fBa

a

\\

¡JJ

+

0,25/ba2

+ e

2b (2f a--2-j\ ■

ib«2 V a

(1 -

er/

0,5'4fB

a

2b - 2f

/■ 2 J

b« a

\W

ja

¡¡¡

(13)

Амплитудный вейвлетный спектр затухающей гармонической функции с учетом выражений (12) и (13) может быть записан в виде

\С(а,Ь\ = |В| • |/| ,

(14)

где B

Ae

yf/B

yl—fBa

QU

I = i e

0

Для вычисления интеграла I восполь

2 I 2b 2 j—fc ,

т*2+\-\t

На рис. 4 представлено сечение поверхности амплитуды вейвлет-коэффициентов (при параметрах анализа и модельного сигнала /в = /с = 1; а = 0,2; а = 0,003) плос-стаШ. (12) костями Ь = сОП^ . Рис. 4,а построен при Ь = 50; рис. 4,б - при Ь = 200. Видно, что рис. 3 и 4 во многом похожи, т. е. форма сече-

a

0

2

2

2

2

b

2

a

a

a

1

а) ~ б)

Рис. 4. Зависимость амплитуды вейвлет-коэффициентов затухающего гармонического сигнала от масштаба преобразования а при b=const. а - Ь=50; б - Ь=200

ния поверхности амплитуды вейвлет-коэффи-циентов гармонических незатухающего и затухающего сигналов плоскостью Ь = СОт( аналогична. Однако с ростом коэффициента Ь амплитуда всех точек последней поверхности монотонно уменьшается (в пределе -до нуля). Сечение поверхности амплитуды вейвлет-коэффициентов плоскостью а = СОт1 при указанных выше параметрах анализа и модельного сигнала представлено на рис. 5. Данная кривая имеет максимум. При масштабе а, определяемом из выражения (11), он соответствует Ь = 45. С уменьшением а величина Ь , при которой имеет место максимум амплитуды вейвлет-коэффи-циентов рассматриваемого сечения, также уменьшается, и, наоборот, при увеличении а - увеличивается значение Ь . Возрастающий характер рассматриваемой кривой при малых Ь объясняется тем, что при уменьшении Ь вейвлетная функция, приближаясь к границе 7 = 0 (что соответствует Ь = 0 ), начинает выходить за пределы области 7 > 0, где модельный сигнал равен 0. Чем меньше масштаб а , тем меньше временной промежуток, на котором локализована вейвлетная функция, что и объясняет приближение максимума коэффициентов |С(а,Ь) рассматриваемого сечения к Ь = 0 с уменьшением а . В дальнейшем анализе будем рассматривать только

Рис. 5. Зависимость амплитуды вейвлет-коэффициентов модельного затухающего гармонического сигнала от параметра Ь при а=а0.

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

Ь > Ь0, где Ь0 - значение временного масштаба, при котором достигается максимум значений коэффициентов |с(а,Ь) сечения

а = СОт1, соотношение (14) может быть записано приближенно в упрощенном виде:

с (а,Ь)

1 — ег/

Ае

/Ва0

—ва0

• 0,25агА жЖве

0,25/ва20 I а-

2Ь ^ /ва2о )

0,5а^/в

2Ьп

\\

а —

а2

Ва0)))

/ва

= Ze ~аЬ.

где а0

2/ а

(15)

к (/в, /с),

Z = 0,25А^а~е'

0,25/ва'а

1 — ег/

0,5а„4Гв

/ва

а2

ва0 )))

2 Ж с

ак 0

и

0,25/ва2

2Ь (2—с ^ ва2 V а )

(1 —

ег/

0,5а 41/в

а

2Ь /

— — 1

Iва

2ж/с

\\

+ а

V

)))

«

«

0,25/ва-

2Ь (2ж/с

/ва2 \ а

(1 —

— ег/

0,5ал/в

а—

2Ь

/ва

— ]

2—с

—а

)))

|С (а0,Ь1) = Ze

С {а0,Ь2 ) = ^

—аЬ1 —аЬ2

(16)

При этом предполагается, что Ь1 > Ь0 и

Ь2 > Ь0. Решение системы уравнений (16) имеет вид:

С(а0,Ь1)

1

а =

1н\

(17)

Следует обратить внимание, что при а0 = сОтг Z = сОтг , т. е. z не зависит от Ь .

При выводе уравнения (15) учтено, что при а = а0 и Ь > Ь0

Анализ соотношения (15) показывает, что при а = а0 = сош1 сечение поверхности амплитуды вейвлет коэффициентов при Ь > Ь0 представляет собой экспоненциально убывающую кривую. Показатель экспоненты представляет собой произведение логарифмического декремента и временного масштаба, взятое со знаком "минус". Для определения логарифмического декремента а по амплитуде вейвлет-коэффициентов составим систему 2 уравнений:

Ь2 — Ь1 \С (а0,Ь2 ) .

Соотношение (17) не изменится и в том случае, если анализируемый сигнал описывается выражением вида

х(г)=1(г—Ь')Ае -а('—Ь,) (г—Ь')),

т. е. импульсное воздействие прикладывается к системе в момент времени г = Ь'. Однако при этом следует учитывать, что соотношение (17) будет справедливо при Ь > Ь0 + Ь'.

Полученное решение может быть распространено и на дискретное вейвлет-преоб-разование.

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

1. Осуществить дискретное комплексное вейвлет-преобразование с использованием вейвлета Морлета.

2. По виду поверхности амплитуды вей-влет-коэффициентов определить величину

масштаба а0 вейвлет-преобразования, при котором амплитуды достигают максимальных значений.

3. Определить величину временного

масштаба Ь0 + Ь', при котором перестают влиять на поверхность амплитуд вейвлет-ко-эффициентов анализируемого сигнала эффекты на его границе.

4. Для двух временных масштабов Ь1 и Ь2 определить амплитуды вейвлет коэффициентов С(а0,Ь1) и С(а0,Ь2 ) и произвести

расчет логарифмического декремента а по соотношению (17).

Ситуация не изменяется в случае, если

2

Ь

а

а—

е

а

2

е

а

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

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

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

плоскостью а = а0 = 32 ;

1 - аппроксимация сечения поверхности амплитуды вейвлет-коэффициентов при а = а0 = 32 и Ь>>200

рованным с помощью генератора случайных чисел, описываемая выражением вида

x\t\ = 1(t - b'X1 + rnd (1))e-a{t-b'] cos[a(t - b')\,

где b' = 100, a = 0,2, a = 0,003, rnd (a) -

случайное число в диапазоне 0...a .

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

a = a0 = 32 при b > 200 экспоненциальной зависимостью (рис. 6,в, кривая 1) по методу наименьших квадратов позволяет найти значение логарифмического декремента a = 0,003074 . Погрешность найденного значения не превышает 2,5%.

Для перехода от вычисленного по приведенной выше методике логарифмического

декремента aebn к его действительному значению a , свободному от предположения, что отсчеты анализируемого сигнала отстоят друг от друга на единичном интервале времени, необходимо воспользоваться соотношением:

a = a еыч ■ /диск .

Для экспериментального определения характеристик демпфирования по разработанной методике использовался трубопровод, описанный в [9]. Схема исследуемого участка представлена на рис. 7. Результаты обработки экспериментальных данных изображены на рис. 8. На рис. 8,а представлена проекция поверхности амплитуды вейвлет-коэффи-циентов. Видно, что амплитуды имеют максимум при a = 29, что при используемом вей-

влете Морлета с /в = /с = 1 и частоте дискретизации /диск = 20 кГц соответствует собственной частоте трубопровода / = 707 Гц .

Рис. 7. Схема исследуемого трубопровода

а)

|С(29,Ь)|

0.6 ■

0.4 ■

0.2 ■

0.0 ■

100

200

300

400

500

б)

600 b

Рис. 8. Вейвлет-преобразование собственных затухающих колебаний трубопровода: а) проекция поверхности амплитуды вейвлет-коэффициентов; б) сечение поверхности амплитуды вейвлет-

коэффициентов плоскостью а = а0 = 29 ; 1 - аппроксимация сечения при Ь>>Ь0.

Аналогичная собственная частота ( f = 703 Гц ) была получена с использованием спектрального анализа в [9]. Она соответствует третьей форме колебаний рассматриваемого участка. Сечение поверхности амплитуды вейвлет-коэффициентов плоскостью

а = а0 = 29 представлено на рис. 8,б. Аппроксимация данного графика экспоненциальной кривой по методу наименьших квадратов позволяет определить логарифмический декремент, который составляет а = 121.

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

Работа выполнена при поддержке Ми-

нистерства науки и образования РФ, Администрации Самарской области и Американского фонда гражданских исследований и развития (CRDF Project SA-014-02) в рамках российско-американской программы "Фундаментальные исследования и высшее образование" (BRHE).

СПИСОК ЛИТЕРАТУРЫ

1. Старцев Н.И. Трубопроводы газотурбинных двигателей. М.: Машиностроение, 1976.

2. Методическое и программное обеспечение автоматизированного эксперимента в динамике машин / М.Б. Левин, А.Б. Оду-ло, Д.Е. Розенберг и др. М.: Наука, 1989.

3. Бабаков И.М. Теория колебаний. М.: Наука, 1965.

4. Микишев Г.Н., Рабинович Б.И. Динамика тонкостенных конструкций с отсеками, содержащими жидкость. М.: Машиностроение, 1971.

5. Короновский А.А., Храмов А.Е. Непрерывный вейвлетный анализ и его приложения. М.: Физматлит, 2003.

6. Дьяконов В.П. Вейвлеты. От теории к практи-

ке. М.: Солон-Р, 2002.

7. M. Prek. Wavelet analysis of sound signal in

fluid-filled viscoelastic pipes // Journal of Fluids and Structures. 2004. Vol. 19. Iss. 1. 8. Интегралы и ряды / Прудников А.П., Брычков Ю.А., Маричев О.И. М.: Наука, 1981.

9. Прокофьев А.Б. Расчет собственных частот и форм колебаний трубопроводов с помощью программного комплекса // Известия Самарского научного центра РАН. 1999. №2.

USAGE OF WAVELET ANALYSIS OF PIPELINE SYSTEM FREE OSCILLATIONS FOR ESTIMATION OF DAMPING CHARACTERISTICS

© 2004 A.B. Prokofiev, E.V. Shakhmatov Samara State Aerospace University

One of the characteristics of vibration suppression of a pipeline system is the logarithmic decrement. The technique of definition of a logarithmic decrement on eigentones oscillogram of a system with the help of the wavelet-analysis using as a mother function complex Morlet wavelet is offered. The results of model and full-scale experiments are adduced.