Метод конечных элементов (МКЭ) является сеточным методом, предназначенным для решения задач микроуровня, для которого модель объекта задается системой дифференциальных уравнений в частных производных с заданными краевыми условиями.
МКЭ основывается на методе взвешенных невязок, суть которого заключается в следующем: подбирается функция, удовлетворяющая дифференциальным уравнениям и краевым условиям, но подбирается не произвольно, поскольку такой подбор вряд ли возможен уже в двухмерном пространстве, а с использованием специальных методов.
Пусть состояние некоторой среды описывается следующим дифференциальным оператором, с заданным граничным условием:
LV + P =0
V(Г) = Vг
Здесь L - дифференциальный оператор (например оператор Лапласа), V - фазовая переменная -неизвестная функция, которую следует найти, P - величина, независящая от V. V(Г)=Vг - граничное условие первого рода (Дирихле), то есть на границе задано значение фазовой переменной. Решение задачи с граничными условиями Неймана будет рассмотрено ниже.
Будем искать решение с помощью функции, имеющей следующий вид:
Здесь V* - приближенное решение, F - функция, удовлетворяющая граничным условиям, Nm - пробные функции, которые на границе области должны быть равны нулю, Am - неизвестные коэффициенты, которые необходимо отыскать из условия наилучшего удовлетворения дифференциальному оператору, M- количество пробных функций.
Если подставить V* в исходный дифференциальный оператор, то получим невязку, принимающую в различных точках области разное значение.
R=LV*+P
Необходимо сформулировать условие, позволяющее минимизировать эту невязку по всей области. Одним из вариантов такого условия может быть следующее уравнение:
Здесь Wn - некоторые весовые
функции, в зависимости от выбора
которых различают варианты метода
взвешенных невязок, S- область
пространства, в которой ищется
решение. При выборе в качестве
весовых функций дельта-фукций
будем иметь метод, который получил
название метод поточечной
коллокации, для кусочно-постоянных
функций - метод коллокации по
подобластям, но наиболее
распространенным является метод
Галеркина, в котором в качестве
весовых функций выбираются пробные
функции N. В этом случае, если
количество пробных функций равно
количеству весовых функций, после
раскрытия определенных интегралов
приходим к замкнутой системе
алгебраических уравнений
относительно коэффициентов A.
KA + Q = 0
где коэффициенты матрицы K и вектора Q вычисляются по формулам:
После нахождения коэффициентов 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' можно обновить результаты расчета при изменении параметров балки, во всех остальных случаях результат обновляется автоматически, при этом справа от рисунка показывается максимальный прогиб балки.