Развитие метода конечных элементов в исследованиях линейного и нелинейного деформирования оболочек как двумерных и трехмерных упругих тел

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

Киселёв Анатолий Петрович РАЗВИТИЕ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ В ИССЛЕДОВАНИЯХ ЛИНЕЙНОГО И НЕЛИНЕЙНОГО ДЕФОРМИРОВАНИЯ ОБОЛОЧЕК КАК ДВУМЕРНЫХ И ТРЕХМЕРНЫХ УПРУГИХ ТЕЛ 05.23.17 – Строительная механика

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

Волгоград 2008 2

Работа выполнена в Федеральном государственном образовательном учреждении высшего профессионального образования Волгоградской государственной сельскохозяйственной академии Научный консультант – доктор технических наук, профессор Николаев Анатолий Петрович, ФГОУ ВПО Волгоградская государственная сельскохозяйственная академия Официальные оппоненты – доктор технических наук, профессор Шапошников Николай Николаевич, ГОУ ВПО Московский государственный университет путей сообщений – доктор технических наук, профессор Овчинников Игорь Георгиевич, ГОУ ВПО Саратовский государственный технический университет – доктор технических наук, профессор Ким Алексей Юрьевич, ФГОУ ВПО Саратовский государственный аграрный университет им. Н. И. Вавилова Ведущая организация – ГОУ ВПО Московский государственный строительный университет.

Защита состоится 22 мая 2008 года в 10.00 часов в аудитории Б-203 на заседании диссертационного совета Д 212.026.01 при Волгоградском государственном архитектурно-строительном университете по адресу:

400074, г. Волгоград, ул. Академическая, 1.

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

Автореферат разослан _ апреля 2008 года.

Ученый секретарь диссертационного совета Кукса Л. В.

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

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

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

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

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

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

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

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

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

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

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

6. На основе разработанной теории реализован алгоритм формирования матрицы жесткости объемного элемента шестигранной формы в геометрически нелинейной постановке;

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

Практическая ценность диссертационной работы заключается:

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

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

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

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

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

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

- векторный способ аппроксимации полей перемещений объемного шестигранного конечного элемента;

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

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

- методика определения напряженно-деформированного состояния пересекающихся оболочек с использованием объемных конечных элементов;

- алгоритм учета смещения конструкции как жесткого целого с использованием восьмиузлового шестигранного конечного элемента;

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

Достоверность научных положений обеспечивалась: корректной математической постановкой задач при использовании теории упругости, методов вычислительной математики и векторного анализа;

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

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

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

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

- на региональных конференциях молодых исследователей (г.

Волгоград, 1999-2001 г.г.);

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

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

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

- на международной научно-технической конференции «Информационные технологии в образовании, технике и медицине» (г.

Волгоград, 2000 г.);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Срединная поверхность произвольной непологой оболочки в исходном состоянии задавалась в декартовой системе координат х1, х2, х3 векторным уравнением R o = x k ( ) ik, (k=1,2,3;

=1,2), (1) где R o -радиус-вектор описывающий срединную поверхность оболочки;

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

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

Поочередным дифференцированием (1) по криволинейным координатам можно получить ковариантные векторы локального базиса xk a0 = R,0 = ik. (2) В (2) нижние символы после запятой обозначают дифференцирование по соответствующей координате, а индексы 1,2 соответствуют направлениям 1 и 2.

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

Ковариантные компоненты метрического тензора определяются скалярными произведениями базисных векторов a0 = a0 = a0 a0. (3) Контравариантные компоненты метрического тензора и контравариантные векторы взаимного базиса a0 0 a = 0 ;

a = a a.

(4) a Ковариантные и смешанные компоненты тензора кривизны b0 = a 0 b0.

b0 = R,0 a 0 ;

(5) Дифференцированием (2) по криволинейным координатам можно получить производные векторов локального базиса.

Положение точки, расположенной на расстоянии от срединной поверхностиоболочки, определяется радиус-вектором R 0 = R 0 + a 0, (6) где - является переменной величиной, определяющей положение точки в пределах толщины оболочки.

В результате дифференцирования (6) можно получить ковариантные векторы базиса g 0 = R,0 = ( R 0 + a 0 ), = a0 b0 a0.

(7) Ко- и контравариантные компоненты метрического тензора в исходном состоянии равны g0 = g0 g 0 ;

g 0 = g 0 / g 0 ;

g 0 = g 0 = g 0 / g 0 (при ), (8) g 0 = g11 g 22 g12 g 21.

0 0 0 где Радиус-вектор, определяющий положение точки срединной поверхности после деформации будет представлен суммой векторов R = R0 + V, (9) где V - вектор, характеризующий перемещение точки срединной поверхности, определяемый в исходном базисе векторов 0 V = ua10 + va2 + wa 0. (10) Величины u, v, w являются компонентами вектора перемещений вдоль координатных линий 1, 2 и нормали, соответственно.

В результате дифференцирования (10) можно получить производные вектора перемещений V по криволинейным координатам 1, 2.

Ковариантные вектора локального базиса в деформированном состоянии определяются из выражения a = R, = a0 + V,. (11) Положение точки, расположенной на расстоянии от срединной поверхности деформированной оболочки, с учетом гипотезы Кирхгофа о прямолинейности нормалей в исходном и деформированном состояниях, определяется радиус-вектором R = R + a, (12) Дифференцированием (12) по криволинейным координатам можно получить выражения для базисных векторов g в деформированном состоянии g = R, = a0 + V, + a,, (13) где a, производная орта нормали.

Ковариантные компоненты тензора деформаций для точки срединной поверхности определяются с использованием соотношений механики сплошной среды = ( a a0 ).

(14) Деформации в точке, расположенной на расстоянии от срединной поверхности оболочки будут равны = ( g g 0 ) / 2. (15) Входящие в (15) ковариантные компоненты метрического тензора в деформированном g и исходном g состояниях определяются скалярными произведениями соответствующих базисных векторов.

Принимая во внимание (14) и (15) можно записать = +, (16) где ( ) = a0 V, + V, a 0 / 2;

= ( a0 a,0 + V, a,0 + a, a 0 + a, V, + (17) + a0 b0 a0 + b0 a0 a 0 ) / 2.

Контравариантные компоненты тензора напряжений в любой точке оболочки выражаются через ковариантные компоненты тензора деформаций соотношением = I1 ( ) g + 2µ g g, (18) где I1 ( ) = g – первый инвариант тензора деформаций;

, µ - параметры Ламе.

На основе соотношений (18) можно записать матричное выражение { } = [ C ]{ }, (19) где { } { } = 12 T,, 11 - вектор напряжений;

{ } = { } T, 22,2 12 - вектор деформаций произвольного слоя оболочки;

[C] – матрица упругости, размером 33.

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

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

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

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

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

Положение произвольной точки сплошной трехмерной среды в исходном недеформированном состоянии в декартовой системе координат может быть определено векторным уравнением (рис. 1) R 0 = x k ( m ) ik, (k,m =1,2,3), (20) где R - радиус-вектор;

x k, ik - соответственно переменные и орты декартовой системы координат;

m – переменные криволинейной системы координат, отличной от декартовой.

M V R M x3 R i x i i x Рис. 1. Перемещение произвольной точки сплошной среды в результате деформации Из (20) путем поочередного дифференцирования по криволинейным координатам m определяются ковариантные векторы локального базиса в исходном состоянии 0 0 R ai = R,i =. (21) i Ковариантные компоненты метрического тензора определяются скалярными произведениями базисных векторов aij = ai0 a 0.

(22) j Производные векторов локального базиса (21) определяются из выражения a 0 0 0i a j,k = = ai Г jk.

j (23) k Символы Кристоффеля второго рода из (23) определяются через символы Кристоффеля первого рода Г 0i = a 0il Г l0, jk = Г kjl, (24) jk где 1 alj alk akj 0 = k+ = Г i0,kj - символы Кристоффеля первого рода.

Г 2 x x j xl l, jk Положение точки сплошной среды после деформации будет определяться радиус-вектором R, определяемым как сумма векторов R = R0 + V. (25) Ковариантные векторы локального базиса в деформированном состоянии определяются путем дифференцирования (25) ak = R,0 + V,k = ak0 + V,k. (26) k В результате скалярного произведения (26) определяются компоненты метрического тензора ( )( ) aij = ai a j = ai0 + V,i a 0 + V, j = aij + ai0V, j + a 0V,i + V,iV, j.

(27) j j Вектор перемещений произвольной точки сплошной среды может быть определен 0 V = a10u + a2 v + a30 w, или V = u m am, (28) m где u - проекции вектора перемещений на направления криволинейной системы координат m (m=1,2,3).

Производные вектора перемещений V в криволинейной системе координат m am um 0 1 2 = f k a1 + f k a2 + f k a3, V,k = am + u (k=1,2,3), (29) k k 1 2 где f k, f k, f k многочлены, содержащие компоненты вектора перемещений и их первые производные.

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

1 (30) Контравариантные компоненты тензора напряжений для упругого тела определяются через ковариантные компоненты тензора деформаций из уравнения механики сплошной среды ij = Aijmn mn, (i,j,m,n=1,2,3) (31) Aijmn = a 0ij a 0 mn + µ ( a 0im a 0 jn + a 0in a 0 jm ) ;

где (32) 0 ij a - контравариантные компоненты метрического тензора полученные для исходного состояния.

На основе (31) можно записать матричную зависимость между вектором контравариантных компонент тензора напряжений и вектором ковариантных компонент тензора деформаций { } = [ C ]{ }, (33) где { } = { 11, 22, 33, 12, 23, 31 } T вектор напряжений;

{ } = {,,,2,2,2,} вектор деформаций;

T 11 22 33 12 23 [ C ] матрица упругости размером 66.

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

Объемный конечный элемент в виде произвольного тетраэдра с четырьмя узловыми точками i,j,k,l или 1,2,3,4 в глобальной системе координат х1, х2, х3 (рис.2) представляется в локальной системе координат,,, изменяющимися в пределах от 0 до +1 (рис.3).

Соотношения между глобальной системой координат х1, х2, х конечного элемента и локальной системой координат определяется выражением x = xi + ( xj xi ) + ( xk xi ) + ( xl xi ), (34) где xi, xj, x, xl - координаты узловых точек в системе х1, х2, х3;

k - поочередно принимает значения 1,2,3.

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

Рис. 2. Объемный конечный элемент в виде произвольного тетраэдра в глобальной системе координат х1, х2, х l (0;

0;

1) i (0;

0;

0) k (0;

1;

0) j (1;

0;

0) Рис. 3. Тетраэдр в локальной системе координат,, В каждом узле конечного элемента неизвестными величинами выбирались компоненты вектора перемещений u, v и w в направлениях осей координат х1, х2 и х3 их первые производные u,1, u,2, u,3, v,1, v,2, v,3, w,1, w,2, w,3.

Вектор перемещений узловых точек конечного элемента в глобальной системе координат имел вид {V у }Т = {u i, u j, u k, u l, u,i1, u,1j, u,k1, u,l1,u,i2, u,j2, u,k2, u,l2, (35) 1 X u,i3, u,3, u,k3, u,l3, v i,..., wi,...} j Для аппроксимации поля перемещений точек внутри конечного элемента использовался полный трехмерный полином третьего порядка q = k1 + k 2 + k 3 + k 4 + k 5 2 + k 6 + k 7 + k8 2 + k9 + k10 + k11 3 + k12 2 + k13 2 + k14 2 + k15 + k16 + k17 3 + (36) k18 + k19 + k 20, 2 2 где k1,k2,…,k20 – искомые коэффициенты.

При использовании конечного элемента в виде тетраэдра с первыми производными число неизвестных узловых перемещений будет равно шестнадцати. Поэтому для корректного определения неизвестных коэффициентов аппроксимирующих полиномов в вектор перемещений, кроме узловых компонент вектора перемещений и их первых производных были дополнительно включены четыре смешанные производные второго 1 2 3 порядка h,n, h,n, h,n, h,n по направлению нормали к грани 2-3-4, определяемые выражением h,in = h,i,n + h,i,n + h,i,n, (37) где h i = q,i + q,i + q,i ;

q – любая из трех компонент вектора перемещений u, v и w;

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

i=1,2,3,4 – номера узлов конечного элемента.

Вектор узловых неизвестных в локальной системе координат будет иметь вид {q } = {q, q, q k, q l, q,i, q,j, q,k, q,l, q,i, q,j, q,k, q,l, q,i, q,j, q,k, q,l, h,in, h,n, h,kn, h,ln }.

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

Производные h,, h,, h, из (37), пользуясь методом конечных разностей, можно представить выражениями h,1 = h 2 h1 ;

h,1 = h 3 h1 ;

h,1 = h 4 h1 ;

h,2 = h 2 h1 ;

h,1 = h 3 h 2 ;

h,1 = h 4 h 2 ;

(38) h,3 = h 2 h 3 ;

h,1 = h 3 h1 ;

h,1 = h 4 h 3 ;

h,4 = h 2 h 4 ;

h,1 = h 3 h 4 ;

h,1 = h 4 h1.

{} * Вектор узловых неизвестных q у содержащий двадцать компонент с учетом (39) можно выразить через вектор с шестнадцатью компонентами {q } = [ ]{ } * M qу, (39) у 20 x16 16 x 20 x где {q } = { q, q } Т, q k, q l, q,i, q,j, q,k, q,l, q,i, q,, q,k, q,l, q,i, q,j, q,k, q,l.

i j j у Используя (39) выражение для любой компоненты вектора перемещения можно представить q = { } T { q у }, (40) 1x16 16 х где { } T - вектор-строка шестнадцати аппроксимирующих функций.

1 x При представлении трехмерного континуума дискретной моделью часто возникает необходимость использования элементов в виде треугольной призмы (рис. 4). Для выполнения численного интегрирования по объему конечного элемента произвольная треугольная призма в глобальной системе координат х1, х2, х3 отображается на треугольную призму в локальной системе координат,, (рис. 5) с координатами узлов i(0;

0;

-1), j(1;

0;

-1), k(0;

1;

-1), l(0;

0;

1), m(1;

0;

1), n(0;

1;

1). Зависимость между координатами х1, х2, х3 и локальными,, конечного элемента определяется соотношением ) (1 ) + (1 ) + (1 ) + x = xi (1 x j xk 2 2 (41) (1 + ) + x ( 1 + ) + x (1 + ), + xl (1 ) m n 2 2 хi,..., xn - координаты узловых точек в системе х1, х2, х3;

где =1,2,3.

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

За узловые неизвестные конечного элемента выбирались компоненты вектора перемещений u, v и w в направлениях осей координат х1, х2 и х3 их первые производные u,1, u,2, u,3, v,1, v,2, v,3, w,1, w,2, w,3 в угловых точках. Таким образом, вектор перемещений узловых точек конечного элемента в глобальной системе координат будет иметь вид {V } = {u, u Т i j, u k, u l, u m, u n, u,i1, u,1, u,k, u,l1, u,m, u,n, j у 1 1 (42) 1 X u, u, u, u, u, u, u, u, u, u, u, u, v,..., w,...}.

i j k l m n i j k l m n i i,2,2,2,2,2,2,3,3,3,3,3, l n i m k x x2 j x Рис. 4. Конечный элемент в виде треугольной призмы в системе координат х1, х2, х Рис. 5. Конечный элемент в виде треугольной призмы в системе координат,, l (0;

0;

1) n m (0;

1;

0) (1;

0;

0) i (0;

0;

1) k j Для аппроксимации поля перемещений внутренней точки конечного элемента через его узловые значения в треугольной призме используются два типа функций формы.

Для аппроксимации полей перемещений внутренних точек в координатах, конечного элемента использовался полный алгебраический полином третьей степени, который в этом случае будет содержать десять членов g (, ) = k1 + k 2 + k3 + k 4 2 + k5 + k 6 2 + k 7 3 + (43) + k8 2 + k9 2 + k10 3, где k1,…,k10 – неизвестные коэффициенты.

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

Вектор узловых неизвестных в этом случае будет иметь следующую структуру {q } = {q, q, q k, qi, q j, qk, qi, qj, qk, qi }, i j (44) у 1 x В результате количество неизвестных оказывалось равным числу уравнений полученной системы.

Перемещение внутренней точки элемента можно выразить через узловые неизвестные матричным соотношением q = { g} { q * }.

T (45) у 1 х10 1 х i Смешанную производную перемещения первого узла q, с использованием метода конечных разностей в локальной системе координат можно выразить через первые производные узловых перемещений q,i = ( q,j q,i + q,k q,i ).

(46) С использованием (46) вектор-столбец узловых перемещений с десятью компонентами {q у } можно определить через вектор-столбец T 1ч {q } = {q, q, q, q, q, q, q, qj, qk }, имеющий только девять составляющих T i j k i j k i у 1x компонент {q } = [ L]{q }, (47) у у 10 x 9 9 x 1 x где [L] – матрица преобразования.

Принимая во внимание (47) выражение (45) примет вид q = { } { q у }, T (48) 1x9 9 x = { g} [ L].

{ } Т T где 1 x10 10 x Перемещение внутренней точки (48) в развернутом виде будет определяться выражением q = G1 (, ) q i + G2 (, ) q j + G3 (, ) q k + + G4 (, ) qi + G5 (, ) q j + G6 (, ) qk + (49) + G7 (, ) q + G8 (, ) q + G9 (, ) q, i j k где G1, …, G9 - аппроксимирующие функции для треугольной области конечного элемента.

С учетом (49) перемещение внутренних точек призматического конечного элемента с треугольным основанием в локальной системе координат будет определяться через перемещения узловых точек q = G1 (, ) h1 ( ) q i + G2 (, ) h1 ( ) q j + G3 (, ) h1 ( ) q k + + G1 (, ) h2 ( ) q l + G2 (, ) h2 ( ) q m + G3 (, ) h2 ( ) q n + + G4 (, ) h1 ( ) q,i + G5 (, ) h1 ( ) q,j + G6 (, ) h1 ( ) q,k + (50) + G4 (, ) h2 ( ) q,l + G5 (, ) h2 ( ) q,m + G6 (, ) h2 ( ) q,n + + G7 (, ) h1 ( ) q,i + G8 (, ) h1 ( ) q,j + G9 (, ) h1 ( ) q,k + + G7 (, ) h2 ( ) q,l + G8 (, ) h2 ( ) q,m + G9 (, ) h2 ( ) q,n + + G1 (, ) h3 ( ) q,i + G2 (, ) h3 ( ) q,j + G3 (, ) h3 ( ) q,k + + G1 (, ) h4 ( ) q,l + G2 (, ) h4 ( ) q,m + G3 (, ) h4 ( ) q,n, где h1 ( ), h2 ( ), h3 ( ), h4 ( ) полиномы Эрмита третьей степени h1 ( ) = ( 3 3 + 2);

h2 ( ) = ( 3 3 2 );

1 4 (51) h3 ( ) = ( + 1);

h4 ( ) = ( + 1).

13 2 4 В качестве объемного шестигранного конечного элемента выбирался произвольный восьмиузловой элемент с узловыми точками 1(i),2(j),3(k),4(l),5(m),6(n),7(p),8(h) (рис.6). Положение узловых точек конечного элемента определялось тремя координатами – х1, х2, х3.

Для выполнения численного интегрирования по объему конечного элемента произвольный восьмиузловой элемент отображался на куб в системе координат,, (рис.7). Локальные координаты кубического элемента изменялись в пределах -1,, 1.

Зависимость между глобальными координатами х1,х2,х3 и локальными,, конечного элемента определяется выражением (1 ) (1 ) (1 ) + x (1 + ) (1 ) (1 ) + x = xi j 8 (1 + ) (1 + ) (1 ) (1 ) (1 + ) (1 ) + + xk + xl 8 (52) (1 ) (1 ) (1 + ) (1 ) (1 ) (1 + ) + xm + xn + 8 (1 ) (1 ) (1 + ) + x (1 ) (1 ) (1 + ), + xp h 8 где =1,2,3;

xi, xj, xk, xl, xm, xn, x p, xh - координаты узловых точек конечного элемента в глобальной системе координат.

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

h m p l n i x k j x x Рис. 6. Произвольный восьмиузловой конечный элемент в системе координат х1, х2, х m h n p l i j k Рис. 7. Восьмиузловой конечный элемент в системе координат,, Вектор перемещений узловых точек конечного элемента в глобальной системе координат будет иметь вид {V уГ } = {u1, u 2, u 3, u 4, u 5, u 6, u 7, u 8, u,11, u,21, u,31, u,41, u,51, u,61, u,71, u,81, (53) 96 X }.

1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 u, 2, u, 2, u, 2, u, 2, u, 2, u, 2, u, 2, u, 2, u,3, u,3, u,3, u,3, u,3, u,3, u,3, u,3, v,..., w,...

Перемещения внутренних точек конечного элемента определяются через перемещения узловых точек с помощью аппроксимирующих функций q = { } { q у }, T (54) 1 x 32 32 х { } T где - вектор-строка тридцати двух аппроксимирующих функций, 1 x основанных на полиномах Эрмита третьей степени.

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

В качестве примера расчета исследовалось напряженно деформированное состояние объемного тела при действии на него сосредоточенной силы (рис. 8). Были приняты следующие исходные данные:

величина сосредоточенной силы Р=39.2104 Н;

сторона куба а=3.05 м. Расчет выполнялся в трех вариантах: в первом расчет выполнялся с использованием конечного элемента в виде тетраэдра, во втором использовался конечный элемент в виде треугольной призмы, в третьем результаты расчета получены при использовании восьмиузлового конечного элемента.

x E H F G A D x B C a x Рис. 8. Задача о действии сосредоточенной силы на объемное тело на примере использования конечного элемента в виде треугольной призмы Результаты расчетов внесены в таблицу № 1. В таблице приводятся величины вертикального перемещения куба под силой в метрах при различном количестве конечных элементов дискретной сетки. При дискретизации объема отдельными тетраэдрами и треугольными призмами вначале пространство разбивалось на отдельные объемные прямоугольники, которые затем разбивались на пять тетраэдров или две призмы. В 3-й колонке таблицы приводятся результаты, полученные с использованием тетраэдальных конечных элементов, в 4-й с использованием элемента в виде треугольной призмы, в 5-й с использованием шестигранных конечных элементов. Во всех вариантах проверялась сходимость вычислительного процесса с различным числом точек интегрирования. Наблюдается сходимость вычислительного процесса при увеличении числа дискретных элементов, также видно, что результаты расчетов, полученные различными конечными элементами, практически совпадают, что говорит о достаточной точности разработанных алгоритмов.

Таблица № Размер Величина вертикального перемещения куба под силой w10-6м п.п. конечно Тетраэдр Треугольная призма Восьмиугольный элементной сетки Элемент 1 1х1х1 3.8 3.3 3. 2 2х2х2 4.8 6.8 6. 3 3х3х3 6.9 10.5 10. 4 4х4х4 9.2 14.2 13. 5 5х5х5 11.5 17.9 17. 6 5х5х6 11.9 18.0 17. 7 5х5х10 12.3 18.4 18. 8 6х6х6 13.8 21.5 21. Рассчитывалась цилиндрическая тонкостенная оболочка, представленная объемными восьмиузловыми конечными элементами.

Оболочка загружалась двумя равными по величине и противоположно направленными сосредоточенными силами, приложенными по концам одного диаметра (рис. 9). За исходные данные приняты: радиус срединной поверхности R=4.953 м;

длина цилиндра L=10.35 м;

модуль упругости материала Е=10.5 106 кПа;

коэффициент Пуассона =0.31;

толщина стенки цилиндра h=0.094 м;

величина сжимающей силы Р=100 кН. С учетом трех плоскостей симметрии, расчет сводился к исследованию восьмой части оболочки. Результаты расчета представлены в таблице № 2. В таблице приводятся величины прогиба оболочки в точке приложения сосредоточенных сил в зависимости от количества конечных элементов дискретной сетки. Результаты расчета сравнивались со значениями, полученными в главе 2 с использованием двумерных дискретных элементов.

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

Рис. 9. Цилиндрическая оболочка, сжатая в середине пролета сосредоточенными силами Таблица Количество Прогиб оболочки Количество Прогиб оболочки под № п.п. двумерных под силой, м объемных силой, м конечных конечных элементов элементов 1 1x1 0,11046 1x1x1 0, 2 2x2 0,11765 3x3x3 0, 3 3x6 0,11817 1x3x3 0, 4 5x14 0,11810 1x10x10 0, 5 6x20 0,11807 1x12x12 0, 6 1x14x14 0, В главе также приводятся примеры расчета других тонкостенных оболочек, представленных объемными конечными элементами. Приведен пример расчета защемленной панели, загруженной линейной нагрузкой результаты расчета сравнивались с аналитическим решением. Приведены примеры практических задач по расчету на прочность корпуса шарового крана по заказу ОАО «Волгограднефтемаш» и шарового резервуара с несовершенством геометрии по заказу Самарского филиала ОАО «Оргэнергонефть».

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

Для описания геометрии двух пересекающихся под прямым углом цилиндрических оболочек вводятся две системы координат x,y,z для основной оболочки с радиусом r и толщиной t и x|,y|,z| для примыкающего патрубка с радиусом r| и толщиной t| (rr|) (рис. 10). Здесь и далее символы без штриха будут относиться к параметрам основной оболочки, со штрихом к патрубку.

Рис. 10. Пересечение цилиндрических оболочек Положение произвольных точек основной оболочки и примыкающей, будет определяться в криволинейных системах координат x,, r и x/, /, r/ радиус- векторами, соответственно для основной оболочки R = xi + r sin j + r cos k (55) и для примыкающей оболочки R / = x / i / + r / sin / j / + r / cos / k /. (56) Дифференцированием (55) и (56) по криволинейным координатам можно получить ковариантные векторы локального базиса основной и примыкающей оболочек.

В зоне примыкания оболочек можно получить матричное выражение векторов локального базиса основной оболочки через векторы локального базиса примыкающей {a 0 } = [ A]{a 0 / } (57) где [ A] - матрица преобразования.

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

гT (58) у,1,2,3,1,2,3,1,2, { u } = { u,u,u,u, v, v,v, v, w, w, w, w }, /г T / / / / / / / / / / / / (59) у,1,2,3,1,2,3,1,2, где цифрам 1,2,3 соответствуют направления координат x,, r и x/, /, r/.

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

Вектор перемещений точки, расположенной на грани сопряжения можно представить компонентами, отнесенными к базисам векторов { a } и {a 0 / } V = ua10 + va2 + wa3 ;

V / = u / a10 / + v / a2 0 / + w / a3 0 /, 0 (60) На основе (60) компоненты вектора примыкающей оболочки выражаются через компоненты вектора основной оболочки u / = vr sin w cos ;

1 1 v / = u / sin / + vr cos / cos / + w sin / cos / ;

(61) r r r w / = u cos / + vr cos sin / + w sin sin /.

Производные вектора перемещений точки, расположенной на грани пересечения оболочек по направлениям x,, r примыкающей оболочки в / / / системе координат x,, r основной оболочки определяются по формуле V, /j = a 0 / gradV, (62) j 01 02 где gradV = a V, x + a V, + a V,r - вектор градиент векторного поля перемещения;

(j=1,2,3).

Из (62) получим ( ).

/ V, j = a 0 / a 01V, x + a 02V, + a 03V,r (63) j На основе (63) получаются зависимости производных компонент вектора перемещений на грани пересечения основной оболочки с примыкающей.

Вектор узловых неизвестных примыкающей оболочки с использованием (61), (63) можно выразить через вектор узловых неизвестных основной оболочки {V } }.

= [ z ]{V уГ /Г T (64) у Матрица жесткости конечного элемента примыкающей оболочки на границе примыкания к основной получается преобразованием [ K ] = [ z] [ K ][ z].

T */ / (65) В качестве примера расчета исследовалось напряженно деформированное состояние в зоне сочленения двух ортогонально пересекающихся цилиндрических оболочек, находящихся под действием внутреннего давления. Были приняты следующие исходные данные: r=0. м, r/=0.064 м, t=0.00254 м, t/=0.00127 м, модуль упругости материала Е=2108 кПа, коэффициент Пуассона v=0.3, величина внутреннего давления q=3.445 кПа. Вследствие наличия двух плоскостей симметрии в расчете рассматривалась четвертая часть конструкции. Изменением количества конечных элементов дискретной модели проверялась сходимость вычислительного процесса.

, даН/см, даН/см 1200 х/, см 0,0 1,0 2,0 3, - - - - х/, см - 0, 0 0, 5 1, 0 1, 5 2, - - Рис. 11 Рис., даН/см, даН/см х, см 0,0 1,0 2,0 3,0 4,0 5, - х,0см 0, 0 1, 0 2, 0 3, 0 4, 0 5, - - Рис.13 Рис. На рисунках 11-14 представлены расчетные кривые изменения меридиональных (сплошные линии) и кольцевых (пунктирные линии) напряжений в продольном сечении рассчитываемой конструкции. На рисунках 11 и 12 изображены графики изменения напряжений во внутренних и наружных волокнах примыкающего цилиндра, соответственно, а на рисунках 13 и 14 – то же самое для основного цилиндра. Экспериментальные значения кольцевых напряжений обозначены светлыми кружками, меридиональных – темными кружками.

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

В пятой главе разработан алгоритм формирования матрицы жесткости шестигранного конечного элемента с векторной аппроксимацией полей перемещений для исследований напряженно-деформированного состояния трехмерных конструкций при значительном смещении их как жесткого целого. За узловые неизвестные конечного элемента в этом случае выбирались векторы перемещений узловых точек и их первые производные {}{ } T V уЛ = V i, V j, V k, V l, V m, V n, V p, V h, V i,..., V i,..., V i,..., ;

1 (66) {V } = {V, V } Г i j k l m n p h i i i T, V, V, V, V, V, V, V,1,..., V, 2,..., V,3,....

у 1 Вектор перемещений внутренней точки объемного восьмиузлового конечного элемента с использованием полиномов Эрмита третьей степени аппроксимируется через вектор-столбец узловых перемещений { }, V = { } V уЛ T (67) 1 32 32 где { } T - вектор-строка аппроксимирующих функций.

Компоненты вектора перемещений в локальной системе координат,, выражаются через компоненты вектора перемещений в глобальной системе координат,, 1 2 {} { }, V уЛ = [ L] V уГ (68) 32 32 32 32 где [ L] - матрица преобразования координат.

С использованием (29), (30) можно получить зависимости для векторов узловых точек конечного элемента, которые в матричном виде представляются следующим образом {V } = [ A] { f }, Г Г (69) у у 32 32 1 96 [] где A – матрица, ненулевыми элементами которой являются векторы базисов узловых точек конечного элемента;

{ f } = {u, v, w, u,..., u k,..., u l,..., u m,..., u n,..., u p,..., u h,..., f11i,..., f12i,..., f 33i,...}. (70) Г i i i j у Учитывая равенства (68) и (69), соотношение (67) можно привести к виду [ A] [ G] { f }, V = { }T Г (71) у 32 96 96 96 96 1 где [ A] = [ A] [ G].

[ L] (72) [] 32 32 32 96 32 96 96 Векторы базисов узловых точек, определяющие матрицу A, можно выразить через базисные векторы внутренней точки конечного элемента, в [] результате матрицу A можно представить [] 0 A = [ A1 ] a10 + [ A2 ] a2 + [ A3 ] a3, (73) С учетом (73) зависимость (71) может быть записана в виде V = { } ( [ A1 ] a10 + [ A2 ] a2 + [ A3 ] a30 )[ G ]{ f уГ }.

0 T (74) Приравнивая правые части выражений (74) и (28) можно получить формулы для каждой компоненты вектора перемещений произвольной точки } T [ A1 ][ G ]{ f уГ };

} T [ A2 ][ G ]{ f уГ };

} T [ A3 ][ G ]{ f уГ }.

u = { v = { w = { (75) Последовательным дифференцированием (75) по криволинейным координатам и приравниванием правых частей полученных равенств и (29), определяются интерполяционные зависимости производных компонент вектора перемещений f mk = {, m } [ Am ][ G ]{ f уГ };

T (76) где,m =,,m +,,m +,,m ;

k, m=1,2,3.

{f } Г Столбец неизвестных можно представить выражением у { f уГ } = [ N ]{u уГ }, (77) {u } = {u, u, u k, u l, u m, u n, u p, u h, u,i1,..., u,i2,..., u,i3,..., v i,..., wi,...};

- вектор ГТ i j где у строка компонент вектора узловых перемещений;

[ N ] - матрица перехода, размером 9696.

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

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

В качестве одного из них рассматривалась оболочка, срединная поверхность которой представляет собой форму усеченного эллипсоида вращения. Оболочка нагружалась равномерно распределенным внутренним давлением интенсивности q=4.0 мПа (рис. 15). Были приняты следующие исходные данные: А=1.0 м;

В=0.5 м;

D=0.9 м;

Е=2105 мПа;

=0.3;

h=0.01 м, (А и В являются главными полуосями эллипса).

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

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

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

Рис. 15. Усеченный эллипсоид вращения, загруженный внутренним давлением Таблица Кольцевые напряжения к, мПа Величина жесткого смещения, м Вариант расчета I II 0.00 126.2 125. 0.01 99.7 125. 0.03 42.1 125. 0.05 -6.02 125. 4.80 -2056.5 124. Аналитическое решение к =125,3 мПа Таким образом, использование векторной аппроксимации полей перемещений конечного элемента позволяет эффективно исследовать напряженно-деформированное состояние конструкций при смещении их как жесткого целого.

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

При выводе основных соотношений рассматриваются три состояния системы. Первое – начальное или исходное, в котором положение произвольной точки сплошной среды М0 определяется радиус вектором R o (рис. 16). Второе – деформированное состояние после j – шагов нагружения.

В результате чего точка М0 займет новое положение М, определяемое радиус вектором R. И третье состояние системы, близкое ко второму – деформированное после (j+1)-го шага нагружения. Точка М в третьем состоянии занимает положение М*, определяемое радиус-вектором R *.

Перемещение точки М0 из первого состояния во второе – в положение точки М - будет характеризоваться суммарным вектором перемещения за j шагов нагружения V, а из второго состояния в третье – шаговым вектором перемещения V.

Векторы a1, a2, a3 ;

a1, a2, a3 и a1, a2, a3 представляют собой 0 0 0 * * * ковариантные векторы локальных базисов произвольной точки сплошной среды в трех рассматриваемых состояниях соответственно.

M V R M x3 R V M* i3 R x i i x Рис. 16. Перемещение произвольной точки сплошной среды в результате деформации при шаговом нагружении конструкции Радиус-вектор R o точки в исходном состоянии, определяется векторным уравнением (21). Поочередным дифференцированием радиуса вектора R o выражаются ковариантные векторы локального базиса (22).

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

Выражения суммарного тензора деформаций за j шагов нагружения определяются по формуле ( ) 1 0 0 mn = amV,n + an V,m + V,mV,n, ( m, n = 1,2,3) (78) где V,i ( i = 1,2,3) производные суммарного вектора перемещений за j шагов нагружения в глобальной системе координат.

Радиус-векторы, определяющие положение точек М и М* после j и (j+1)-ом шагов нагружения, будут равны (рис. 16) R = R0 + V ;

R* = R + V. R* = R 0 + V + V. (79) Вектор перемещений точки на шаге нагружения V = a10 u 1 + a20 u 2 + a30 u 3, (80) где u приращения компонентов вектора перемещения на (j+1)-ом шаге нагружения ( = 1,2,3).

Дифференцированием (79) по глобальным координатам можно определить ковариантные векторы локального базиса после (j+1)-го шага нагружения a* = R,* = a + V,, (81) 0 где V, = a, u + a u, производные вектора перемещений на шаге нагружения.

Ковариантные компоненты тензора приращения деформаций на (j+1) ом шаге нагружения определяются на основе соотношения ij = ( aij aij ), 1* (82) где aij = ai a j, aij = ai* a *.

* (83) j Ковариантные векторы базисов после j шагов и (j+1) –го шага нагружения a = R, = a0 + V, и a* = a0 + V, + V,, =1,2,3. (84) Из (82) можно получить [( ] ) ( ) 1 0 mn = am + V,m V,n + an0 + V,n V,m + V,m V,n. (85) Приращения напряжений на шаге нагружения определяются матричным выражением { } = [ C ]{ }, (86) где { } = { 11, 22, 33, 12, 23, 31 } T вектор приращений напряжений;

{ } = { 11, 22, 33,2 12,2 23,2 31,} T вектор приращений деформаций;

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

Закон сохранения энергии статически нагружаемого тела в приращениях может быть записан в виде равенства [ ] { } + k { } { } dV = { v} [{ P} + k { P} ] ds, T T T (87) V S где { P}, { P} - вектор нагрузок и их приращений;

{ V } T = { u, v, w} - вектор приращений перемещений.

Значение коэффициента k принимается равным 1/2.

Вектор приращений перемещений { V } внутренней точки конечного элемента определяется через вектор приращений узловых перемещений { V } = [ A]{ VУГ }, (88) где [A] – диагональная матрица, элементами которой являются векторы строки аппроксимирующих функций;

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

Выделяя в приращениях деформаций линейную и нелинейную части вектор приращений { } можно представить суммой { }{ } { } = л + н, (89) {} н где элементы матрицы зависят от компонент вектора перемещений { V } и не зависят от компонент вектора перемещений {V } за предыдущие {} шаги нагружения, элементы матрицы л зависят от всех указанных компонент.

{} Матрицу л можно представить в виде { л } = [ B]{ V уГ }. (90) Используя (89), (90), из (87) можно получить { VуГ } Т k [ B] Т [ C ][ B] dV { V уГ } + { Н } T { } dV = V V (91) = { V} [ A] [{ P} + k { P} ] ds { } { } dV ГТ ЛТ T у S V Здесь { V } - вектор приращений узловых смещений на (j+1)-м шаге Г у нагружения.

В качестве примера рассчитывалась тонкостенная цилиндрическая панель, загруженная в середине пролета линейной нагрузкой интенсивностью q (рис. 17). Были использованы следующие исходные данные: внутренний радиус r=4.953 м;

ширина панели L=1.0 м;

модуль упругости материала Е=1.05107 кПа;

коэффициент Пуассона =0.27;

толщина стенки цилиндра h=0.094 м.

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

Рис. 17. Цилиндрическая панель, загруженная линейной нагрузкой Рис. 18. График зависимости между прогибом оболочки и величиной нагрузки Таблица № п.п., Величина нагрузки на Вертикальный прогиб, м Кол-во шаге нагружения За один шаг Суммарный w w шагов Р, кН/м 1 2 3 1 20 0.0050 0. 2 40 0.0055 0. 3 60 0.0056 0. 1 2 3 4 80 0.0061 0. 5 100 0.0064 0. 6 120 0.0069 0. 7 140 0.0073 0. 8 160 0.0079 0. 9 180 0.0085 0. 10 200 0.0092 0. 11 220 0.0099 0. 12 240 0.0109 0. 13 260 0.0119 0. 14 280 0.0132 0. 15 300 0.0146 0. 16 320 0.0164 0. 17 340 0.0187 0. 18 360 0.0216 0. 19 380 0.0256 0. 20 400 0.0312 0. 21 420 0.0399 0. 22 440 0.0553 0. 23 460 0.0875 0. 24 480 0.1827 0. 25 500 0.7451 1. В примере расчета цилиндрической панели с защемленными концами загруженной линейной нагрузкой, результаты расчета сравнивались с экспериментальными данными других авторов. Анализ полученных графиков показывает удовлетворительное соответствие значений перемещений, определенных с использованием описанного алгоритма и полученных экспериментально. Здесь же приводится алгоритм дискретного продолжения решения по параметру в окрестности предельной точки деформирования оболочки на основе метода конечных элементов в геометрически нелинейной постановке. Приводится пример расчета.

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

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

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

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

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

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

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

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

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

1. Клочков, Ю. В. Особенности формирования матрицы жесткости треугольного конечного элемента размером 54х54 / Ю. В. Клочков, А. П.

Николаев, А. П. Киселёв // Изв. вузов, сер. Строительство. – 1998. - № 2. – С.

32-37.

2. Клочков, Ю. В. Применение четырехугольного конечного элемента с матрицей жесткости 72х72 для расчета оболочечных конструкций / Ю. В.

Клочков, А. П. Николаев, А. П. Киселёв // Изв. вузов, сер. Строительство. – 1998. - № 4-5. – С. 36-41.

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

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

5. Николаев, А. П. О функциях формы в алгоритмах формирования матрицы жесткости в треугольном конечном элементе / А. П. Николаев, Ю.

В. Клочков, А. П. Киселёв // Изв. вузов, сер. Строительство. – 1999. - № 10. – С. 23-27.

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

Николаев, А. П. Киселёв, С. Н. Дейнего // Труды междунар. научн. конф.

«Актуальные проблемы механики оболочек». - Казань, 2000. – С. 218-219.

7. Киселёв, А. П. Восьмиугольный конечный элемент для расчета толстостенных оболочек вращения / А. П. Киселёв, А. П. Николаев, В. Н.

Юшкин // Сб. трудов междунар. научно-техн. конф. «Актуальные проблемы механики оболочек». - Казань, 2000. – С. 27-30.

8. Киселёв, А. П. Функции формы объемных конечных элементов / А.

П. Киселёв, А. П. Николаев // Сб. междунар. научно-техн. конф.

«Информационные технологии в образовании, технике и медицине». – Волгоград. - ч. 2 - 2000. – С. 110-113.

9. Николаев, А. П. Использование теории упругости трехмерного тела в расчетах оболочек / А. П. Николаев, А. П. Киселёв // Сб. трудов междунар.

научной конф. «Архитектура и прочностной расчет тонкостенных строительных и машиностроительных конструкций сложной формы». Москва, 2001. – С. 58-59.

10. Николаев, А. П. Расчет оболочек с использованием трехмерных конечных элементов в виде треугольной призмы и восьмиугольника / А. П.

Николаев, А. П. Киселёв // Сб. трудов междунар. научной конф.

«Архитектура и прочностной расчет тонкостенных строительных и машиностроительных конструкций сложной формы». - 2001. – С. 319-323.

11. Киселёв, А. П. Восьмиугольный конечный элемент для расчета на прочность участков соединения труб водохозяйственного назначения / А. П.

Киселёв, В. Н. Юшкин // Материалы V регион. конф. молодых исследователей Волгогр. обл. – Волгоград, 2001.

12. Николаев, А. П. Расчет на прочность грунтовых плотин на основе метода конечных элементов / А. П. Николаев, А. П. Киселёв, С. С. Марченко // Матер. междунар. научно-практ. конф. «Проблемы научного обеспечения экономической эффективности орошаемого земледелия в рыночных условиях». - г. Волгоград, 2001. – С. 28-29.

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

1, техн. науки. - Москва, 2001. – С. 54-55.

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

15. Киселёв, А. П. Объемный конечный элемент в форме треугольной призмы / А. П. Киселёв, А. П. Николаев // Межвуз. сб. научных трудов «Концептуальное проектирование в образовании, технике и технологии». Волгоград, 2001. – С. 137-141.

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

107-111.

17. Николаев, А. П. Решение задачи нелинейного деформирования оболочки на основе МКЭ при наличии особой точки / А. П. Николаев, А. П.

Киселёв // Деп. в ВИНИТИ 26.12.2002, - № 2262- В2002.

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

19. Николаев, А. П. Напряженно-деформированное состояние в зоне сочленения двух цилиндрических оболочек / А. П. Николаев, А. П. Киселёв, В. Н. Юшкин // Научный вестник, сер. Инж. науки. - №3.- ВГСХА, 2002. С.

167-170.

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

21. Киселёв, А. П. Сравнительный анализ результатов прочностного расчета тонкостенных оболочек с использованием двумерных и трехмерных конечных элементов / А. П. Киселёв // Матер. междунар. научно-практ. конф.

«Проблемы АПК», сер. Инж. науки. – ВГСХА, 2003. – С. 175-176.

22. Киселёв, А. П. Инженерное обследование и выполнение прочностных расчетов эксплуатируемых гидросооружений / А. П. Киселёв, А. П. Николаев // Научный вестник, сер. Инж. науки. - № 4. - ВГСХА. – Волгоград, 2003. – С. 86-88.

23. Киселёв, А. П. Формирование матрицы жесткости конечного элемента в геометрически нелинейной постановке / А. П. Киселёв // Матер.

III междунар. научно-техн. конф. «Эффективные строительные конструкции:

теория и практика». - Пенза, 2004. – С. 367-370.

24. Киселёв, А. П. Деформационная теория пластичности при использовании трехмерных конечных элементов / А. П. Киселёв // Матер.

междунар. научн.-практ. конф. «Современные оросительные мелиорации – состояние и перспективы». - Волгоград, 2004. – С. 129-132.

25. Николаев, А. П. Расчет многослойных оболочек на основе метода конечных элементов / А. П. Николаев, А. П. Киселёв, Р. З. Киселёва // Научные сообщ. клуба докторов наук. – №14. - Волгоград, 2005. – С. 21-23.

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

27. Николаев, А. П. Расчет оболочек в трехмерной постановке с учетом геометрической нелинейности на основе МКЭ / А. П. Николаев, А. П.

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

28. Киселёв, А. П. Объемный конечный элемент в виде треугольной призмы с первыми производными узловых перемещений / А. П. Киселёв, А.

П. Николаев // Изв. Вузов, сер. «Строительство». -2006.- № 1. – С. 13-18.

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

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

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

32. Киселёв, А. П. Расчет тонких оболочек на прочность в трехмерной постановке без упрощающих гипотез / А. П. Киселёв // Изв. Вузов, сер.

«Строительство». -2008.- № 1.

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

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

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

Киселёв Анатолий Петрович Развитие метода конечных элементов в исследованиях линейного и нелинейного деформирования оболочек как двумерных и трехмерных упругих тел Автореферат диссертации на соискание ученой степени доктора технических наук Подписано в печать 25.01.08. Формат 60х84 1/16.

Усл. печ. л. 2,5. Тираж 100. Заказ № 143.