Лекция N 4

§ 2. 5 Метод конечных элементов

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

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

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

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

Сама задача вычисления узловых значений формулируется как задача минимизации функционала, характеризующего качество аппроксимации решения v функцией u. В конечном счете задача сводится к решению СЛАУ, неизвестными которой являются значения в искомой функции в узлах сетки КЭ.

Рассмотрим основные положения МКЭ несколько подробнее.

Решение задачи МКЭ начинают с разбиения исследуемой области на конечные элементы.

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

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

Если задача осесимметричная, то КЭ представляет собой тело вращения, образуемое поворотом на 360^о относительно оси симметрии треугольного или четырехугольного элемента:

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

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

При разбиении плоской области на треугольные КЭ вручную соблюдают следующие правила:

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

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

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

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

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

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

Коэффициенты a можно выразить через координаты узлов (x, y) и значения функции в узлах (T_i, T_j, T_k) следующим образом

Решая систему получим значения коэффициентов a₁ a₂ a₃ то интерполяционный полином для аппроксимации скалярной величины (в данном случае температуры T) будет иметь вид

где:

Формулы для N_j, N_k получают из выражений для N_i круговой перестановкой индексов. Функции N получили названия функций формы. Их отличительная особенность - они равны 1 в том узле, индекс которого носят и равны 0 во всех остальных узлах элемента. Например для узла i

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

В матричном виде зависимость имеет вид

где [N]=[N_i, N_j, N_k] - матрица-строка функций формы, {T^e}={T_i, T_j, T_k}^Т - вектор столбец узловых значений температур для произвольного элемента.

В приведенном примере интерполировалась скалярная величина. Векторную величину (например перемещение) представляют через компоненты - проекции на оси координат, которые рассматривают как неизвестные скалярные величины.

Например, в случае плоской деформации вектор перемещений {q} произвольной точки будет характеризоваться компонентами {u, v}^т, направленным соответственно по осям x, y. В матричном виде соотношения между узловыми перемещениями и перемещениями точек элемента будет иметь вид:

Здесь функции формы N те же, что и в предыдущем случае.

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

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

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

где П - потенциальная энергия, L - энергия деформации, W - работа внешних сил. При конечно-элементной аппроксимации составляющие полной энергии получаются суммированием по элементам:

Для вычисления интегралов используем следующие известные соотношения:

- соотношения Коши, связывающие перемещения и деформации, где [¶] - матрица дифференциального оператора (например для плоской задачи e_x=¶ u/ ¶ x, e_y=¶v / ¶y, g_xy=¶ u/ ¶y+¶ v/ ¶x).

Поскольку окончательно можно записать

- где [B] - матрица градиентов -частных производных функций формы по пространственным координатам

- Закон Гука, где [D] - матрица упругости.

С учетом полученных выражений функционал полной энергии принимает вид, принимая во внимание, что ([A][B])^T= [B]^T[A]^T

Минимизацию функционала выполним продифференцировав его по перемещениям и приравняв дифференциал нулю.

Здесь учтено правило дифференцирования матричных соотношений:

Это выражение преобразовывают к виду

Здесь

[K] - глобальная матрица жесткости системы, представляющая собой квадратную матрицу n ´ n, где n - количество узловых неизвестных. Для плоской задачи n=2m - где m - количество узлов. Глобальная матрица жесткости получают суммированием матриц жесткости каждого элемента

{w} - вектор столбец узловых перемещений. Для плоской задачи

{F} - вектор столбец узловых нагрузок, представляет собой сумму векторов внешних сил, приложенных к узлам сетки и сумму поверхностных и объемных сил, приведенных к узлам сетки.

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

Для плоской задачи (плоское напряженное состояние)

здесь коэффициенты b, c - определены ранее, h - толщина элемента, D - площадь элемента, E - модуль Юнга (модуль упругости), m - коэффициент Пуассона.

После вычисления матрицы глобальной матрицы жесткости (алгоритм будет рассмотрен ниже) уравнение [K][w]=[F] превращается в СЛАУ, которая решается либо прямым (обычно методом Гаусса), либо итерационным (метод простой итерации, Зейделя, Ньютона, и т.д.) методами. После получения вектора узловых неизвестных {w}, используя выражения {e}=[B]{w} и {s}=[D]{e} определяют напряженно-деформированное состояние области.

Этапы решения задач МКЭ:

1. Выбор расчетной схемы (объемная, осесимметричная, плоская, одномерная, смешанная);

2. Выбор типа конечных элементов;

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

4. Описание элементов, включающее нумерацию узлов и самих элементов, определение координат узлов, механических характеристик материала;

5. Составление системы уравнений путем минимизации функционала, соответствующего данной задаче;

6. Решение полученной системы уравнений относительно неизвестных в узлах;

7. Вычисление предусмотренных постановкой задачи выходных параметров (напр. деформаций, напряжений) в теле по значениям неизвестных в узлах;

В настоящее время разработано значительное количество комплексов МКЭ, решающих самые различные задачи. Как правило, все программные продукты имеют пре, - и постпроцессоры для автоматизации ввода данных и разбиения области на КЭ, а также графического представления результатов расчета. Пользователь выполняет только 1 и 2 пункт, - т.е. формирует расчетную модель, все остальное выполняет ЭВМ (см лаб. работу N3, - программы ANSYS и LS-DYNA3D.).

Типы анализа, выполняемые с помощью программы ANSYS:

1. Статический и динамический анализ конструкций с учетом геометрической и физической нелинейности, ползучести и пластичности. Анализ усталостных разрушений.

2. Задачи линейной и нелинейной устойчивости конструкций.

3. Контактные задачи (например, 2-D и 3-D - штамповка). (ANSYS или/и LS-DYNA3D)

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

5. Задачи гидравлики и гидродинамики, ламинарное и турбулентное течение сжимаемой и несжимаемой жидкости (с учетом вязкости). Задачи обтекания тел произвольной формы, дозвуковой и сверхзвуковой режимы. (FLOTRAN)

6. Задачи акустики. (ANSYS and/or COMMET/ACOUSTIC)

7. Смешанные типы анализа (термо-механический, гидротепловой, магнитопрочностной и т.д.)

8. Задачи оптимизации.

9. Задачи электромагнитных полей и электростатики, с учетом анизотропии материалов.

Среди отечественных разработок можно выделить FORM2D предназначенный для анализа процессов обработки давлением плоских и осесимметричных задач с учетом температурного поля обрабатываемого материала. Этот комплекс ориентирован на персональные ЭВМ.

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

[1] Норма в данном случае вычисляется как корень квадратный из суммы квадратов

[2] Функционалом Ф, зависящем от функции f называется такая переменная величина, которая принимает конкретное числовой значение, при подстановке в нее каждой функции f из некоторого класса функций.