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

МАТЕМАТИКА И МЕХАНИКА

УДК 539.3, 534.1

В.А. Крысько, М.В. Жигалов, Т.В. Яковлева, Е.Ю. Крылова, И.В. Папкова МЕТОД УСТАНОВЛЕНИЯ В НЕЛИНЕЙНЫХ ЗАДАЧАХ БАЛОК И ПЛАСТИН С УЧЕТОМ ЛОКАЛЬНОСТИ НАГРУЖЕНИЯ*

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

Метод установления, нелинейная динамика, балки, пластины, оболочки

V.A. Krysko, M.V. Zhigalov, T.V. Yakovleva, E.Y. Krylova, I.V. Papkova THE METHOD OF ESTABLISHING IN NON-LINEAR PROBLEMS OF BEAMS AND PLATES CONSIDERING FOR LOCAL LOADING

The application of the method of establishing problems of algebra and differential equations in partial derivatives are considers. The convergence of a method of establishing in nonlinear problems of beams and plates under local loading was investigated. Results are obtained for different types of local loads.

Method of establishing, nonlinear dynamics, beams, plates, shells

Введение. Идея рассмотрения решений стационарных задач как предела решений нестационарных при возрастании времени впервые использовалась в 30 годах прошлого века А.Н. Тихоновым.

Задача стационарного сверхзвукового обтекания тел газом методом установления была решена С.К. Годуновым, А.В. Забродиным и Г.П. Прохоровым в 1961 году [1]. Решение уравнения Пуассона указанным методом с использованием разностной аппроксимации было предложено в 1947 г. Л. А. Люстерником [2]. Различные модификации этого подхода использовались для решения различных задач Г.И. Марчуком и Ю.А. Кузнецовой, Р.П. Федоренко и многими другими авторами.

Для нелинейных задач теории оболочек этот метод впервые был применен В.И. Феодосьевым [3]. Широкое использование его для различных задач теории оболочек и его модификации было сделано В.А. Крысько и его учениками [4-6].

Изложим коротко идею метода установления, для этого запишем исходную систему нелинейных уравнений в виде

A[x]+[q ] = 0 (1)

где A - нелинейный оператор, который может быть как алгебраическим, так и дифференциальным; [x ] - вектор искомых переменных; [q ] - заданная вектор-функция.

Введем члены, содержащие производную по времени, т.е. искусственно построим динамический процесс:

A[x ]+[q ] = [x] + £[x ] (2)

где x, x - производные по времени; £ - коэффициент диссипации.

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

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

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

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

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

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

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

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

f (л) = 0. (3)

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

c(л + £х) = f (x) (4)

или

- c(X + £к) = f (x). (5)

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

На практике этот алгоритм можно реализовать следующим образом. Задавая произвольным образом начальное приближение л0 и решая дифференциальное уравнение (4), получим решение л1. Далее, задавая начальное значение в виде Л0 = л1 + 5 и решая уравнение (5), можно получить решение л2

. Затем, вновь решая уравнение (4), получаем третий корень уравнения (3) и т.д. В качестве примера приведем решение уравнения [7]

x4 - 12x3 + 47x2 - 60x = 0, (6)

имеющего четыре действительных корня л1 = 0, л2 = 3, л3 = 4, л4 = 5 методом установления (МУ) и

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

Начальное приближение 10 4.5 3.5 1

Метод Ньютона (МН) 5 3 5 5

Метод установления (МУ) 5 4 3 0

Кол-во итераций МН 10 9 12 7

Кол-во итераций МУ 60 87 54 68

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

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

А) балка Бернулли - Эйлера.

Разрешающие уравнения в перемещениях после введения безразмерных параметров по формулам:

А = а/(2к), й = иа/(2^)2, х = х/а, Л = а/(2Л), д = (да4)/(2Л)4Е (7)

имеют вид

дх

2 + Ц (м, м) = 0.

-11 Ь2 (w, м) + Ц (и, м)- -1^ і + д = 0-

здесь w(х, ,)- прогиб балки; и(х, ,) - перемещение вдоль оси ОХ; /^(и, А) =

(8)

д 2и дм ди д2 м

дх2 дх дх дх2

, £ - коэффициент диссипации; д = д(х, г) - попереч-

^ А) = I!х?(!т) • ^ А) аг дх

ная нагрузка; Е - модуль Юнга; V - коэффициент Пуассона; р - плотность материала; б - удельный вес материала балки; g - ускорение свободного падения.

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

и(0,ґ) = и(1,ґ) = 0; м(0,ґ) = м(1,ґ) = 0;

д 2 м(0,ґ) д 2 м(1, ґ)

дх

2

дх

2

= 0.

(9)

Вводя динамические члены, получим следующую систему уравнений:

д2 и д2

дґ2

м

д2

дх2 дм

+ Ь3 (м, м),

1 д4 м

(10)

где безразмерные параметры для временных членов введены по формулам

г = г/т, т = а/с, с = 7Е^б, £ = ш/с. (11)

Для задачи (10) кроме граничных условий (9) необходимо ввести начальные условия:

и,=0 = А0 ; й\г=0 = 0; А,=0 = А0 ; А|г=о = 0 (12)

Заменим далее дифференциальные операторы по пространственной переменной х, разностными

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

Рассмотрены следующие типы приложения нагрузки (рис. 1): а) равномерная - нагрузка задана по всей длине балки; б) локальная симметричная - нагрузка распределена на 4 отрезках разбиения балки (1/10 часть длины) симметрично расположенных, относительно середины; в) локальная несимметричная - нагрузка распределена на 4 отрезках разбиения балки, расположенных слева от центра балки (8-12 отрезки).

и

и

а)равномерная

в) локальная несимметричная

б) локальная симметричная Рис. 1. Типы нагружения

В результате численных экспериментов, было выявлено, что не способ приложения нагрузки (симметрично / несимметрично относительно центра балки) ни ширина полосы приложения нагрузки не влияют на выбор Єкр. Во всех случаях єкр = 1.5.

Рис. 2. График а(,) для различных значения коэффициента диссипации £

Полученные численные результаты проиллюстрированы на рис. 2 для случая, когда нагрузка приложена несимметрично относительно центра балки к 4 отрезкам. Как видно из графика, при £ < £кр

процесс устанавливается на более далеких временах - , > 35, при £ < £кр процесс установления также сходится медленно.

Рассмотрим графики нагрузка - прогиб для различного вида нагрузок при £кр = 1.5 .

Рис. 3

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

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

Б) Г ибкая пластинка.

Рассмотрим пластинку произвольного плана, из изотропного материала, с учетом геометрической нелинейности по модели Т. фон Кармана [8]. Введем в систему дифференциальных уравнений [9] динамические члены в результате приходим к следующей системе уравнений:

+ Є = д-А2 м + Ь(м, Е), А2 Е +1 ь(м, м) = 0. 2

(13)

2 1 д4(-) п2 д4(0 д4(-) т( \ д2и д2V д2и д2V д2V д2и

где А О = тт+ А + 2-, 2^ 2 , 1(и’v)^ТТЧТ - \ а а а + ТГЧТ Безразмер-А2 дх ду дх ду дх ду дх ду дх ду дх ду

ные величины введены по формулам:

м = Км; Е = ЕН2Е ; г = г0г; є = є/ т ; х = ах; у = Ьу ;

- 12(1 -V )ЕН аЬ

д=д-------^2—; Т=~тл

а Ь Н \

р . А=а (14)

Ея ’ Ь ’

где а, Ь - максимальные размеры пластинки в плане по х и у соответственно, к - толщина, g - ускорение силы тяжести, р = }к , У - объемный вес материала, V - коэффициент Пуассона.

В качестве граничных условий к системе (13) необходимо присоединить граничное условие:

д 2а р д2 ^ п

А=э^ = ? = э^=0' (15)

где п - нормаль к границе области.

К уравнениям (13) помимо краевых условий (15) присоединяем начальные условия

4=0 = А0 ,=0 = 0 (16)

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

Аппроксимация уравнений по пространственным координатам проводится МКЭ на основе процедуры метода Бубнова - Г алеркина. В связи с тем, что область не является прямоугольной, используются шестиузловые треугольные элементы [10].

В качестве примера рассмотрим применение М.У. для пластинки произвольного плана (рис. 4а)), под действием локальной нагрузки (место приложения обозначено другим цветом). Здесь И=1.0; г=0.5;а=1.0; Ь=0.5; с=0.366; ё=0.368.

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

Рассмотрим результаты для точки (0; 0), полученные по методу установления.

Рис. 5. Графики нагрузка-прогиб и прогиб-время

В левой части рис. 5 приведен график нагрузка-прогиб, в правой части - график время-прогиб. Как видно из правого графика, при величине прогиба, превышающей 1, график не является «убывающей синусоидой», как это происходит при равномерно распределенной нагрузке.

В) Пакет пластина - балка, соединенные через краевые условия

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

r ~ Wl - hk - w2)] (17)

Верхняя пластинка описывается линейным уравнением Жермен-Лагранжа [11], а балка - линейным уравнением Эйлера-Бернулли.

Для учета контакта между пластинкой и балкой, используется винклерова связь между обжатием и контактным давлением между двумя балками [12]:

1

¥= ^ [1 + ^П(м1 - Нк - м2 )]

Функция ¥ = 1, если w1 > м2 + Нк то есть имеется контакт между пластиной и балкой, иначе

¥ = 0.

Отметим два, на наш взгляд, существенных обстоятельства:

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

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

Вводя динамические члены указанным выше приемом, получаем следующую систему уравнений:

1

W9 + єА&2 = q(t) +

2 2 12(1 -v2)

f 1 д4wl +л2 д4wl + 2- — 4Wl ^

Л2 дx4 ду4 —x ду2

У&2 + ЄА 2 = '

2(1 -V ) V Л2 —x4 —y 4 —x 2 —y 2

1 д4 w

12 —x

v J J

4

w

+ K(Wl - W2 - hk )^

+ K(Wl - W2 - hk )^,

(18)

здесь w1• ^2 - функции прогибов пластины и балки соответственно, К - коэффициент жесткости транс -версального обжатия структуры в зоне контакта, кк - зазор между пластиной и балкой.

Система уравнений (18) приведена к безразмерному виду. Для простоты записи черточки над безразмерными параметрами в системе уравнений (18) опущены.

К уравнениям (18) следует присоединить граничные условия:

—2 д2

wt = 0; —Wl = 0; при x = 0;1; wt = 0; —Wt = 0; при y = 0;1 ox dy

—2 a2

w2 = 0; —2^ = 0; при x = 0;1; w2 = 0; —^ = 0 при y = 0;1, Ox —у

(19)

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

Wm |t=0 = 0; Wm |t=0 = 0^ m=1, 2. (20)

Полученные системы нелинейных дифференциальных уравнений в частных производных (18) сводятся к системе обыкновенных дифференциальных уравнений второго порядка методом Бубнова-Галеркина в высших приближениях в форме В.З. Власова по пространственным переменным x и у , которая решается методом Рунге-Кутта четвертого порядка точности.

Рассмотрим характер колебаний пластины и балки под действием поперечной равномерно распределенной нагрузки на пластину q = q0 = const.

На рис. 7 приведены зависимости q(w ) для квадратной в плане пластины (Л = 1) с коэффициентом диссипации среды £ = 5 .

Графики, построенные для величины зазора hk = 0.2 (сплошная линия) и hk = 0.1 (пунктирная линия), полностью совпали до момента контакта пластины и балки. Следует отметить, что при прогибе w = hk = 0.1 происходит соприкосновение пластины и балки и система становится жестче. Расхождение значения происходит при наличии касания пластины и балки.

Рис. 7. Графики нагрузка-прогиб д(^уст) и прогиб-время ) для пластинки

Также приведены графики зависимостей w(0,5;0,5;t) и w(0,5;t) для пластины при различных значениях нагрузки и величине зазора Нк = 0,2. При маленьких нагрузках колебания пластины с течением

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

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

Для исследования колебательных процессов в механических системах (балки, пластины и оболочки) под действием переменных нагрузок необходимо использовать другой аппарат, а именно методы нелинейной динамики, которые включают в себя преобразование Фурье и вейвлет-анализ [13]. Ниже приведем пример применения этих методов при исследовании сценария перехода в хаос при локальной несимметричной нагрузке для геометрически нелинейной балки Бернулли-Эйлера.

Локальная нагрузка, приложенная к 4 отрезкам, симметрично расположенным относительно центра балки (балка разбита на 40 отрезков), имеет вид

д = А • $,тЩрг), (21)

где Щр = 5,3 - частота возбуждения, А - амплитуда.

При увеличении амплитуды до 56250, на спектре мощности наблюдается появление линейно независимой от Щр частоты С01 = 1,54 и линейно зависимой частоты Щ = Щр — Щ = 3,76. С ростом амплитуды внешней нагрузки до д0 = 56710 появляется пара Щ = 2,27 со4 = 3,07 и при д0 = 56750 пара Щ = 0,75 Щ = 4,55, связанные следующим образом. Щ+1 = Щр — Щг,г = 3,5. (Частоты с нечетным индексом линейно не зависят от частоты возбуждения Щр) Частоты Щ и С02 немного меняют свои

значения, вышеуказанная связь между ними остается неизменной (их значения 1,5 и 3,8, соответственно). Наблюдается усемерение периода колебаний системы. Разница между соседними частотами составляет 0,75.

При А = 58000 наблюдаем уодиннадцатирение периода колебаний системы. Описанные выше частоты немного «плывут» (щ = 1,44, Щ = 3,86 и Щ = 2,41, щ4 = 2,89 и Щ = 0,96, Щ = 4,35) к

ним присоединяется еще две пары Щ = 0,48, Щ = 4,82 и Щ = 1,93, Щ0 = 3,37 вида

Щ+1 =Щр —Щ, г = 7Я

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

Спектр мощности Фурье

2-Р вейвлет спектр Морле

500

56250

СО = 5,3 щ = 1,54, а>2 = 3,76

56750

-ю

-15

5 ' 1

сор _

а>2

(Оі

- 1 СОз 0>4 1 (&б ~

1 і ®

со = 5,3 щ = 1,5, со2 =3,8 щ = 2,27,со4 = 3,07 щ = 0,75, со6 = 4,55

63000

Авторами были построены карты типа колебаний размером 300x300 точек для управляющих параметров {А,Щр} для различных случаев приложения нагрузки. Амплитуда внешнего воздействия менялась при этом на интервале (0; 60000), частота на интервале (3,45; 10,3) для всех вариантов.

5353484848484853232348535323485353484848

Подобные карты позволяют изучить все многообразие поведения конструкции и выявить области с благоприятными типами колебаний (гармонические, затухающие, квазипериодические), а также с колебаниями, приводящими к негативным последствиям, вплоть до разрушения (хаос и бифуркации). На основании построенных карт, можно путем управления параметрами нагрузки {А, С0р } удерживать режимы колебаний балки в благоприятных зонах.

Идентификация типа колебаний при построении данных карт {А, С0р } для каждого сигнала )

проводилась с помощью анализа спектра мощности Б (®р ) , фазового портрета, сечения Пуанкаре и показателей Ляпунова. Таким образом, при построении таких карт необходимо решить и проанализировать 90000 задач.

еоооо

SS000 50000 4S000 40000 35000 3000 0 25000 20000 1S000 10000 5000

1 I

Periodic vibrations with cop Bifurcations

Periodic vibrations with cjpj 2 Chaotic vibrations

Periodic vibrations with (apj 3 Quasi-periodic vibrations

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

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

G)p < 6.9), лишь при небольших значениях амплитуды возбуждающей силы 0,5 < A < 10000 наблюдаются неустойчивое поведение системы (появляются зоны независимых частот - Quasi-periodic vibrations). В остальном колебания системы носят гармонический характер, для всех значений амплитуды внешнего воздействия. С увеличением частоты, появляются зоны независимых частот, и зоны бифуркаций Андронова - Хопфа. С ростом амплитуды в данном диапазоне частот возбуждающей силы обнаруживаются зоны хаотических колебаний, но они заметно меньше хаотических зон при нагружении балки по всей длине.

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

ЛИТЕРАТУРА

1. Годунов С. К. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной I С. К. Годунов, А. В. Забродин, Г. П. Прокопов II Ж. вычисл. матем. и матем. физ., 1:б (19б1), 1020-1050.

2. Люстерник Л. А. Замечания к численному решению краевых задач уравнения Лапласа и вычислению собственяых значений методом сеток I Л. А. Люстерник II Тр. Матем. ин-та им. В. А. Стеклова, 20 (1947), 49-б4.

3. Феодосьев В.И. Об одном способе решения нелинейных задач устойчивости деформируемых систем I В.И. Феодосьев II ПММ, 19б3. Т.27. Вып. 2. С. 2б5-274.

4. Крысько В.А. Определение неустойчивых решений при расчете пластин и оболочек I В.А.Крысько, С.А.Комаров II Труды XVII международной конференции по теории оболочек и пластин. Казан. гос. ун-т, 199б. Т. 2. С. 19-24.

5. Крысько В.А. Выпучивание гибких пластин под действием продольных и поперечных нагрузок I В.А.Крысько, С.А.Комаров, Н.В. Егурнов II Прикл. Механика. 199б. Т.32. № 9. С.80-87.

6. Крысько В.А. О применении метода В.И. Феодосьева к решению нелинейных задач статики пологих оболочек I В.А. Крысько, Н.В. Егурнов, А.А. Сопенко I СПИ. Саратов, 1985. 8 с. Деп. в ВИНИТИ, № 8174-В85.

7. Деннис Дж. Численные методы безусловной оптимизации и решения нелинейных уравнений I Дж. Деннис, Р. Шнабель. М.: Мир, 1988. 440 с.

8. Karman Th. Festigkeitsprobleme in Maschinenbau I Th. Karman II Encykle. D. Math. Wiss. 1910. Vol. 4, №4. P. 311-385

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

с.

10. Сагерлинд Л. Применение метода конечных элементов I Л. Сагерлинд. М.: Мир, 1979. 393 с.

11. Тимошенко С.П. Теория упругости I С.П. Тимошенко, Дж. Гудьер. М.: Наука, 1975. 57б с.

12. Кантор Б.Я. Метод решения контактных задач нелинейной теории оболочек I Б.Я. Кантор, Т.Л. Богатыренко II Докл. АН УССР. Сер. А. 198б. №1. С. 18-21.

13. Awrejcewicz J. Chaos in Structural Mechanics I J. Awrejcewicz, V.A. Krysko. Berlin, London: Springer, 2008. 424 p.

* - Работа выполнена при поддержке Министерства образования и науки РФ в рамках ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы проект 2012-1.4-12-000-1004-00б.

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

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

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

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

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

Vadim A. Krysko -

Dr. Sc., Professor,

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

Maksim V. Zhigalov -

Ph. D., Associated Professor,

Department of Mathematics and Modeling, Gagarin Saratov State Technical University

Tatyana V. Yakovleva -

Postgraduate,

Department of Mathematics and Modeling, Gagarin Saratov State Technical University

Крылова Екатерина Юрьевна -

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

Папкова Ирина Владиславовна -

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

Ekaterina Yu. Krylova -

Postgraduate,

Department of Mathematics and Modeling, Gagarin Saratov State Technical University

Irina V. Papkova -

Ph. D., Associated Professor,

Department of Mathematics and Modeling, Gagarin Saratov State Technical University

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