Хаотическая динамика гибких криволинейных балок Бернулли-Эйлера (часть 1) Текст научной статьи по специальности «Физика»

УДК 534.1

Н.А. Загниборода, В.В. Добриян, М.В. Жигалов, А.В. Крысько, В.А. Крысько

ХАОТИЧЕСКАЯ ДИНАМИКА ГИБКИХ КРИВОЛИНЕЙНЫХ БАЛОК БЕРНУЛЛИ-ЭЙЛЕРА (ЧАСТЬ 1)

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

Хаотические колебания криволинейной балки Эйлера-Бернулли, аттракторы, бифуркации, фазовые портреты, показатели Ляпунова, хаос, хаос-гиперхаос, пространственно-временной хаос

N.A. Zagniboroda, V.V. Dobriyan, M.V. Zhigalov, A.V. Krysko, V.A. Krysko CHAOTIC DYNAMICS OF FLEXIBLE CURVILINEAR BERNOULLI-EULER BEAM (PART 1)

The work is devoted to the theory of chaotic dynamics offlexible curvilinear Eu-ler-Bernoulli beam. A mathematical model was constructed. Differential equations for the area, boundaries and initial conditions have been formulated. The basis of the mathematical model is based on the Euler-Bernoulli hypothesis, the geometric nonlinearity in the

form of Karman is taken into account, and the conditions of flatness curvilinear beams in the form of V.Z.Vlasov are accepted. The numerical methods for converting partial differential equations to ordinary differential equations (finite difference method of second-order accuracy and finite elements method) have been developed. The Cauchy problem is solved by the Runge-Kutta fourth and second order accuracy. The maps of the character of oscillations for a range of values kx, the appearance of elastic-plastic deformation maps, and the maps of dynamic instability, depending on the amplitude andfrequency ofcompelling vibration have been build.

Chaotic oscillations curvilinear beams, attractors, bifurcation, phase portraits, Lyapunov exponents, chaos, chaos-hyperchaos, spatio-temporal chaos

Введение

Анализ хаотической динамики гибких криволинейных балок Бернулли-Эйлера с позиции распределенных систем с бесконечным числом степеней свободы проводится впервые. Достоверность получаемых результатов обеспечивается решением задач несколькими методами: конечных разностей и конечных элементов для пространственной аппроксимации (для сведения к задаче Коши), и решением задачи Коши двумя методами Рунге-Кутта 4-го и 6-го порядков точности. Статья разбита на две части. В первой части осуществляется постановка задачи, приводятся методы её решения, а также осуществляется количественный анализ (исследуются различные виды карт). Во второй части производится качественный анализ, исследуется сценарий перехода колебаний тонкой криволинейной балки от гармонических к хаотическим, производится анализ колебаний с помощью показателей Ляпунова.

1. Математическая модель и метод решения

В работе рассматриваются гибкие однослойные, тонкие криволинейные балки с длиной а, высотой И и геометрической кривизной кх = —, где Ях - радиус кривизны. Балка нагружается распре-

их

деленной по ее поверхности нагрузкой д(х, £), действующей в направлении нормали к серединной поверхности балки (рис. 1).

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

— любое поперечное сечение, нормальное к серединной поверхности до деформации, остается после деформации прямым и нормальным к серединной поверхности, вместе с тем высота сечения не изменяется;

— инерция вращения элементов балки не учитывается, однако учитываются силы инерции, отвечающие за перемещения вдоль нормали к серединной поверхности и вдоль балки;

— внешние силы не меняют своего направления при деформации балки;

— геометрическая нелинейность учитывается в форме Т. Кармана [1];

— тонкая криволинейная балка обладает свойствами пологости [2].

Приведенные выше гипотезы основаны на идеях Бернулли-Эйлера и считаются математической моделью первого приближения, но она является достаточной точной для анализа тонких криволинейных балок [3, 4].

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

1 Э V

1

1 12 Эх4

Эи

Эх

-км-

1 Г Эм

21 Эх

-м'-

ЭV

"ЭХ^

+ Д (и. м>) + Х2(^. м") > + д-

Э V Э?2"'

Эм

= "■

(1.1)

Ц(и. м) =

Э2и Эм + Эи Э2м Эх2 Эх Эх Эх2

; £2(м.м) = ■

з Г Эм ^ э 2м

21 Эх ^ Эх

■; Х3(м. м) =

Эм Э 2м Эх Эх2

где ^1(и. м). Ь2(м. м). £3(м. м) - нелинейные операторы. м( х. г) - прогиб элемента в направлении нормали; и( х. г) - перемещение элемента в продольном направлении; е1 - коэффициент диссипации окружающей среды; Е - модуль Юнга; Н - высота поперечного сечения панели; у - удельный вес материала; g - ускорение свободного падения; кх - геометрическая кривизна балки; г - время; д = д0 8т(йр0 - внешняя нагрузка. д0 - амплитуда внешней нагрузки. О р - частота воздействия

внешней нагрузки.

Безразмерные параметры введены следующим образом:

„а_м_иа_х-г а Л = -; м = —; и = —; х = -; Г = —; Т = —; р НИН

а

Р

'V

Е^ - ^ ^ да - ка

; *=_; д =т^; кх=1. (12)

у р НЕ 1

Черточка над безразмерными параметрами в уравнениях (1.1) для простоты опущена.

К уравнениям системы (1.1) следует присоединить одно из граничных условий

1. Оба конца балки имеют жесткую заделку (х = 0, х = а)

м(0. г) = м(а. г) = и(0. г) = и(а. г) = м\ (0. г) = м\ (а. г) = 0; (1.3)

2. Оба конца балки закреплены шарнирно-неподвижно (х = 0. х = а)

м(0.0 = м(а. 0 = и(0.0 = и(а. г) = м 'хх (0. о = м' 'хх (а. г) = 0; (1.4)

3. Один конец шарнирно-неподвижен (х = 0). а другой имеет жесткую заделку (х = а )

м(0. г) = м(а. г)=и(0. г)=и(а. г) = м\ (а. г) = м 'хх (0. г) = 0; (1.5)

4. Один край имеет жесткую заделку ( х = 0 ). а другой свободен ( х = а )

м(0.г) = м’х(0.г) = и(0.г) = 0; мх(а.г) = N.(а.г) = Ох(а.г) = 0; (1.6)

Мх (а. 0 = N. (а. 0 = Ох (а. 0 ,

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

м( х.0) = м (х.0) = м(х.0) = й(х.0) = 0. (1.7)

Далее приведены численные экспериментальные данные для граничных условий типа (1.5) и начальных условий (1.7) при воздействии внешней знакопеременной нагрузки д = д0 8т(&рг).

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

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

Оы = {0 < х1 < 1. х1 = .1 = 0.^}.

Введем следующие разностные операторы при аппроксимации 0(с2). где с - шаг по пространственной координате [5]:

Л ( \ = ()г +1 — ()г-1 . Л ( \ = ()г +1 — 2()г + ()г-1 .

Л (') ------2С-----■ Л '■(') "-----------------С5-•

2

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

Полученная система обыкновенных дифференциальных уравнений второго порядка (1.8) с соответствующими граничными и начальными условиями, записанными после применения конечноразностной аппроксимации второго порядка точности, методом замены переменных сводим к системе обыкновенных дифференциальных уравнений первого порядка, которую решаем методом Рунге-Кутта четвертого порядка точности - классический метод Рунге-Кутта, или методом Бутчера - метод Рунге-Кутта 6-го порядка [6].

Проверки достоверности результатов численного решения задачи была реализована двумя методами - методом конечных разностей (МКР) и методом конечных элементов (МКЭ). Сравнение сигналов, спектров мощности Фурье и материнских вейвлетов Морле для результатов полученных МКР и МКЭ показало, что они качественно сходятся. На практике МКР показал себя значительно экономичнее МКЭ, поэтому для дальнейших исследований использовался он. При использовании МКР для решения системы уравнений (1.8) применялся классический метод Рунге-Кутта 4-го порядка точности. Метод Бутчера 6-го порядка точности давал те же результаты, при вдвое больших затратах по времени. Подробнее о сходимости результатов можно узнать в [7].

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

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

2. Количественный анализ (карты) характера колебаний криволинейной балки

Рассмотрим криволинейную балку с параметрами п = 120, 1 = 100, кх = 48, е1 = 1.

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

предшествующей работе [7].

В табл. 1 приведены карта режимов колебаний для данной криволинейной балки, а также условные обозначения режимов колебаний криволинейной балки для построенных карт зависимости характера колебаний от управляющих параметров \а>р, q0}.

Кроме режима колебаний балки, нас также интересует ряд количественных характеристик колебательного процесса. Так, существенную ценность имеет карта предельных прогибов. Для её построения у каждого сигнала определяется максимальный по модулю прогиб. Ниже (табл. 2) приведены две карты предельных прогибов. Первая карта содержит абсолютные прогибы, полученные при моделировании, вторая карта является производной от первой.

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

Другой важной характеристикой является предельное растяжение, которому подвергается криволинейная балка. Для построения карты (табл. 3) предельных растяжений в каждом эксперименте определяется максимальное растяжение криволинейной балки.

С

Лх2(^і) кхАх('№і') + Лх(и'()Лх2(и'(),

V

(1.8)

Таблица 1

Обозначении режимов ка-лева-мй

И

и

2 независимые частоты супер поз и-иня лез о ВИСИМЬІИ частот режим не опр^лелен. затухающие колебания.

| лай т и чес кие ко ле Дан и я ■ гармонически* колебания.

| бифуркации

Таблица 2

Карта с асболютными значениями прогибов: ,

.Стоит отметить, что на больших областях прогибы выходят за рамки гипотез. Данная карта носит исключительно теоретический характер и объясняет происхождение карты, приведенной справа

Карта со срезом по прогибу в 5 единиц: min = 0, max = 5 (области, отмеченные белым цветом, выходят за рамки гипотез)

Таблица 3

Карта предельных растяжений: min (є) = 1, max(£) = 1,2339.

Карта с зонами упругопластических деформаций для термически упрочненного алюминия: т1п (г) = 1, тах(г) = 1,002785, а тах = 195 (МПа)

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

ввести относительное удлинение балки --------. то закон Гука в относительных единицах примет вид

а = Ее. где Е - модуль упругости. Зная предельные деформации для конкретных материалов. можно построить карты с зонами упруго-пластических деформаций. В табл. 3 приведены карты предельных растяжений и карта с зонами упругопластических деформаций для термически упрочненного алюминия 1915Т [8] (данный материал был выбран из-за высокого порога пластических деформаций). Черный цвет соответствует минимальному растяжению. белый - максимальному. градации серого отра-

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

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

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

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

Пусть а и 7 і - последовательности прогибов для сигналов а и Ь соответственно. Тогда можно определить расстояние между сигналами по формуле 8 = Егіїк - 7і|, где М- количество точек в сигнале. Результатом будет скалярная величина, характеризующая расхождение или различие рядов между собой, фактически это метрика, заданная на множестве сигналов с одной частотой вынуждающих колебаний. Чтобы не работать с абстрактной величиной (расстояние между рядами), перейдем к формуле 8 = (2?<1Іаі — 7г|)/М. Таким образом, будет получено среднее расстояние (в прогибах) между точками рядов а и Ь. Если расстояние стремится к нулю, то сигнал а сходится к сигналу Ь, если расстояние равно нулю - это один и тот же сигнал.

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

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

В табл. 4 приведены карты расхождения сигналов в абсолютном выражении и со срезами по максимуму в 2, 0.5, и 0.25 единиц прогиба (белый цвет - расхождение выше максимума, градации серого - расхождение ниже максимума, вкрапления серого цвета - результат моделирования не определен).

В табл. 5 для наглядности приведены вместе карта режимов колебаний и карта расхождения сигналов со срезом по максимуму в 0,25 единиц.

На основе результатов, приведенных в табл. 4, 5, можно сделать следующие выводы:

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

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

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

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

хЮ3

Рис. 2. Карта средних прогибов, тіп = — 0.22,тах = 10.69

5. Для эффективного отображения карты подобия срез определяется эмпирически.

Таблица 4

Карта расхождения сигналов (абсолютная), min = 0, max = 9,95

Карта расхождения сигналов (срез по 1 единице среднего прогиба расхождения), min = 0, max = 1.

Карта расхождения сигналов (срез по 0,5 единицы среднего прогиба расхождения), min = 0, max = 0,5.

Карта расхождения сигналов (срез по 0,25 единицы среднего прогиба расхождения), min = 0, max = 0,25.

Таблица 5

Карта режимов колебаний

Карта расхождения сигналов (срез по 0,25 единицы среднего прогиба расхождения)

Ниже приведены результаты исследований влияния кривизны балки кх на режимы колебаний. величины предельных прогибов. величины расхождения сигналов. В табл. 6 для сравнения приведены карты режимов колебаний с параметрами криволинейной балки п = 120. Л = 100. е1 = 1. с кривизной кх = 0. 12. 48 соответственно (при кх = 0 речь идет об обычной балке).

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

В табл. 7 приведены карты предельных прогибов (срез по 5 единицам).

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

В табл. 8 приведены карты расхождения сигналов (срез по 0.25 единиицы среднего прогиба расхождения). Данные карты позволяют обнаружить зоны динамической потери устойчивости.

Таблица 6

Таблица 7

Таблица 8

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

ЛИТЕРАТУРА

1. Karman Th. Festigkeitsprobleme in Maschinebau / Th. Karman // Encykle D. Math. Wiss. 1910. Vol. 4. № 4. P. 311-385.

2. Власов В.З. Общая теория оболочек / В.З. Власов. М.: Гостехиздат, 1949. 784 с.

3. Вольмир А. С. Нелинейная динамика пластинок и оболочек / А. С. Вольмир. М.: Наука, 1972. 432 с.

4. Euler L. Sur la force des colones / L. Euler // Memories de L’Academie de Berlin. 1757. Vol. B. P. 252-282.

5. Самарский А.А. Введение в численные методы / А.А. Самарский. М.: Наука, 1987. 459 с.

6. Grehard E. Hairer Solving ordinary differential equations I: Nonstiff problems? Second edition / E. Hairer Grehard, Wanner Syvert, Paul Norsett. Berlin: Springer Verlag, 1993.

7. Analysis of chaotic vibrations of flexible plates using fast Fourier transforms and wavelets / J. Awrejcewicz, A.V. Krysko, I.E. Kutepov, N.A. Zagniboroda, M.V. Zhigalov, V.A Krysko // Int. J. Str. Stab. Dyn., DOI: 10.1142/S0219455413400051.

8. СНиП 2.03.06-85. М.: Госстрой СССР, 1988.

Работа выполнена при финансовой поддержке РФФИ № НК 12-01-31204.

Загниборода Николай Александрович -

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

Nikolai A. Zagniboroda -

Postgraduate

Department of Applied Mathematics and System Analysis

Yuri Gagarin State Technical University of Saratov

Добриян Виталий Вячеславович -

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

Vitaly V. Dobriyan -

Postgraduate

Department of Applied Mathematics and System Analysis

Yuri Gagarin State Technical University of Saratov

Жигалов Максим Викторович -

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

Maksim V. Zhigalov -

Ph. D., Associated Professor,

Department of Mathematics and Modeling,

Yuri Gagarin State Technical University of Saratov

Крысько Антон Вадимович -

доктор физико-математических наук. профессор. заведующий кафедрой «Высшая математика и механика» Энгельского технологического института (филиала) Саратовского государственного

Anton V. Krysko -

Dr. Sc., Professor,

Head: Department of Higher Mathematics and Mechanics

Engels Institute of Technology

технического университета имени Гагарина Ю.А. Yuri Gagarin State Technical University of Saratov

Крысько Вадим Анатольевич -

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

Vadim A. Krysko -

Dr. Sc., Professor,

Head: Department of Mathematics and Modeling, Yuri Gagarin State Technical University of Saratov

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