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

В следующем разделе основы МКЭ будут изложены в объеме, необходимом пользователям программ, реализующих данный метод. Более подробное изложение этой темы можно найти в многочисленных учебниках, например, в работах [14,18-22].

Общая формулировка уравнений метода конечных элементов

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

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

Рассмотрим дифференциальное уравнение (DE) для некоторой величины и(х) в зависимости от ее расположения

DE(u) = 0.

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

На первом шаге производится аппроксимация формы функции и для конечного элемента. Эта аппроксимация имеет следующий вид:

u(x) - ЕфдК*) • uN, N

где N — количество узлов в элементе; фл, — интерполяционные функции; uN — их значения (постоянные) в узловых точках.

Интерполяционные функции фд, называются функциями формы. Как будет пока­зано далее, в методе конечных элементов они имеют ключевое значение.

Значения uNинтерполяционных функций (4.53) в точках интерполяции являют­ся постоянными; от координат зависят только функции формы. Таким образом, пер­вая производная функциям, например, по координатехбудет выглядеть следующим образом:

(4.54)

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

При подстановке уравнений (4.53) и (4.54) в дифференциальное уравнение (уравнение (4.52)), в правой части уравнения, как правило, получается ненулевое значение е, называемое невязкой:

Г _ зф;У(х)

DE фдг(х),—^—,uN =c.

Чтобы минимизировать эту невязку, используется скалярное произведение не­вязки и так называемых весовых функций 1Т;, где (г = 1,2,3,..., N). Затем это произве­дение интегрируется по площади конечного элемента О., и результат приравнивается к нулю:

(4.56)

n

Таким образом, ошибка приближенных вычислений для дифференциального уравнения принудительным образом устанавливается на нуль путем осреднения по элементу (так называемая слабая формулировка). В этой обобщенной форме данная процедура называется методом взвешенных невязок {Methodof Weighted Residuals).

В зависимости от выбора весовых функций ИТ можно получать различные выраже­ния, одно из которых, известное как метод Галеркина, приобрело наибольшее значе­ние вследствие того, что данный метод обеспечивает наилучшую аппроксимацию [19].

При использовании метода Галеркина N весовых функций ИТ приравниваются к функциям формы ф(, и уравнение (4.56) приобретает вид

(4.57)

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

Если в дифференциальном уравнении DE (уравнение (4.52)) появляются произ­водные высших порядков, то с помощью теоремы Грина-Гаусса уравнение (4.57) сле­дует преобразовать в один объемный интеграл с производной порядка, уменьшенного на 1, и один поверхностный интеграл. Это важно, поскольку, как будет показано в следующем разделе [ 18,19], производные порядков выше единицы для функций фор­мы, как правило, не определены.

Функции формы и точность решения

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

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

Наибольшее значение имеет выбор функций формы для аппроксимации желае­мой величины. Они могут совпадать с функциями, аппроксимирующими форму рас­четной области, хотя это условие и не является обязательным. Невязка, определяе­мая уравнением (4.55), тем меньше, чем выше точность интерполяционного уравнения (4.53), описывающего желаемую функцию в элементе. Эта зависимость отчетливо проявляется на рис. 4.4 [23].

На рис. 4.4 показан профиль скоростей при течении псевдопластичной жидкости в трубе. Решение методом конечных элементов с выбором одного элемента по радиусу трубы и линейной интерполирующей функции привело к результату, помеченному цифрой 1. Хотя площадь, ограниченная отрезком прямой и границами расчетной об­ласти (объемный расход) и очень близка к действительному объемному расходу (так как осредненная невязка принудительным образом устанавливается на нуль, см. урав­нение (4.57)), в отдельных точках расчетный профиль скоростей может сильно

Рис. 4.4. Зависимость точности аппроксимации от порядка полинома и количества элементов

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

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

Наконец, необходимо обсудить вопрос описания первой производной (здесь и да­лее по тексту под производной подразумевается локальная производная) уравнения­ми (4.54) и (4.55). Этот вопрос важен, поскольку в дифференциальных уравнениях встречаются производные как минимум первого порядка; кроме того, параметры, по­лученные на основании поля скоростей (например, скорости сдвига), необходимы для расчета течения. Функции формы гарантируют непрерывность интерполируе­мых величин не только в пределах элемента, но и на его границах [ 18-22]. Как прави­ло, это условие не соблюдается для первых производных [23], как показано на рис. 4.5 на примере двух элементов с линейными функциями формы. Вследствие кривизны

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

Тем не менее, проводя последовательные вычисления, можно разработать приемы улучшения качества расчетных величин, получаемых на основе локальных производ­ных, и обеспечить сглаженные формы кривой во всей расчетной области (это так называемые процедуры сглаживания [20,23-25]).