Спектральный анализ свободных колебаний гибкой оболочки Текст научной статьи по специальности «Физика»

УДК 539.1

А.Н. Куцемако

СПЕКТРАЛЬНЫЙ АНАЛИЗ СВОБОДНЫХ КОЛЕБАНИЙ ГИБКОЙ ОБОЛОЧКИ

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

Динамика геометрически нелинейной оболочки, скрытые периодичности, амплитудно -частотная характер истика

A.N. Kutsemako

FREE VIBRATIONS SPECTROLOGY OF THE FLEXIBLE SHELL

A new technique dealing with the numerical solution of the problem was offered for finding hidden periodicities. Amplitude-frequency characteristics were developed for free vibrations of the geometrically nonlinear plate and shell.

Dynamics of the geometrically nonlinear shell, hidden to periodicity, gain-frequency characteristics

1. Свободные нелинейные колебания пластин и оболочек

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

Методика решения задачи. Для решения поставленной задачи в рамках геометрически нелинейной модели Кирхгофа-Лява динамические уравнения оболочки в смешанной форме в безразмерном виде запишем [1], положив в них q = 0, £ = 0 :

д2 w 2 / ч 1 г 1 д4 w .2д4 w д4 w

—=V2F+L(w, F )-12^ [-—+2дд71 • (11)

£l£ + 2-dJ— = -V2w-1 L(w, w). (1.2)

l2 д x4 д y' д xд у2 ' 2 V '

Граничные условия - шарнирное опирание на гибкие несжимаемые (нерастяжимые) в касательной плоскости ребра, которые могут быть представлены в виде

д2 w д2 F

w =-= F =-= 0 при x = 0;1,

1 2 1 2 F ' '

дх д Х (1.3)

д2 w д2 F w =-= F =-= 0 при y = 0;1.

д У2 д у2

Начальные условия

д w

=t0 д "=к где w0 берем из решения соответствующей статической задачи.

Исходные уравнения (1.1), (1.2) решаем методом Бубнова-Галеркина в высших приближениях. Для этого искомые функции w,F , удовлетворяющие граничным условиям (1.3), представим в виде:

w = ^ Aij (t )sin(i px)sin( jpy),

i, j

F = ^ Atj (t )sin(ipx)sin( jpy), ( 1 . 5 )

i, j

W = w0, —I = 0, (1.4)

li=to 0' ■} s Ii=t„

(/ = 1,2,...Мх; ] = 1,2,...Му; и = Мх *Му ).

После применения к (1.1) по пространственным координатам процедуры Бубнова-Галеркина, получим систему обыкновенных дифференциальных уравнений по времени, а к (1.2) - систему алгебраических уравнений.

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

Рис. 1.1. Изменение во времени прогиба в центре квадратной пластинки при различных начальных отклонениях

На рис. 1.1 показано движение центральной точки квадратной пластинки во времени при и = Мх *Му = 25 для w0 = 0,5; 1; 2.

Колебания носят, хотя и слабо выраженный, нелинейный характер и только при малых начальных отклонениях близки к гармоническим.

На рис. 1.2 - то же, но для оболочки с геометрическими параметрами кх = ку = 12 для

w0 = 1; 3; 5 . Здесь даже при малых начальных отклонениях колебания носят ярко выраженный нелинейный характер.

Рис. 1.2. Изменение во времени прогиба в центре квадратной оболочки (кх=ку=12) при различных начальных отклонениях

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

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

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

2. Спектральный анализ решения

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

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

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

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

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

Пусть число ординат, если процесс задан в виде графика, есть нечетное число 2Л^+1. Примем за начальный момент отсчета времени 1 = 0 момент, соответствующий средней точке наших ординат. Пусть ординаты получены с интервалом времени т , тогда, с переходом к новой переменной

1

t =—11 t

ординаты будут соответствовать моментам времени

(2.1)

tk = 0,± 1,± 2,...± N.

Для ординат введем обозначение

Л = / (1к ). (2.2)

Функция /(1), вообще говоря, имеет следующую форму:

f (t) = Е ( Aa cos eat + Ba sin в J ) .

(2.3)

a=1

Число периодических компонент, обозначенных здесь через 7, обычно заранее не известно.

ва „

Точно так же мы ничего не можем сказать о величинах угловых частот (Оа = — . Однако мы заранее

Т

знаем, что в новой переменной 1 угловые частоты 0а должны ограничиваться неравенством 0а < К, так как, еслива превышает этот предел, то две частоты К + ¡5 иК-¡5 не могут быть различимы. Далее, положим

а Ж

ва= ^Ра, (2.4)

и пусть ра принимает значения между 0 и N .

Отделим синусы от косинусов, образовав суммы и разности

f (t) + f (-t) = 2Е Aa cos eat, f (t) - f (-t) = 2Е Ba sin в J .

a=1 a=1

Ординаты разделятся соответственно на две группы:

(2.5)

uk = fk + f-k, Vk = fk-f_k, (k=0, 1, 2,..., N). (2.6)

Воспользуемся методом преобразования Фурье, приспособленным для суммирования вместо интегрирования. Первоначальный ряд данных ик преобразуем в новый ряд N +1 амплитуд ak при

косинусах, а данные Vk — в ряд N — 1 амплитуд bk при синусах. Эти амплитуды получаются умножением исходных данных на заранее составленную матрицу, содержащую косинусы и синусы углов, Ж

кратных —:

N

N, Ж N—1

ak = Z и a eos—di, bk = ^Va sin —ak. (2.7)

a=0 N d=í N

Штрих в символе отмечает тот факт, что первое и последнее данные значения функции и0 и uN входят в сумму с множителем 0,5. Полученные коэффициенты йк и bk можно представить себе как «линейный спектр», соответствующий целым значениям pa непрерывного параметра p . Весь дальнейший анализ будет применять эти две последовательности коэффициентов.

Прежде всего, рассмотрим случай, когда заданные частоты ва таковы, что все или некоторые pa в формуле (2.4) оказываются целыми числами. В этих местах получается одиночный пик без каких - либо «хвостов». Отдельный максимум с обеих сторон сопровождается нулями.

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

в = Nk. (2.8)

Соответствующая амплитуда при косинусе и, соответственно, при синусе будет

A = Nak, 5 = . (2.9)

Если, однако, максимум не попадает на целое значение, то максимальная амплитуда сопровождается с обеих сторон меньшими амплитудами. Наименьший спад получается, если максимум находится в точности на полпути между двумя целыми значениями p . Тогда схема амплитуд будет

1 — 1 1 — 1 1 1 — 1 1 — 1

9, 7, 5, 3, , , 3, 5, 7, .. (2.10)

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

^k = Ук—1 + 2Ук + Ук+1. (2.11)

Предыдущая схема изменяется следующим образом:

—_8__— 8 8 —_8__L

693,3í5, í05, í5, 3, 3, í5, 105' 315'...' (212)

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

Точное положение максимума, вычисленное с помощью вторых разностей (фактически вторых сумм, ввиду чередующихся знаков схемы), может быть получено следующим образом. Проверяем последовательность yk и обращаем внимание особенно на правильное чередование знаков ±. В некоторых точках отмечаем, что чередование нарушается следованием знаков ++ или — . Подчеркиваем эти неправильности следования и устанавливаем, что пик должен находиться между двумя такими значениями p :

p = к и p=k+1. (2.13)

В общем случае точное значение p , соответствующее некоторой частоте ea, будет лежать

между двумя целыми числами к и к +1. Далее полагаем

р = к + е. (2.14)

и получаем е при помощи интерполирования. Для этой цели пользуемся методом «вторых сумм». Составляем отношение двух последовательных г , соответствующих р = к и р=к+1:

p

Z,,

Но zk пропорционально

а zk+1 пропорционально

qk =

Zk+1

2 1 2 ■ +--

(2.15)

e-1 e e +1 (e + 1)e(e-1)' (216)

1 2 1 2 ■ + ■

e e-1 e-2 (e-2)(e-1)e' (217) Составляя отношение, получаем

zk 2 -e

qk = тг=7+7 ■ (218)

откуда

2 - qk

e = -

1 + Чк

Наконец, после определения е получаем частоту, используя соотношение ва

(2.19)

Р

qa= N(k + e) . (2.20)

Коэффициенты Aa,Ba, соответствующие этой частоте, выражаются через амплитуды ряда Фурье йк и Ьк:

Aa = , Ba= . (2.21)

N sin pe N sin pe

Описанный метод выделения периодических составляющих весьма удовлетворителен, если общее число наблюдений 2N +1 достаточно велико для того, чтобы изолировать два соседних «пиковых» значения. Необходимо, чтобы два последовательных значения pa в соотношении (2.4) были

отделены друг от друга, по меньшей мере, четырьмя единицами.

Если они расположены ближе друг к другу, то их взаимное влияние возрастает; в этом случае трудно выделить их надлежащим образом.

Этот метод отыскания «скрытых периодичностей» суммы периодических функций может быть назван спектроскопическим методом, так как он математически имитирует действие спектроскопа. Спектроскоп вскрывает частоты, из которых состоит световое излучение возбужденного атома. Эти частоты могут быть вычислены из положения «спектральных линий». Спектральные линии не являются линиями в математическом смысле, но имеют некоторую конечную ширину. Исключительная точность спектроскопических измерений обусловливается большим постоянством оптических колебаний, которое соответствует весьма большому значению величины N . Существует возможность уменьшить число N и все же сохранить высокую точность ввиду большого числа значащих цифр, которыми мы располагаем. Хотя пики теперь расположены значительно теснее, можно все же скорректировать наши предварительные результаты, рассчитав интерференцию соседних пиков и вычтя их влияние.

2. Оценим необходимую длину интервала записи (наблюдения) T непрерывного сигнала и шаг отсчета. В общем случае отсчеты большинства непрерывных сигналов s(t) производятся через фиксированный интервал А, и полученные таким образом дискретизованные сигналы (ординаты) используются для дальнейших вычислений.

Дискретизованный сигнал можно рассматривать как результат умножения первоначального непрерывного сигнала на сигнал i(t), состоящий из бесконечного ряда единичных импульсов, или дельта-функций [3]:

1

¡(г) = £ ё(г — пк). (2.22)

П= — ¥

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

5. (г) = )1(г). (2.23)

Следовательно, воспользовавшись теоремой о свертке, находим

(/) = | ^ (/ — % )1 (ё ,

(2.24)

где I (ё) является преобразованием Фурье от ¡(г). Используя для I (ё) выражение

1 +¥ п 5 (/)Л),

ЛП^ Л

(2.25)

преобразуем

+¥ 1 +¥ П 1 +¥ П (1) = / 5(/ — ёЛ £ — = Л £ ^(I— Л). (2.26)

— ¥ Л П=—¥ Л Л П=—¥ Л

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

сигнал я. (г) имеет периодическое преобразование Фурье с периодом — и, если 5 (/) обращается в

Л

нуль при 1/1 >-, то Sj (/) является просто периодически повторяемой функцией 5(/). Это

2Л

означает, что можно восстановить 5(/) по (/) , умножив (/) на Н (/)

Л / ^

Н (/) =

2Л (2.27)

0,1/1 >—. 2Л

Так как умножение в частотной области соответствует свертке во временной области, то отсюда следует, что

5(г) = г $т{т1Л) ^ (г — и. ^^

—¥ ши / Л

Функция ^п(ши / Л) служит идеальным фильтром для восстановления непрерывного сигнала ши / Л

) из дискретизованного сигнала (г) .

, $лп(ши / Л) „ „ ,

Иначе говоря, функция - является идеальной интерполирующей функцией для рав-

ши / Л

ноотстоящих ординат. Следовательно, длина записи Т определяет степень различимости пиков в преобразовании Фурье, а интервал отсчета Л определяет максимальную частоту, которую можно различать. Частота

/• 1

/ы = — (2.29)

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

Таким образом, для разделения двух пиков на частотах /1 и / 2 необходимо использовать запись Т порядка

Т

/,— / (2.30)

для прямоугольного окна [3].

3. Сходимость метода

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

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

Дано:

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

Требуется:

Исходя из получаемой таблицы как исходных данных с помощью предложенной методики:

1) выделить гармоническую составляющую, т.е. численно установить ее частоту и амплитуду,

2) численно оценить получаемую при этом погрешность.

Пусть такая функция будет

Р

у = 5* соз(— ^)

(3.1)

причем значение частоты гармоники намеренно выбрано так, чтобы применительно к нашим основным объектам исследования (оболочкам и пластинкам), это была низкочастотная гармоническая составляющая и соответствующие косинусы на основном периоде (—Р,+Р) были достаточно произвольными числами.

Переменная / принимает значения, начиная с нуля с шагом отсчета

Р

—. (3.2)

N

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

Исходя из этой методики следует ожидать точное решение задачи, если выбирать N кратным 7, т.е. N = 7, 14,..., что на самом деле и происходит.

В общем случае, когда мы не знаем или не можем выполнить условие кратности, будем задавать N произвольно (не учитывая кратность), наращивая его с постоянным шагом.

В нашей задаче, очевидно, начальное значение N полагать меньше 7 не имеет смысла, так как в этом случае не обеспечивается основное условие, чтобы два соседних пика были разделены, по меньшей мере, четырьмя единицами.

Результаты проведенного исследования представлены на рис. 3.1, 3.2.

На рис. 3.1а показана зависимость частоты гармоники во всем диапазоне изменения параметра N . На рис. 3.1б в увеличенном масштабе представлен тот же график, из которого исключены малые значения параметра N.

0.455 0.45 0.445 0.44 0.435 0.43

ча< яот а—

V

0.44975 0.44955 0.44935 0.44915 0.44895 0.44875

час тота

7 14 21 28

35

а

42 49 56 63 N

14 21 28

35 б

42 49 56 63 N

Рис. 3.1. Зависимость частоты гармоники от N числа наблюдений, сделанных через равные промежутки времени

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

5,5 5,25 5

4,75 4,5

Л / /1 ампл итудг

/ / /1 л л л, Л Л Л Л Л Л -Л -Л -л -л

/ 1 1/ V V V V V \У

1/

21

35

49

63

77

91

105 119 133 147 N

Рис. 3.2. Зависимость амплитуды гармоники от N числа наблюдений, сделанных через равные промежутки времени

4. Спектральный анализ свободных колебаний

Свободные нелинейные колебания были рассмотрены в п. 1. Возвратимся к этому вопросу вновь и на основе предложенной методики проведем спектральный анализ получаемого решения. Несмотря на то, что пластинка или оболочка как колебательная система имеет бесконечное число степеней свободы, мы не можем ставить перед собой задачу определения всего спектра частот системы.

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

Краевая задача (1.1)-(1.3) решалась при п = Мх * Му = 25, а соответствующая задача Коши -

с шагом по времени Л1 = 0,001. Следовательно, помня о том, что два соседних пика могут быть удовлетворительно отделены друг от друга в том случае, если между ними находится не менее четырех единиц в (2.4), реально на таких данных можно рассчитывать на выделение гармоник с частотами не выше 80-100.

Чтобы обеспечить удовлетворительную степень различимости пиков в преобразовании Фурье, длина записи исходных данных (ординат) должна быть достаточно большой, т.е. охватывать несколько периодов исследуемого колебательного процесса. Руководствуясь сказанным, расчет по времени в безразмерных единицах проводим до 6-10.

2,5

1,5 1

0,5 0

^

М 1 ^ / к / / / , // / /

а / / /

1 /

1 частота

10

20

30

40

50

60

70

80

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

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

На рис. 4.1, 4.2 построены амплитудно-частотные характеристики для однородной квадратной пластинки и оболочки (кх = к = 12) соответственно.

7

2

0

7,5

5

2,5

0

0 10 20 30 40 50 60 70

Рис. 4.2. Амплитудно-частотная характеристика свободных нелинейных колебаний оболочки (кх=ку=12)

Из полученных результатов можно сделать следующие выводы:

1. Частота всех детектированных гармоник как для оболочки, так и для пластинки зависит от амплитуды свободных колебаний системы.

2. При малых амплитудах свободных колебаний частота низшей гармоники практически совпадает с первой частотой линейных колебаний системы.

3. Качественно пластинка - система с жесткой характеристикой, оболочка - система с мягкой характеристикой.

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

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

ЛИТЕРАТУРА

1. Куцемако А.Н. Устойчивость решения нелинейных дифференциальных уравнений в частных производных / А.Н. Куцемако. Саратов: СГТУ, 2007. 88 с.

2. Ланцош К. Практические методы прикладного анализа / К. Ланцош. М.: Физматгиз, 1961.524с.

3. Дженкинс Г. Спектральный анализ и 1972. Вып. 2. 356с.

Куцемако Анатолий Николаевич -

доктор технических наук, профессор кафедры «Математика и моделирование» Саратовского государственного технического университета имени Гагарина Ю.А.

Статья поступила в редакцию 11.11.13, принята к опубликованию

его применение / Г. Дженкинс, Д. Ваттс. М.: Мир,

Anatoly N. Kutsemako -

Dr. Sc., Professor

Department of Mathematics and Modeling Yuri Gagarin State Technical University of Saratov