Портал | Содержание | О нас | Авторам | Новости | Первая десятка | Дискуссионный клуб | Чат Научный форум --> Первая десятка "Русского переплета" Темы дня:

Президенту Путину о создании Института Истории Русского Народа. |Нас посетило 40 млн. человек | Чем занимались русские 4000 лет назад?

| Кому давать гранты или сколько в России молодых ученых?

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

Л. А. РОЗИН

Санкт-Петербургский государственный технический университет

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

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

Допустим, требуется построить такое семейство функций u(x) при a # x # b. Интервал ab разбивается на конечное число частей (элементов), соединяющихся между собой и с концами интервала в узловых точках (узлах) xi (рис. 1). В пределах каждого элемента задается функция, например в виде линейного полинома. Она определяется своими значениями u(xi) в узлах на концах элемента. Если отыскиваемая функция является непрерывной, то значения ее в каждом узле для соседних элементов совпадают. В результате имеем семейство кусочно-линейных непрерывных функций, которые изображаются в виде ломаных и определяются конечным числом параметров - своими узловыми значениями. На рис. 1 показана одна из функций такого семейства. Здесь 5 элементов, 6 узлов и 6 узловых параметров u(xi) = ui . В случае нескольких переменных схема метода конечных элементов в принципе не меняется. Таким образом, метод конечных элементов заменяет задачу отыскания функции на задачу отыскания конечного числа ее приближенных значений в отдельных точках-узлах. При этом если исходная задача относительно функции состоит из функционального уравнения, например дифференциального уравнения с соответствующими граничными условиями, то задача метода конечных элементов относительно ее значений в узлах представляет собой систему алгебраических уравнений.

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

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

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

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

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

Первая работа, где рассматривалась схема типа метода конечных элементов, принадлежит известному математику Р. Куранту [2]. Построение метода с использованием физических соображений и его название "метод конечных элементов" содержатся в статье, написанной инженерами [3]. Такое сочетание специальностей авторов характерно для работ по методу конечных элементов. В последующем было опубликовано много статей и книг, посвященных этому методу и его различным модификациям. Некоторое представление об этом можно получить из списка литературы, например в [4]. Метод конечных элементов реализован в больших универсальных компьютерных пакетах программ, которые имеют широкое применение.

ПРИМЕР ДИСКРЕТНОЙ ЗАДАЧИ

Обратимся к дискретной задаче, состояние которой точно определяется конечным числом параметров. Рассмотрим упругий стержень в виде прямого кругового цилиндра, длина которого значительно больше его диаметра. Это позволяет отождествить стержень с его осью. Пусть три таких стержня расположены на оси x и соединены между собой в точках 1 и 2 (рис. 2, а). Точки a и b закреплены, что условно изображено на рисунке. К осям стержней вдоль x приложим внешнюю нагрузку. Очевидно, точки на осях стержней перемещаются вдоль x. Силы и перемещения считаются положительными, если они направлены в положительном направлении x. Задача состоит в определении перемещений точек, принадлежащих осям стержней и продольных внутренних сил в поперечных сечениях стержней.

Согласно методу конечных элементов, представим стержневую систему в виде элементов, соединенных в узлах. В качестве элементов примем отдельные стержни, а узлов - точки 1 и 2. На рис. 2, а в скобках указаны номера элементов. Обратимся к типовому для данной системы элементу r. На элемент r с узлами i, j (рис. 2, б ) может действовать распределенная нагрузка интенсивности q(r)(x) и перемещения его узлов Примем x = 0 в узле i и обозначим длину элемента l (r). Получим задачу для функции u(r)(x) перемещений точек оси r. Бесконечно малая часть элемента dx находится в равновесии под действием нагрузки q(r)(x)dx и продольных внутренних сил N (r)(x), действующих так, как показано на рис. 2, в. Из условия равновесия dx имеем

Согласно закону Гука, для упругого стержня

где c(r) > 0 носит название продольной жесткости стержня и определяется из опыта. Пусть c(r) = const для элемента r. Подставляя (2) в (1) получим задачу относительно u(r)(x) в виде дифференциального уравнения и граничных условий

В нашем примере для случая дискретной задачи положим q(r) = 0 и будем считать, что на стержневую систему действуют только сосредоточенные силы P1 и P2 соответственно в узлах 1 и 2. Тогда решение (3) примет вид

На основании (4) можно заключить, что состояние типового элемента r, то есть u(r)(x), N(r), точно определяется двумя параметрами - перемещениями его узлов что делает задачу дискретной.

Рассмотрим внутренние силы действующие в узлах i, j на элемент r. Поскольку имеет место линейная задача, то они линейно зависят от

Здесь (l = i, j; t = i, j) есть внутренняя сила действующая на элемент r в узле l и возникающая от единичного перемещения узла t. При этом перемещение другого узла равно нулю. На рис. 3, а показана сила для сжатого элемента при

Соотношения (5) можно представить в матричной форме. Введем столбцы f (r), u(r) и матрицу K(r)

Тогда (5) можно записать в виде

f (r) = K(r)u(r).

Для упругой пружины коэффициент пропорциональности между силой и перемещением называется коэффициентом жесткости пружины. Аналогично K(r) носит название матрицы жесткости элемента r.

От типового элемента перейдем к отдельным элементам данной системы. Для элемента 2 с двумя узлами 1, 2 справедливы все зависимости (4)-(7), где следует положить r = 2, i = 1, j = 2. Поскольку точки a, b неподвижны, то состояние элемента 1 определяется перемещением узла 1, а элемента 3 - перемещением узла 2. На основании (4) будем иметь для элементов 1 и 3 соотношения

В (8) координата x для элемента 1 равна нулю в точке a, а для элемента 3 - в узле 2. Зависимости (5) для элементов 2 и 3 примут соответственно вид

Сравнивая (5) или (7) для элемента 2 с (9) для элементов 1 и 3, можно заключить, что вместо матрицы жесткости для двуузлового элемента 2 фигурируют коэффициенты жесткости для одноузловых элементов 1 и 3. Теперь все известно о каждом отдельном элементе системы. Следующим шагом является соединение элементов в узлах на основе условий

где u1 , u2 - перемещения узлов системы 1, 2. Отсюда следует, что состояние соединенных элементов или системы в целом определяется двумя узловыми перемещениями u1 , u2 и рассматриваемая задача является дискретной.

Для всей системы можно записать соотношения типа (5) относительно суммарных для смежных элементов внутренних сил в узлах 1, 2. Обозначим их f1 , f2 . Очевидно, как и в случае типового элемента, они должны линейно зависеть от u1 , u2 :

f1 = f11u1 + f12u2 ,

f2 = f21u1 + f22u2 .

Введем столбцы f, u и матрицу жесткости всей системы K по формулам

Тогда матричное соотношение типа (7) для всей системы будет

f = Ku.

Здесь flt (l = 1, 2; t = 1, 2) есть суммарная для смежных элементов в узле l внутренняя сила fl , возникающая от единичного перемещения узла t. При этом перемещение другого узла равно нулю. Эти суммарные силы определяются через узловые силы в смежных элементах

На рис. 3, б показаны силы при u1 = 1, u2 = 0. При этом элемент 1 растянут, а элемент 2 сжат. В (13) поскольку узел 2 не принадлежит элементу 1, а узел 1 не принадлежит элементу 3. Из (13) следует, что матрица жесткости системы строится на основе коэффициентов жесткости для отдельных элементов. Алгоритмически выполнить это можно по-разному [5]. Например, можно для всех элементов строить матрицы жесткости одинаковой размерности равной размерности матрицы K, основываясь на столбце u перемещений всех узлов системы. Это возможно, поскольку если по крайней мере один из узлов l или t не принадлежит элементу r. В данном примере будем иметь

и, согласно (13), K = K(1) + K(2) + K(3). Из условия равновесия элемента r следует при где N (r) > 0 при растяжении и N (r) < 0 при сжатии. В результате на основании (4) и (8) будем иметь

а матрица K примет вид

Из условия равновесия узлов 1, 2 следует f1 = P1 , f2 = P2 или для столбцов f = P, где P - столбец из P1 , P2 . Подставляя сюда вместо f его выражение согласно (12), окончательно получим систему алгебраических уравнений относительно u1 , u2

После определения u1 , u2 в результате решения (15) находятся u(r)(x), N (r) во всех элементах системы при помощи (4) и (8).

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

СИСТЕМА ЛИНЕЙНЫХ

АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ

Метод конечных элементов сводит решение линейной задачи к решению системы линейных алгебраических уравнений

f11u1 + f12u2 + _ + f1nun = P1 ,

.............................................

fn1u1 + fn2u2 + _ + fnnun = Pn .

Здесь ui (i = 1, 2, _, n) - неизвестные, Pi (i = 1, 2, _, n) - заданные свободные члены, fi j (i, j = 1, 2, _, n) - коэффициенты при неизвестных. По аналогии с (11) коэффициенты fi j образуют квадратную матрицу, состоящую из n строк и n столбцов:

Если обозначить столбец неизвестных u, а столбец свободных членов P, то (16) принимает матричную форму (15). Система алгебраических уравнений должна быть невырожденной, то есть иметь единственное решение. Казалось бы, дальнейшее ясно. Можно воспользоваться для решения (16), например, методом исключения Гаусса. Однако при применении приближенных методов обычно приходится иметь дело с системами большого порядка n, и матрица, вообще говоря, может иметь такую структуру, которая затрудняет получение решения. При этом на точности результата в той или иной степени сказываются неизбежные в процессе вычислений ошибки округления. Одним из важных достоинств метода конечных элементов является то, что он обычно приводит к таким системам алгебраических уравнений, матрицы K которых позволяют эффективно строить решение.

Выясним, какой желательно иметь матрицу K в (16). Пределом мечты была бы система (16) с диагональной матрицей K, когда все fi j = 0 при i ? j. В этом случае (16) распадается на отдельные уравнения fi iui = Pi . Такое может быть, только если в физической системе, рассчитываемой методом конечных элементов, узлы между собой не связаны, то есть по существу системы не существует. Однако теперь уже ясно, к чему надо стремиться: следует так выполнять процесс построения алгебраической системы уравнений, чтобы матрица по возможности содержала больше нулевых коэффициентов и была близка к диагональной, другими словами, желательно, чтобы в каждое уравнение входило относительно небольшое число неизвестных в соседних узлах.

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

где знак * заменяет коэффициенты, отличные от нуля, а главная диагональ и параллельные ей линии указаны пунктиром. Ленточную матрицу характеризует ширина ленты t = t1 + t2 + 1, равная наибольшему числу коэффициентов в строке в пределах ленты. В данном случае t1 = t2 = 2 и t = 5. Для диагональной матрицы t = 1. При решении системы уравнений с ленточной матрицей участвуют только те коэффициенты, которые расположены в пределах ленты. Число арифметических операций, необходимых для решения системы алгебраических уравнений с полностью заполненной матрицей методом Гаусса, при больших n имеет порядок n3. В то же время для ленточной матрицы при t1 = t2 и t1 ! n он составляет

Для примера ленточной матрицы обратимся к задачам предыдущего раздела, но с пятью узлами и шестью элементами на рис. 2, г. Аналогично (11) матрица K будет иметь коэффициенты flt . По смыслу flt они отличны от нуля только для тех узлов l, где перемещение узла t вызывает отличную от нуля силу при условии, что остальные узлы, кроме t, неподвижны. Отсюда при нумерации узлов, показанной на рис. 2, г, слева от оси x имеем

Здесь t = 3 и матрица K является трехдиагональной.

При применении метода конечных элементов ширина полосы ленточной матрицы зависит от нумерации узлов. Например, если пронумеровать узлы так, как показано на рис. 2, г справа от оси x, то K примет вид (17). Вообще если элементы имеют несколько узлов, то при t1 = t2 величина t1 равна максимальной по элементам величине наибольшей разности между номерами узлов в отдельном элементе. В первом случае нумерации узлов слева на рис. 2, г t1 = 1, а при нумерации справа t1 = 2.

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

В прямоугольной системе координат u1 , u2 на рис. 4 уравнение прямой будет u2 = u1 tg a + g, где a - угол между прямой и положительным направлением оси u1 , g - отрезок отсекаемый прямой на оси u2 . Уравнения (18) описывают две прямые на рис. 4, а решение (18) представляет собой координаты точки пересечения этих прямых. Здесь

Если a1 = a2 и прямые параллельны, то решение системы (18) не существует и она является вырожденной. Если a1 и a2 различаются мало, то система близка к вырожденной. При этом незначительные изменения углов a1 , a2 сильно скажутся на координатах точки пересечения прямых, то есть на решении. Таким образом, плохая обусловленность объясняется тем, что система является почти вырожденной.

В качестве примера обратимся к (18). Пусть +(1) = +(3) = + и +(2) @ +, то есть элемент 2 на рис. 2, а значительно более жесткий, чем элементы 1 и 3. При этом tg a1 © tg a2 и система (18) почти вырожденная. В данном случае разумно изменить постановку задачи и считать элемент 2 абсолютно жестким по сравнению с элементами 1 и 3. Это позволяет объединить узлы 1 и 2 в один узел, который обозначим 12, и приложить к нему суммарную силу P12 = P1 + P2 . Если в (18) положить u1 = u2 = u12 , вычесть из первого уравнения второе и после преобразований пренебречь + по сравнению с +(2), то задача сведется к одному уравнению 2+u12 = P12 .

ПРИМЕР КОНТИНУАЛЬНОЙ ЗАДАЧИ

Обратимся к задаче (3) для одного элемента. В общем случае задания q(r)(x) она является континуальной задачей. Для простоты положим с(r) = 1, l (r) = 1, и опустим индекс r, тогда задача (3) будет

- u" = q(x), u(0) = u(1) = 0,

где штрихи означают дифференцирование по x. Согласно схеме метода конечных элементов разобьем интервал 0, 1 на элементы, соединенные в узлах xi (i = 0, 1, _, n + 1) (рис. 5). Будем разыскивать приближенное решение (19) среди функций семейства с конечным числом параметров в виде

u(x) = u0j0(x) + u1j1(x) + _ + unjn(x) + un + 1jn + 1(x).

Здесь u(x) приближенно представлена линейной комбинацией некоторых функций ji(x) с коэффициентами (параметрами) ui = u(xi) - неизвестными значениями искомой функции в узлах xi . Для того чтобы в (20) u(xi) = ui во всех узлах xi , функции ji(x) должны удовлетворять условиям ji(xi) = 1 и ji(xj) = 0 для всех узлов xj при j ? i. Кроме того, чтобы выполнялись граничные условия (19), следует в (20) положить u0 = un + 1 = 0. В остальном функции ji(x), которые носят название пробных, можно выбирать в довольно широких пределах. Общие требования к ним состоят в возможности выполнить процесс построения приближенного решения и на основе (20) при n ? осуществить сколь угодно точно соответствующую аппроксимацию любой функции, среди которых разыскивается решение задачи. Очевидно, выбор ji(x) играет важнейшую роль как в отношении трудоемкости расчета, так и точности результата. Метод конечных элементов оперирует в качестве ji(x) кусочно-полиномиальными функциями, отличными от нуля в пределах небольшого числа элементов вблизи узла xi . Именно это делает метод максимально эффективным. Поскольку u(x) по своему физическому смыслу должна быть непрерывной функцией, выберем ji(x) в виде кусочно-линейных функций-"домиков", отличных от нуля на двух элементах (см. рис. 5). Каждая такая функция ji(x), i = 1, 2, _, n, равна единице в xi и нулю во всех остальных узлах. При этом набор функций u(x) в (20) будет состоять из непрерывных функций линейных в пределах элементов с изломами в узлах и определяемых своими узловыми значениями ui , i = 1, 2, _, n. На концах интервала 0, 1 они обращаются в нуль. Каждую из таких функций можно изобразить в виде ломаной линии.

Остается определить ui в (20). Это можно сделать по-разному путем приближенного удовлетворения уравнению в (19). Однако, поскольку уравнение в (19) содержит u", а уже u' в (20) терпит разрывы непрерывности в узлах, воспользуемся следующим приемом. Обозначим R(x) = u"(x) + q(x) невязку уравнения в (19). Точное решение дает R(x) = 0, и, следовательно,

для любых функций j(x), которые носят название тестовых. Поскольку разыскивается приближенное решение в форме (20) и для него, как правило, R(x) ? 0, то выполнение тестового условия (21) на базе (20) невозможно. Смягчим выполнение условия (21), потребовав, чтобы оно выполнялось только для n функций jj(x), которые совпадают с пробными. Такой прием носит название метода Гал╦ркина. Выполним в (21) интегрирование по частям при условии j(x) = jj(x) и jj(0) = = jj(1) = 0, тогда вместо (21) получим

Теперь уже в задачу (22) входит u' и можно подставить u из (20) в (22), что дает систему линейных алгебраических уравнений относительно ui вида (16) с коэффициентами fi j и свободными членами Pi

Здесь fi j = fj i и матрица K симметричная, что характерно для метода Гал╦ркина. Для простоты примем длину элементов одинаковой и равной h. Согласно рис. 5, наклон функции ji равен 1/ h на интервале xi - 1 , xi и -1/ h на интервале xi , xi + 1 .

Кроме того, произведение отлично от нуля только при j = i, j = i 7 1, когда соответствующие два элемента, которые несут на себе функции ji и jj , перекрываются (см. рис. 5). В противном случае Если i = j, то

Аналогично fi j = -1/ h при j = i 7 1. Следовательно, матрица K в данном случае оказывается трехдиагональной:

Согласно (23), интегрирование в Pi совершается только на двух соседних элементах. Решение полученной системы алгебраических уравнений дает ui и позволяет представить приближенное решение в форме (20).

В данном примере непрерывность пробных функций ji позволила воспользоваться (22). Кроме того, все функции ji отличны от нуля на разных интервалах 2h, что делает их существенно различными и построенную на их основе при помощи (23) систему линейных алгебраических уравнений невырожденной. Более того, матрица (24) оказалась ленточной и каждое уравнение связывает не более трех неизвестных в соседних узлах. Полученное приближенное решение (20) в виде ломаной линии хорошо аппроксимирует решение задачи при достаточно больших n.

Таким образом, для континуальной задачи метод конечных элементов осуществляет приближенный переход к дискретной задаче на основе (20) и соответствующих кусочно-полиномиальных функций ji , отличных от нуля на нескольких соседних элементах, содержащих узел xi . Дальнейшие процедуры метода конечных элементов для континуальной и дискретной задач в основном совпадают. Здесь, так же как и в случае дискретной задачи, можно выполнить построение матрицы K(r) для типового элемента и из них в процессе соединения элементов в систему сформировать матрицу K для всей системы. Аналогично формируются и свободные члены уравнений. Алгоритм метода конечных элементов особенно эффективен для решения двух- и трехмерных задач, где проявляются основные преимущества этого метода.

ЛИТЕРАТУРА

1. Стренг Г., Фикс Дж. Теория метода конечных элементов. М.: Мир, 1977. 349 с.

2. Courant R. // Bull. Amer. Math. Soc. 1943. Vol. 49. P. 1-43.

3. Turner M., Clough R., Martin H., Topp L. // J. Aeronaut Sci. 1956. Vol. 23, ╧ 9. P. 805-823.

4. Зинкевич О., Морган К. Конечные элементы и аппроксимация. М.: Мир, 1986. 318 с.

5. Розин Л.А. Стержневые системы как системы конечных элементов. Л.: Изд-во ЛГУ, 1976. 232 с.

Рецензент статьи Ю.Г. Мартыненко

* * *

Леонид Александрович Розин, доктор физико-математических наук, профессор, зав. кафедрой строительной механики и теории упругости Санкт-Петербургского государственного технического университета, заслуженный деятель науки и техники РФ. Область научных интересов - численные методы решения задач механики деформируемых систем. Автор более 160 статей и семи монографий.