Лекция 10. Метод конечных элементов. Плоская задача.

10.1. Содержание метода. Метод конечных элементов (МКЭ) представляет

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

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

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

Разделим воображаемыми линиями диск в условиях плоской задачи на некоторое количество элементов конечных размеров, например, треугольной формы и примем за узловые точки их вершины. Очевидно, что если диск находится в равновесии то и его элемент, определенный узлами i, j, k, под воздействием напряжений (усилий) от смежных элементов, также уравновешен. Приложим затем к этому элементу вместо фактических усилий, действующих вдоль его граней, статически эквивалентные узловые силы, т.е. силы, вызывающие внутри элемента действительное напряженно- деформированное состояние.

Поставив в соответствие каждому узловому усилию узловое перемещение, представим сплошной диск набором конечных элементов (КЭ), взаимодействующих между собой в конечном числе узловых точек.

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

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

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

Алгоритм решения задач по МКЭ содержит следующие этапы:

1. Дискретизация - разбиение заданной области на КЭ; нумерация узлов и КЭ.

2. Аппроксимация перемещений узлов в КЭ.

3. Построение матриц жесткости (МЖ) конечных элементов.

4. Построение глобальной матрицы жесткости общей системы уравнений.

5. Сведение нагрузок и воздействий, приложенных к КЭ, к узловым силам, учет условий закрепления.

6. Решение общей системы уравнений.

7. Определение напряжений и (при необходимости) деформаций в КЭ.

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

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

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

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

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

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

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

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

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

Ширина полосы B вычисляется по формуле

B=(R+1)Q (10.1)

где R – максимальная по элементам величина наибольшей разности между номерами узлов в отдельном элементе; Q - число неизвестных (число степеней свободы в каждом узле). Минимизация величины В связана с минимизацией R, что может быть осуществлено последовательной нумерацией узлов при движении в направлении наименьшего размера тела против часовой стрелки. Правильная нумерация узлов экономит машинную память более чем на 60%, хотя не влияет на вычислительные аспекты задачи.

Построение градиентной матрицы. Для плоского КЭ «e » в форме треугольника с вершинами i, j, k с узловыми координатами соответственно (, ), (, ), (, ), перемещения каждого узла в направлении осей x и y имеют две компоненты u и v.

 

Так, для узла i

= , (10.2)

а для трех узлов i, j, k шесть компонент образуют вектор перемещений внутри КЭ «e »:

= = (10.3)

Простейшими интерполяционными функциями являются линейные полиномы

u = +

v = + (10.4)

Значения шести постоянных ÷ можно найти их двух систем, каждая из которых состоит из трех уравнений, полученных в результате подстановки в (10.4) узловых координат и приравнивания перемещений соответствующим узловым перемещениям.

Система уравнений для перемещений u узлов i, j, k:

= + ;

= + ; (10.5)

= + ;

позволяет определить коэффициенты , , по правилу Крамера через определители:

; ; ,

где ∆ - площадь треугольника ijk, удвоенное значение которой равно определителю

2∆=det = . (10.6)

После подстановки , , в первое уравнение (10.4), получим

u = , (10.7)

где коэффициенты определяются через координаты узлов КЭ:

; ; = . (10.8)

Остальные коэффициенты получаются циклической перестановкой индексов i, j, k.

Аналогично можно представить перемещение v, пользуясь вторым уравнением (10.4):

v = . (10.9)

Вектор перемещений внутри КЭ «e »принимает вид

= = , (10.10)

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

 

 

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

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

= = . (10.11)

Подставляя (10.10) в (10.11) получим

 

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

Матрица напряжений. Для упругого изотропного материала физические соотношения т.е. зависимости между напряжениями и деформациями, в данном случае закон Гука, линейны и определяются уравнениями:

где ̶ матрица упругости, содержащая упругие постоянные материала. Для изотропного материала в случае обобщенного плоского напряженного состояния ОПНС:

где 𝐸 ‒модуль упругости; 𝜇 ‒ коэффициент Пуассона.

Используя (10.12) и (10.13), напряжения в КЭ 𝑒 можно выразить через перемещения его узлов

где матрица называется матрицей напряжений:

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

Матрица жесткости конечного элемента 𝑒 (𝑖𝑗𝑘) устанавливает соответствие между узловыми силами и перемещениями по направлению этих сил. Матрица осредненно характеризует жесткостные свойства среды в объеме -го конечного элемента, а вектор ‒ приходящееся на этот элемент внешнее воздействие. Известно, что решение плоской задачи эквивалентно отысканию функций и , удовлетворяющих на контуре краевым условиям и минимизирующих функционал

Π =𝑈‒𝐴;

где энергия деформации 𝑈 определяется выражением (в соответствии с формулой Клапейрона):

Энергия деформации 𝑈 для КЭ 𝑒, с учетом матричного представления и полученных ранее соотношений

Напомним, что транспонирование произведения двух матриц выполняется по формуле (лк.1 ЧМ ч.1):

Работа 𝐴, совершаемая приложенными силами, может быть разделена на три части:

‒ работа совершаемая внешними сосредоточенными силами в узлах;

‒работа , совершаемая распределенными поверхностными силами приложенными на гранях КЭ;

‒работа , выполняемая объемными силами 𝑋, 𝑌.

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

Поверхностные и объемные силы можно также привести к узловым. Тогда

𝐴 = = .

При этом выражение для полной потенциальной энергии для КЭ принимает вид:

Чтобы согласно вариационному принципу Лагранжа минимизировать величину , необходимо продифференцировать выражение (10.18) по и результат приравнять нулю:

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

 

 

Теперь матрица жесткости элемента 𝑖𝑗𝑘 определяется так:

Отметим единообразие выражений для матриц жесткости стержневых и плоских (10.20) КЭ.

Если толщина элемента 𝑒 постоянна, то:

= h∆,(10.21)

где: Δ ‒ площадь КЭ; ‒ его объем.

Матрица жесткости является симметричной и имеет блочную структуру:

=

где ‒ квадратная симметричная подматрица порядка 𝑚 (𝑚 ‒ число компонент перемещений в узле), элементами которой являются реактивные усилия в связях, наложенных на узел 𝑟, вызванные единичными перемещениями и узла 𝑠.

Реактивные усилия в связях, наложенных на узлы 𝑖, 𝑗, 𝑘, вызванные перемещениями всех узлов, равны:

так, например, - реактивные усилия в связях, наложенных на узел 𝑖, вызванные перемещениями всех узлов

При ОПНС подматрицы имеют второй порядок и определяются следующим образом:

где и ‒ блоки градиентной матрицы :

;

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

Узловые силы, связанные с объемными силами X и Y (проекции их интенсивностей на оси x и y), определяются выражениями:

dxdy,

= - , и т.д.

…………………………

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

= - ∆/3 (10.26)

Узловые силы, вызванные поверхностными нагрузками и ,равны:

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

 

10.3. Построение глобальной матрицы жесткости. Поскольку полная потенциальная энергия (функционал) для рассматриваемой области может быть представлен суммой всех КЭ

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

,

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

 

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

Рис. 1

Пусть – типовой узел расчетной схемы области .

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

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

Из условий равновесия узла следует:

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

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

 

где – число узлов конечных элементов ;

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

Из подстановки (10.32) в (10.31) следует

где – матрица жесткости всего диска.

– вектор узловых сил, эквивалентных внешней нагрузке.

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

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

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

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