МЕХАНИКА ТВЕРДОГО ТЕЛА < 6 • 2008
УДК 539.3;534.1
© 2008 г. А.С. ДЕСЯТОВА, М.В. ЖИГАЛОВ, В.А. КРЫСЬКО, О.А. САЛТЫКОВА
ДИССИПАТИВНАЯ ДИНАМИКА ГЕОМЕТРИЧЕСКИ НЕЛИНЕЙНЫХ БАЛОК БЕРНУЛЛИ-ЭЙЛЕРА
В статье рассмотрены нелинейные диссипативные колебания балки Бер-нулли-Эйлера. Выявлено, что при действии поперечной знакопеременной нагрузки колебания могут переходить в хаотические. Исследован сценарий перехода гармонических колебаний в хаотические - сценарий Фейгенбаума. Определена константа Фейгенбаума. В работе уделено значительное внимание достоверности получаемых результатов. Для этого используются два метода - конечных разностей O(h2) и конечных элементов в форме Бубнова-Галеркина, в работе исследована сходимость этих методов.
1. Введение. Исследованию сложных колебаний распределенных систем в настоящее время уделяется должное внимание, но имеющиеся результаты получены для математических моделей с одной степенью свободы. Так, например, зачастую колебания распределенных систем заменяются исследованием по одномодовому приближению в методе Бубнова-Галеркина [1]. В работах [2, 3] тем же методом исследованы нелинейные колебания балок Бернулли-Эйлера в одномодовом приближении с учетом геометрической нелинейности. Результаты этих работ следует рассматривать как качественные, так как увеличение количества мод зачастую приводит к существенно другим режимам колебаний. В работе [4] представлены исследования глобальных бифуркаций и хаотической динамики в нелинейных неплоских колебаниях консольной балки под действием осевых гармонических возбуждений и поперечных возбуждений на свободном конце балки. Применение метода конечных элементов для исследования вынужденных нелинейных колебаний балки рассмотрено в статье [5]. Изучению явления хаоса в задачах колебаний пластин и оболочек с геометрической нелинейностью посвящены работы [6-14]. В данной работе изучаются колебания геометрически нелинейной балки Бернулли-Эйлера как системы с множеством степеней свободы.
2. Постановка задачи. Объектом исследования является однослойная балка, представляющая собой двумерную область пространства R2 с декартовой системой координат, введенной следующим образом: в теле балки фиксируется линия приведения, называемая срединной линией г = 0, ось OX направлена слева направо вдоль срединной линии, ось OZ - вниз, перпендикулярно OX. В указанной системе координат балка, как двумерная область Q определяется следующим образом: Q = {x е [0, a]; -h < г < h}, 0 < t < го. Здесь и в дальнейшем будем использовать обозначения: 2h - высота, a - длина балки.
Математическая модель строится на основе гипотезы Бернулли-Эйлера и с учетом нелинейной зависимости между деформациями и перемещениями в форме Кармана [15]. Система дифференциальных уравнений в перемещениях, описывающих движения балки с учетом диссипации энергии, получена из системы уравнения для нелинейной пластинки [15] в случае одной пространственной координаты:
Е2Н[^ + Ь3(*, *)}-гн&^-гмМ = 0
[Эх2 ] 8 дI2 8д г
Е2Н[Ь,(х, *) + Ь2(*) - ^^ 1 + 9 -2й*^-2Ле,*^ = 0 (2.1)
[ 12 дх ) 8 дг 8 дг
, , , д2мд* дмд2* , , , 3д2*^д*^2 , , , д2*д*
Ь (»• *> + дхдх?- Ь2(*• *) = 3д?IжЬз(*•*) -где *(х, г) - прогиб балки, х(х, г) - перемещение вдоль оси ОХ, е1; е2 - коэффициенты диссипации, 9 = 9(х, г) - поперечная нагрузка, Е - модуль Юнга, V - коэффициент Пуассона, р, у - плотность и объемный вес материала балки, 8 - ускорение свободного падения. Введены безразмерные переменные
4
ха - х - а а
* - -Л-, х - -х - -, Л - —, 9 - 9-4-
2Н (2Н)2 а 2Н (2 Н )4 Е (2.2)
г - г/т, т - а/с, с - л/е8/у, е - еа/с
С учетом (2.2) система (2.1) записывается в виде:
д2х/дх2 + Ь3(*) - д2х/дг2- е2дх/дг - 0
1 ( ) г( ) 1 д4* 1 д2 * д * „ (2.3)
где черточки над безразмерными параметрами для простоты опущены.
К уравнениям (2.3) следует присоединить один из следующих вариантов краевых условий:
*(0, г) - *(а• г) - х(0, г) - х(а• г) - д*(0, г)/дх - д*(а• г)/дх - 0 (2.4)
*(0, г) - *(а, г) - х(0, г) - х(а, г) - д2■*(0, г)/дх2 - д2*(а, г)/дх2 - 0 (2.5)
22
*(0, г) - х(0, г) - *(а, г) - х(а, г) - д*(0, г)/дх - д2*(а, г)/дх - 0 (2.6)
*(0, г) - д2*(0, г)/дх2 - х(0, г) - 0; Мх(а, г) - (а, г) - 0х(а, г) - 0 (2.7)
Имеем также начальные условия:
*(х,г)|г-0 - д*(х,г)/дг|г-0 - х(х,г)|г-0 - дх(х,г)/дг|г-0 - 0 (2.8)
3. Методы решения. Дискретизация функций *(х, г), х(х, г) проведена с помощью метода конечных элементов (МКЭ) в форме метода Бубнова-Галеркина и метода конечных разностей (МКР) с погрешностью О(Н2). Таким образом, уравнение в частных производных (2.3) сводится к решению системы обыкновенных дифференциальных уравнений по времени. Использование МКР и МКЭ позволило решить вопрос о достоверности получаемых результатов в задачах сложных колебаний распределенных систем, который является очень важным. Рассмотрим эти два подхода более подробно.
Согласно теории МКЭ, для построения балочного элемента необходимо ввести функции формы. На элементе рассматривается четыре степени свободы (*1, *2, 01, 02) и
5 Механика твердого тела, № 6
129
для аппроксимации используется кубический полином w = а1 + а2х + а3х2 + а4х3 и -й*/йх = -(а2 + 2а3х + 3а4х2). После определения коэффициентов а^О, а2(0, а3(0, а4(0 в зависимости от ^(0, *2(0, 0Х(?), 02(О для аппроксимации функции *(х, I) получим:
* = ]{ W }
[^] = (1-3^ + 2- № -1 )2; 3^2-2£3; - I^2- £))
где [N^1 - матрица формы; {W} = 0:, *2, 02)г - матрица узловых смещений; £ = хЦ -безразмерная величина (локальная координата).
Аппроксимация перемещения и(х, 0 выглядит следующим образом: и = [N„1^}, где [N„1 = (1 - £} - матрица формы; {и} = (и1, и2)Т - матрица узловых значений.
Применяя процедуру метода Бубнова-Галеркина с учетом введенных аппроксимаций получаем следующие разрешающие уравнения МКЭ:
Г МД У ] + С1[ У] + К 1[ Ж] = д, и)
] .. . (3.1)
IМ2[и] + С2[и] + к2[и] = р, Ж)
где М,, С, К, - матрицы масс, демпфирования и жесткости соответственно. К системе (3.1) присоединено одно из краевых условий (2.4)-(2.7) и начальное условие (2.8).
Заменим далее дифференциальные операторы по пространственной переменной х, разностными операторами для функций *(х, 0, „(х, 0 с помощью МКР с аппроксимацией 0(Н2). Таким образом, система (2.3) сводится к системе обыкновенных дифференциальных уравнений:
* ,+ £1 * г = ч Н(„ ,0. и + £2„г = ¿2.Н(„ ,0. 1 = 0,..., Л (3.2)
где л - число разбиений по пространственным координатам, а разностные операторы ¿1 н, ¿2 н имеют стандартный вид. К системе (3.2), так же как и для МКЭ присоединены одно из краевых условий (2.4)-(2.7) и начальное условие (2.8), которые записаны в разностном виде.
Системы обыкновенных дифференциальных уравнений, полученных с помощью аппроксимации по пространственной переменной с помощью МКЭ или МКР, решаются методом Рунге-Кутта четвертого порядка точности.
4. Анализ результатов и примеры. Для исследования поведения нелинейных балок под действием знакопеременной поперечной нагрузки д = д08т(юрО разработан пакет программ, позволяющий строить карты характера колебаний в зависимости от управляющих параметров {д0, юр}. Для построения указанной выше карты решается задача динамики, строится и анализируется спектр мощности и ляпуновские показатели для каждого набора значений {д0, юр}. Алгоритм позволяет выделять на картах зоны гармонических колебаний, зоны бифуркаций удвоения периода, зоны колебания балки на независимых частотах и зоны хаоса.
На основе полученных результатов, более подробный анализ которых приведен ниже, сделан вывод о зависимости сходимости решения по сетке для МКР от вида граничных условий. Поэтому для апробации процедуры исследования колебаний балки с помощью МКЭ и МКР рассмотрена задача с краевыми условиями (2.7), дающими наихудшую сходимость. Приняты следующие данные: коэффициенты демпфирования равны £1 = 1; £2 = 0; юр = 5.1 - частота вынуждающей силы, д0 = 55000 - амплитуда вынуждающей силы; отношение длины к толщине равно X = 50.
Ниже приведен график зависимости прогиба от времени *(() в середине балки при разбиениях 40 (МКР - линия 1, МКЭ - на графике линия 2) и 80 (МКР - линия 3, МКЭ -
Фиг. 1
линия 4). Сигнал отражает временной интервал t е [1002; 1005] для граничного условия (2.7) и соответствует хаотическим колебаниям системы (фиг. 1).
Как видно из графиков, при увеличении количества разбиений в два раза значения прогиба, полученные по МКР и МКЭ, отличаются на 15% и 16% соответственно. При этом время вычислений возрастает в несколько раз. Исходя из этого, при дальнейшем исследовании колебаний балки с граничным условием (2.7) использовалось разбиение, равное 40. Кроме этого, для получения карты характера колебаний в зависимости от управляющих параметров размерностью 300 х 300 необходимо просчитать 9 ■ 104 вариантов. На персональном компьютере с процессором Celeron 1700 с использованием метода конечных разностей, при количестве разбиений по пространственной координате n = 40 расчет занимает примерно 70 суток, при n = 80 - около 140 суток, а при n = 160 - около 280 суток. В случае метода конечных элементов время расчета при тех же количествах разбиений возрастает примерно в 2 раза для n = 40 и в 2.7 раза для n = 80. Сравнение результатов, полученных по МКР и МКЭ показывает, что при одних и тех же параметрах нагрузки и краевых условиях значения прогиба близки, спектры мощности, полученные с помощью МКР и МКЭ при одних и тех же разбиениях практически совпадают, т.е. имеется сходимость методов в среднем, что говорит о достоверности получаемых результатов.
Данные, полученные для граничных условий (2.4)-(2.6) с помощью метода конечных разностей и метода конечных элементов, также подтверждают сделанные выше выводы. На основании этого и учитывая увеличение временных затрат при использовании МКЭ основные расчеты проводились с помощью МКР.
Далее приведены результаты исследования колебаний балки Бернулли-Эйлера методом конечных разностей для краевых условий вида (2.4), (2.5), (2.6), (2.7) и начальных условий (2.8). В результате численного эксперимента для симметричных граничных условий (2.4), при X = 50, были построены карты зависимости х
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.