Метод конечных элементов

Метод конечных элементов возник из необходимости новых путей решения задач строительной механики и теории упругости. Одним из первых предшественников МКЭ считается метод «ферменной аналогии» (framework method) Александра Хренникова, опубликованный в 1941 году: сплошная упругая среда заменялась эквивалентной системой стержней, допускающей расчёт обычными методами строительной механики^[2].

В 1943 году Рихард Курант в работе по вариационным методам предложил использовать кусочно-линейные функции на треугольных элементах для приближённого решения задач равновесия и колебаний: этот подход, по существу, и является современным МКЭ^[3].

В 1950-х годах Иоаннис Аргирис разработал энергетические теоремы матричного анализа конструкций, а в 1956 году М. Тёрнер, Р. Клаф, Х. Мартин и Л. Топп опубликовали статью о расчёте жёсткости и деформаций сложных конструкций, в которой был описан прямой метод жёсткости — непосредственный предшественник МКЭ в его современной матричной форме. Сам термин «конечный элемент» (finite element) был предложен Рэем Клафом в 1960 году^[4][5].

Хотя перечисленные подходы различаются между собой в деталях, они имеют одну общую черту: дискретизация непрерывной области сеткой на набор дискретных поддоменов, обычно называемых элементами. Дальнейшее развитие метода конечных элементов связано с решением задач авиастроения и космических исследований в 1950-х годах. В Китае в 1950-х годах Фэн Кан предложил численный метод решения дифференциальных уравнений в частных производных для расчёта конструкций плотин; этот метод, названный методом конечных разностей на основе вариационного принципа, рассматривается как ещё одна независимая реализация идеи метода конечных элементов^[6].

В СССР распространение и практическая реализация МКЭ в 1960-х годах связана с именем Леонарда Оганесяна — одного из пионеров конечноэлементных расчётов в нашей стране^[7].

Существенный толчок в своём развитии МКЭ получил после того, как было доказано, что его можно рассматривать как один из вариантов распространённого в строительной механике метода Рэлея — Ритца, который путём минимизации потенциальной энергии сводит задачу к системе линейных уравнений равновесия. После того, как была установлена связь МКЭ с процедурой минимизации, он стал применяться к задачам, описываемым уравнением Лапласа или уравнением Пуассона. Область применения МКЭ значительно расширилась, когда было установлено, что уравнения, определяющие элементы в задачах, могут быть легко получены с помощью вариантов метода взвешенных невязок, таких как метод Галёркина или метод наименьших квадратов. Это сыграло важную роль в теоретическом обосновании МКЭ, так как позволило применять его при решении многих типов дифференциальных уравнений. Таким образом, метод конечных элементов превратился в общий метод численного решения дифференциальных уравнений и систем дифференциальных уравнений^[4][6][8].

Развитие матричных методов анализа конструкций и их превращение в метод конечных элементов подробно прослежено в обзоре К. Фелиппы (2001), выделяющем три периода: становление до Второй мировой войны, период усложнения матричных методов и эпоху прямого метода жёсткости^[9].

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

Область, в которой ищется решение дифференциальных уравнений, разбивается на конечное количество элементов (дискретизация области). В каждом из элементов выбирается определённый вид аппроксимирующей функции (обычно полиномиальной), обеспечивающий требуемую гладкость и сходимость метода; в простейшем случае это полином первой степени, и элементы называются линейными. Вне своего элемента аппроксимирующая функция равна нулю. Значения функций на границах элементов (в узлах) являются решением задачи и заранее неизвестны. Коэффициенты аппроксимирующих функций обычно ищутся из условия равенства значений соседних функций на границах между элементами (в узлах); затем эти коэффициенты выражаются через значения функций в узлах элементов. Составляется система линейных алгебраических уравнений (СЛАУ). Количество уравнений равно количеству неизвестных значений в узлах, прямо пропорционально количеству элементов и ограничивается только возможностями вычислительной техники. Так как каждый из элементов связан с ограниченным количеством соседних, система линейных алгебраических уравнений имеет разрежённый вид, что существенно упрощает её решение^[12].

Если говорить в матричных терминах, то собираются так называемые матрицы жёсткости (или матрицы Дирихле) и матрицы масс. Далее на эти матрицы накладываются граничные условия: например, при условиях Неймана в матрицах не меняется ничего (они формируют так называемые естественные граничные условия), а при условиях Дирихле из матриц вычёркиваются строки и столбцы, соответствующие граничным узлам, так как в силу краевых условий значения соответствующих компонент решения известны. Затем собирается система линейных уравнений и решается одним из известных методов (например, методом сопряжённых градиентов или прямыми методами для разрежённых матриц)^[13][14].

Пусть в одномерном пространстве необходимо решить следующее дифференциальное уравнение для нахождения функции ${\displaystyle u}$ на промежутке от 0 до 1; на границах области значение функции ${\displaystyle u}$ равно 0 (стандартная модельная задача)^[19]:

${\displaystyle {\mbox{ P1 }}:{\begin{cases}-u''(x)=f(x){\mbox{ in }}(0,1),\\u(0)=u(1)=0,\end{cases}}}$

где ${\displaystyle f}$ — известная функция, ${\displaystyle u}$ — неизвестная функция от ${\displaystyle x}$, ${\displaystyle u''}$ — вторая производная от ${\displaystyle u}$ по ${\displaystyle x}$. Знак «минус» перед второй производной обеспечивает положительную определённость итоговой матрицы жёсткости. Решение поставленной задачи методом конечных элементов разбивается на два этапа:

• переформулируем граничную задачу в так называемую слабую (вариационную) форму — на этом этапе вычислений почти не требуется;
• на втором этапе разобьём слабую форму на конечные отрезки-элементы.

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

Если ${\displaystyle u}$ есть решение, то для любой гладкой функции ${\displaystyle v}$, которая удовлетворяет граничным условиям ${\displaystyle v=0}$ в точках ${\displaystyle x=0}$ и ${\displaystyle x=1}$, можно записать следующее выражение:

(1) ${\displaystyle \int _{0}^{1}f(x)v(x)\,\mathrm {d} x=\int _{0}^{1}(-u''(x))v(x)\,\mathrm {d} x.}$

С помощью интегрирования по частям преобразуем выражение (1) к следующей форме:

(2) ${\displaystyle \int _{0}^{1}f(x)v(x)\,\mathrm {d} x=-u'(x)v(x){\bigg |}_{0}^{1}+\int _{0}^{1}u'(x)v'(x)\,\mathrm {d} x=\int _{0}^{1}u'(x)v'(x)\,\mathrm {d} x=\phi (u,v).}$

Оно получено с учётом того, что ${\displaystyle v(0)=v(1)=0}$, поэтому граничный член зануляется.

Разобьём область, в которой ищется решение, и перейдём к конечномерному подпространству ${\displaystyle V_{h}\subset H_{0}^{1}}$: найти ${\displaystyle u_{h}\in V_{h}}$ такое, что для всех ${\displaystyle v_{h}\in V_{h}}$ выполняется ${\displaystyle \phi (u_{h},v_{h})=\int _{0}^{1}f(x)v_{h}(x)\,\mathrm {d} x}$. Пространство ${\displaystyle V_{h}}$ состоит из кусочно-полиномиальных функций, непрерывных на всей области.

Есть много способов для выбора базиса ${\displaystyle V_{h}}$. Выберем в качестве базисных функций такие ${\displaystyle v_{k}}$, чтобы они представлялись прямыми линиями (полиномами первой степени):

${\displaystyle v_{k}(x)={\begin{cases}{x-x_{k-1} \over x_{k}\,-x_{k-1}}&{\mbox{, }}x\in [x_{k-1},x_{k}],\\{x_{k+1}\,-x \over x_{k+1}\,-x_{k}}&{\mbox{, }}x\in [x_{k},x_{k+1}],\\0&{\mbox{, }}x\not \in [x_{k-1},x_{k+1}],\end{cases}}}$

для ${\displaystyle k=1,\ldots ,n-1}$ (в данном примере ${\displaystyle n=5}$).

Если теперь искомое приближённое решение представить в виде ${\displaystyle u_{h}(x)=\sum _{k=1}^{n-1}u_{k}v_{k}(x)}$, а функцию ${\displaystyle f(x)}$ аппроксимировать как ${\displaystyle f(x)=\sum _{k=0}^{n}f_{k}v_{k}(x)}$, то с помощью слабой постановки можно получить следующую систему уравнений относительно искомых коэффициентов ${\displaystyle u_{k}}$:

${\displaystyle \sum _{k=1}^{n-1}u_{k}\phi (v_{k},v_{j})=\sum _{k=0}^{n}f_{k}\int _{0}^{1}v_{k}v_{j}\,\mathrm {d} x,}$

где ${\displaystyle j=1,\ldots ,n-1}$. Матрица этой системы — матрица жёсткости задачи; она трёхдиагональна и положительно определена, поскольку базисные функции с номерами, отличающимися более чем на единицу, имеют непересекающиеся носители^[20][21].

Выбор типа элемента зависит от размерности задачи и требуемой точности. Наиболее распространённые типы^[22]:

• Одномерные: линейные отрезки (2 узла), квадратичные отрезки (3 узла).
• Двумерные:
  • Треугольные элементы (линейные с 3 узлами, квадратичные с 6 узлами). Обеспечивают высокую гибкость при аппроксимации сложных криволинейных границ.
  • Четырёхугольные элементы (билинейные с 4 узлами, биквадратические с 8 или 9 узлами). Обычно обеспечивают более высокую точность при том же числе узлов, но сложнее в генерации сетки для произвольных областей.
• Трёхмерные: тетраэдры (4 или 10 узлов) и гексаэдры (шестигранники, 8 или 20 узлов).

Повышение степени полинома (${\displaystyle p}$-версия МКЭ) позволяет увеличить точность без измельчения сетки (${\displaystyle h}$-версия МКЭ), что особенно эффективно при наличии гладких решений и позволяет достичь экспоненциальной скорости сходимости (${\displaystyle hp}$-версия МКЭ).

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

• Метод Петрова — Галёркина: пространство пробных функций и пространство функций веса различаются. Широко используется для задач с преобладанием конвекции (например, уравнение переноса или Навье — Стокса при высоких числах Рейнольдса), где стандартный метод Галёркина даёт нефизичные осцилляции решения^[23].
• Смешанный метод конечных элементов (Mixed FEM): независимая аппроксимация нескольких физических величин (например, давления и скорости в задачах гидродинамики или напряжений и перемещений в механике). Это позволяет обходить ограничения, накладываемые условием inf-sup (условием Бабушки — Брецци), и избегать эффекта «объёмного защемления» (volumetric locking)^[24].

При уменьшении характерного размера элементов (измельчении сетки) приближённое решение сходится к точному. Для эллиптических краевых задач получены строгие оценки погрешности: при использовании кусочно-полиномиальных элементов степени ${\displaystyle p}$ погрешность в энергетической норме имеет порядок ${\displaystyle O(h^{p})}$, где ${\displaystyle h}$ — максимальный диаметр элемента; в норме пространства ${\displaystyle L_{2}}$ — порядок ${\displaystyle O(h^{p+1})}$. Первые общие оценки погрешности МКЭ принадлежат И. Бабушке (1971); систематическое изложение теории сходимости дано в монографиях Ф. Сьярле и С. Бреннер — Л. Скотта^[18][25][26].

Качество сетки существенно влияет на точность: вырожденные (почти плоские) элементы приводят к потере точности и ухудшению обусловленности матрицы жёсткости. Поэтому для построения сеток применяются алгоритмы на основе триангуляции Делоне, обеспечивающие «почти равносторонние» треугольники^[27].

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

Долгое время широкому распространению МКЭ мешало отсутствие алгоритмов автоматического разбиения области на «почти равносторонние» треугольники (погрешность, в зависимости от вариации метода, обратно пропорциональна синусу самого острого или самого тупого угла в разбиении). Эту задачу удалось успешно решить с помощью алгоритмов, основанных на триангуляции Делоне, что дало возможность создавать полностью автоматические конечноэлементные системы автоматизированного проектирования^[13].

К другим ограничениям относятся: чувствительность решения к качеству сетки, необходимость контроля численных погрешностей при решении плохо обусловленных систем, значительные требования к памяти и вычислительным ресурсам при моделировании трёхмерных задач.

МКЭ является вычислительным ядром большинства современных инженерных CAE-систем (Computer-Aided Engineering).

• Коммерческие пакеты: ANSYS, Abaqus, COMSOL Multiphysics, Nastran, LS-DYNA, SolidWorks Simulation.
• Открытое ПО и библиотеки: FEniCS, FreeFEM, CalculiX, Code_Aster, deal.II, MOOSE.

Эти системы автоматизируют все этапы решения: от импорта геометрической модели и генерации расчётной сетки до сборки глобальных матриц, решения СЛАУ и постпроцессинга (визуализации полей напряжений, температур, скоростей и т. д.).