Естественные координаты элемента

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

Рис. 10.1. Четырехсторонние двумерные элементы в криволинейных косоугольных координатах.

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

При разработке алгоритмов МКЭ часто бывает удобно визуализировать произвольный конечный элемент в декартовых координатах на плоскости (Рис. 10.2). Любой четырехсторонний конечный элемент в этом случае может быть представлен квадратом со стороной 2 в безразмерных единицах в пределах . Взаимооднозначное преобразование между локальной системой координат и глобальной описывается математическими уравнениями общего вида:

,

(10.1)

В случае произвольного двумерного конечного элемента первая часть формул (10.1) могут быть представлена следующим образом:

,

(10.2)

где n – число узлов на элементе.

Соотношения типа (10.1) или (10.2) носят название изопараметрического отображения. Вопрос корректности отображения может быть достаточно острым при расчетах методом конечных элементов. Необходимо правильное построение сетки конечных элементов, чтобы не было нарушения взаимооднозначности отображения (Рис. 10.2).

Рис. 10.2. Отображение элемента, заданного в естественных локальных координатах , на тот же элемент в глобальной системе координат .

В случае корректного отображения имеет место взаимооднозначное соответствие точек элемента в плоскостях и . В случае некорректного отображения взаимооднозначное соответствие нарушается. Так например, видно, что любой одной точке на линии 1-2 в плоскости соответствует две точки в плоскости .

Четырехсторонний билинейный элемент

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

(10.3)

Рис. 10.3. Четырех-узловой четырехсторонний билинейный элемент.

Функции формы элемента выражены следующими соотношениями через локальные естественные координаты:

(10.4)

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

Рис. 10.4. Четырехсторонний билинейный элемент: (а) геометрия элемента; (б) уравнения сторон, противоположных первому узлу; (в) вид функции формы .

Четырехсторонние элементы высшего порядка

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

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

(10.5)

Рис. 10.5. Четырехсторонние конечные элементы высшего порядка: (а) девяти-узловой биквадратичный элемент; (б) восьми-узловой «серендипов» элемент.

Приведем формулы для некоторых характерных функций формы биквадратичного конечного элемента (Рис. 10.6):

(10.6)

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

,

(10.7)

где функции формы могут быть представлены следующим образом:

,

(10.8)

где - локальные координаты i -го узла на элементе.

Рис. 10.6. Графики характерных функций формы биквадратичного конечного элемента.

Заметим, что соотношения (10.4) линейного элемента также могут быть представлены в компактной форме, удобной для программирования:

(10.9)

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

(10.10)

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

Свойство полноты

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

,

(10.11)

где - произвольные константы.

Проверим, удовлетворяют ли рассмотренные функции формы и соответствующие приближенные решения (9.4, лекция 9) данному условию. Вычислим перемещения в узлах согласно соотношениям (10.11):

(10.12)

Подставим в первое интерполяционное соотношение (9.4, лекция 9). Получим:

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

Контрольные вопросы

1. Дать понятие естественных координат четырехстороннего конечного элемента.

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

3. Дать понятие полноты аппроксимирующей функции.

Лекция 11

Вычисление матрицы градиентовизопараметрического элемента

Введение

В настоящей и следующей лекции будет подробно рассмотрен алгоритм формирования матриц жесткости конечных элементов (6.23, лекция 6):

(11.1)

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

В общем случае алгоритм формирования ключевой характеристики механического элемента – его матрицы жесткости – включает следующие шаги:

1. Построение функций формы элемента;

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

3. Численное интегрирование по площади (или объему в трехмерном случае) для окончательного вычисления матрицы жесткости.

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

Матрица градиентов

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

Итак, рассмотрим произвольный четырехсторонний изопараметрический конечный элемент, заданный n узлами. Матрица градиентов имеет вид (5.8, лекция 5):

(11.2)

Введем следующие обозначения для сокращения записи. Пусть матрица градиентов представлена в виде матричной строки, образованной блоками (3x2):

,

(11.3)

где произвольный блок, связанный с узлом i имеет вид:

,

(11.4)

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

Матрицы Якоби

Пусть функции формы представлены сложной функцией многих переменных:

(11.5)

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

(11.6)

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

,

(11.7)

где J – матрица Якоби ( называемая также якобианом преобразования), состоящая из частных производных и определяющая однозначное отображение элемента, заданного в глобальных координатах , на «единичный» элемент в локальной плоскости :

Рассмотрим вычисление матрицы Якоби:

(11.8)

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

(11.9)

После того как определены компоненты матрицы Якоби, система (11.7) может быть разрешена относительно неизвестных функций :

,

(11.10)

где представляет собой обратную матрицу Якоби.

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

(11.11)

где - определитель матрицы Якоби, также иногда называемый якобианом:

(11.12)

Окончательно, можно записать выражения 2 n искомых компонент матрицы градиентов в развернутом виде:

,

(11.13)

где - общее количество элементов в конечно-элементной сетке, называемое глобальным числом конечных элементов.

Обратим внимание на особенности формул (11.13):

1. В общем случае все входящие в формулы переменные являются функциями локальных естественных координат элемента , включая компоненты матрицы Якоби и якобиан, что следует из правил их вычисления (11.9);

2. Компоненты матрицы Якоби и якобиан зависят только от номера элемента (геометрии элемента), но не связаны с номером узла на элементе. В то же время, функции формы естественным образом связаны с номером узла на элементе, что приводит к однозначной зависимости компонент матрицы градиентов от локального номера узла на элементе.

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

Таким образом, определены выражения 2 n искомых компонент матрицы градиентов и, следовательно, могут быть сформированы матрицы градиентов любого элемента согласно формулам (11.3).

Замечание

Заметим, что обратная матрица Якоби определяет обратное отображение «единичного» элемента в локальной плоскости на элемент, определенный в глобальных координатах:

По определению якобиана преобразования, можно записать:

,

(11.14)

где – матрица Якоби обратного преобразования.

Покажем, что совпадает с . Для этого достаточно перемножить две матрицы:

,

(11.15)

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

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

Ограничения на геометрию элементов

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

и

,

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

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

Рис. 10.1. Четырехсторонний конечный элемент с требуемыми узловыми углами.

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

Рис. 10.2. Примеры четырехсторонних конечных элементов с некорректными значениями узловых углов.

Контрольные вопросы

1. Дать понятие матрицы Якоби и якобиана.

2. Вычислить матрицу градиентов изопараметрического четырехстороннего конечного элемента.

3. Перечислить ограничения на геометрию изопараметрического четырехстороннего конечного элемента.

Лекция 12

Формирование матрицы жесткости изопараметрического элемента

Введение

В настоящей лекции будет подробно рассмотрен алгоритм формирования матрицы жесткости конечных элементов (6.23, лекция 6):

,

(12.1)

где - матрица градиентов элемента е;

- матрица упругих модулей элемента е;

- толщина элемента е;

- площадь элемента е.

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

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

2. Численное интегрирование по области или объему (для трехмерных элементов).

Рассмотрим первую из задач.

Структура матрицы жесткости

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

Обозначим:

(12.2)

Тогда:

(12.3)

Раскроем выражение (12.2), перемножив матрицы градиентов и упругих модулей, пользуясь блочным представлением (11.3, лекция 11):

,

(12.4)

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

,

(12.5)

Соответственно, и матрица жесткости имеет аналогичную матричную структуру, состоящую из блоков:

(12.6)

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

Матрица упругих модулей

Пусть двумерное тело находится в плоско-напряженном или плоско-деформированном состоянии (см. лекцию 4). Тогда определяющие соотношения имеют вид:

Плоско-напряженное состояние

(12.7а)

Плоско-деформированное состояние

(12.7б)

Запишем оба выражения (12.7а и б) в едином виде для повышения универсальности алгоритма:

,

(12.8)

где определяются соотношениями (12.7а) и (12.7б).

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

Выражение компонент матрицы жесткости элемента

Раскроем формулу (12.5), перемножив матрицы:

(12.9)

Таким образом, видно, что элементарный блок матрицы жесткости двумерного изопараметрического элемента (12.6) есть матрица размером 2x2. Поскольку таких блоков в элементной матрице жесткости n x n, то окончательный размер матрицы жесткости есть 2 n x2 n.

Численное интегрирование

Полученные выражения (12.9) должны быть проинтегрированы согласно формуле (12.6). Очевидно, что подынтегральные функции могут быть представлены как функции локальных естественных координат элемента согласно формулам (9.3, лекция 9) и (11.13, лекция 11):

(12.10)

Интеграл (12.10) есть интеграл по плоской замкнутой области. Для его вычисления применим известную математическую процедуру. Рассмотрим произвольный элемент е одновременно в локальной и глобальной системах координат. Выделим бесконечно малый элемент площади конечного элемента (Рис. 12.1).

Рис. 12.1. Бесконечно малый элемент площади конечного элемента в локальной и глобальной системах координат.

Элемент площади образован сторонами ОА и ОВ параллелограмма ОАСВ (Рис. 12.1), которые в свою очередь совпадают с координатными линиями и (см. рис. 10.1, лекция 10). Пусть радиус-вектор вершины параллелограмма точки О имеет радиус вектор:

(12.11)

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

(12.12)

Заметим также, что имеет место очевидное геометрическое соотношение:

(12.13)

Разложим радиус-вектор точки А (12.12) в ряд Тейлора относительно точки О:

(12.14)

Сравнивая соотношения (12.13) и (12.14), найдем с точностью до величин второго и высшего порядка малости вектор, определяющий сторону параллелограмма ОА:

(12.15)

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

(12.16)

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

(12.17)

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

,

(12.18)

где - определитель матрицы Якоби – якобиан, рассмотренный в предыдущей лекции (11.12, лекция 11).

Таким образом, окончательно находим:

(12.19)

С помощью полученного выражения (12.19) можно перейти от интеграла по области к повторному интегралу по локальным координатам в выражении блочных компонент матрицы жесткости (12.10):

(12.20)

Применим известные формулы интегрирования Гаусса для вычисления полученных интегралов. Согласно квадратурам Гаусса, одномерный интеграл от произвольной функции на интервале от -1 до +1 вычисляется следующим образом:

,

(12.21)

где - координаты гауссовых точек (узлов квадратуры);

- весовые коэффициенты квадратуры;

р – порядок квадратуры.

На рис. 12.2 схематично показаны четыре квадратурные схемы Гаусса соответственно 1-го, 2-го, 3-го и 4-го порядков.

Рис. 12.2. Четыре квадратурные схемы Гаусса порядка р =1, 2, 3, 4. Гауссовы точки показаны черным цветов. Радиусы точек пропорциональны весовым коэффициентам квадратур.

Формула (12.21) легко распространяется на двумерный интеграл. Ее применение к вычислению интегралов (12.20) дает следующее окончательное алгоритмическое выражение для вычисления компонент матрицы жесткости:

(12.22)

Расположение точек Гаусса на элементе схематично показано на рис. 12.3 для четырех низших квадратур.

Рис. 12.3. Низшие квадратурные схемы Гаусса порядка р = 1, 2, 3, 4 для четырехстороннего элемента. Гауссовы точки показаны черным цветов. Радиусы точек пропорциональны весовым коэффициентам квадратур.

Контрольные вопросы

1. Структура матрицы жесткости изопараметрического четырехстороннего конечного элемента.

2. Вычисление элемента площади изопараметрического четырехстороннего конечного элемента.

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

Лекция 13

Формирование векторов узловых сил изопараметрического элемента

Введение

В настоящей лекции будет подробно рассмотрен алгоритм формирования элементного вектора узловых сил, статически эквивалентного заданным объемным силам (6.27, лекция 6), и элементного вектора узловых сил, статически эквивалентного заданным поверхностным силам (6.28, лекция 6):

	(13.1)
	(13.2)

где - заданный вектор объемных сил;

- заданный вектор поверхностных сил;

- толщина элемента е;

- площадь элемента е;

- граница элемента е;

- матрица интерполирующих функций, или функций формы элемента.

Структура элементных векторов сил

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

Раскроем подынтегральное произведение (формулы (5.6, лекции 5) и (4.5в, лекции 4)):

(13.3)

Подставляя полученный матричный вектор (13.3), в формулу (13.1) получим выражение элементного вектора узловых сил, статически эквивалентного заданным объемным силам:

,

(13.4)

где компоненты вектора имеют вид:

(13.5)

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

,

(13.6)

где компоненты вектора имеют вид:

(13.7)

Вычисление компонент элементного вектора объемных сил

Компоненты элементного вектора объемных сил (13.5) вычисляются путем перехода к повторному интегралу аналогично вычислению элементных матриц жесткости, подробно рассмотренному в лекции 12:

(13.8)

Применения к (13.8) квадратуры Гаусса, в результате получим окончательное выражение для вычисления компонент элементного вектора объемных сил:

(13.9)

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

,

(13.9а)

где - значение заданной силы в узле i конечного элемента;

- значение толщины элемента е в узле i.

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

Вычисление компонент элементного вектора поверхностных сил

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

,

(13.10)

где - стороны элемента е.

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

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

Как уже было отмечено, стороны элемента, по определению, совпадают с координатными линиями локальной естественной системы координат (см. лекцию 10): и . Следовательно, параметрическое представление кривой 2-3, совпадающей с первой стороной элемента и описываемой уравнением в локальной системе координат, имеет вид:

(13.11)

Принимая во внимание изопараметрические соотношения (10.2, лекция 10), уравнение первой стороны элемента виде может быть представлено следующим образом:

(13.12)

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

,

(13.13)

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

Дифференцируя уравнения (13.12) и подставляя найденные значения в (13.13), получим выражение дифференциал дуги :

(13.14)

Ведем обозначение функции:

,

(13.15)

имеющей смысл якобиана преобразования между координатами точек кривой, заданной в локальной системе координат уравнением , и в глобальной системе координат уравнением (13.12). Другая, но близкая по сути, интерпретация функции (13.15) – якобиан отображения кривой 2-3 на отрезок прямой .

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

(13.16)

Численное интегрирование

Применения к (13.16) квадратуры Гаусса, в результате получим окончательное выражение для вычисления компонент элементного вектора поверхностных сил, заданных на первой стороне конечного элемента:

(13.17)

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

,

(13.17а)

где - значение заданной силы в узле i конечного элемента;

- одномерные функции формы;

n^* - число узлов на стороне конечного элемента.

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

Сторона II, узлы 4-1, уравнение :

(13.18)

Сторона III, узлы 3-4, уравнение :

(13.19)

Сторона IV, узлы 1-2, уравнение :

(13.20)

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

Контрольные вопросы

1. Как вычисляется формирования элементный вектор узловых сил, статически эквивалентного заданным объемным силам.

2. Как вычисляется дифференциал дуги.

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

Лекция 14

Формирование и решение глобальной системы конечно-элементных уравнений

Введение

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

Матричное уравнение (6.34, лекция 6) представляет собой стандартную форму записи системы линейных алгебраических уравнений метода конечных элементов (СЛАУ МКЭ):

(14.1)

где K – глобальная матрица жесткости конечно-элементной модели;

U – глобальный вектор степеней свободы (узловых перемещений) модели;

F _V – глобальный вектор узловых сил, статически эквивалентный заданным объемным силам;

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

Структура глобальной матрицы жесткости

Матрица жесткости конечно-элементной модели формируется путем сложения по определенным правилам элементных матриц жесткости:

,

(14.2)

где - матрицы кинематических связей.

Эта процедура называется сборкой или ансамблированием.

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

Пусть сетка образована двумя четырех-узловыми линейными изопараметрическими элементами (Рис. 14.1). Глобальное число узлов равно шести, при этом каждый узел имеет свой уникальный глобальный номер от 1 до 6. Введем локальную естественную систему координат на каждом элементе. Пронумеруем узлы на каждом элементе в соответствии с ориентацией локальной системы координат от 1 до 4. Таким образом, получаем, что каждому узлу ставится в соответствие два номера: локальный и глобальный , где - глобальное число узлов конечно-элементной сетки.

Данное соответствие нумерации узлов помогает легко сформировать матрицы кинематических связей , которые, по определению (6.18, лекция 6), связывают глобальный вектор узловых перемещений U с каждым элементным вектором узловых перемещений:

(14.3)

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

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

В данном случае имеем:

,	(14.4)
,	(14.5а)
.	(14.5б)

Сравнивая выражения (14.5а) и (14.5б) с (14.4) и принимая во внимание, что должны выполняться уравнения (13.3), получим для каждого элемента выражения матриц кинематических связей:

,	(14.6а)
,	(14.6б)

где 1 - единичная матрица размерности 2x2;

0 – нулевая матрица размерности 2x2:

,

(14.7)

Вычислим первое слагаемое в сумме (14.2):

,

(14.8)

где - матричный блок размерности 2x2, определяемый формулами (12.9 и 12.22, лекция 12).

Перемножив матрицы в формуле (14.8), получим «отображение» элементной матрицы первого элемента на «плоскость» глобальной матрицы:

(14.9)

Сразу заметим, что положение блока в матрице определяется соответствующими глобальными номерами. Так, например, блок имеет локальные номера i =1 и j =2. Узлы с локальными номерами i =1 и j =2 на первом элементе имеют глобальные номера I =4 и J =5. В результате блок перемещается на 4-ю строку, 5-й столбец глобальной матрицы .

Аналогично рассуждая, вычислим второе слагаемое в сумме (14.2):

(14.10)

После перемножения получим:

(14.11)

Остается сложить матрицы (14.9) и (14.11), чтобы получить искомую глобальную матрицу жесткости рассматриваемой конструкции:

(14.12)

Структура глобального вектора узловых сил

Глобальный вектор узловых сил формируется аналогичным образом:

(14.13)

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

,	(14.14а)
,	(14.14б)

где - матричный блок размерности 2x1, состоящий из компонент вектора силы в узле i элемента е (13.5, 13.9, лекция 13).

Подставляя выражения (14.14а, б) в (14.13) с учетом вида матриц кинематических связей (14.6а, б), получим глобальный вектор узловых сил:

(14.15)

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

Решение глобальной системы конечно-элементных уравнений

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

Однако при решении и выборе метода расчета систем конечно-элементных уравнений необходимо принимать во внимание ряд свойств СЛАУ и, в частности, свойств матрицы жесткости K:

1. При расчете реальных механических конструкций СЛАУ МКЭ может достигать десятков и сотен тысяч уравнений. Современный рекорд формирования и решения СЛАУ МКЭ равняется 50 миллионам уравнений. Заметим здесь же, что если задача отлична от статической линейной, то необходимое число решений системы может быть сколько угодно большим, например, при расчете динамики конструкции. К счастью, сложности, связанные с большим количеством уравнений, частично компенсируются другими свойствами матрицы жесткости, что и делает МКЭ таким привлекательным с точки зрения инженерных расчетов.

2. Матрица жесткости симметрична:

(14.16)

3. Хотя мы специально не доказывали свойство (14.16), его легко проверить на примере (14.12), если учесть, что элементные матрицы жесткости симметричны, т.е.:

(14.17)

4. Матрица жесткости является положительно-определенной, что означает, что квадратичная форма:

(14.18)

5. при любых .

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

7. Матрица я

Дата добавления: 2014-10-15