МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ

Хиллиг В.Б.

Источник: Хиллиг В.Б. Сопротивление материалов. М- Металлургия, 1972.

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

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

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

LV + P =0

V(Г) = Vг

     Здесь L - дифференциальный оператор (например оператор Лапласа), V - фазовая переменная -неизвестная функция, которую следует найти, P - величина, независящая от V. V(Г)=Vг - граничное условие первого рода (Дирихле), то есть на границе задано значение фазовой переменной. Решение задачи с граничными условиями Неймана будет рассмотрено ниже.

R=LV^*+P

KA + Q = 0

     После нахождения коэффициентов A и подстановки их в (1), получаем решение исходной задачи.

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

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

1) в узле аппроксимации функции имеют значение равное единице;

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

     В результате на каждом КЭ действует строго определенное число ненулевых глобальных базисных функций ( в данном примере две) и вместо вычисления интеграла по всей области можно вычислить интегралы по КЭ и сложить их. Процедура сложения получила название ансамблирование. Использование глобальных базисных функций приводит к тому, что процедура вычисления интегралов по КЭ становиться достаточно простой и, поскольку в узлах аппроксимации Nm=1, коэффициенты A приобретают физический смысл, они становятся равными значению фазовой переменной в узлах. В аппроксимации (1) теперь можно отказаться от использования функции F, поскольку удовлетворить граничные условия можно естественным образом, задавая значения V в узлах, расположенных на границе.

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

V = (1-X/L)Vi + X/L Vj = N _eV

     Здесь X- текущая координата, отсчитываемая от начала КЭ, L- его длина, Vi и Vj - значения фазовых переменных в узлах КЭ. Компоненты вектора Ne получили название функций формы конечного элемента. Функции формы можно получить и из других соображений. Зададимся полиномом, аппроксимирующим решение внутри конечного элемента, например:

V = A0 + A1* X.

при X=0 Vi = A0, при X=L Vj=A0 + A1*L,

     Находим коэффициенты A0 и A1 и подставляем их в аппроксимацию

V= Vi + (Vj-Vi)*X/L = (1-X/L)Vi + X/L Vj

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

     В уравнении: Y - прогиб балки, M - момент силы, E - модуль Юнга, I - момент инерции сечения.

     Этапы решения этой задачи будем рассматривать с общих позиций.

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

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

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

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

где L - длина КЭ.

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

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

     Этап 5. Ансамблирование или получение глобальных матрицы жесткости и вектора нагрузок. Вектор неизвестных составляют прогибы в узлах, то есть матрица коэффициентов будет иметь четвертый порядок. Рассматривая поочередно каждый КЭ, располагаем элементы локальной матрицы жесткости в глобальной в соотвествии с номерами улов подключения КЭ. Проделать эту процедуру можно с помощью приведенного ниже апплета. Задать нумерацию узлов можно последовательно нажимая круглые кнопки, при этом последний узел нумеруется автоматически. Получить аддитивный вклад от отдельного КЭ в глобальную матрицу жесткости можно с помощью кнопок, расположенных над КЭ. Повторять нумерацию и получение аддитивного вклада можно многократно, пользуясь кнопкой 'reset' и кнопками, расположенными над конечными элементами. После ансамблирования получаем математическую модель объекта.

     Этап 6. Учет граничных условий. Количество неизвестных в полной математической моделе равно шести при четырех уранениях, но учитывая, что Y(0) =0 и первая производная в нулевом узле также равна нулю, получаем замкнутую систему уравнений.

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

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

     В качестве исходных параметров можно задать: 'number FE' - число КЭ, выбираемых из меню, в том числе можно задать разбиение при помощи 'мыши' в режиме 'Manual'; 'diametr' - диаметр балки (в миллиметрах), по нему расчитывается момент инерции сечения; 'Mmax'-максимальный по модулю момент силы (в килоньютонах на метр), прилагаемый к балке, ему на графике соответствуют крайние по высоте светлые линии; 'length' - длина балки (в метрах); Yong's modulus - модуль Юнга (в гигапаскалях); Fmax- максимальное значение силы по модулю (в килоньютонах), ему на графике соотвествуют крайние по высоте светлые линии. Режимы работы с апплетом задаются кнопкамим,расположенными под рисунком. Кнопка 'Diskret' - поволяет разбить балку на КЭ, как автоматически, так и с помощью 'мыши', кнопка 'Set/Reset M' - позволяет задать момент силы, воздействующий на балку с помощью 'мыши', в пределах светлых линий, кнопка 'Set/Reset F' - позволяет задавать силы, воздействующие на балку, в пределах светлых линий. С помощью кнопок 'Clear' производиться очистка соответствующих значений. С помощью кнопки 'Show result' можно обновить результаты расчета при изменении параметров балки, во всех остальных случаях результат обновляется автоматически, при этом справа от рисунка показывается максимальный прогиб балки.