CC BY
65
УДК 539.3
РАСЧЕТ ПЛАСТИН И ПОЛОГИХ ОБОЛОЧЕК МЕТОДАМИ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ
С.И. Трушин
Кафедра информатики и прикладной математики Московского государственного строительного университета 129337 Москва, Ярославское шоссе, 26
Представлены алгоритмы расчета пластин и пологих оболочек, основанные на вариационно-разностном подходе в сочетании с методами прямой минимизации дискретного аналога исходного функционала. Используются процедуры метода сопряженных градиентов Флетчера-Ривса и Полака-Рибьера, квазиньютоновского метода Дэви-дона-Флетчера-Пауэлла. С помощью разработанных численных методик решен ряд тестовых задач. Даны примеры расчета квадратной пластинки, пологой сферической и цилиндрической оболочек.
Задача статического расчета пластин и оболочек формулируется как задача отыскания минимума функционала Лагранжа, представляющего собой полную потенциальную энергию системы относительно функций перемещений. Исходные геометрические соотношения линейной теории пологих оболочек с учетом деформаций поперечного сдвига записываются в известной форме:
ди м> ду м' ду ди
в,,——I—; ; вп— +-
дх Я1 ду дх ду
30, 30, 30, 30, дм> дм>
к.к„=—к,,= + ; еп=—+0.; е,,=—+0,.
" дх 22 ду 12 дх ду 13 дх 1 23 ду 2
(1)
где и и V - тангенциальные составляющие перемещения; \\> - нормальная составляющая; 01 и 02 - углы поворота поперечных сечений оболочки; Я\, Яг- радиусы кривизн в направлении координатных осей.
В задачах теории оболочек исходный трехмерный функционал Лагранжа после интегрирования по толщине преобразуется к двумерному:
J(u) - есЬс1(Ьс2 - иск,(Ьс2, (2)
где 1) - матрица упругости; е = ( ви е22 Кц к22 К12 £¡3 )Т - вектор, компонентами которого являются составляющие тензора деформаций (1 );# = (#/ (¡2 п?! ^ )т - вектор
внешней нагрузки, компоненты которого имеют направления, соответствующие компонентам вектора перемещений; и - вектор функций перемещений.
Предположим, что на область О, занимаемую конструкцией, наложена сетка
х, = х0 + /Д; у] = Уо+^'г /' = 0,1,...,от,; ] = ОД,...,от2; к = шах(/г, ,1ц).
Рассмотрим ячейку области £2 с вершинами в точках (хь у}), (*¡+1, У}), (хь У)+\), 0*ч+ь _У)+1) и обозначим значения сеточных функций в этих точках и°, И1, и , и . Функционал энергии (2) может быть представлен в виде суммы интегралов по ячейкам сетки
У = ^ JІJ . Для единственности решения дискретной задачи необходимо потребовать,
'.У
чтобы аппроксимационная схема сохраняла свойство строгой выпуклости исходного функционала. Кроме того, в целях построения наиболее экономичных схем, желательно, чтобы число подсчетов подынтегральной функции было минимальным.
Для вычисления интеграла ./у по ячейке используется его дискретный аналог .Ру. Наиболее простая аппроксимационная схема имеет следующий вид:
=*Л/
(3)
йЧн' + йЧи3 и'+и1-и°-и2 (<*3+и2-и1-и°)
’ V ' ~ 2/г, ’ 2/г,
где Й1 , /?2 - линейные размеры ячейки в направлении координатных осей; / - подынтегральная функция; иш - значения искомых функций перемещений в узлах. Однако схема (3) обеспечивает решение задачи только в том случае, когда функция И задана своими значениями на границе области О., что выполняется не всегда.
В результате дискретизации исходная вариационная задача сводится к задаче минимизации функции /?(«) = ^ Ру , где и—(Ы] Мг ••• Мы)Г - вектор узловых перемещений. Данная
и
функция является дискретным аналогом полной потенциальной энергии системы (2) и представляет собой сумму потенциальной энергии деформации И^(и) и потенциала внешних силЛ(и):
Г(и) = Щи) + А(и). (4)
Для решения задачи минимизации как квадратичной, так и неквадратичной функции Р\и), представленной в виде (4), весьма эффективными могут быть методы нелинейного программирования, в частности метод Флетчера-Ривса, метод Полака-Рибьера и метод Дэ-видона-Флетчера-Пауэлла. Вопросам численной реализации методов нелинейного программирования и построения эффективных алгоритмов посвящены, в частности, работы [1,2,3].
В методе сопряженных градиентов Флетчера-Ривса на каждом шаге по ведущему параметру строится последовательность направлений поиска 5 , являющихся линейными комбинациями текущего градиента и предыдущих направлений поиска. При этом весовые коэффициенты выбираются из условия сопряженности направлений. Схема этого метода имеет вид:
ик+1=ик+Х5к.
Алгоритм метода Флетчера-Ривса включает в себя следующие этапы:
1. Выбирается начальная точка (начальное решение) и .
2. В и° вычисляется 5° =
3. На к-м шаге с помощью одномерного поиска по А>0 в направлении 5^ находится минимум функции Р(ик +Хк 5*). Для одномерного поиска используется алгоритм Девиса, Свена и Кэмпи [1]. Согласно этому алгоритму при одномерном поиске проводятся возрастающие по величине шаги до тех пор, пока не будет пройден минимум, а затем выполняется одноразовая квадратичная интерполяция.
4. Полагается I/**1 = ик + Хк 5*.
* к м
5. Проверяется условие сходимости. Алгоритм заканчивается, когда || £ ||<8 , где £ - малое число.
6. Проверяется, выполнен ли поиск во всех п+1 направлениях, где п - размерность задачи. Если цикл поиска по п+1 направлениям закончен, то осуществляется переход на шаг 2.
7. Вычисляется ¥Р(ик+1^у
8. Направление поиска 5 определяется из соотношения
5*+' =
1 * (V р(ик))^ ¿У) '
9. Полагается к=к+1 и осуществляется переход на шаг 3.
Существенным достоинством метода Флетчера-Ривса является то, что здесь не требуется обращать матрицу и не требуется большой объем машинной памяти по сравнению с методами, использующими матрицы размера пхп.
Алгоритм метода Полака-Рибьера аналогичен алгоритму метода Флетчера-Ривса, отличие заключается лишь в том, что на 8-м шаге используется несколько иная формула для оп-
ск+1
ределения 5
5к+1
(ЧР(ик)) ^(«*)
Метод Дэвидона-Флетчера-Пауэлла относится к группе квазиньютоновских методов. В нем матрица, обратная к матрице Гессе, аппроксимируется с помощью первых производных. Схема этого метода имеет вид:
ик+,=ик-Хк |/(и* )УР(«*),
где матрица Л(и*) представляет собой аппроксимацию матрицы [V2И*)]'1 . Для квадратичной целевой функции Ди) матрица Г) перевычисляется таким образом, что после п ша-
2 | 0 гов она становится равной матрице [V . На первом шаге полагается Т] = /. В этом
случае исходное направление минимизации совпадает с направлением наискорейшего спуска.
Алгоритм метода Дэвидона-Флетчера-Пауэлла состоит из следующих этапов:
1. В точке и° вычисляется направление
5° = -П°УР(и0) = -Г7Р{и°У
2. На к-м шаге с помощью одномерного поиска по X* находится ттР(ик — X)) и определяется величина X* .
3. Вычисляется Аик = - Х*т\(ик)^ Р(и^) .
4. Определяется ик+1 = ик +Аик .
5. Проверяются условия:
і М«‘)
<є7;
Аик
, <Е,, где Єї , Єг - заданные малые числа.
и
В случае выполнения записанных выше условий решение заканчивается.
6. Вычисляется вектор ЧР(ик+і).
7. Определяются вектор и матрицы:
Л пг/ ь/\ А Лик{Лик)Г к пк^кШ)(гікї
^8к= УР\и )~УР\и ) * = 77;
[Аик) V (V) і/к^к
цк+І =і1к+Ак+Вк
Осуществляется переход на шаг 2.
Достаточно универсальный алгоритм вычисления компонент вектора градиента Р(и) и матрицы вторых производных (матрицы Гессе) дискретного аналога Щи) потенциальной энергии деформации строится на основе приближенных разностных формул:
+ Ни + Ье )]; г=1, 2,..., п\ _/-1, 2,..., п, (5) ди. 26
где п - размерность задачи; 5 = (10 3 -г 10 5) II и II - некоторое малое отклонение; || и 1 -сферическая норма вектора перемещений; Є,, е, - единичные векторы.
В качестве тестовой задачи взято решение задачи об изгибе квадратной пластинки, подверженной действию поперечной равномерно распределенной нагрузки. Решение представлено в монографии С.П.Тимошенко и С.Войновского-Кригера [4], где оно получено методом Мориса Леви. Принята следующая расчетная схема: нагрузка - равномерно распреде-
ленная; граничные условия - свободное опирание по контуру. Максимальный прогиб в центре пластинки, согласно [4], определяется по формуле: м>тах = 0,00406да4 / £), где 0 -интенсивность равномерно распределенной нагрузки; а - размер пластинки в плане; £> = £й3 /(і2(і-V2)) - цилиндрическая жесткость. В расчетах взяты следующие исходные данные: Е=2,1 -1011 Па; у=0,3; й=0,01 м; 0=104 Па.
Данная тестовая задача об изгибе квадратной пластинки решалась, кроме того, с помощью ВК «Лира». Ввиду симметрии рассматривалась четверть пластинки, которая разбивалась на 36 конечных элементов (6x6). Все исходные данные принимались такими же, как и для приведенного выше аналитического решения.
Полученные решения сравнивались с результатами расчетов по методу сопряженных градиентов Флетчера-Ривса. Значения прогиба в центре пластинки даны в табл. 1. Как видно из таблицы максимальная относительная погрешность решения не превышает 1,6% даже при грубой сетке 2x2.
Таблица 1
Прогиб в центре пластины, м
Метод Флетчера-Ривса ВК «Лира» (6x6 КЭ) Аналитическое решение [4]
Сетка 2x2 Сетка 4x4 Сетка 8x8
0,2145-10'2 0,2121-Ю'2 0,2115-Ю"2 0,2128-10'2 0,2111-Ю'2
В качестве другой тестовой задачи взята задача расчета пологой сферической оболочки на прямоугольном плане, находящейся под действием поперечной равномерно распределенной нагрузки (рис.1). Граничные условия соответствуют опиранию по контуру на жесткие в плоскости и абсолютно гибкие из плоскости диафрагмы. Приняты следующие исходные данные: размеры оболочки в плане а=Ь—1ы; толщина оболочки /г=0,01 М; кривизны к]=к2 —0,2 м *; модуль упругости материала Е—2,1 * 1011 Па; коэффициент Пуассона У=0,3; интенсивность поперечной равномерно распределенной нагрузки 0=0,1289-106 Па. Ввиду симметрии задачи рассматривалась четверть оболочки. Результаты решения методом Флетчера-Ривса при различной густоте разностной сетки даны в таблице 2. Здесь же приводится решение, взятое из монографии М.С. Корнишина [5].
Таблица 2
Прогиб в центре оболочки, м
Метод Флетчера-Ривса Метод Ритца [5]
Сетка 4x4 Сетка 5x5 Сетка 6x6 Сетка 7x7 Сетка 13x13
0,2099-102 0,1971-Ю'2 0,1922-10'2 0,1898-10'2 0,1859-10'2 0,221-10’2
Применение алгоритма метода Девидона-Флетчера-Пауэлла, относящегося к группе квазиньютоновских методов, иллюстрируется на примерах расчета пологих сферической и цилиндрической оболочек, опертых по контуру на жесткие в плоскости и абсолютно гибкие из плоскости диафрагмы (рис. 1,2).
Рис. 1. Пологая сферическая оболочка
Рис. 2. Пологая цилиндрическая оболочка
Исходные данные для расчета сферической оболочки: размеры плане 0=6=6 М; толщина й=0,1 м; кривизны ^1=^2=0,0555 м"1; модуль упругости материала £=2-Ю10 Па; коэффициент Пуассона У=0,2; интенсивность поперечной равномерно распределенной нагрузки <7=0,3 • 104 Па. На четверть оболочки накладывалась сетка 6x6, в разностных фор-мулах (5) принимался шаг 8=0,1 • 10 .
Исходные данные для расчета цилиндрической оболочки: размеры плане а—10 м, 6=20 М; толщина И—0,06 м; кривизны К1=0,1 м'1, Кг=0; модуль упругости материала £=3,06-Ю10 Па; коэффициент Пуассона У=0,2; интенсивность поперечной равномерно распределенной нагрузки #=0,3 104 Па. На четверть оболочки накладывалась сетка 8x8, в разностных формулах (5) принимался шаг 6=0,1 -10".
На рис. 3, 4 показаны изолинии нормальных перемещений для четверти оболочки.
Возможности решения геометрически нелинейных задач теории оболочек методом Дэ-видона-Флетчера-Пауэлла иллюстрируются на примере расчета пологой шарнирнонеподвижно закрепленной цилиндрической панели с параметрами к = Ь2/{ЯИ) = 10, £ = (1 -у)(Ь/И)2 /2 = 3500 при различных уровнях нагружения
р - дЪ* (1 — V2 )/(£Уг4). Дуга окружности разбивалась на 8 элементов (9 узловых точек и 27 неизвестных). Результаты расчетов приведены в табл. 3, где даны значения минимизируемой функции Р* (и) - 2Ь3(Т 2) И пР0ГИ®а = ~ ■
_______________________________________________________________________________Таблица 3
Номер шага нагружения Текущее значение нагрузки Количество итераций на шаге Значение минимизируемой функции Прогиб в центре панели
1 4,333 27 -0,1783 0,06697
2 8,667 27 -0,7535 0,1467
3 13,00 27 -1,815 0,2489
4 17,33 28 -3,547 0,4075
10
I
М I и \
I \\ ' *\\\
\ \ \ \ '
\ Vх 4
\
- -0.04--------
Рис. 3. Изолинии нормальных перемещений г10'2мв сферической оболочке
1 3 5
Рис.4. Изолинии нормальных перемещений уи- 10"2 м в цилиндрической оболочке
При использовании метода сопряженных градиентов квадратичная функция п переменных минимизируется за п шагов, то есть необходимо выполнить по одному шагу в каждом из сопряженных направлений [1].
Алгоритмы методов Флетчера-Ривса, Полака-Рибьера и Дэвидона-Флетчера-Пауэлла реализованы в виде пакета прикладных программ, написанного на алгоритмическом языке Фортран с использованием Fortran Power Station.
ЛИТЕРАТУРА
1. Хчммельблау Д. Прикладное нелинейное программирование. - М.: Мир, 1975. -534 с.
2. Моисеев Н.Н., Иванилов Ю.П., Столярова Е.М. Методы оптимизации. - М.: Наука, 1978. -352 с.
3. Вержбицкий В.М. Численные методы (линейная алгебра и нелинейные уравнения) - М.: Высшая школа, 2000. -266 с.
4. Тимошенко С.П., Войновский-Кригер С. Пластинки и оболочки. - М.: Физматгиз, 1963.-636 с.
5. Корнишин М.С. Нелинейные задачи теории пластин и пологих оболочек и методы их решения. - М.: Наука, 1964. - 192 с.
ANALYSIS OF PLATES AND SHELLS BY NONLINEAR PROGRAMMING METHODS
S.I.Trushin
Department of Computer Science and Applied Mathematics Moscow State University of Civil Engineering Jaroslavscoe Shosse, 26, Moscow 129337, RUSSIA
Some affective algorithms for analysis of plates and shells based on finite-difference energy procedure and methods for nonlinear programming are described. To solve the problem of minimisation of quadratic function the most effective methods may be conjugate-gradient method by Fletcher-Reeves and Quasi-Newton method by Davidon-Fletcher-Powell. Some test problems were solved with the help of the proposed numerical algorithms. These numerical procedures are used for analysis of plates and shallow shells.
Трушин Сергей Иванович, родился в 1951 году, окончил в 1974 году Саратовский политехнический институт. Доктор технических наук, профессор кафедры информатики и прикладной математики Московского государственного строительного университета. Автор более 60 публикаций, в том числе 2-х монографий и 5-ти учебных пособий. Основные научные интересы: численные методы решения линейных и нелинейных статических и динамических задач расчета тонкостенных пространственных конструкций.
Trushin Sergey Ivanovich (b. 1951) graduated from Saratov Politechnic Institute in 1974, doctor of engineering science, professor of Department of Computer Science and Applied Mathematics of the Moscow State University of Civil Engineering. Author of 60 scientific publications including 2 monographs and 5 manuals for students. His primary research interests are numerical methods for nonlinear stability and dynamic analysis of thin-walled structures.
