CC BY
58
Вычислительные технологии
Том 8, № 6, 2003
ВЫБОР ОТСЧЕТНОЙ ПОВЕРХНОСТИ В УРАВНЕНИЯХ ПЛАСТИН И ОБОЛОЧЕК*
С. Н. КОРОБЕЙНИКОВ Институт гидродинамики им. М.А. Лаврентьева СО РАН,
Новосибирск, Россия
А. В. Шутов Новосибирский государственный университет, Россия e-mail: [email protected], [email protected]
An influence of the reference surface placement on the solution of plates and shells deformation problem is investigated. A problem on a beam bending by the central concentrated force is considered. For some types of boundary conditions the reference surface placement strongly influences the value of a deflection in the center of the beam. The C0-formulation of isoparametric shell element with a possibility of the reference surface placement control is developed. The use of such element allows to receive a correct pattern of deformation for various types of boundary conditions. It is underlined that some problems on plates and shells deformation cannot be solved using the elements, which formulation does not stipulate the reference surface placement control.
Введение
В практике решения задач механики деформируемого твердого тела (МДТТ) возникают ситуации, когда приходится иметь дело с областями, один размер которых существенно меньше двух других (далее такие области называем пластинами или оболочками). Задачи о деформировании оболочек в классическом подходе решаются заменой исходной трехмерной области в евклидовом пространстве на двумерную, так что основные геометрические особенности этих областей воспроизводятся отсчетной поверхностью и толщиной оболочки. "Платой" за уменьшение размерности является формулировка уравнений в общем случае (кроме уравнений пластин) в двумерном неевклидовом пространстве. Уравнения теории пластин и оболочек (ТПО) получаются наложением некоторых дополнительных кинематических и статических гипотез либо прямой формулировкой [1], либо последовательным сведением трехмерных уравнений к двумерным [2]. Второй путь вывода уравнений ТПО более предпочтителен. В том и в другом случае получается довольно громоздкая система уравнений, использующая аппарат тензорного анализа в двумерных
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 02-01-00195) и Межвузовской научно-технической программы "Фундаментальные исследования высшей школы в области естественных и гуманитарных наук. Университеты России" (тема УР.04.01.005).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2003.
неевклидовых пространствах. При этом теряется простота (по крайней мере в геометрически линейной постановке) исходных трехмерных уравнений.
Первые работы (60-е годы прошлого века) по применению метода конечных элементов (МКЭ) к решению задач о деформировании оболочек были основаны на дискретизации уравнений ТПО (обзор работ см, например, в [3]). При этом обнаружился неприятный эффект, связанный с решением методом конечных элементов задач ТПО, основанных на кинематической гипотезе Кирхгофа — Лява (волокно, направленное по нормали к от-счетной поверхности оболочки, в процессе деформирования не изменяет своей длины и остается нормальным к этой поверхности): аппроксимируемые поля перемещений должны обладать свойством непрерывной дифференцируемости при переходе через границы элементов, в то время как для решения трехмерных уравнений МДТТ, записанных в слабой форме, такое требование не накладывается. Для решения таких задач ТПО требуется использовать усложненные (по сравнению с решениями трехмерных задач) эрмитовы базисные функции [4, 5], обеспечивающие сходимость численных решений к точным (С^подход к формированию конечных элементов оболочки [3, 5]).
Начиная с 70-х годов прошлого века исследователи при решении задач о деформировании оболочек МКЭ стали использовать уравнения ТПО, основанные на кинематической гипотезе Тимошенко: волокно, направленное по нормали к отсчетной поверхности оболочки в недеформированном состоянии, в процессе деформирования остается прямым, но не обязательно нормальным к отсчетной поверхности, при этом его длина не изменяется. Развитие формулировок элементов такого типа представлено в [5, 6]. При дискретизации таких уравнений ТПО можно использовать те же базисные функции, что и при решении трехмерных задач МДТТ (аппроксимация полей перемещений с пониженным требованием на гладкость этих полей). Такой подход к формулировке конечных элементов оболочки принято называть С 0 -подходом [3, 5].
Также в 70-х годах прошлого века возник более практичный подход к решению задач о деформировании оболочек МКЭ, основанный на прямом введении в трехмерные уравнения статической гипотезы Кирхгофа — Лява (нормальные напряжения пренебрежимо малы) и кинематической гипотезы Тимошенко, минуя этап формирования уравнений ТПО [7]. Этот вариант С0-подхода в последнее время стал широко распространенным (см., например, монографии [3, 5, 8]), поскольку давал более простой путь дискретизации геометрически и физически нелинейных уравнений МДТТ по сравнению с первым вариантом С0
-подхода.
Результаты исследований по разработке и применению МКЭ к решению ряда задач деформирования пластин и оболочек в рамках второго варианта
С0
-подхода приведены
также в работах авторов [9-17]. В формулировку разработанного конечного элемента, представленного в этих работах, заложена, следуя [18, 19], возможность управления от-счетной поверхностью оболочки: она может как совпадать со срединной поверхностью оболочки (стандартный подход при выводе уравнений ТПО), так и помещаться на любом уровне от срединной поверхности в трехмерной области, занимаемой оболочкой, совпадая, в частности, с одной из лицевых поверхностей. В [19] рекомендуется при решении контактных задач отсчетную поверхность всегда совмещать с лицевой контактируемой поверхностью оболочки. При этом подразумевается, что расположение контактной поверхности не должно сильно влиять на решение задачи о деформировании оболочки.
Цель настоящей работы состоит в исследовании влияния расположения отсчетной поверхности и связанных с ней линий, на которых ставятся граничные условия, на решение задач о деформировании оболочек. Численные решения получены МКЭ как с исполь-
зованием второго варианта C 0 -подхода к решению задач о деформировании оболочек, так и решением трехмерных задач без использования каких-либо гипотез, принятых в теории оболочек. Численные расчеты проведены с помощью вычислительного комплекса PIONER [12]. В целях сравнения и подтверждения достоверности полученных решений приведены решения тех же задач с использованием широко распространенного пакета MSC/NASTRAN for Windows 4.6, разработанного корпорацией MSC software, демо-версия которого любезно предоставлена авторам московским представительством этой корпорации.
Получено аналитическое решение задачи об изгибе балки сосредоточенной силой, приложенной в центре. Показано, что при решении этой задачи для граничных условий шарнирного опирания торцов, стесненных в продольном направлении, положение отсчетной поверхности оказывает существенное влияние на величину прогиба балки. Численные решения этой задачи получены с использованием как трехмерных изопараметрических элементов, так и элементов пластин и оболочек. Сравнение численных решений с аналитическими показывает, что использование трехмерных элементов и элементов пластин и оболочек, в формулировке которых заложена возможность управления положением отсчетной поверхности, позволяет получить правильную картину деформирования для любого типа граничных условий.
Главный вывод работы состоит в том, что если в формулировке конечного элемента оболочки не предусмотрено управление положением отсчетной поверхности, то некоторые задачи деформирования пластин и оболочек в принципе не решаются с использованием таких элементов.
1. Применение МКЭ к решению задач о деформировании пластин и оболочек
В настоящем разделе рассматриваются формулировки статических уравнений линейной теории упругости. Приводится схема дискретизации этих уравнений МКЭ, рассматриваемого как вариант метода Ритца со специальным набором базисных функций, определенных на локальных подобластях исследуемой области. Описывается технология сведения трехмерных элементов к элементам оболочек, причем в последних предусматривается возможность управления отсчетной поверхностью.
1.1. Статические уравнения линейной теории упругости
Приведем уравнения, описывающие деформирование тел из упругих материалов в предположении малости деформаций, поворотов и перемещений [20].
1. Дифференциальные уравнения равновесия (без учета действия объемных сил) и граничные условия
V ■ а = 0 в V, и = и* на Би, N ■ а = Т* на Бт. (1)
Здесь и далее а — тензор напряжений Коши; V ■ а — дивергенция тензора напряжений (вектор); и — вектор перемещений; V — область, занимаемая телом; Б — замкнутая поверхность, ограничивающая V; Би, Бт — части поверхности Б = Би и Бт, на которых заданы векторы перемещений и и напряжений Т = N ■ а соответственно; N — единичный
вектор внешней нормали к поверхности Бт; знак * обозначает заданную величину; точка между тензорами и/или векторами обозначает операцию их скалярного (внутреннего) произведения.
2. Кинематические соотношения
с = + Уит), (2)
где с — тензор деформаций Коши; Уи — тензор градиента перемещений; здесь и далее индекс "т" обозначает операцию транспонирования.
3. Определяющие соотношения закона Гука имеют вид линейной однородной связи тензора напряжений с тензором деформаций
^ = С : с, (3)
где С — тензор четвертого ранга констант материала (для изотропного материала компоненты тензора С образуются с помощью двух констант: модуля Юнга Е и коэффициента Пуассона V). Компоненты тензора С обладают симметриями вида С^ы, = Сш = Сыщ (Сцш — компоненты тензора в декартовой системе координат, здесь и далее индексы компонент тензоров и векторов пробегают значения 1,2,3). Здесь и далее две точки между тензорами обозначают операцию их двойного скалярного произведения (свертки по двум индексам, например, в декартовой системе координат соотношение (3) в компонентном виде имеет следующую запись с учетом использования правила суммирования по немому индексу: сщ = С"Ы6Ы).
Рассмотрим альтернативную к (1)-(3) (ослабленную) вариационную формулировку уравнений. Введем функционал
I (и) = [ и (и) дУ - [ Т* ■ и ¿Б (и = и* на Би), (4)
иу ¿вт
представляющий полную потенциальную энергию упругого тела, являющуюся суммой потенциальной энергии внутренних сил (интеграл по области У) и потенциальной энергии поверхностных (внешних) сил (интеграл по поверхности Бт). Поля вектора перемещений и, удовлетворяющие главному граничному условию (на Би) и имеющие конечное значение интеграла по области У в правой части (4) (т. е. принадлежащие пространству функций Н1 [21]), называются допустимыми. В (4) введена скалярнозначная однородная функция второй степени компонент тензора деформаций
и = 1 с : С : с, (5)
называемая удельной энергией деформаций. Учитывая кинематические соотношения (2), получаем и(и).
Пользуясь (5), определяющие соотношения (3) можно переписать в виде
. . ди(и) ди
ст(и) = — ~ а« = ■ (6)
Допустимые поля перемещений удовлетворяют главным граничным условиям (на Би), но могут не удовлетворять естественным граничным условиям (на Бт). Подставляя (6) в (1),
получим уравнения равновесия и естественные граничные условия, записанные в перемещениях
V ■ а(и) = 0 в V, N ■ а(и) = Т* на Бт. (7)
Вариационная формулировка уравнений равновесия основана на минимизации функционала (4). Необходимое условие минимума приводит к вариационному уравнению
51 (и,5и) = 0, (8)
где 5и — виртуальные перемещения (вариации): произвольные векторные поля из пространства функций Н1, удовлетворяющие однородным граничным условиям на Би. Отметим, что уравнение (8) можно трактовать как слабую форму уравнения (7). Принцип минимума полной потенциальной энергии утверждает, что для допустимых полей перемещений из пространства С2 (дважды непрерывно дифференцируемых) вариационная формулировка уравнений равновесия (8) эквивалентна дифференциальной формулировке (7). В то же время вариационная формулировка может использоваться для более широкого пространства допустимых функций (из Н1) по сравнению с дифференциальной (минимальное требование — принадлежность функций к пространству С2).
1.2. Пространственная дискретизация трехмерных уравнений
Идея метода Ритца заключается в поиске минимума функционала (4) не на исходном бесконечномерном пространстве допустимых функций, а на его проекциях на некоторые конечномерные пространства и поиске минимума в этих проекционных пространствах, последовательность которых может сколь угодно близко аппроксимировать исходное бесконечномерное пространство. Если в классическом подходе к методу Ритца в качестве базисных функций конечномерных пространств используются полиномиальные или тригонометрические функции, определенные на всей области, то специфика МКЭ состоит в выборе "почти ортогональных" базисных функций, отличных от нуля на локальных подобластях — конечных элементах. В С0-формулировке МКЭ базисные функции принадлежат пространству Н1 : на границах элементов накладывается только условие их непрерывности, а внутри элементов — непрерывной дифференцируемости.
В последнее время наибольшее распространение получили элементы, основанные на отображениях подобластей (элементов) на стандартный куб с длиной ребра, равной двум при определении матриц и векторов элементов. При этом в каждом элементе вводится локальная (естественная) система координат г, в, Ь (рис. 1)
-1 ^ г ^ +1, -1 ^ в ^ +1, -1 ^ Ь ^ +1.
В качестве базисных функций, определенных в некотором элементе, используем полиномы Нк(г,в,Ь), обладающие свойствами [3, 4, 5, 8]
, , , г 7 / ч п /, , —т-рч г (1, если к = I,
Нк (Г1,в1 ,Ь) = 5к1, > Нк (г,в,Ь) = 1 (к, I = 1,М), 5 к1 =< ,
к=1 10, если к = I.
Здесь N ^ 8 — число узловых точек в элементе (на рис. 1 приведен криволинейный конечный элемент с N = 20); г, вг, ¿г — локальные координаты 1-й узловой точки элемента. Рассматриваем изопараметрические конечные элементы, радиусы-векторы материальных
71
-2-
Рис. 1. Моделирование деформируемого тела конечными элементами (приведен конечный элемент с квадратичной аппроксимацией геометрии и перемещений).
г
2
г
точек которых х = [жьж2,ж3]т и вектор допустимых перемещений и = [п1,п2,и3]т аппроксимируются одним и тем же набором полиномов кк:
N N
*) = £ кк (т,8,г)ик. (9)
к=1 к=1
Здесь верхний правый индекс к указывает на то, что данные величины определяются в к-й узловой точке элемента:
хк = х(гк ,вк ,Ьк), ик = и(гк ,вк ,Ьк). Используя стандартную технику МКЭ, получаем дискретный аналог функционала (4)
4(И) = 2итки - ИтИ, (10)
где и — вектор глобальных степеней свободы ансамбля элементов
и = [и ,и2,...,^вд]т,
и (г = 1,NEQ) — одна из компонент вектора перемещений в некоторой узловой точке; NEQ — общее число независимых компонент вектора перемещений ансамбля узловых точек; К — глобальная матрица жесткости ансамбля элементов (размер NEQ х NEQ); И, — глобальный вектор внешних сил (размер NEQ) представляет собой сумму векторов поверхностных Иу и сосредоточенных Ис сил
И = Иу + Ис,
так что скаляр — Ит Иу является дискретным аналогом потенциальной энергии поверхностных сил (второй член в правой части (4)), а Ис — вектор сосредоточенных сил, который получается в результате добавления в правую часть (4) потенциальной энергии этих сил. Отметим, что сосредоточенные силы вводятся в дифференциальные уравнения или естественные граничные условия в (7) (т. е. они могут действовать как в области V, так и на границе Бу) с помощью ^-функции Дирака. Необходимым условием минимума функционала определенного в (10), является удовлетворение вектором И равенства
ки = И. (11)
1.3. Пространственная дискретизация трехмерных уравнений, описывающих деформирование пластин и оболочек
Приведенную в подразд. 1.2 пространственную дискретизацию трехмерными конечными элементами можно использовать для решения задач о деформировании тонкостенных конструкций (например, для решения задач об изгибе пластинки можно использовать дискретизацию, приведенную на рис. 2, а). Однако при решении системы (11) с помощью дискретизации трехмерными конечными элементами, приведенной на рис. 2, б, встречаются трудности, связанные с тем, что если один линейный размер элемента существенно превосходит другой (другие), то матрица K становится плохо обусловленной [3]. Плохой обусловленности можно избежать, вводя статическую гипотезу Кирхгофа — Лява: напряжения, действующие вдоль волокна, нормального к срединной поверхности оболочки, пренебрежимо малы [1].
Техника введения этой гипотезы в формулировку криволинейных трехмерных конечных элементов приведена в [13, 15] и состоит в следующем. Локальная система координат в элементе всегда ориентируется таким образом, что координата t направлена по "нормали" к срединной поверхности оболочки (рис. 3). В каждой точке интегрирования вводится собственная система координат слоя. Слой определяется поверхностью t = tint, где tint — значение локальной координаты t в точке интегрирования. В этой точке вводятся два касательных к слою вектора единичной длины e1, e2, ортогональных друг другу. Третий вектор единичной длины e3 вводится как нормальный к слою вектор, ортогональный первым двум. Локальную декартову систему координат с ортонормальными базисными векторами e1, e2, e3 назовем системой координат слоя. Компоненты тензоров в этой системе координат обозначаем крышкой сверху. Введение статической гипотезы Кирхгофа — Лява соответствует предположению ст33 = 0. Матрица жесткости K модифицируется путем учета последнего ограничения.
В направлении t в этом элементе используется линейная интерполяция. При наложении некоторых граничных условий на оболочку использование такого элемента позволяет избежать плохой обусловленности касательной матрицы жесткости и получить решение задачи. Тем не менее в общем случае применение этого элемента не приводит к желаемому результату вследствие наличия нулевой паразитической собственной моды (при некоторых граничных условиях ранг матрицы K для ансамбля тонкостенных элементов может
Рис. 2. Варианты дискретизации пластинки трехмерными конечными элементами: а — дискретизация, которую можно использовать для решения задачи о деформировании пластинки; б — дискретизация, которая может привести к недостоверному решению.
Рис. 3. Естественная система координат: а — тонкостенный элемент; элемент оболочки с£ = 0 (б), £ = +1 (в), £ = —1 (г) (отсчетная поверхность заштрихована).
оказаться меньше, чем для обычного трехмерного элемента). Эта ситуация обусловлена слишком большой "мягкостью" элемента в нормальном направлении.
Для увеличения "жесткости" элемента в нормальном направлении кроме статической гипотезы введем кинематическую гипотезу Тимошенко: волокно, направленное вдоль вектора нормали к срединной поверхности оболочки, не искривляется и не изменяет своей длины в процессе деформирования1. В силу того, что в направлении координаты Ь используются два узла и вводится линейная интерполяция геометрии и неизвестных по "толщине" элементов, первая часть гипотезы Тимошенко (в процессе деформирования волокно, направленное в отсчетной конфигурации по нормали, остается прямым) выполняется автоматически. Более того, эта часть гипотезы Тимошенко ослабляется: вместо волокна, направленного по нормали к срединной поверхности, используем волокно, соединяющее два узла: на нижней (Ь = —1) и верхней (Ь = +1) поверхностях элемента. Такое ослабление позволяет существенно упростить конечно-элементное моделирование тонкостенных конструкций, например, при моделировании оболочек с изломом срединной поверхности (рис. 4).
Интерполяцию геометрии элемента, выраженную первой формулой в (9), можно представить в виде
N
х(г,М) = (^)Хк (*), (12)
к=1
где N — число пар "верхних" и "нижних" узлов элемента; к — текущий номер пары узлов;
1 Иногда накладывается и более сильная кинематическая гипотеза Кирхгофа — Лява: в дополнение к гипотезе Тимошенко принимается, что волокно, направленное по нормали к срединной поверхности оболочки, в процессе деформирования остается направленным по нормали. Это предположение накладывает более жесткие условия на гладкость исследуемых функций и создает математические трудности, связанные с обеспечением сходимости решения при применении МКЭ [3, 4].
Рис. 4. Моделирование тонкого тела с изломом срединной поверхности: а — тонкостенными элементами; б — элементами оболочки.
/1к (г, 5) — интерполяционные функции слоя; хк (¿) — интерполяция радиуса-вектора х для к-й пары узлов в ¿-направлении
-к / 1+ £ к 1 £ хк (¿) = "Г" хи + —
х
(13)
Здесь х^, хгк — радиусы-векторы "верхнего" и "нижнего" узлов к-й пары узлов:
хи = х(Гк, 5к, +1), хк = х(Гк, 5к, -1).
Пусть // — некоторая заданная величина (-1 ^ // ^ +1). Поверхность £ = определенную в элементе, назовем отсчетной поверхностью элемента оболочки. В частности, при // = 0 отсчетная поверхность совпадает со срединной поверхностью оболочки (рис. 3, б), при £ = +1 отсчетной поверхностью служит "верхняя" поверхность оболочки (рис. 3, в), а при // = —1 — "нижняя" поверхность оболочки (рис. 3, г). Определим радиус-вектор к-й узловой точки отсчетной поверхности соответствующей к-й пары, полагая в (13) £ = //,
хк =
1 +/ * —хи+
1-/
х
(14)
Вектор "нормали" единичной длины Vk и "толщина" оболочки ак в к-м узле определяются следующим образом:
хи — хк = ак V, ||2 = 1. (15)
Здесь и далее ||у||2 = \/v2 + + v2 — евклидова норма вектора V = [г>1, Подста-
вляя (14) и (15) в (13), получаем
хк (¿) = хк + —^—ак V.
(16)
Элемент с аппроксимацией по толщине функцией хк (¿), представленной формулой (13),
называем тонкостенным элементом, а формулой (16) — элементом оболочки. В первом
случае для описания геометрии используются координаты радиусов-векторов "верхнего" кк
хи и нижнего хк узлов элемента, а во втором — параметр / и координаты радиуса-вектора отсчетной поверхности оболочки хк, вектора "нормали" единичной длины и "толщины" ак оболочки в узле с локальным номером к. Слова "нормаль" и "толщина" взяты в кавычки, так как из приведенного вывода следует, что на самом деле требуется лишь их близость к настоящим нормали и толщине.
2
Вектор перемещений изопараметрического конечного элемента определяется из (12) по формуле
N
и = X - 0х = ^ кк(г, в) икЦ),
к=1
где
-Ли\ — 0-к
ик(г) — хк(г) - 0хк(г). (17)
Здесь х, Xк и 0х, 0Хк — текущие и начальные положения радиусов-векторов материальных точек элемента. Для тонкостенного элемента вектор ик (г) определяется из (13), (17):
ик (г) — ^ ии + ^ ик, (18)
а для элемента оболочки — из (16), (17) с учетом неизменности "толщины" оболочки ак (составляющая кинематической гипотезы Тимошенко):
г — и
ик (г) - ик + —ак ДУП- (19)
Здесь введены обозначения
ик — хк 0хк ик — хк 0хк ик — хк 0хк ДУк — Ук 0Ук
ии — хи — хи иг — х1 — х1 , и — х — х , Д Уп — Уп — УП-
Из (19) следует, что вектор перемещения любой материальной точки волокна, совпадающий с вектором "нормали" к отсчетной поверхности оболочки в узловой точке к элемента оболочки, можно представить в виде суммы вектора перемещений ик этой узловой точки, находящейся на отсчетной поверхности оболочки, и вклада вектора ДУП, характеризующего поворот вектора 0УП в новое положение УП.
Обозначим через к!, к2, к3 ортонормальные базисные векторы декартовой системы координат. Для описания поворота вектора "нормали" рассмотрим две возможные системы координат.
1. В каждой узловой точке к вводится местная система координат с ортонормальным базисом 0Ук, 0Ук, 0УП, где векторы 0Ук, 0У2к определяются, следуя [3, 8]:
0ук =) к2 х 0УП/\\к2 х 0УП¡2, если ||к2 х 0УП¡2 = 0, 0Ук = 0Ук х 0Ук 1 ^ кз, если \\к2 х 0УП\\2 = 0, 2 П
Знак х обозначает операцию векторного произведения.
2. Используется общая декартова система координат с базисными векторами к1, к2, к3. В линейном приближении (малые углы поворота) для каждой из введенных выше систем координат получаем формулы, описывающие поворот вектора "нормали".
Местная система координат
дуп = вк 0Ук - а 0Ук -
Здесь ак, вк — углы поворотов вектора "нормали" относительно осей, совпадающих по
0т тк 0т тк
направлению с векторами °Ук, 0У2. Глобальная система координат
ДУк = 0Ук
ПП
Здесь 8к — матрица углов поворота
Бк =
0 —ак ак "
ак 0 —ак
— а 2к ак 0
где а, , — углы поворотов вектора "нормали" относительно осей, совпадающих по направлению с векторами к, к2, к3 соответственно.
Таким образом, для элемента оболочки получаем две альтернативных записи вектора
ик(г).
В местной системе координат, служащей для описания поворота вектора "нормали", получаем
ик (г) = ик + ак (—ак °Ук + вк X).
Число степеней свободы в узловой точке к равно пяти: «к, «2, «3, ак, вк. В общей декартовой системе координат имеем
-к (г) = ик +
г — -
ак °Ук
Число степеней свободы в узловой точке к равно шести: «к, п^к, «3, ак, акк, .
Напомним, что для тонкостенного элемента вектор —к(г) имеет вид (18) и число степеней свободы в паре узловых точек к равно шести: 1, п^2, п^3, «к, «к2, «к3.
Выше приведены формулы, служащие основой для получения матриц и элементов оболочки с пятью и шестью степенями свободы на узел. У каждой из этих формулировок есть свои преимущества и недостатки. Так как поворот вектора "нормали" полностью описывается двумя углами поворота, шестая степень свободы в общем случае лишняя. С этой точки зрения использование пяти степеней свободы в узле представляется более логичным по сравнению с использованием шести степеней свободы. Недостатком же формулировки матриц элементов с пятью степенями свободы в узле является использование поворотных степеней свободы ак и вк относительно векторов °У к и , которые определяются через вектор "нормали" . Так как при восстановлении векторов "нормали" в узлах, где стыкуются элементы, происходит разрыв, эти векторы должны вводиться независимо от геометрии элементов (или осредняться в узлах стыковки элементов). Это создает определенные трудности при автоматизации входной информации для описания элементов оболочки в практических расчетах. Для элементов с шестью степенями свободы на узел поворотные степени свободы ак, а^к, ак рассматриваются относительно осей глобальной системы координат х 1, х2, х3 с базисными векторами к 1, к2, к3. Поэтому в местах соединения элементов независимо для каждого элемента можно восстановить свой собственный вектор "нормали". Повороты всех векторов "нормали", приписанных к одной и той же узловой точке, описываются одними и теми же степенями свободы. Это позволяет ограничиться только вводом координат узловых точек срединной поверхности элементов. Недостатком этих элементов является одна лишняя поворотная степень свободы на гладких участках оболочки (соответствующая повороту вектора нормали вокруг себя). Возникающая при этом сингулярность в матрице жесткости обычно регуляризуется введением штрафа [3, 22].
2
2. Задача о деформировании балки сосредоточенной силой
В данном разделе приведены аналитические решения задачи о деформировании балки сосредоточенной силой, приложенной в центре, с различными вариантами задания условия шарнирного опирания на ее торцах. Эти решения полагаются эталонными при тестировании разного типа конечных элементов для оценки полученных с их помощью численных решений.
2.1. Уравнения деформирования балки
Рассмотрим балку с прямоугольным поперечным сечением (рис. 5). Длину балки обозначим через L, ширину и высоту поперечного сечения — через b и h соответственно. Наряду с координатной осью z введем координатную ось t, так что t = —1 соответствует пересечению этой оси с нижней, а t = +1 — с верхней лицевой поверхностью балки. Значение i £ [—1, +1] этой координаты соответствует значению z = 0 (при i = —1 точка z = 0 расположена на нижней лицевой поверхности, а при t = +1 — на верхней лицевой поверхности). Координаты z и t связаны аффинным преобразованием
z = h(t — t)/2. (20)
Начало координатной оси y равноудалено от боковых, а оси x — от торцевых поверхностей балки. Пересечение оси балки (оси x) с торцевыми поверхностями соответствует расположению шарниров, на которые опирается балка. Рассмотрим условия шарнирного опирания двух типов:
— SS1: стесненные в продольном направлении торцы;
— SS2: свободные в продольном направлении торцы.
Стандартное положение оси балки соответствует ее равноудалению от верхней и нижней лицевых поверхностей (i = 0). В настоящей работе рассматривается более общий случай положения этой оси. Варьируя параметр ¿, можно расположить шарниры на любом расстоянии от верхней (нижней) лицевой поверхности балки. Например, значения i = 0; +1; —1 свидетельствуют о том, что шарниры находятся на пересечении торцевых поверхностей с плоскостью срединной поверхности, верхней и нижней лицевой поверхностью балки соответственно. Задавая значение координаты t =1 в (20), получаем расстояние от оси балки до верхней лицевой поверхности hi = h(1 — t)/2 (рис. 5).
Рассмотрим деформирование балки при действии поперечной силы P, приложенной к ее центру (рис. 5). Под центром здесь понимается точка, соответствующая началу системы координат (x = z = y = 0).
Рис. 5. Шарнирно опертая балка под действием сосредоточенной силы Р.
Приведем уравнения плоского изгиба балки, учитывающие сдвиг поперечного сечения к оси балки (т. е. уравнения, основанные на кинематической гипотезе Тимошенко) [23]. В силу симметрии решения относительно плоскости x = 0 рассматриваем деформирование половины балки (x G [0,L/2]).
Уравнения равновесия имеют вид
N'(x) = 0, M'(x) = Q(x), Q'(x) = 0 для x G (0,L/2). (21)
Здесь введены продольная сила N (x), изгибающий момент M (x) и перерезывающая сила Q(x). Здесь и далее штрих у некоторой величины обозначает ее производную по x. Рассматриваем граничные условия
— симметрии (x = 0)
u = а = 0, Q = P/2; (22)
— шарнирного опирания (x = L/2)
551 : u = M = w = 0, (23a)
552 : N = M = w = 0. (23b)
Здесь u, w, a — соответственно x- и z-компоненты векторов перемещений точек, расположенных на оси балки, и угол поворота поперечного сечения вокруг оси y (см. рис. 5).
При использовании гипотезы Тимошенко ненулевые компоненты тензора деформаций записываются в виде
ежж = u' — za', Yxz = 2exz = w' — a. (24)
Запишем определяющие соотношения, принимая закон Гука,
^жж ^"xz GTxz, (25)
где E — модуль Юнга; G = kE/(2(1 + v)) — модуль сдвига (в соответствии с кинематической гипотезой Тимошенко v = 0); k = 5/6 — корректирующий множитель [23]. Систему уравнений замыкает определение сил N, Q и момента M :
N = ажж dS, Q = dS, M = axxz dS. (26)
./s JS JS
Здесь S — область, занимаемая поперечным сечением балки (x = const):
Г Г hi z-b/2
/ (•) dS = / (■) dydz, ./S Jh2 J-b/2
где h2 = z11=—i = — h(1 + i)/2.
2.2. Аналитическое решение задачи
Из (24)-(26) получаем выражения обобщенных сил через обобщенные перемещения оси балки
N = E (Su' — loa'), M = E (/gu' — Да'), Q = SG(w' — a). (27)
Здесь введены обозначения
S = hb, /0 = z dS, /1 = z2 dS. ./S JS
Из (21)-(23) находим
Р Р Ь
N(х) = с, д(я) = -, М(я) = -(я - (28)
где с — константа, которая определяется далее в зависимости от вида граничных условий шарнирного опирания (881 или 882).
2.2.1. Решение задачи с граничными условиями SS1
Из второго соотношения в (27) и третьего в (28) находим а'. Подставляя выражение для а' в первое соотношение (27), получаем
N = 4S - Й u' + ^ (x-2)- (29)
Интегрируя левую и правую части (29) по x в пределах от 0 до L/2, получаем с учетом N(x) = c и того, что и(0) = u(L/2) = 0,
N = c = - 8L P. (30)
Исключая из первых двух соотношений (27) и' и учитывая выражение для M(x) в (28) и N(x) в (30), получаем уравнение для определения а'
а' = (ai + a2x)P, (31)
где
L 2/iS - /02 S
ai = т^ /т n-t2\°T , a2 = -
E (liS - /2)8/i' 2E(IiS - /о2)-
Интегрируя уравнение (31) с учетом граничного условия а(0) = 0 из (22), находим
x2
а^) = ( aix + a2 — ) P- (32)
Из третьего соотношения в (27), второго в (28) и (32) находим выражение для и'. Используя граничные условия и(Ь/2) = 0 из (23а), получаем
L/2
' L2e 1
w(0) = - w'(x) dx = -LP
+
16E(/iS - /02) 4GS
(33)
где
в = T - 6- (34)
2.2.2. Решение задачи с граничными условиями SS2
Из (23b) находим N(x) = 0 (так как N(x) = c = const). Исключая u' из первых двух соотношений в (27) с учетом выражения для M(x) в (28), получаем уравнение (31) с новыми значениями:
= L S = S
ai = 4E (/iS - /2)' a2 = - 2E (/iS - /2)
Интегрируя уравнение (31) с учетом граничного условия а(0) = 0, получаем решение (32) с a1 и a2 из (35). Определяя w(0) так же, как и при решении задачи с граничными условиями SS1, получаем решение (33) со значением
в = S/3.
(35)
2.2.3. Определение прогиба в центре балки
Найдем значение прогиба балки И = —и осевой силы N для задач с граничными условиями шарнирного опирания и и значениями параметра ¿" = 0; —1. Определим интегральные характеристики поперечного сечения балки:
— £ = 0 : 1о = 0, /1 = 6Л,3/12;
— £ = — 1 : /о = 6й2/2, /1 = 6й3/3.
Подставляя эти значения в (30), (33), (34), получаем
— 88!(0). 882(0). 882(-1) : N = О, И, = + I) ;
— 831(-1) : N = — = Ж( 16Ш + I) .
Здесь и далее используются следующие обозначения задач:
— ЯЯ1(0) — задача с граничными условиями и параметром £ = 0;
— Я81(-1) — задача с граничными условиями и параметром £ = —1;
— 8Б2(0) — задача с граничными условиями и параметром £ = 0;
— 882(-1) — задача с граничными условиями и параметром £ = —1.
Проведем анализ полученных решений. В общепринятых уравнениях изгиба балок ось балки равноудалена от верхней и нижней лицевых поверхностей = 0). В этом случае оба типа граничных условий шарнирного опирания и приводят к одному и
тому же решению задачи с нулевым значением осевой силы N. В том случае, когда ось балки смещена с этого положения = 0), в решении с граничными условиями
882(-1)
(свободными для перемещений в осевом направлении торцами) по-прежнему осевая сила N нулевая и решение для прогиба И не изменяется. Однако при задании граничных условий 8Я1(-1) (стесненными для перемещений в осевом направлении торцами) появляется сжимающая в осевом направлении сила N< 0 (Р > 0) и решения для прогиба И могут значительно различаться для двух типов граничных условий и ББ2). Максимальные различия проявляются при £ = +1 или £ = — 1 (числовые значения для одной конкретной задачи приведены в разд. 3).
Обычно при постановке граничных условий на расположение осевой линии балки не обращается внимание и предполагается, что достаточно получить решение задачи об изгибе балки для значения £ = 0 и согласно принципу Сен-Венана все другие решения будут отличаться от полученного незначительно, и это различие проявляется только в области,
расположенной около торцов балки (см., например, [24]). Исследование, проведенное в настоящем разделе, показывает, что следует проявлять осторожность в интерпретации полученных решений и в некоторых случаях надо тщательно моделировать граничные условия. В качестве примера служит решение с граничными условиями SS1, которое сильно изменяется далеко от торцов балки в зависимости от значения параметра i.
3. Результаты численных расчетов
Управление положением отсчетной поверхности пластин и оболочек через задание параметра i £ [—1, +1] введено в формулировку конечного элемента оболочки вычислительного комплекса PIONER [12]. В этом элементе объединены как собственно элементы оболочки с пятью или шестью степенями свободы в узлах, так и тонкостенные элементы с шестью степенями для пары узлов. Более того, часть узлов одного элемента можно рассматривать как узлы оболочки (расположены на отсчетной поверхности), а другую часть — тонкостенного элемента (расположены на "нижней" и "верхней" поверхностях элемента). Такие элементы называются переходными [13, 15]. Они позволяют достаточно гибко моделировать геометрию составных тонкостенных конструкций. Число узлов/пар узлов в элементе (см. разд. 1) варьируется от 4 до 16 (допускается билинейная, биквадратич-ная и бикубическая аппроксимации радиуса-вектора и вектора перемещений оболочки в (r, й)-поверхности).
Продемонстрируем особенности численных решений задач МКЭ, аналитические решения которых приведены в разд. 2. Численные решения получим как с помощью трехмерных изопараметрических конечных элементов, так и с помощью элементов пластин и оболочек. Исходная геометрическая модель с граничными условиями SS1 приведена на рис. 6, а, ас граничными условиями SS2 — на рис. 6, б. Для того чтобы максимально приблизиться к постановке задачи о плоском изгибе балки, сосредоточенная сила P, действующая на балку (см. рис. 5), заменяется погонной силой F (см. рис. 6) постоянной амплитуды, равнодействующая которой равна P.
В силу симметрии решений каждой из задач проводим конечно-элементное моделирование левой половины балки. На левом торце этой модели ставятся граничные условия шарнирного опирания SS1 или SS2, а на правом торце — условия симметрии. Геометрические параметры балки и приложенная нагрузка приведены на рис. 7, в расчетах используется значение силы F = 1 Н. Материал балки — упругий со следующими механическими характеристиками: E = 103 ГПа, v = 0 (такое значение коэффициента Пуассона используется для максимального приближения решений задачи в трехмерной постановке и по теории пластин и оболочек к решению тестовой задачи о деформировании балки, исполь-
Рис. 6. Шарнирно опертая балка под действием погонной силы F: а — граничные условия SS1, б — граничные условия SS2.
а бе
Рис. 7. Половина балки, рассматриваемая в численном моделировании: а — геометрические параметры балки и условия симметрии деформирования; б, в — расположение отсчетной поверхности при £ = 0 и £ = —1 соответственно.
Рис. 8. Конечно-элементные модели балки: трехмерное моделирование линейными 8-узловыми (а) и квадратичными 27-узловыми (б) элементами; моделирование линейными 4-узловыми (в) и квадратичными 8-узловыми (г) элементами оболочек (пластин).
зующей кинематическую гипотезу Тимошенко).
Расчеты проводились с использованием вычислительных комплексов PIONER [12] и MSC/NASTRAN for Windows 4.6 (далее для краткости называется NW4). Для получения входного файла пакета PIONER и графического отображения результатов расчетов, полученных этим пакетом, используется пакет FEMAP 7.1, входящий в качестве пре/пост-процессора в пакет NW4.
Первая серия расчетов проведена с использованием трехмерных элементов с линейной (рис. 8, а) и квадратичной (рис. 8, б) аппроксимациями радиуса-вектора геометрии и компонент вектора перемещений. Приведем некоторые детали численного моделирования. Смещения всех узлов в направлении оси y полагаются равными нулю. Условия симметрии реализуются запретом перемещения по оси x для всех узлов правого торца. Граничным условиям SS1 соответствует равенство нулю всех компонент вектора перемещений для узлов, лежащих на пересечении срединной поверхности (задача SS1(0)) или нижней лицевой поверхности (задача SS1(-1)) с левой торцевой поверхностью. Граничные условия SS2 отличаются от соответствующих условий SS1 тем, что в направлении оси x компонента вектора перемещений считается неизвестной. Узловые точки, соответствующие шарнирному опиранию по линии, проходящей через срединную поверхность, показаны для левого торца на рис. 8 треугольниками. Линией приложения силы для задач SS1( ) и SS2(0) является линия пересечения срединной поверхности балки с правым торцом модели, а для задач SS1(-1) и SS2(-1) — линия пересечения нижней лицевой поверхности с тем же тор-
Значения нормального прогиба Ш (10 4см) центра шарнирно опертой балки при ее изгибе сосредоточенной силой Р=1 Н
Трехмерные элементы Элементы пластин и оболочек
Задача Аналит. РЮШИ РЮШИ
решение 8 27 8 20 4 8 4 8
881(0) 3.887 3.677 3.880 3.881 3.881 3.885 3.887 3.887 3.887
8Я1(-1) 1.730 1.713 1.819 1.784 1.799 1.727 1.730 1.730 1.730
ЯЯ2(0) 3.887 3.677 3.880 3.881 3.881 3.885 3.887 3.887 3.887
8Я2(-1) 3.887 3.686 3.891 3.891 3.891 3.885 3.887 3.887 3.887
цом. При решении задачи трехмерными линейными элементами использовалась одна и та же конечно-элементная модель в расчетах по комплексам РЮМЕИ и МШ4 с 80 восьмиуз-ловыми элементами (рис. 8, а). В то же время расчеты с применением элементов второго порядка аппроксимации проводились по пакету МШ4 с использованием 20-узловых элементов, а по пакету РЮМЕИ — 27-узловых. При этом балка разбивается на 10 элементов (на рис. 8, б приведена конечно-элементная модель, используемая в расчете по пакету РЮМЕИ). Отметим, что общее число узлов конечно-элементных моделей с линейными и квадратичными элементами в расчетах по пакету РЮМЕИ — одно и то же (189), а в расчетах по пакету МШ4 число узлов конечно-элементной модели с линейными элементами равно 189, а с квадратичными элементами — 138.
Вторая серия расчетов проведена с использованием изопараметрических элементов оболочки пакета РЮМЕИ и элементов пластин пакета МШ42. В обоих случаях узловые точки располагаются на отсчетной поверхности пластины, которая моделирует балку. Пластина моделируется или сорока четырехузловыми линейными элементами (рис. 8, в) или десятью восьмиузловыми квадратичными элементами (рис. 8, г) по обоим пакетам. Так же, как и при моделировании балки трехмерными элементами, перемещения всех узлов в направлении оси у полагаем равными нулю. Для удовлетворения условиям симметрии полагаем перемещения по оси х и повороты вокруг оси у узловых точек правого торца равными нулю. Отсчетная поверхность для задач и совпадает со срединной
поверхностью балки, а для задач и — с нижней лицевой поверхностью
балки. Условия шарнирного опирания левого торца реализованы аналогично тому, как описано выше для трехмерных моделей. Сила Р прикладывается по линии, проходящей через узлы правого торца. Общее число узлов в моделях с линейными элементами равно 63, а в моделях с квадратичными элементами — 53. Так как для элементов пластин в пакете МШ4 нельзя использовать значение коэффициента Пуассона V = 0, в расчетах по этому пакету задавалось значение V = 0.0001.
В таблице приведены значения для прогибов И, полученные как из аналитических (см. подразд. 2.2), так и из численных решений задачи об изгибе балки для различных типов граничных условий шарнирного опирания. Напомним, что в аналитическом решении под прогибом понимаем вертикальную компоненту вектора перемещений центра оси балки, положительное значение которой совпадает с направлением действия силы Р. В численном же решении прогиб И определяется (с тем же соглашением о знаках) как вертикальная компонента вектора перемещений в узловых точках, где действует сила Р. Из приведенных значений прогиба И видно его существенное отличие, полученное в решении
2В пакете NW4 имеется возможность управления положением отсчетной поверхности элементов пластины.
задачи 1), от остальных трех задач. Все полученные численные решения качественно воспроизводят эту особенность.
Проведем анализ полученных численных решений. Отметим, что формулировки элементов в комплексах РЮМЕИ и МШ4 различаются, поэтому для одинаковых конечно-элементных моделей результаты расчетов различны. В особенности это различие проявляется в решениях задач линейными трехмерными элементами. В пакете РЮМЕИ доступны только стандартные изопараметрические трехмерные конечные элементы с трилинейной аппроксимацией геометрии и вектора перемещений. При этом численные интегрирования матриц жесткости элементов проводились по квадратурным формулам Гаусса — Лежан-дра с полным порядком 2 х 2 х 2 по локальным координатам г, Уменьшенное интегрирование (с одной точкой) приводит к сингулярной матрице жесткости, т. е. решение этой задачи с такими элементами невозможно. В пакете МШ4 используется более сложная формулировка 8-узловых элементов, когда к трилинейной добавляется триквадра-тичная неполная аппроксимация вектора перемещений с использованием ЬиЬЬ1е-функции [5]. Кроме того, используется раздельное интегрирование матрицы жесткости: численное интегрирование проводится с уменьшенным порядком по сдвиговым компонентам деформации (одна точка интегрирования) и с полным порядком (2 х 2 х 2) по нормальным компонентам деформации. Применение раздельного интегрирования или интегрирования с уменьшенным порядком позволяет ускорить сходимость численного решения (детали использования уменьшенных порядков интегрирования для улучшения сходимости решения см., например, в [3, 5, 8]). Эти усложнения в формулировке восьмиузлового конечного элемента позволяют получить более точное решение задачи линейными элементами пакета МШ4 по сравнению с элементами аналогичного типа пакета РЮМЕИ. В то же время формулировки элементов с триквадратичной аппроксимацией радиуса-вектора геометрии и вектора перемещений основаны на стандартной изопараметрической концепции построения матриц жесткости элемента как в РЮМЕИ, так и в МШ4. Поэтому результаты расчетов с использованием квадратичных элементов по обоим комплексам близки. При интегрировании матриц жесткости 27-узловых элементов в расчетах, полученных по комплексу РЮМЕИ, использовался полный порядок 3 х 3 х 3 интегрирования по квадратурным формулам Гаусса — Лежандра. Расчеты с использованием элементов оболочки вычислительного комплекса РЮМЕИ получены с применением правила раздельного интегрирования как для четырехузловых, так и для восьмиузловых конечных элементов. Значения прогибов Ш, полученные в численных решениях по обоим комплексам, близки и находятся в хорошем соответствии с аналитическим решением.
Численное решение с использованием трехмерных конечных элементов воспроизводит существенное отличие решений для неклассической (в теории изгиба балок) постановки граничных условий шарнирного опирания (задача 8Б(-1)) по сравнению с остальными типами. Эту же особенность воспроизводят расчеты с использованием элементов оболочек/пластин, в формулировку которых включена возможность управления положением отсчетной поверхности.
Типичные деформированные конфигурации приведены на рис. 9. Чтобы отличать элементы пластин со смещенной отсчетной поверхностью (£ = —1) от аналогичных элементов со стандартной отсчетной поверхностью (£ = 0), в пре/пост-процессоре ЕЕМАР 7.1 используется особое графическое изображение последних (рис. 9, г). По-видимому, такое отображение отсчетной поверхности нельзя признать удачным, так как из рисунка следует, что узлы с шарнирным опиранием торца и узлы приложения силы расположены на срединной поверхности пластины, что противоречит постановке граничных условий и
Рис. 9. Деформированные конфигурации (вектор перемещений увеличен в 1000 раз): моделирование двадцатиузловыми квадратичными трехмерными элементами задачи SS1(0) (а) и задачи SS1(-1) (б); моделирование квадратичными восьмиузловыми элементами пластин SS1(0) (в) и задачи SS1(-1) (г).
определению линии приложения силы в математической модели.
Заключение
Приведена C0-формулировка изопараметрического конечного элемента оболочки с возможностью управления положением отсчетной поверхности. Такая возможность позволяет решать задачи с неклассической постановкой граничных условий. Управление положением отсчетной поверхности введено в формулировку конечного элемента оболочки вычислительного комплекса PIONER. Сравнением численных расчетов с аналитическим решением, полученным авторами, показано, что наличие такого элемента позволяет решать задачи с неклассической постановкой граничных условий шарнирного опирания. Показано следующее:
1. Существует класс задач теории пластин и оболочек, когда от изменения положения отсчетной поверхности сильно изменяется решение; это надо учитывать при численном моделировании.
2. Если в формулировке конечного элемента пластины/оболочки не предусмотрено управление отсчетной поверхностью, то некоторые задачи деформирования пластин и оболочек в принципе не решаются с использованием таких элементов.
Авторы благодарят сотрудников московского представительства корпорации MSC software Ю.Р. Мартыненко и А.Н. Солдаткина за оказанные консультации по работе с пакетом MSC/NASTRAN for Windows 4.6.
Список литературы
[1] Новожилов В.В. Теория тонких оболочек. Л.: Гос. союзное изд-во судостроительной промышленности, 1962.
[2] NAGHDI P.M. Foundations of elastic shell theory // Progress in Solid Mechanics. Vol. 4 / L.N. Sneddon, R. Hill (Eds). Amsterdam: North-Holland Publishing Company, 1963. P. 1-90.
[3] ZiENKiEwicz O.C., Taylor R.L. The Finite Element Method. L. et al.: McGraw Hill, 1991.
[4] Strang G., Fix G.J. An Analysis of the Finite Element Method. Englewood Cliffs. New Jersey: Prentice Hall, 1973; Стренг Г., Фикс Дж. Теория метода конечных элементов. М.: Мир, 1977.
[5] Hughes T.J.R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Englewood Cliffs: Prentice-Hall, 1987.
[6] РикАрдс Р.Б. Метод конечных элементов в теории оболочек и пластин. Рига: Зи-натне, 1988.
[7] Ahmad S., Irons B.M., Zienkiewicz O.C. Analysis of thick and thin shell structures by curved elements // Intern. J. Num. Meth. Eng. 1970. Vol. 2. P. 419-451.
[8] BATHE K.-J. Finite Element Procedures in Engineering Analysis. Englewood Cliffs: Prentice-Hall, 1982.
[9] Коробейников С.Н. Расчет подкрепленных оболочек методом конечных элементов // Динамика сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Ин-т гидродинамики. Новосибирск, 1985. Вып. 71. С. 55-64.
[10] Коробейников С.Н. Многоцелевая вычислительная программа по решению задач линейной теории упругости // Динамика сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Ин-т гидродинамики. Новосибирск, 1986. Вып. 75. С. 78-89.
[11] Коробейников С.Н. Применение МКЭ для статического и динамического анализа прямоугольной пластинки // Динамика и прочность элементов авиационных конструкций: Сб. науч. тр. / Новосибирск: Новосибирский электротехн. ин-т, 1987. С. 43-49.
[12] Korobeinikov S.N., Agapov V.P., Bondarenko M.I., Soldatkin A.N. The general purpose nonlinear finite element structural analysis program PIONER // Proc. of the Intern. Conf. on Numerical Methods and Applications / B. Sendov et al. (Eds). Sofia: Publ. House of the Bulgarian Acad. of Sci., 1989. P. 228-233.
[13] Коробейников С.Н. Геометрически нелинейный анализ оболочек с учетом больших приращений поворотов // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Вычисл. центр, Ин-т теор. и прикл. механики. Новосибирск, 1990. Т. 4, № 4. С. 119-126.
[14] Korobeinikov S.N., Alyokhin V.V., Bondarenko M.I. Application of a finite element method for the solution of three dimensional contact problems // Advances in Simulation and Interaction Techniques: Proc. 2nd Intern. Conf. on Computational Structures Technology / M. Papadrakakis, B.H.V. Topping (Eds). Edinburgh: Civil-Comp Press, 1994. P. 165-175.
[15] Korobeinikov S.N, Bondarenko M.I. A material and geometric nonlinear analysis of shells including large rotation increments // Numerical Methods in Engineering' 96 : Proc of the 2nd ECCOMAS Conf. / J.-A. Desideri et al. (Eds). Chichester: John Wiley & Sons, 1996. P. 754-762.
[16] Курзин В.Б., Коробейников С.Н., Рябченко В.П., Ткачева Л.А. Собственные колебания лопастей однородной решетки гидротурбин в жидкости // ПМТФ. 1997. Т. 38, № 2. С. 80-90.
[17] Коробейников С.Н., Худяков Ю.С., Шутов А.В. Математическое моделирование хрупкого разрушения тонких тел // Вычисл. методы и программирование. 2002. Т. 3, № 2. С. 94-117.
[18] Hughes T.J.R., Liu W.K. Nonlinear finite element analysis of shells: Pt 1. Three-dimensional shells // Comp. Meth. Appl. Mech. Eng. 1981. Vol. 26. P. 331-362.
[19] HALLQUIST J.O., Benson D.J. A comparison of an implicit and explicit implementation of the Hughes-Liu shell // Finite Element Metdods for Plate and Shell Structures, Vol. 1. / Eds T.J.R. Hughes, E. Hinton. Swansea: Pineridge Press, 1986. P. 394-431.
[20] Коробейников С.Н. Нелинейное деформирование твердых тел. Новосибирск: Изд-во СО РАН, 2000.
[21] Curnier A. Computational Methods in Solid Mechanics. Dordrecht et. al.: Kluwer Acad. Publ., 1994.
[22] Onate E., Oliver J. A finite element formulation for the geometrically nonlinear analysis of shell // Finite Element Metdods for Plate and Shell Structures. Vol. 2. / T.J.R. Hughes, E. Hinton (Eds). Swansea: Pineridge Press, 1986. P. 83-101.
[23] Washizu K. Variational Methods in Elasticity and Plasticity. Oxford et al.: Pergamon Press, 1982; Васидзу К. Вариационные методы в теории упругости и пластичности. М.: Мир, 1987.
[24] Гастев В.А. Краткий курс сопротивления материалов. М.: Физматгиз, 1959.
Поступила в редакцию 10 апреля 2003 г., в переработанном виде —16 июня 2003 г.
