Линейное и нелинейное деформирование упругих тел на основе трехмерных КЭ при вариативной интерполяции перемещений тема автореферата и диссертации по механике, 01.02.04 ВАК РФ

На права* рукописи

¿А"

Киселёв Анатолий Петрович

ЛИНЕЙНОЕ И НЕЛИНЕ ЙНОЕ ДЕФОРМИРОВАНИЕ УПРУГИХ ТЕЛ НА ОСНОВЕ ТРЕХМЕРНЫХ КЭ ПРИ ВАРИАТИВНОЙ ИНТЕРПОЛЯЦИИ ПЕРЕМЕЩЕНИЙ

01.02.04 - Механика деформируемого твердого тела

Автореферат диссертации на соискание ученой степени доктора технических наук

4 АПР 2013

Волгоград 2013

005051377

005051377

Работа выполнена на кафедре «Водохозяйственное строительство» в Федеральном государственном бюджетном образовательном учреждении высшего профессионального образования «Волгоградский государственный аграрный университет»

Научный консультант - доктор технических наук, профессор

Николаев Анатолий Петрович

Официальные оппоненты: - Серазутдинов Мурат Нуриевич

доктор физико-математических наук, профессор Казанский национальный исследовательский технологический университет, кафедра «Сопротивление материалов», заведующий;

- Козлов Владимир Анатольевич доктор физико-математических наук, профессор Воронежский государственный архитектурно-строительный университет, профессор кафедры «Теоретическая механика»;

- Зверяев Евгений Михайлович доктор технических наук, профессор Институт прикладной математики им. М. В. Келдыша РАН, ведущий научный сотрудник.

Ведущая организация - ФГБОУ ВПО Южно-Российский

государственный технический университет (Новочеркасский политехнический институт)

Защита состоится « 29 » апреля 2013 года в 10.00 часов на заседании диссертационного совета Д.212.028.04 при Волгоградском государственном техническом университете по адресу: 400005, г. Волгоград, пр. Ленина, 28, ауд. 209.

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

Автореферат разослан « 20 » марта 2013 года. Ученый секретарь /п

диссертационного совета ¿ТЗйАг«-*^ Водопьянов Валентин Иванович

Общая характеристик» работы

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

Одним из наиболее эффективных численных методов расчета конструкций является метод конечных элементов (МКЭ). Значительный вклад в развитие метода внесли следующие зарубежные и отечественные исследователи - Р. Клаф, Дж. Аргирис, Дж. Оден, О. К. Зенкевич, К.- Ю. Бате, В. А. Постнов, Н. Н. Шапошников, Н. Г. Бандурин, А. П. Николаев, Ю. В. Клочков и другие. Анализ современных научных публикаций, посвященных вопросу исследования процессов деформирования различных констру кций на основе МКЭ, привел к выводу, что существует еще ряд В1ЖНЫХ проблем, которые требуют принципиально нового решения.

Центральным разделом МКЭ являегся аппроксимация перемещений внутренних точек конечного элемента через узловые величины, подлежащие определению. Одной из основных проблем аппроксимации полей перемещений в МКЭ является учет смещения КЭ как жесткого целого. О существовании этой проблемы неоднократно упоминается в работах многих исследователей МКЭ, в том числе в работах Голованова А. И., Скопинского В. Н. и др. Были также попытки ее решения, которые имели частный характер и не приводили к решению проблемы в целом. Поэтому задача дальнейшего развития теории линейного и нелинейного деформирования инженерных конструкций на основе МКЭ является достаточно актуальной и представляет собой как теоретический, так и практический интерес.

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

Для достижения этой цели были поставлены и решены следующие задачи:

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

• Разработать КЭ четырех-, пяти- и шестигранной формы, за узловые неизвестные которых выбираются компоненты век-

тора перемещений и их первые производные (при размерах матриц жесткости 48x48, 72x72, 96x96 соответственно), с векторной аппроксимацией полей перемещений.

• Разработать аппроксимирующие функции для КЭ в форме тетраэдра с использованием полных трехмерных полиномов и функции формы треугольной призмы на основе использования полиномов Эрмита в комбинации с полными двумерными полиномами.

• Разработать трехмерные КЭ для исследования НДС оболочек произвольной толщины с учетом поперечных и сдвиговых деформаций без привлечения упрощающих гипотез.

• Разработать условия сочленения КЭ с узловыми неизвестными в виде компонент вектора перемещений и их первых производных для двух пересекающихся оболочек вращения.

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

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

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

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

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

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

Методы проведения исследований. Методы линейной и нелинейной теории упругости, вариационные методы механики сплошной среды, методы теории аппроксимации полей перемещений КЭ, чис-

ленные методы интегрирования по объему КЭ, реализованные в алгоритмах формирования матриц жесткости КЭ с использованием функционала Лагранжа.

Основные научные положения, выносимые на защиту:

• новый вариант получения функций формы, основанный на векторном способе аппроксимации полей перемещений трехмерных конечных элементов, для учета смещения конструкции как жесткого целого;

• вариант получения функций формы и матриц жесткости четырех* и пятигранных КЭ, за узловые неизвестные которых приняты компоненты вектора перемещений и их первые производные;

•алгоритмы формирования матриц жесткости четырех-, пяти- и шестигранных объемных КЭ на основе разработанного векторного способа аппроксимации полей перемещений и их первых производных;

• алгоритмы деформирования оболочки как трехмерного тела в геометрически нелинейной постановке при шаговом способе нагру-жения;

• методика определения НДС пересекающихся оболочек и оболочек слоистой структуры с использованием трехмерных конечных элементов;

• алгоритм дискретного продолжения по параметру нагружения в окрестности предельной точки деформирования оболочек с использованием шестигранного конечного элемента в геометрически нелинейной постановке.

Научная новизна диссертационной работы заключается в следующем:

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

2. Разработаны на основе векторного способа аппроксимации полей перемещений алгоритмы формирования матриц жесткости КЭ четырех-, пяти- и шестигранной формы с узловыми неизвестными в виде компонент вектора перемещений и их первых производных для определения НДС однородных и разнородных конструкций, в том числе и в зонах их сочленения.

3. Разработаны на шаге нагружения матрицы жесткости трехмерных КЭ различных конфигураций с узловыми неизвестными в виде перемещений и их первых производных для определения НДС конструкций с учетом геометрической нелинейности.

4. Разработан и реализован алгоритм дискретного продолжения решения по параметру нагружения на основе КЭ шестигранной формы в окрестности предельной точки деформирования оболочки в геометрически нелинейной постановке.

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

Программные модули разработанных КЭ использовались в расчетах на прочность конструктивных элементов нефтегазового и химического оборудования: в ОАО «ОРГЭНЕРГОНЕФТЬ», ОАО «ВОЛГОГРАДНЕФТЕМАШ», ОАО «ЛУКОЙЛ-Волгограднефте-переработка».

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

Апробация работы. Материалы диссертационной работы докладывались и обсуждались:

- на ежегодных научно-практических конференциях Волгоградского государственного аграрного университета 1999-2012 г.;

- на международной конференции «Актуальные проблемы механики оболочек» (г. Казань, 2000 г.);

- на международной научной конференции «Архитектура оболочек и прочностной расчет тонкостенных строительных и машиностроительных конструкций сложной формы (РУДН, г. Москва, 2001 г.);

- на международной научно-технической конференции «Эффективные строительные конструкции: теория и практика» (г. Пенза, 2004 г.);

- в Волгоградском государственном архитектурно-строительном университете на расширенном заседании кафедры «Строительная механика» (г. Волгоград, 2005 г.);

- на международной конференции «Актуальные проблемы нелинейной механики оболочек» (г. Качань, 2008 г.);

- на международной научно-практической конференции «Использование инновационных технологий для решения проблем АПК в современных условиях» (г. Волгоград, 2009 г.);

- на ежегодных международных научно-технических конференциях «Инженерные системы - 2008-2012» (РУДН, г. Москва, 2008-12 г.).

Публикации. Основные результаты исследований, выполненных по теме диссертационной работы, опубликованы в 34 научных: статьях, из которых 20 включены в перечень периодических изданий., рекомендованных ВАК РФ для публикации материалов диссертаций. Из совместных публикаций в диссертацию включены разработки, принадлежащие автору. Список опубликованных работ приводится в конце данного реферата.

Струю-ура и объем диссертации. Диссертация изложена на 288 страницах машинописного текста, содержит титульный лист, оглавление, введение, семь глав основного текста, заключение;, список литературы, приложения, содержит 22 таблицы и 79 рисунков, библиографический список из 367 литературных источников.

Основное содержание работы

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

Отмечается вклад, который внесли в развитие метода конечных элементов отечественные и зарубежные исследователи.

В первой главе изложен краткий исторический обзор развития численного метода КЭ в задачах исследования НДС инженерных конструкций по опубликованным материалам отечественных и зарубежных авторов.

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

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

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

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

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

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

Имеется незначительное число опубликованных работ по исследованию НДС достаточно тонких оболочек в трехмерной постановке, оболочек со ступенчатым изменением толщины стенок, с учетом на-

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

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

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

Во второй главе на основе соотношений теории тонких оболочек с использованием гипотезы Кирхгофа о прямолинейности нормалей разработаны алгоритмы формирования матриц жесткости двумерных треугольных и четырехугольных КЭ, за узловые неизвестные которых принимались компоненты вектора перемещений их первые и вторые производные с размерами матриц жесткости 54x54 и 72x72 соответственно.

Срединная поверхность произвольной непологой оболочки в исходном состоянии задавалась в декартовой системе координат* , х , х3 радиус-вектором

R° =хк(ва)\к, (к=1,2,3; а~1,2), (1)

где

хк,\к - соответственно координаты и орты декартовой системы координат;

в" - криволинейные координаты, определяющие положение точек срединной поверхности.

Ковариантные компоненты тензора деформаций для точки срединной поверхности и точки, расположенной на расстоянии £от срединной поверхности оболочки, будут равны

£aß=\aaß-aüaß)12' 4ß = {^ß ~ slß)'^ (cc,ß = \,l). (2)

Входящие в (2) ковариантные компоненты метрического тензора в деформированном aaß, gaß и исходном a°aß, g"aß состояниях для

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

Принимая во внимание (2), можно записать

:= £аД+СХар, (3)

где £ар, хар - компоненты тензора деформаций и изменения кривизн срединной поверхности.

Контравариантные компоненты тензора напряжений в любой точке оболочки выражаются через ковариантные компоненты тензора деформаций соотношением

= + ^ («>Р,Г,Р = 1.2), (4)

где /, (г ) = - первый инвариант тензора деформаций;

¿^ - контравариантные компоненты метрического тензора;

к,ц - параметры Ламе.

На основе соотношений (4) можно записать матричное выражение

НЦс] И, (5)

где {<т} = {сг11 ,сг22,сг12- вектор-строка напряжений;

{г^ ]■ = , е[г ,1е\г }Г - вектор-строка деформаций произвольного слоя оболочки;

[с]-матрица упругости.

3*3

В качестве примера на основе МКЭ в двумерной постановке исследовалось НДС цилиндрической оболочки конечной длины при воздействии в вертикальной плоскости в середине пролета двух равных по величине и противоположных по направлению сосредоточенных сил, приложенных по концам одного диаметра. Рассчитывалась достаточно длинная цилиндрическая оболочка с круговым неподкреп-ленным отверстием на боковой поверхности, растянутая равномерно распределенной по ее торцам нагрузкой. Результаты расчетов сравнивались с экспериментальными данными других авторов.

Приводится пример расчета напряженно деформированного состояния осесимметрично нагруженного объемного тела вращения с использованием кольцевых конечных элементов четырехугольного сечения. Работа выполнялась по заказу ОАО «Волгограднефтемаш».

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

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

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

Я0 = дгЛ(б'от)іА, (к, т =1,2,3), (6)

где Я0 - радиус-вектор;

Xі,ік - соответственно переменные и орты декартовой системы координат;

вт - переменные криволинейной системы координат.

Дифференцированием (6) по криволинейным координатам в т определяются ковариантные векторы локального базиса в исходном состоянии

дв'

00

В (7) нижние символы после запятой обозначают дифференцирование по соответствующей координате, а индекс 1=1,2,3 соответствуют направлениям в1, в2 и в3.

Рис. 1. Перемещение произвольной точки сплошной среды в результате деформации

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

(8)

а°=а°.а°.

Дифференцированием (7) определяются производные векторов локального базиса а°к, (],к =1,2,3).

Положение точки сплошной среды после деформации определится радиус-вектором

Я = Я°+У. (9)

Ковариантные векторы локального базиса ак в деформированном состоянии определяются дифференцированием (9), в результате скалярного произведения которых определяются компоненты метрического тензора

%=ага,=(а°+У(}(а«+У>4+а°-Уу+а°.Х/+УгХг (10) Вектор перемещений произвольной точки сплошной среды может быть определен

Ч = ита.°т, (11)

где и" - составляющие контравариантные компоненты вектора перемещений в криволинейной системе координат в" (т=1,2,3).

Выражение производных вектора перемещений V в криволинейной системе координат имеет вид

V, = + /*2а° + Л3а°, (і,т,к=1,2,3), (12)

со

где Г%к - символы Кристоффеля второго рода;

/*>/*>/* " многочлены, содержащие компоненты вектора перемещений и их первые производные.

Ковариантные компоненты тензора деформаций определяются на основе соотношений механики сплошной среды

(13)

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

а" = (Цт,»=1,2,3), (14)

где А1""" = Хаща*тп + ц + аша^я)\ (15)

а"1 - контравариантные компоненты метрического тензора, полученные для исходного состояния.

На основе (14) записывается матричная зависимость между вектором конгравариантных компонент тензора напряжений и вектором ковариантных компонент тензора деформаций

И=[сК4 (іб)

где {сг}= {сгп,а22,а-13,ст12,сг2„сг,1}Т - вектор напряжений; {є}={є11,є22,£„,2єа,2є]3,2є,1,}т -вектор деформаций; [с]- матрица упругости.

бхб

а) б)

Рис. 2. Объемный конечный элемент в виде тетраэдра в глобальной х', х2, х и в локальной 4, т}, С, системах координат

Для выполнения численного интегрирования по объему КЭ (с узловыми точками Ц,к,1 или 1,2,3,4) в глобальной системе координат х', х2, X3, (рис.2а) отображается на КЭ в локальной системе координат <£ т], £ изменяющимися в пределах от 0 до +1 (рис.2б).

Координаты х', х2, х3 внутренней точки конечного элемента аппроксимируются через координаты узловых точек линейными соотношениями

+(*; + + & (17)

где х",х",х",х° - координаты узловых точек; а - поочередно принимает значения 1,2,3.

Из (17) можно получить производные глобальных координат в локальной системе и производные локальных координат в глобальной системе.

В каждом узле конечного элемента неизветными величинами выбирались компоненты вектора перемещений и', и2, и или а, V и » в направлениях осей координат х', х2, х3 и их первые производные ил, иь и}, V/, V:, Vз, гс,и н'^. Вектор перемещений узловых точек конечного элемента в глобальной системе координат имел вид

1*48

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

Ч = кх+к1£ + къц + + к£г +...+ к„т)г + к1ят]2 + к19т)С2 + к20С\ (19)

где к1,к2,...,к2о - искомые коэффициенты.

При использовании конечного элемента в виде тетраэдра с первыми производными число неизвестных узловых перемещений будет равно шестнадцати. Поэтому для корректного определения неизвестных коэффициентов аппроксимирующих полиномов в столбец узловых неизвестных, кроме узловых компонент вектора перемещений и их первых производных, были дополнительно включены четыре смешанные производные второго порядка Л,1,,А2,А3 ,А4 по направлению нормали к грани 2-3-4, определяемые выражением

К = + Кп, + . (¡=1,2,3,4), (20)

где

¡7 - любая из трех компонент вектора перемещений и, V или

и — направление нормали к наклонной плоскости тетраэдра 2-3-4.

Вектор узловых неизвестных в локальной системе координат будет иметь вид

Решением системы из двадцати уравнений определяются выражения двадцати аппроксимирующих функций.

Производные , используя метод конечных разностей,

можно представить в виде

А2-А1; = А3 --А1; А^=А4-А';

А2=А2-А'; Ад =А3 - А2; А2=А4-А2; (21)

А3 -А2-А3; А3 = А3 --А1; /£=А4-А3;

А4=А2-А4; А4 = А3 - А4; А4=А4-А'.

Вектор узловых неизвестных содержащий двадцать компонент, с использованием (21) можно выразить через вектор с шестнадцатью компонентами

20x1 201,6 16x1

где \с,у}Т = У ,Чк ^ ^, , , , ^, , }.

Используя (22), выражение для любой компоненты вектора перемещения можно представить в виде

(23)

1х1< 16x1

где - вектор-строка аппроксимирующих функций.

1x16

При представлении трехмерного континуума дискретной моделью часто возникает необходимость использования элементов в виде треугольной призмы (рис. За). Для выполнения численного интегри-

(24)

рования по объему конечного элемента произвольная треугольная призма в глобальной системе коор динат х', х2, х3 отображается на треугольную призму в локальной системе координат £ 77, £ (рис. 36) с координатами узлов і(0;0;-1), ](!;0;-1), к(0;1;-1), 1(0;0; 1), т(1;0;1), п(0;1;1). Зависимость между координатами х', х2, х3 и локальными £ г/, £ конечного элемента определяется соотношением

где х°,...,х°- координаты узловых точек в системе*', х2, х3; а= 1,2,3.

Из (24) определяются производные х1, х2, х3- координат в локальной системе координат и производные локальных координат в координатной системе .х1,*2,*3.

За узловые неизвестные конечного элемента выбирались компоненты вектора перемещений и, V и те в направлениях осей координат* ,Х2 ИХ3 ИХ первые производные II), иа, 11,3, V/, Уг, V Л \і>і, IV. 3 В узловых точках. Таким образом, вектор перемещений узловых точек

конечного элемента в глобальной системе координат будет иметь вид

..... (25)

И,2 » и!гУг> ",2.«"."". «З,«і, и,з ,и'3, а™, м"3, у',..., и-',...}

4

седо)

>(ОД-1)

6)

а)

Рис. 3. Конечный элемент в виде треугольной призмы в системах коордииатд/, х , х3 и ¿¡, г], С

Для аппроксимации перемещений внутренней точки КЭ через узловые неизвестные в треугольной призме используются функции формы с двумя типами полиномов: алгебраические полиномы третьей степени в плоскости £,77 и полиномы Эрмита третьей степени в направлении , перпендикулярном

Для аппроксимации полей перемещений внутренних точек в координатах £ г} конечного элемента использовался полный алгебраический полином третьей степени, который в этом случае содержит десять членов

ЖШ = *1 + ^ + М + К? + к54ч + к6Л2+ *7£3 + (26)

+ к^2т] + к94712+к10г}\ где к,, ...,кю - коэффициенты, подлежащие определению.

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

1x10

В результате количество неизвестных коэффициентов kj оказывалось равным числу компонент вектора узловых неизвестных, что позволило для каждой аппроксимирующей функции составить

систему десяти уравнений с десятью неизвестными к, (/=/...,/0), решением которой определялись коэффициенты А:, . Таким способом получены функции gm

Смешанную производную перемещения первого узла q'Sn, с использованием метода конечных разностей в локальной системе координат, можно выразить через первые производные узловых перемещений

С использованием (28) вектор-столбец узловых перемещений с десятью компонентами можно определить через вектор-столбец = У,q',qkИМеЮЩИЙ ТОЛЬКО девять составляющих компонент

кШкЬ <29>

1,10 9x1

где [£] - матрица преобразования.

В результате любая из трех компонент вектора перемещения будет выражена в виде

где {pf _ вектор строка аппроксимирующих функций.

1*9

Перемещение внутренней точки (30) в развернутом виде будет определяться выражением Я = С1(^г?)Я' +С2(£/7У + + (31)

где (7/, ..., Од - аппроксимирующие функции для треугольной области КЭ.

С учетом (31) перемещение внутренней точки призматического конечного элемента с треугольным основанием в локальной системе координат будет определяться через перемещения узловых точек

д = ОМ,пШУ +С2(£;77)й1(<ГУ +С3&пЫСУ +•••+ (32)

+

+о2(£,пМе)<12,

где

4 4 - полиномы Эрмита. (33)

В качестве объемного шестигранного КЭ выбран произвольный восьмиузловой элемент с узловыми точками 1 (¡),2ф,3(к),4(1),5(т),6(п), 7(р),8(Ь). Положение узловых точек КЭ определяется тремя координатами - х', х2, х3.

Для выполнения численного интегрирования по объему конечного элемента произвольный восьмиузловой элемент отображался на куб в системе координат £,, 77,^ (рис.4). Локальные координаты кубического элемента изменялись в пределах Т],^ < I.

т

/И 'С /

1 у р Г? У

хг ]

Рис. 4. Восьмиузловой КЭ в системах координат х', х2, х1 и

Зависимость между глобальными координатами х1, х2, х3 и локальными ^ Ц, С конечного элемента определяется выражением

8 8

+ Хк Г" 8 (34)

8 8

+ г * 8

где а=1,2,3;

х° ,х" ,х" ,х° ,х"„,К >х">х1 - координаты узловых точек конечного

элемента в глобальной системе коорд инат.

Дифференцированием (34) определялись производные х',х2,х3 -координат в локальной системе координат и производные локальных координат в глобальной систе ме.

Вектор перемещений узловых точек конечного элемента в гло-б!1пьной системе координат будет иметь вид

«XI ' ' ' (35)

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

<7=Игк}> (36)

1*32 32*1

где - вектор-строка аппроксимирующих функций, основанных на

1*32

полиномах Эрмита третьей степени.

Матрицы жесткости объемных конечных элементов формировались на основе условия равенства работ внутренних и внешних сил. Интегрирование, необходимое при формировании матриц жесткости КЭ выполнялось численно.

В качестве примера расчета исследовалось НДС объемного тела при действии на него сосредоточенной силы (рис. 5). Были приняты сле-' дующие исходные данные: величина сосредоточенной силы Р=39.2-10 Н; сторона куба а=3.05 м. Расчет выполнялся в трех вариантах: в первом расчет выполнялся с использованием конечного элемента в виде тетраэдра, во «тором использовался конечный элемент в виде треугольной призмы, в третьем результаты расч<гга получены при использовании

восьмиузлового конечного элемента. Результаты расчетов внесены в таблицу № 1. В таблице приводятся величины вертикального перемещения куба под силой в метрах при различием количестве КЭ дискретной сетки. В 3-й колонке таблицы приводятся результаты, полученные с использованием тетраэдальных КЭ, в 4-й с использованием элемента в виде треугольной призмы, в 5-й с использованием шестигранных КЭ. Во всех вариантах проверялась сходимость вычислительного прюцесса с различным числом точек интегрирования.

1

Рис. 5. Задача о действии сосредоточенной силы на объемное тело на примере использования конечного элемента в віще треугольной призмы

Таблица 1

№ Размер Величина вертикального перемещения

п. п. КЭ сетки куба под силои \ух 1 О м

Тетраэдр Треугольная призма Восьмиузловой КЭ

1 1x1x1 3.8 3.3 3.5

2 2x2x2 4.8 6.8 6.5

3 3x3x3 8.9 10.5 10 2

4 4x4x4 11.2 14.2 13 9

5 5x5x5 13.5 17.9 17 5

6 5x5x6 15.9 18.0 178

7 5x5x10 17.3 18.4 18.1

8 6x6x6 18.8 21.5 212

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

Рассчитывался тонкостенный цилиндр, дискретизация которого выполнялась восьмиузловыми КЗ1. Оболочка загружалась ДЕ;умя рав-

ными по величине и противоположно направленными сосредоточенными силами, приложенными по концам одного диаметра (рис. 6). За исходные данные приняты: радиус срединной поверхности Я=-4.953 л<; длина цилиндра ¿=10.35 м; модуль упругости материала Е= 10.5x10 кПа; коэффициент Пуассона у=0.31; толщина стенки цилиндра ¡1=0.094 м; величина сжимающей силы Р=100 кН. С учетом трех плоскостей симметрии расчет сводился к исследованию восьмой части оболочки. Результаты расчета представлены в таблице 2. В таблице приводятся величины прогиба оболочки в точке приложения сосредоточенных сил в зависимости от количества КЭ дискретной сетки. Результаты расчета сравнивались со значениями, полученными в главе 2 с использованием двумерных дискретных элементов.

т.

Рис. 6. Цилиндрическая оболочка, сжатая в середине пролета сосредоточенными силами

Таблица 2

№ п.п.

Размер сетки дискретизации двумерными КЭ

1x1

2x2

3x6

5x14

6x20

Прогиб оболочки в точке приложения силы, м

0,11046 0,11765 0,11817 0,11810 0,11807

Размер сетки дискретизации трехмерными КЭ

1x1x1 3x3x3 1x3x3 1x10x10 1x12x12 1x14x14

Прогиб оболочки в точке приложения силы, м

0,00064 0,05812 0,05846 0,11207 0,11339 0,11431

В таблице приводятся величины прогиба оболочки в точке приложения сосредоточенных сил в зависимости от количества КЭ дискретной сетки. Результаты расчета сравнивались со значениями, полученными в главе 2 с использованием двумерных дискретных элементов. Как видно из примера результаты, полученные с использованием объемного вось-

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

В главе также приводятся примеры расчета других тонкостенных оболочек, представленных объемными КЭ. Приведен пример рас-, чета защемленной панели, загруженной линейной нагрузкой. Для верификации разработанных алгоритмов, полученные результаты расчета сравнивались с аналитическим решением. Приведены примеры практических задач по расчету на прочность корпуса шарового крана по заказу ОАО «Волгограднефтемаш» и шарового резервуара диаметром 10,5м с несовершенством геометрии по заказу Самарского филиала ОАО «Оргэнергонефть».

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

Для описания геометрии двух пересекающихся под прямым углом цилиндрических оболочек вводятся две системы координат Х,у,2 для основной оболочки с радиусом г и толщиной I и для примыкающего патрубка с радиусом г1 и толщиной (г>г) (рис. 7). Здесь и далее символы без штриха будут относиться к параметрам основной оболочки, со штрихом - к патрубку.

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

К = лт + г5П1# ] + гсо50к (37)

и для примыкающей оболочки

К' + г' йав' У +г'сОБв' к7. (38)

Дифференцированием (37) и (38) по криволинейным координатам можно получить ковариантные векторы локального базиса основной {а°)Г = {а® а.^ а°} и примыкающей (а/0}Г = {а®' а®7 а"7} оболочек.

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

{ Г = { > ".2 > м,з >' ^' ъ, IV, И',, }; (39)

{и'у = {«'У.^Уу'.у'р^.у^ и^У^'з), (40) где цифрам 1 Д,3 соответствуют направления координат*, в, г и х, ё, г.

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

Векторы перемещений точки, расположенной на грани сопряжения, представляются компонентами, отнесенными к базисам векто-ров{а°}и{а0/}

,, 0,0 О тг/ /01 / О/ , / О/ п

V = ма, + га2 + \ча3; V =ы а, +у а2 + ус а3 , (41)

На основе равенства векторов V = V' компоненты вектора примыкающей оболочки выражаются через компоненты вектора основной оболочки в виде

и' =уг ьтв-мсовв;

V1 + ч^тв-усозб?';

г г г

м' =ы соъв' + угсов^втб' + и'зтбзтб'.

Производные вектора перемещений точки, расположенной на грани пересечения оболочек по направлениям х',в',г' примыкающей оболочки в системе координат х,в,г основной оболочки определяются по формуле

V', = ухам(а? • а?)+ Увсо3(а° • а°)+ Угсо5(а° - а?), 0=1.2.3) (43)

Из (43) получим выражения производных компонент вектора перемещения примыкающей оболочки через производные компонент вектора перемещения основной оболочки

",/= "/("д. И',/ V,е. И",, )■ (44)

На основе (44) получаются зависимости производных компонент вектора перемещений на грани пересечения основной оболочки с примыкающей.

Вектор узловых неизвестных примыкающей оболочки с использованием (42), (44) можно выразить через вектор узловых неизвестных основной оболочки

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

НЧгГИ^Н^ГИ- (46)

В качестве примера расчета исследовалось НДС в зоне сочленения двух ортогонально пересекающихся цилиндрических оболочек, находящихся под действием внутреннего давления. Были приняты следующие исходные данные: г=0.127 м, /=0.064 м, /=0.00254 м, /=0.00127 м, модуль упругости материала £=2><108 кПа, коэффициент Пуассона у=0.3, величина внутреннего давления <7=3.445 кПа. Вследствие наличия двух плоскостей симметрии в расчете рассматривалась "Л часть конструкции. Изменением количества конечных элементов дискретной модели проверялась сходимость вычислительного процесса.

На рисунках 8-11 представлены расчетные кривые изменения меридиональных (сплошные линии) и кольцевых (пунктирные линии)

напряжений в продольном сечении рассчитываемой конструкции. На рисунках 8 и 9 изображены графики изменения напряжений во внутренних и наружных волокнах примыкающего цилиндра, соответственно, а на рисунках 10 и 11 — то же самое для основного цилиндра. Экспериментальные значения кольцевых напряжений обозначены светлыми кружками, меридиональных - темными кружками.

а, даН/см2

1ЮО

г

0 •/ ( » 2, » 3,

/

2А СО -1--

Рис. 8.

(Т, даН/см2 1200

%

О, >| 1,0 2.0 З, I)

сг, даН/см2

2403

Xі, СМ <«»

\

\

V

V \п

А ■ \ \

О.

0 0 5 ^^ 2

х'.с

Рис. 9.

Рис. 10.

Рис. 11.

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

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

1x32 (47)

1x32

Вектор перемещений внутренней точки объемного восьмиузлового конечного элемента с использованием полиномов Эрмита третьей степени аппроксимируется через вектор-столбец узловых перемещений

V = (48)

1*32 32x1

где {у}г - вектор-строка аппроксимирующих функций.

Векторы перемещений в локальной системе координат С выражаются через векторы перемещений в глобальной криволинейной системе координат в\вгв виде

32x1 32x32 32x1

где [¿]- матрица преобразования координат.

С использованием (11), (12) можно получить зависимости для векторов узловых точек конечного элемента, которые в матричном виде представляются следующим образом

{<}=[А] {/Л (50)

32x1 32x96 96x1

где [а] - матрица, ненулевыми элементами которой являются векторы базисов узловых точек конечного элемента;

Учитывая равенства (49) и (50), соотношение (48) можно привести к виду

У = МГ[А]И{/;}, (52)

1x32 32x96 96x96 %)<1

где использовано полученное без принципиальных затруднений выражение

[¿][А]=[А][С]. (53)

32x3232x96 32x9696x96

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

+ (54)

С учетом (54) зависимость (52) может быть записана в виде

У = {у,}т (№ (55)

Приравнивая правые части выражений (55) и (11), можно получить выражения для каждой компоненты вектора перемещений произвольной точки

и=У}т1ЛоЬУг\ №]{/-;} ™=ЫТ[АФУГ\ (56)

Последовательным дифференцированием (56) по криволинейным координатам и приравниванием правых частей полученных равенств и выражения (12), определяются интерполяционные зависимости производных компонент вектора перемещений

ЛМ^ИсШ (57)

где к, т=1,2,3.

Столбец неизвестных можно представить выражением

{//ЬММ (58)

где {к/}г={м^м^г/^м^мm,м^м^м^«;1,...,^,...,м;3,...,v^...,wí,4 -векторстрока компонент вектора узловых перемещений; [*] - матрица преобразования.

96x96

Выражения (56), (57) с использованием (58) можно представить через столбец узловых неизвестных (35)

И1

И ИММ И

(58')

Таким образом, в отличие от интерполяционных зависимостей, в которых каждая компонента вектора перемещений внутренней точки конечного элемента аппроксимируется через узловые значения этой же компоненты и ее производные (36), в соотношении (58;) лю-

бая из трех компонент вектора перемещений зависит от всех составляющих компонент узлового вектора.

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

В качестве одного из примеров рассматривалась оболочка, срединная поверхность которой представляет собой форму усеченного эллипсоида вращения. Оболочка нагружалась равномерно распределенным внутренним давлением интенсивности <7=4.0 МПа (рис. 12). Были приняты следующие исходные данные: А= 1.0 м; £=0.5 м; £>-0.9 м; £=2-105 МПа; 1^=0.3; /1=0.01 м, (А и В являются главными полуосями эллипса).

Расчеты выполнялись в двух вариантах: в первом реализовыва-лась общепринятая скалярная аппроксимация полей перемещений; во втором - предложенная векторная аппроксимация.

Величины кольцевых напряжений оболочки при х=0,9 м в зависимости от величины смещения конструкции А как жесткого целого приводятся в таблице 3. Анализ табличного материала показывает, что результаты вычислений в первом и втором вариантах расчета при жестком закреплении оболочки (Л=0.0 м) практически одинаковы и совпадают с аналитическим решением. А при незначительном смещении оболочки как жесткого целого значения напряжений, полученные с использованием традиционной независимой аппроксимации, отличаются от первоначальных значений и с увеличением А становятся неприемлемыми. В то время как напряжения, вычисленные с использованием векторного способа аппроксимации, остаются постоянными независимо от А, и лишь при величине смещения 4,8 м, что больше чем в пять раз превышает размеры самой конструкции, несколько отличаются от полученных аналитически.

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

Рис. 12. Усеченный эллгясонд вращения, загруженный внутренним давлением

Таблица 3

Величина жесткого смещения Д, м

Кольцевые напряжения СГК, МПа

Вариант расчета

I

II

0.00 0.01 0.03 0.05 4.80

126.2 99.7 42.1 -6.02 -2056.5

125.6 125.6 25.6

125.6

124.7

Аналитическое решение ак=125,3 МПа

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

В шестой главе изложен вывод основных соотношений теории деформирования оболочек как трехмерных тел в геометрически нелинейной постановке при шаговом нагружении. При выводе основных соотношений рассматриваются три состояния системы. Первое — начальное или исходное, второе — деформированное состояние после] -шагов нагружения и третье состояние системы, близкое ко второму — ■деформированное после 0+1)-го шага нагружения (рис. 13). Перемещение прои звольной точки среды А4° из первого состояния во второе -в положение точки М - будет характеризоваться суммарным вектором перемещения за _/ шагов нагружения V, а из второго состояния в третье - шаговым вектором перемещения АУ.

м'

Рис. 13. Перемещение произвольной точки сплошной среды в результате деформации при шаговом натр ужении конструкции

Векторы а°,а2,а°; а,,а2,а3 и а^а^а, представляют собой ко-вариантные векторы локальных базисов произвольной точки сплошной среды в трех рассматриваемых состояниях соответственно.

Радиус-вектор И.0 точки в исходном состоянии, определяется выражением (6). Дифференцированием радиус-вектора Я" выражаются ковариантные векторы локального базиса (7). После чего определяются ковариантные и контравариантные компоненты метрического тензора в исходном и деформированном состояниях. Вектор перемещений V произвольной точки сплошной среды в этом случае будет являться суммарным вектором перемещения за ./ шагов нагружения.

Выражения суммарного тензора деформаций за j шагов нагружения определяются по формуле

^ = ~(ат• V, + а» • \'т + V,, • УД (т,п = 1,2,3), (59)

где V, (< = 1,2,3)—производные суммарного вектора перемещений за

] шагов нагружения в глобальной системе координат.

Радиус-векторы, определяющие положение точек М и Л/* послеу и 0+1)-ом шагов нагружения, мо1уг быть представлены выражениями (рис. 13)

И = 110+У; 1Г=К+ДУ; 11*=110 + У + ДУ. (60)

Вектор перемещений точки на шаге нагружения имеет вид

АУ = а® Дм1 + а°Дн2 + а °Ди3 = а° Аиа, (61)

где Ли" - приращения компонент вектора перемещения на. (¡+1)-очл шаге нагружения (а = 1,2,3).

Дифференцированием (60) по глобальным координатам можно определить ковариантные векторы локального базиса после 0+1)-го шага нагружения в виде

а;=Я>аа+ДУа, (62)

где ДУв=а»вДи'+а®Ди£ -производные вектора перемещений

на шаге нагружения.

Ковариантные компоненты тензора приращения деформаций на {/+!)-ом шаге нагружения определяются на основе соотношения

(63)

где

а,;=а,.-ау, 4= а*-а*. (64)

Ковариантные векторы базисов после у шагов и 0+1) -го шага нагружения могут быть представлены в виде

аа = Ё_а = а° + Уа и а*в = а° + Уа + АУа, а= 1,2,3. (65) Из (63) можно записать приращения деформаций д®«. + У^ДУ, + (а; + У^ду,, + ДУ^-ДУ,]. (66)

Приращения напряжений на шаге нагружения определяются матричным выражение

{Лег} = [с]{Дг}, (67)

где {Дсг}= {Дсг",До-22,Дсг33,Дст12,До-23,Дсг3,}г - вектор-строка приращений напряжений;

{Ае}={Аеи,А£21,А£г},2А£12,2Ае1},2А£3„}т -вектор-строка приращений деформаций;

[с]- матрица упругости.

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

|1ИГ +*Мг]{д4^= {{ДуГ [{Р}+*{ДР}]Л, (68)

V 5

где

{;Р},{Д/>} - вектор нагрузок и их приращений; {Д У}т = {Дг/,Ду, Дм»} - вектор приращений перемещений. Значение коэффициента А: принимается равным 1/2. Вектор приращений перемещений {ДК} внутренней точки конечного элемента определяется через вектор приращений узловых перемещений

{ДГИГ]{Л(^}, (69)

где {дгуг} - вектор приращений узловых неизвестных конечного элемента.

Выделяя в приращениях деформаций линейную и нелинейную части вектор приращений {Дг} можно представить суммой

{де}={д£л)+{д£н|, (70)

где элементы матрицы {д^"} зависят от компонент вектора перемещений {А К} и не зависят от компонент вектора перемещений {V} за предыдущие шаги нагружения, элементы матрицы | зависят от всех указанных компонент.

Матрицу |дгл | можно представить в виде

{д^}=[в]{дг;}. (71)

Используя (89), (90), из (87) можно получить

(72)

5 V

Здесь {д Уу}~ вектор приращений узловых смещений на (/+1)-м шаге нагружения.

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

Анализ графиков показывает удовлетворительное соответствие значений перемещений, определенных с использованием описанного алгоритма и полученных экспериментально. В этой главе приводится алгоритм дискретного продолжения решения по параметру в окрестности предельной точки деформирования оболочки на основе метода конечных элементов в геометрически нелинейной постановке. Приводится пример расчета.

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

Г'ис. 14. Графики зависимости вертикального прогиба и> оболочки от нагрузки

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

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

1. Векторы перемещений узлов КЭ расположенных на границе контакта смежных слоев принимаются равными

У'=Р, (73)

откуда получаем равенства компонент г/ = м; у' = V; и>' = м-.

2. Производные векторов перемещений вдоль координатных ли' ний х, у в плоскости сопряжения слоев будут равны

У1,а = и £\ф=£ар> (а, р = 1,2) (74)

на основании чего имеем

и', = и,; у', = V,; у'2=У2; и-', =и»л; = (75)

3. Напряжения в направлении нормали к плоскости сопряжения слоев в совместных узлах КЭ равны

(76)

(77)

сг = и ,

что приводит к соотношению

¿'(4 + + 4)+ 2м'е31 = КЕп + е22 + £п)+ 2№З.

откуда

зз

5/ ../

3 Л' + 2ц

7[(л - Л')(и, + у2)+ (Л + 2ц)м>ъ\ (78)

Л,ц,Л // - параметры Ляме.

4. Будут равны касательные напряжения на границах слоев

_/23 —.732 23 _-,-32 „. /13 _/31_1_1Э „31

<Т = СГ = СГ = СГ И <Т =<Т = <У = <Т ,

следовательно

(79)

(80)

Выражения для производных и'ъ и у', определяются из выражений

Л Є23 — /'£23 > ^ ¿13 ~ /«13 •

V з Ч- і 1 -

./Л

р)

(81)

На основе равенств (75), (78) и (81) формируется матрица преобразования матриц жесткости и векторов узловых нагрузок КЭ, примыкающих к границе раздела слоев конструкции.

Пример. Рассчитывалась трехслойная пластина с защемленным концом и загруженной на свободном крае линейной распределённой нагрузкой интенсивности д (рис. 15).

Рис.15. Трехслойная пластина

В качестве исходных данных были приняты: 6=0.2м; /=0.2м; ¿7=100 Н/см; И/=0.002 м; Л2=0.006 м; £=2*105 МПа;£'=2*104 МПа; // = 0.3; ц' = 0.25.

В таблице 4 приведены результаты численного расчёта, где даны нормальные напряжения для точек 1, 2, ... 4 (рис.15). Штрихами отмечены точки, относящиеся к среднему слою панели. Как видно, в граничных точках имеются скачки в значениях нормальных напряжений.

Таблица 4

Напряжения МПа Точки

1 2 2» 3' 3 4

-151.049 -85.783 -8.32 8.322 85.783 151.049

1.118 1.163 1.175 -1.163 -1.163 -1.118

Для контроля вычислений выполнена проверка условий равновесия по силам(г:з' = 0)( которая выполняется с точностью

6 = 0.0039 % и по моментам =о)точность составила

5 = 0.03

В таблице 5 для сравнения приведены результаты численных расчётов пластины при различных значениях модулей упругости слоев. В строках таблицы в числителе представлены нормальные напряжения (Ууу для точек 1, 2,.. .4 (рис.15), полученные с использованием разработанного алгоритма, а в знаменателе представлены нормальные напряжения СУуу для точек 1, 2, ... 4, полученные на основе конечно - элементного программного комплекса АВАСШБ при модулях упругости слоев материалов, различающихся в 10 раз, в 100 раз, в 1000 раз, в 10000 раз и при различных отношениях длин пролета к толщине панели.

По результатам нормальных напряжений сг^,, полученных с

использованием разработанной программы, и по результатам, полученным на основе конечно - элементного программного комплекса АВА<ЭШ (табл. 5), построен ряд эпюр при £=2*105 МПа, £'=2*102МПа и при отношениях: 1/11=20 (рис.16, рис.16а); 1/Ь=10 (рис.17, рис. 17а); 1/Ь=5 (рис.18, рис.18а).

Таблица 5

Напряжения Отношения длины пролета

СГууМПа, при к толщине панели 1/Ь=20 для точек

1 2 2' 3' 3 4

Е=2*105МПа, -151.05 -85.78 -8.32 8.32 85.78 151.05

£'=2*104МПа. -152.75 -83.69 -8.13 8.13 83.76 152.68

Различие, % 1.1 2.4 2.3 1

Е=2*105 МПа, -176.25 -63.63 -1.77 1.77 63.63 176.25

£'=2*103МПа. -178.23 -59.88 -0.4 0.3 60.57 177.87

Различие, % 1.1 5.9 4.8 0.9

Е=2*105 МПа, -246.02 23.16 0.55 -0.55 -23.16 246.02

£'=2*102МПа. -258.2 33.85 0.36 -0.48 -30.17 255.4

Различие, % 4.7 31.6 23.2 5

Е=2*105 МПа, -488.75 300.6 0.59 -0.59 -300.6 488.75

£'=2*ЮМПа. -521.8 346.89 0.2 -0.7 -330.47 511.85

Различие, % 6.3 13.4 9 4.5

Отношения длины пролета к толщине панели 1/Ь=10

Е=2*105МПа, -155.58 -80.54 -7.52 7.52 80.54 155.58

£"=2*10'МПа. -158.75 -77.22 -7.08 7.08 77.49 158.66

Различие, % 1.8 4 3.8 1.9

Е=2*105 МПа, -198.89 -33.81 0.54 0.54 33.81 198.89

£'=2*10'мПа. -208.24 -25.49 -0.006 -0.007 26.61 207.56

Различие, % 4.5 24.6 21.3 4.2

Е=2* 105 МПа, -341.34 139.98 1.18 -1.18 -139.98 341.34

£'=2*102МПа. -368.44 163.59 1.21 -1.49 -158.16 365.46

Различие, % 7.3 14.4 11.5 6.6

Е=2*103 МПа, -780.68 660.58 1.18 -1.18 -660.1 780.12

£'=2*10МПа. -825.412 718.28 0.89 -1.9 -719.98 852.13

Различие, % 5.4 8 8.3 8.4

Отношения длины пролета к толщине панели 1/11=5

Е=2*103МПа, -164.6 -70.11 -5.94 5.95 70.11 164.6

£'=2*104МПа -164.24 -65.29 -5.76 5.77 65.63 164.02

Различие, % 0.2 6.8 6.3 0.3

Е=2*105 МПа, -247.22 25.55 1.99 -1.98 -25.57 247.25

£'=2*103МПа -254.85 35.03 1.4 -1.5 -33.08 253.56

Различие, % 3 27 22.7 2.5

Е=2* 105 МПа, -543.56 374.49 2.44 -2.44 ■370.11 541.3

£'=2*102МПа -572.6 415.36 2.86 -3.79 -406.65 577.56

Различие, % 5 9.8 9 6.3

Е=2*103 МПа, -1064.27 1025.72 2.39 -2.26 -1038 1076.33

£'=2*10МПа - - - - - -

Различие. % - - - - - -

-23.16 ^ ■^Г 115.!

1

1 ога

-0.83

лза

0.55

___—" •

-

•300 -200 -100

100 200 300

--Н-] . _——■ 255

-30,17^ -- 112,7

-0.48_

, "о'ії-

0^6

33.86

---- *

•"Зьа.2 1 --,-1-»-

-300 -200 -100 0 100 200 300

Рис.16

Рис. 16а

-400 -200 0 200 400 600

Рис.17

Рис. 17а

--

- . ——--.

1.»4

7 о3'

.174

г44

■ЙО.Ь 3

4 -Г—----" 1

♦3415 1 -,--,-а-

-600 -400 -200

200 400 600

--н~ - --___*- ___'--- 577,56

. ———т" 80,45

-3,79

,

"С

______4

*3/2,6 1 --,-,-1—е-

-800 -600 -400 -200 0 200 400 600 800

Рис.18 Рис.18а

По полученным результатам нормальных напряжений сгуу (табл. 5) с использованием эпюр нормальных напряжений в таблице 6 приведены для сравнения результаты проверки условия равновесия по силам (г:_у = 0) и по моментам {?Му = о), где в числителе представ-

лены результаты проверки, полученные с использованием разработанной программы, а в знаменателе представлены результаты проверки, полученные на основе конечно - элементного программного комплекса АВАСЗиБ при модулях упругости слоев материалов, различающихся в 10 раз, в 100 раз, в 1000 раз, в 10000 раз и при различных отношениях длины пролета к толщине панели.

Таблица 6

Условия равновесия Модули упругости слоев материалов при отношении длины пролета к толщине панели 1/Ь=20

Е=2*10 МПа, £'=2*104МПа Е=2*10 МПа, Е'<\* 103МПи Е=2*103 МПа, £'.=2*10^3 Е=2*105МПи, £'=2*10МПа

£Г = 0 6,% 0.0039 0.5 0.004 0.5 0.004 0.5 0.004 0.5

ЪМу =0 <5,% 0.03 0.1 0.08 0.1 0.08 0.7 0.04 1.5

Отношение длины пролета к толщине панели 1/Ь=10

ХУ = 0 д,% 0.0039 3 0.004 2 0.004 2 0.056 3

ХМу =0 8,% 0.03 1 0.09 2 0.05 1 02 1

Отношение длины пролета к толщине панели 1/11=5

17 = 0 6,% 0.0039 0.5 (1.001 1 0.0015 1.3 0.02

ъму= 0 0.03 0.1 0.5 3 0.5 6.4 М

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

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

Основные результаты и выводы

На основании проведенных исследований можно сделать следующие выводы:

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

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

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

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

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

6. Полученные зависимости между перемещениями и их производными на границе слоев материалов с различными физико-механическими свойствами, позволяют преобразование матриц жесткости и узловых усилий КЭ, примыкающих к границам, что приводит к корректному определению НДС слоистых структур.

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

Основные положения диссертации опубликованы в работах.

Публикации в изданиях, рекомендованных Высшей аттестационной комиссией:

1. Киселёв, А. П. Особенности формирования матрицы жестко-ста треугольного конечного элемента размером 54x54/ Ю. В. Клочков, А. П. Николаев, А. П. Киселёв// Изв. вузов, сер. Строительство. -1998.-№2.-С. 32-37.

2. Киселёв, А. П. Применение четырехугольного конечного элемента с матрицей жесткости 72x72 для расчета оболочечных конструкций / Ю. В. Клочков, А. П. Николаев, А. П. Киселёв // Изв. вузов, сер. Строительство. - 1998. - № 4-5. - С. 36-41.

3. Киселёв, А. П. Решение проблемы учета смещения конечного элемента как жесткого целого на основе векторной интерполяции полей перемещений / А. П. Николаев, Ю. В. Клочков, А. П. Киселёв // Изв. вузов, сер. Машиностроение. - 1998. - № 1-3. - С. 3-8.

4. Киселёв, А. П. Треугольный конечный элемент произвольной непологой оболочки с матрицей 54x54 при учете смещений как жесткого целого / А. П. Николаев, Ю. В. Клочков, А. П. Киселёв // Изв. Вузов. Сев.-кавказский регион. Техн. науки. - 1999. — №2. — с. 80-84.

5. Киселёв, А. П. О функциях формы в алгоритмах формирования матрицы жесткости в треугольном конечном элементе / А. П. Николаев, Ю. В. Клочков, А. П. Киселёв // Изв. вузов, сер. Строительство. - 1999. - № 10. - С. 23-27.

6. Киселёв, А. П. К расчету оболочек на основе метода конечных элементов / А. П. Николаев, А. П. Киселёв // Вестник Российского университета дружбы народов, сер. Инж. исследования. - 2002. — № 1.-С. 107-111.

7. Киселёв, А. П. Решение задачи нелинейного деформирования оболочки на основе МКЭ при наличии особой точки / А. П. Николаев, А. П. Киселёв // Деп. в ВИНИТИ 26.12.2002, - № 2262- В2002.

8. Киселёв, А. П. Определение напряжений в стенках изотермического резервуара для транспортировки сжиженного газа в местах действия опор / А. П. Николаев, Н. Г. Бандурин, А. П. Киселёв, А. А. Сизых // Изв. Вузов. Сев.-Кавказский регион. Техн. науки. - 2005. -№2. - С. 54-57.

9. Киселёв, А. П. Расчет оболочек в трехмерной постановке с учетом геометрической нелинейности на основе МКЭ / А. П. Николаев, А. П. Киселёв // Научно-техн. журнал «Строительная механика инженерных конструкций и сооружений» - 2005. - № 1, РУДН. -С. 119-122.

10. Киселёв, А. П. Объемный конечный элемент в виде треугольной призмы с первыми производными узловых перемещений / А. П. Киселёв, А. П. Николаев // Изв. Вузов, сер. «Строительство». -2006. — № 1.-С. 13-18.

11. Киселёв, А. П. Векторная аппроксимация полей перемещений объемного шестигранного конечного элемента / А. П. Киселёв // Научно-техн. журнал «Строительная механика инженерных конструкций и сооружений» - 2007. -№ 1, РУДН. - С. 21-24.

12. Киселёв, А. П. Использование трехмерных конечных элементов в расчетах прочности с учетом геометрической нелинейности / А. П. Киселёв // Изв. Вузов, сер. «Строительство». - 2007. -№11.

13. Киселёв, А. П. Метод конечных элементов в решении трехмерных задач теории упругости /А. П. Киселёв// Научно-техн. журнал «Строительная механика инженерных конструкций и сооружений» -2007. - № 4, РУДН. - С. 11-17.

14. Киселёв, А. П. Расчет тонких оболочек на прочность в трехмерной постановке без упрощающих гипотез /А. П. Киселёв// Изв. Вузов, сер. «Строительство». - 2008. — № 1.

15. Киселёв, А. П. К расчету двух пересекающихся оболочек на основе объемных КЭ /А. П. Киселей'/ Научно-техн. журнал «Строительная механика инженерных конструкций и сооружений» — 2008. -№ 1, РУДН.

16. Киселёв, А. П. Использование трехмерных конечных элементов в расчетах прочности многослойных панелей / Н. А. Гуреева, Р. 3. Киселёва/ Научно-техн. журнал «Строительная механика инженерных конструкций и сооружений» - 2009. - № 4, РУДН. - с. 37-40.

17. Киселёв, А. П. Расчет многослойных оболочек вращения и пластин с использованием объемных конечных элементов /Р. 3. Киселёва, Н. А. Гуреева / Изв. Вузов, сгр. «Строительство». - 2010. -№ 1.-е. 106-112.

18. Киселёв, А. П. Расчет многослойной оболочки с использованием объемного конечного элемента /Р. 3. Киселёва, Н. А. Гуреева / Изв. ВолГТУ. - Волгоград. - 2010.-№ 4. - с. 125-128.

19. Kiselyev, A. The finite elements of a quadrilateral shape for analysis of shells taking into consideration a displacement of a body with rigid body modes / Yu. Klochkov, A. Nikolaev, A. Kiselyev/ Научно-техн. журнал «Строительная механика инженерных конструкций и сооружений» - 2011. - № 3, РУДН. - с. 49-59.

20. Киселёв, А. П. Определение напряжений в зоне пересечения пластин при плоском нагружении на основе МКЭ / Н. А. Гуреева, Р. 3. Киселёва, В. В. Леонтьева / Научно-техн. журнал «Строительная механика инженерных конструкций и сооружений» -2012. - № 2,

. РУДН. - с. 55-62.

21. Киселёв, А. П. Определение напряжений в зоне соединения оболочек вращения на основе МКЭ при осесимметричном нагружении / А. П. Николаев, Н. А. Гуреева, Р. 3. Киселёва, В. В. Леонтьева / Научный журнал «Фундаментальные исследования» - 2012. - № 6.

Публикации в других изданиях:

1. Киселёв, А. П. Вариант получения функций формы тетра-эдального конечного элемента с первыми производными узловых перемещений / А. П. Николаев, А. П. Киселёв, С. Н. Дейнего // Труды междунар. научн. конф. «Аетуальные проблемы механики оболочек». - Казань, 2000. - С. 218-219.

2. Киселёв, А. П. Использо вание теории упругости трехмерного тела в расчетах оболочек / А. П. Николаев, А. П. Киселёв // Сб. трудов междунар. научной конф. «Архитектура и прочностной расчет тонкостенных строительных и машине строительных конструкций сложно й формы». - Москва, 2001. - С. 58-59.

3. Киселёв, А. П. Расчет оболочек с использованием трехмерных конечных элементов в виде треугольной призмы и восьмиугольника / А. П. Николаев, А. П. Киселёв // Сб. трудов междунар. научной конф. «Архитектура и прочностной расчет тонкостенных строительных и машиностроительных конструкций сложной формы». - 2001. - С. 319-323.

4. Киселёв, А. П. Аппроксимация в методе конечных элементов в приложениях к векторным полям / А. II. Николаев, А. П. Киселёв // Материалы международной конф. «Естествознание на рубеже столетий». - т. 1, техн. науки. - Москва, 2001. - С. 54-55.

5. Киселёв, А. П. Уравнения функции формы треугольного и тетраэдального конечных элементов / А. П. Киселёв, А. П. Николаев // Научный вестник, сер. Инж. науки.- № 3.» ВГСХА, 2002. - С. 170-174.

6. Киселёв, А. П. Конечно-элементное представление тензорных полей в криволинейных системах координат/ А. П. Николаев, Ю. В. Клочков, А. П. Киселёв // Успехи современного естествознания -2003. —№ 1—С.7-11.

7. Киселёв, А. П. Сравнительный анализ результатов прочностного расчета тонкостенных оболочек с использованием двумерных и трехмерных конечных элементов / А. П. Киселёв // Матер, междунар. научно-практ. конф. «Проблемы АПК», сер. Инж. науки. - ВГСХА 2003.-С. 175-176.

8. Киселёв, А. П. Формирование матрицы жесткости конечного элемента в геометрически нелинейной постановке / А. П. Киселёв // Матер. III междунар. научно-техн. конф. «Эффективные строительные конструкции: теория и практика». - Пенза, 2004. - С. 367-370.

9. Киселёв, А. П. Деформационная теория пластичности при использовании трехмерных конечных элементов / А. П. Киселёв // Матер. междунар. научн.-практ. конф «Современные оросительные мелиорации - состояние и перспективы». - Волгоград, 2004. - С. 129-132.

10. Киселёв, А. П. Расчет трехслойных оболочек вращения методом конечных элементов в трехмерной постановке У А. П. Николаев, А. П. Киселёв, Р. 3. Киселёва/ Матер, междунар. семинара «Актуальные проблемы нелинейной механики оболочек». - г. Казань. - 2008.

11. Киселёв, А. П. Расчет слоистых плит в трехмерной постановке / А. П. Николаев, А. П. Киселёв, Р. 3. Киселёва/ Труды междунар. научно-практ. конф. «Инженерные системы-2009» - г. Москва,

РУДН. - 2009. - с. 33.

12. Киселёв, А. П. Использование современных программных комплексов, основанных на МКЭ в расчетах на прочность инженерных конструкций /А. П. Киселёв / Матер, междунар. научно-практ. конф. «Новые направления в решении проблем АПК», ВГСХА, 2010.

13. Киселёв, А. П. Расчет конструкций из разнородных материалов на основе МКЭ / А. П. Киселёв, А. П. Николаев/ Сб. трудов междунар. научно-практ. конференции «Инженерные системы - 2011». -

г. Москва, РУДН.

14. Киселёв, А. П. Расчет плосконагруженной оболочки на основе объемного конечного элемента при наличии предельной точки / А. П. Николаев, А. П. Киселёв, С. С. Марченко/ Сб. трудов междунар. научно-практ. конференции «Инженерные системы - 2011». - т. 2. -г. Москва, РУДН.

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

Киселёв Анатолий Петрович

ЛИНЕЙНОЕ И НЕЛИНЕЙНОЕ ДЕФОРМИРОВАНИЕ

УПРУГИХ ТЕЛ НА ОСНОВЕ ТРЕХМЕРНЫХ КЭ ПРИ

ВАРИАТИВНОЙ ИНТЕРПОЛЯЦИИ ПЕРЕМЕЩЕНИЙ

Автореферат диссертации на соискание ученой степени доктора технических наук

Подписано в печать 13.03.13 г. Формат 60x84 1/|6.

Усл.-печ. л. 2,0 Тираж 100. Заказ 88. ИПК ФГБОУ ВПО Волгоградский ГАУ «Нива». 400002, Волгоград, г р. Университетский, 26.

Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования Волгоградский государственный аграрный университет

На правах рукописи

05201351278

Киселёв Анатолий Петрович

ЛИНЕЙНОЕ И НЕЛИНЕЙНОЕ ДЕФОРМИРОВАНИЕ УПРУГИХ ТЕЛ НА ОСНОВЕ ТРЕХМЕРНЫХ КЭ ПРИ ВАРИАТИВНОЙ ИНТЕРПОЛЯЦИИ ПЕРЕМЕЩЕНИЙ

01.02.04 - Механика деформируемого твердого тела диссертация на соискание ученой степени доктора технических наук

научный консультант

доктор технических наук, профессор

Николаев Анатолий Петрович

Волгоград 2013

ОГЛАВЛЕНИЕ

Стр.

ВВЕДЕНИЕ.................................................... 6

1. КРАТКИЙ ОБЗОР РАЗВИТИЯ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ В ИНЖЕНЕРНЫХ РАСЧЕТАХ...................................16

2. ОСНОВНЫЕ ПОЛОЖЕНИЯ РАСЧЕТА ТОНКОСТЕННЫХ ОБОЛОЧЕК НА ОСНОВЕ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ ... .38

2.1. Основные соотношения теории тонких непологих оболочек произвольного очертания.........................................38

2.1.1. Геометрия произвольной оболочки в исходном состоянии.........38

2.1.2. Геометрия оболочки в деформированном состоянии..............42

2.1.3. Физические соотношения упругих произвольных непологих оболочек........................................................48

2.2. Последовательность выполнения основных операций метода конечных элементов .............................................50

2.3. Треугольный криволинейный конечный элемент................52

2.3.1. Геометрия элемента..........................................52

2.3.2. Выбор узловых неизвестных и аппроксимирующих функций.......54

2.3.3. Матрица жесткости..........................................61

2.4. Четырехугольный криволинейный конечный элемент...........64

2.4.1. Геометрия элемента..........................................64

2.4.2. Выбор узловых неизвестных и аппроксимирующих функций.......66

2.5. Матрица жесткости конечно-элементной модели................70

2.6. Примеры расчета.............................................71

2.7. Деформация объемного тела вращения при осесимметричном нагружении.......................................................85

2.7.1. Основные соотношения.......................................86

2.7.2. Матрица жесткости конечного элемента.........................88

2.7.3. Пример расчета..............................................89

Выводы по главе.................................................96

3. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ В РЕШЕНИИ ТРЕХМЕРНЫХ ЗАДАЧ ТЕОРИИ УПРУГОСТИ...................................97

3.1. Основные соотношения теории упругости сплошной среды.......97

3.1.1. Исходное состояние..........................................97

3.1.2. Зависимости между компонентами тензора деформаций и составляющими компонентами вектора перемещения..............99

3.1.3. Соотношения между напряжениями и деформациями для сплошной изотропной среды............................................101

3.2. Объемный конечный элемент в виде тетраэдра с четырьмя узловыми точками..........................................103

3.2.1. Геометрия элемента.........................................103

3.2.2. Выбор узловых неизвестных и аппроксимирующих функций......106

3.2.3. Матрица жесткости конечного элемента........................114

3.3. Объемный конечный элемент в виде треугольной призмы с первыми производными узловых перемещений.........................117

3.3.1. Геометрия элемента.........................................117

3.3.2. Выбор узловых неизвестных и аппроксимирующих функций......119

3.3.3. Матрица жесткости........................................125

3.4. Объемный восьмиузловой конечный элемент..................125

3.4.1. Геометрия элемента.........................................125

3.4.2. Выбор узловых неизвестных..................................128

3.4.3. Матрица жесткости.........................................130

3.5. Примеры расчета............................................131

3.6. Примеры расчета тонкостенных конструкций..................154

Выводы по главе................................................166

4. РАСЧЕТ ПЕРЕСЕКАЮЩИХСЯ ОБОЛОЧЕК НА ОСНОВЕ

ОБЪЕМНЫХ КОНЕЧНЫХ ЭЛЕМЕНТОВ........................168

4.1. Основные соотношения двух пересекающихся цилиндрических оболочек....................................................168

4.1.1. Геометрия оболочек в исходном состоянии на границе пересечения.................................................168

4.1.2. Матрица преобразования компонент вектора перемещения одной оболочки через компоненты вектора перемещения другой оболочки .. 173

4.2. Пример расчета.............................................179

Выводы по главе...............................................184

5. ВЕКТОРНАЯ АППРОКСИМАЦИЯ ПОЛЕЙ ПЕРЕМЕЩЕНИЙ ОБЪЕМНОГО ВОСЬМИУЗЛОВОГО КОНЕЧНОГО ЭЛЕМЕНТА . . .185

5.1. Матрица жесткости восьмиузлового конечного элемента с векторной аппроксимацией полей перемещений...............186

5.2. Примеры расчета............................................191

Выводы по главе................................................199

6. НЕЛИНЕЙНАЯ ТЕОРИЯ УПРУГОСТИ ТРЕХМЕРНОГО ТЕЛА В ИССЛЕДОВАНИИ НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ ИНЖЕНЕРНЫХ КОНСТРУКЦИЙ.................200

6.1. Основные соотношения нелинейной теории упругости..........200

6.1.1. Геометрия тела.............................................200

6.1.2. Суммарные деформации и напряжения после завершения /-шагов нагружения.................................................202

6.1.3. Деформации и напряжения на шаге нагружения.................203

6.2. Формирование матрицы жесткости конечного элемента на шаге нагружения.................................................205

6.3. Примеры расчета...........................................210

6.4. Нелинейное деформирование при наличии особой точки ........ 215

6.4.1. Алгоритм метода дискретного продолжения по параметру в окрестности особой точки.....................................215

6.4.2. Реализация метода дискретного продолжения по параметру в нелинейной конечно-элементной процедуре......................219

6.5. Пример расчета.............................................221

Выводы по главе................................................229

7. ОСОБЕННОСТИ РАСЧЕТА НА ПРОЧНОСТЬ СЛОИСТЫХ КОНСТРУКЦИЙ И КОНСТРУКЦИЙ ИЗ РАЗНОРОДНЫХ МАТЕРИАЛОВ НА ОСНОВЕ МКЭ В ТРЕХМЕРНОЙ ПОСТАНОВКЕ 230 7.1. Матрица жесткости конечного элемента и условия преобразования

векторов узловых неизвестных на границе контакта слоев.......230

7.2. Примеры расчета..........................................234

Выводы по главе................................................246

ЗАКЛЮЧЕНИЕ................................................247

ЛИТЕРАТУРА.................................................249

ПРИЛОЖЕНИЕ................................................287

ВВЕДЕНИЕ

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

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

МКЭ во всех его различных формулировках предусматривает следующие основные этапы расчета: замена исходной конструкции дискретной моделью; выбор узловых неизвестных и аппроксимирующих функций; формирование матрицы жесткости или податливости и вектора нагрузок; определение компонентов напряженно-деформированного состояния (НДС) исследуемой конструкции путем решения полученной СЛАУ.

Для МКЭ характерны - широкий диапазон применимости, инвариантность по отношению к геометрии конструкции и механическим

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

Понятие конечных элементов (КЭ) было впервые введено М. Тернером, Р. Клафом, X. Мартином, JI. Топпом. Дальнейшее развитие метода отражено в работах зарубежных и отечественных исследователей Дж. Аргириса, E.JI. Вильсона, М.Р. Айронса, Р.У. Клафа, У.М. Дженкинса, O.K. Зенкевича, A.B. Александрова, A.M. Масленникова, J1.A. Розина, H.H. Шапошникова, В.А. Постнова, И.Я. Хархурима, Д. В. Вайнберга, A.C. Сахарова и др.

Литература, посвященная теории и реализации МКЭ, довольно обширна (в последние годы изданы книги [15, 31, 204, 206]). Особо следует отметить книги O.K. Зенкевича [54] и В.А. Постнова, И.Я. Хархурима [179], в которых исчерпывающе изложена теория метода и дано ясное представление его реализации на ЭВМ. Однако анализ современных научных публикаций [15, 204], посвященных вопросам исследования процессов деформирования различных конструкций на основе МКЭ, позволяет заключить, что остается ряд весьма важных проблем, которые требуют иного подхода или принципиально нового решения. Наиболее сложными проблемами в МКЭ являются учет смещения конструкции как жесткого целого и использование объемных высокоточных конечных элементов в расчетах геометрически линейных и нелинейных задачах механики деформируемого твердого тела. Поэтому задача дальнейшего развития теории линейного и нелинейного деформирования инженерных конструкций на основе МКЭ является, достаточно актуальной и представляет собой как теоретический, так и практический интерес.

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

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

Для достижения этой цели были поставлены и решены следующие задачи:

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

• Разработать КЭ четырех-, пяти- и шестигранной формы, за узловые неизвестные которых выбираются компоненты вектора перемещений и их первые производные (при размерах матриц жесткости 48x48, 72x72, 96x96 соответственно), с векторной аппроксимацией полей перемещений.

• Разработать аппроксимирующие функции для КЭ в виде тетраэдра с использованием полных трехмерных полиномов и функции формы треугольной призмы на основе использования полиномов Эрмита в комбинации с полными двумерными полиномами.

• Разработать трехмерные КЭ для исследования НДС оболочек произвольной толщины с учетом поперечных и сдвиговых деформаций без привлечения упрощающих гипотез.

• Разработать условия сочленения КЭ с узловыми неизвестными в виде компонент вектора перемещений и их первых производных для двух пересекающихся оболочек вращения.

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

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

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

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

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

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

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

интегрирования по объему КЭ, реализованные в алгоритмах формирования матриц жесткости КЭ с использованием функционала Лагранжа.

Основные научные положения, выносимые на защиту:

•новый вариант получения функций формы, основанный на векторном способе аппроксимации полей перемещений трехмерных конечных элементов, для учета смещения конструкции как жесткого целого;

•вариант получения функций формы и матриц жесткости четырех- и пятигранных КЭ, за узловые неизвестные которых приняты компоненты вектора перемещений и их первые производные;

• алгоритмы формирования матриц жесткости четырех-, пяти- и шестигранных объемных КЭ на основе разработанного векторного способа аппроксимации полей перемещений и их первых производных;

• алгоритмы деформирования оболочки как трехмерного тела в геометрически нелинейной постановке при шаговом способе нагружения;

• методика определения НДС пересекающихся оболочек и оболочек слоистой структуры с использованием трехмерных конечных элементов;

• алгоритм дискретного продолжения по параметру нагружения в окрестности предельной точки деформирования оболочек с использованием шестигранного конечного элемента в геометрически нелинейной постановке.

Научная новизна диссертационной работы заключается в следующем

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

2. Разработаны на основе векторного способа аппроксимации полей перемещений алгоритмы формирования матриц жесткости КЭ четырех-, пяти-и шестигранной формы с узловыми неизвестными в виде компонент вектора перемещений и их первых производных для определения НДС однородных и разнородных конструкций, в том числе и в зонах их сочленения.

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

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

4. Разработан и реализован алгоритм дискретного продолжения решения по параметру нагружения на основе КЭ шестигранной формы в окрестности предельной точки деформирования оболочки в геометрически нелинейной постановке.

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

Программные модули разработанных КЭ использовались в расчетах на прочность конструктивных элементов нефтегазового и химического оборудования в ОАО «ОРГЭНЕРГОНЕФТЬ», ОАО

«ВОЛГОГРАДНЕФТЕМАШ», ОАО «РЕМГАЗКОМПЛЕКТПОСТАВКА», ОАО «ЛУКОЙЛ- Волгограднефтепереработка».

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

Апробация работы. Материалы диссертационной работы докладывались и обсуждались:

- на ежегодных научно-практических конференциях