Метод конечных элементов

Решение методом конечных элементов двухмерной магнитостатической задачи (линии и цвет означают направление и величину магнитной индукции)

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

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

Идея метода[править | править код]

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

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

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

Метод получил широкое применение при проектировании сооружений, а также при моделировании моделей движения, к примеру, грунта. За рубежом метод почти сразу начал повсеместно использоваться, а в России — только в 2000-х годах сменил вариационно-разностные, конечно-разностные методы и пр.[источник не указан 148 дней]

Из недостатков метода стоит отметить влияние размера сетки на конечные результаты.

Иллюстрация метода на одномерном примере[править | править код]

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

Базисные функции vk (голубые) и линейная комбинация из них, которая аппроксимирует искомую функцию (красная).

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

Пусть в одномерном пространстве Р1 необходимо решить следующее одномерное дифференциальное уравнение для нахождения функции на промежутке от 0 до 1. На границах области значение функции равно 0:

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

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

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

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

(1)

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

(2)

Оно получено с учётом того, что .

Разобьём область, в которой ищется решение

такое, что

на конечные промежутки, и получим новое пространство :

(3) такое, что

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

для (в данном примере )

Если теперь искомое приближённое решение представить виде , а функцию аппроксимировать как , то с помощью (3) можно получить следующую систему уравнений относительно искомых :

,

где .

Преимущества и недостатки[править | править код]

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

Долгое время широкому распространению МКЭ мешало отсутствие алгоритмов автоматического разбиения области на «почти равносторонние» треугольники (погрешность, в зависимости от вариации метода, обратно пропорциональна синусу или самого острого, или самого тупого угла в разбиении). Впрочем, эту задачу удалось успешно решить (алгоритмы основаны на триангуляции Делоне), что дало возможность создавать полностью автоматические конечноэлементные САПР.

История развития метода[править | править код]

Метод конечных элементов возник из необходимости новых путей решения задач строительной механики и теории упругости в 1930-х годах. Одними из основоположников идей, лежащих в основе МКЭ, считаются Александр Хренников и Рихард Курант. Их работы опубликованы в 1940-х годах. Впервые эффективность МКЭ была продемонстрирована в 1944 году Иоаннисом Аргирисом, который реализовал метод с применением ЭВМ.

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

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

Дальнейшее развитие метода конечных элементов связано также с решением задач космических исследований в 1950-х годах.

В СССР распространение и практическая реализация МКЭ в 1960-х годах связана с именем Леонарда Оганесяна.

Существенный толчок в своём развитии МКЭ получил в 1963 году после того, как было доказано, что его можно рассматривать как один из вариантов распространённого в строительной механике метода Рэлея — Ритца, который путём минимизации потенциальной энергии сводит задачу к системе линейных уравнений равновесия. После того, как была установлена связь МКЭ с процедурой минимизации, он стал применяться к задачам, описываемым уравнениями Лапласа или Пуассона. Область применения МКЭ значительно расширилась, когда было установлено (в 1968 году), что уравнения, определяющие элементы в задачах, могут быть легко получены с помощью вариантов метода взвешенных невязок, таких как метод Галёркина или метод наименьших квадратов. Это сыграло важную роль в теоретическом обосновании МКЭ, так как позволило применять его при решении многих типов дифференциальных уравнений. Таким образом, метод конечных элементов превратился в общий метод численного решения дифференциальных уравнений или систем дифференциальных уравнений.

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

См. также[править | править код]

• Метод дискретного элемента
• Метод конечных разностей
• Метод конечных объёмов
• Метод подвижных клеточных автоматов
• Метод граничного элемента

Литература[править | править код]

• Галлагер Р. Метод конечных элементов. Основы: Пер. с англ. — М.: Мир, 1984
• Деклу Ж. Метод конечных элементов: Пер. с франц. — М.: Мир, 1976
• Зенкевич О. Метод конечных элементов в технике — М.: Мир, 1975.
• Зенкевич О., Морган К. Конечные элементы и аппроксимация: Пер. с англ. — М.: Мир, 1986
• Сегерлинд Л. Применение метода конечных элементов — М.: Мир, 1979. — 392 С.

Ссылки[править | править код]

• Боровков А.И. и др. Компьютерный инжиниринг. Аналитический обзор — учебное пособие. — СПб.: Изд-во Политехн. ун-та, 2012. — 93 с. — ISBN 978-5-7422-3766-2.

Примечания[править | править код]

Лекция третья. Введение в Метод конечных элементов (МКЭ)

У Ч Р Е Ж Д Е Н И Е ЦЕНТР НЕЗАВИСИМОЙ ЭКСПЕРТИЗЫ НА АВТОМОБИЛЬНОМ ТРАНСПОРТЕ «ЦНЭАТ»

443098 г. Самара, ул. Пугачевская 73А, (АТП-5) тел. (846) 958-87-45 тел/факс. (846) 958-84-09, : [email protected]

Лекция третья. Введение в МКЭ

«Не пора ли, друзья мои, нам замахнуться на Вильяма, понимаете ли, нашего Шекспира?»

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

Итак, пусть есть абстрактная задача найти некоторую зависимость yот xв интервале от aдо b. Можно поступить двояко — 1) искать аналитический вид зависимости в виде функции y=y(x), т.е. в виде некоторой формулы, как например, все делали в школе, когда брали интеграл, или 2) искать функцию в виде набора точек с некоторой нужной или заданной точностью.

Так, например, на рисунке выше показана некоторая искомая функция y=y(x), вместо которой найдена ломаная 1-2-3-4-5-6 без особого ущерба для точности. Понятно, что чем в большем числе точек искать значение функции, тем точнее будет результат ее представления ломаной. А методы, которые вместо аналитической зависимости находят искомую функцию (или много функций) в виде конечного числа точек (чисел), называются численными.

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

Для определенности на рисунке выше наша балка уже разбита (представлена) в виде двух (не будем усложнять) конечных элементов, первый из которых имеет узлы 1 и 2, а второй — узлы 2 и 3. Т.е. узел 2 является общим для соседних конечных элементов.

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

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

1) точка может перемещаться по вертикали на расстояние v,

2) поперечное сечение в точке может развернуться на некоторый угол Ө.

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

.

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

Известно (например, из сопромата), что функция угла разворота сечения есть первая производная от функции перемещения, т.е.

.

Положим, что начало системы координат расположено в узле iконечного элемента, а координатная ось xнаправлена вправо. Тогда выразим коэффициенты полинома через перемещения узлов iи j, т.е. vi, Өi, vjи Өj.

При x=0, т.е. в узле i

В узле j, т.е. при x=xj=l, где l=xj-xi- длина конечного элемента, и с учетом полученных выше значений и

Решая эту систему уравнений, окончательно получаем

Подставляя полученные выражения для в формулу для v, получаем

Выглядит громоздко, но результат получен — мы выразили перемещение в любой точке конечного элемента через перемещения в его узлах — основных неизвестных задачи, так как, разбивая конструкцию на конечные элементы, мы и сбирались вычислять искомую функцию перемещений лишь в некоторых точках, а между ними представлять ее приближенно (аппроксимируя). Теперь, зная длину элемента l, перемещения в узлах конечного элемента vi,Өi,vj,Өjи координату какой-либо точки xмы можем по этой формуле найти перемещение этой точки v.

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

.

Введем вектор формы элемента в виде

.

Тогда выражение для перемещения vможно компактно записать в виде

,

где верхний индекс Тобозначает операцию транспонирования матрицы.

Далее следует договориться, что имеется в виду под словом «деформация». В английском языке русскому слову «деформация» соответствуют два несколько разных понятия — deformationи strain. Первое из них обозначает изменение геометрии чего-либо, т.е. то же самое, что и русское слово «деформация», применяемое в разговорном языке. Это — обобщенное понятие для растяжения, сжатия, скручивания, изгиба, смятия и т.п. Английское же слово strainобозначает меру деформации, или степень деформации.

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

.

Величину и будем называть деформацией, или относительной деформацией. Недостатком относительной деформации является то, что она не аддитивна — при сложении относительных деформаций накапливается ошибка. Например, образец длиной 10см некий экспериментатор растянул на 1см, т.е. получил деформацию =1/10. Теперь наш образец имеет длину 11см. Другой экспериментатор, не зная о прошлой истории образца, растянул его еще на 1см, т.е. получил деформацию =1/11. Но тогда суммарная относительная деформация равна

,

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

Наша конструкция имеет некоторую толщину, пусть постоянную по длине конечного элемента. Малый фрагмент элемента длиной dxдо изгиба показан на рисунке выше слева, а после изгиба — справа. При изгибе некоторое материальное волокно (ось) АВ не изменит своей длины, т.е. АВ=А1В1, а радиус его кривизны составит . Относительная деформация произвольного волокна CDбудет равна

откуда окончательно получаем

.

Из математики известна зависимость для малых величин углов разворота сечения

.

Тогда окончательно имеем

,

или

,

где вектор есть

.

Согласно закону Гука напряжения (силы на единицу сечения) связаны с деформациями через модуль упругости (модуль Юнга) E

.

Известно (из сопромата), что изгибающий момент связан с напряжениями в сечении, с учетом осевого момента инерции сечения, выражением

,

откуда получаем вектор узловых сил

.

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

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

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

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

,

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

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

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

.

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

,

где — подматрицы размером 2х2, а подстановка вместо iи j номеров узлов конечного элемента показывает место каждой подматрицы в глобальной матрице жесткости конечно-элементного аналога.

Аналогично структуру вектора перемещений и вектора узловых сил конечного элемента можно представить в виде двух подвекторов, а подстановка вместо iи j номеров узлов конечного элемента показывает место каждого подвектора в глобальном векторе узловых перемещений и узловых сил конечно-элементного аналога

.

Сейчас, читатель, начнем решать реальную задачу в цифрах, и все станет окончательно понятно. Вновь обратимся к нашему второму рисунку. Пусть конечный элемент №1 имеет узлы 1 и 2 с координатами и мм. Пусть конечный элемент №2 имеет узлы 2 и 3 с координатами мм и мм, т.е. длина l=100мм одинакова для обоих конечных элементов. Для обоих конечных элемента примем модуль упругости кг/мм2, осевой момент инерции сечения мм4.

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

.

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

Получаем глобальную матрицу жесткости размером 6х6 в виде

.

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

.

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

,

тем самым имеем систему линейных уравнений

,

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

.

Теперь из решения системы уравнений всегда будет и .

Зададим перемещение узла №3 вниз на 10мм, или . Для этого сформируем некоторый вектор в виде

,

и получим измененный вектор по алгоритму

,

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

.

Обнулим пятую строку и пятый столбец матрицы , поставим в главной диагонали на пятой строке 1, поставим в векторе на пятой строке -10. Граничное условие внесено, а матрица и вектор теперь выглядят как

,

и

.

Теперь можно решать систему уравнений. Тогда решение есть

,

или , , мм, рад, мм, рад.

Ну и теперь, когда узловые перемещения найдены, можно подставить из в формулу для перемещения v и построить его график — форму деформированной балки.

С аналитическим решением можно не сравнивать — совпадение будет на все 100%, так как реально упругая линия балки описывается уравнением третьей степени.

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

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

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

Часто задают вопрос: — «По какой формуле Вы сделали расчет?». Ответ: , причем всегда и для всех ДТП.

В.Н.Никонов,

ведущий научный сотрудник Института механики Уфимского научного центра РАН,

кандидат технических наук

Статья в формате PDF

Смотри далее, лекция № 4

© 2007, ЦНЭАТ , г. Самара, ссылка на ЦНЭАТ и страницу обязательны (www.cneat.ru)

Литература:

1. ОФС.1.2.1.2.0003.15 Тонкослойная хроматография // Государственная фармакопея, XIII изд.
2. Мирский, «Медицина России X—XX веков» (Москва, РОССПЭН, 2005, 632 с.).
3. ОФС.1.2.1.1.0003.15 Спектрофотометрия в ультрафиолетовой и видимой областях // Государственная фармакопея, XIII изд.
4. https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%BA%D0%BE%D0%BD%D0%B5%D1%87%D0%BD%D1%8B%D1%85_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2.
5. https://cneat.ru/lex3.html.
6. З.С. Смирнова, Л.М. Борисова, М.П. Киселева и др. Противоопухолевая эффективность прототипа лекарственной формы соединения ЛХС-1208 для внутривенного введения // Российский биотерапевтический журнал. 2012. № 2. С. 49.
7. А.В. Ланцова, Е.В. Санарова, Н.А. Оборотова и др. Разработка технологии получения инъекционной лекарственной формы на основе отечественной субстанции производной индолокарбазола ЛХС-1208 // Российский биотерапевтический журнал. 2014. Т. 13. № 3. С. 25-32.
8. Bangun H., Aulia F., Arianto A., Nainggolan M. Preparation of mucoadhesive gastroretentive drug delivery system of alginate beads containing turmeric extract and anti-gastric ulcer activity. Asian Journal of Pharmaceutical and Clinical Research. 2019; 12(1):316–320. DOI: 10.22159/ajpcr.2019.v12i1.29715.