Ansys, Abaqus

В наши дни, МКЭ — это наверное самый распространенный метод для решения широкого спектра прикладных инженерных задач. Исторически, он появился из механики, однако впоследствии был применен к всевозможным не механическим задачам.

Сегодня имеется большое разнообразие программных пакетов, таких как ANSYS, Abaqus, Patran, Cosmos, и т.д. Эти программные пакеты позволяют решать задачи строительной механики, механики жидкости, термодинамики, электродинамики и многие другие. Сама реализация метода, как правило считается достаточно сложной и объемной.

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

Картинка для привлечения внимания:


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

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

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

Входные данные

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



На рисунке мы видим девять узлов и восемь элементов. Для описания сетки, нужно два списка:

• 
  Список узлов, который определяет координаты каждого узла.
  
• 
  Список элементов.
  

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

Список узлов:

0.000000 31.500000

15.516667 31.500000
0.000000 19.489759
18.804134 23.248226
0.000000 0.000000
20.479981 11.365822
27.516667 19.500000
27.516667 11.365822
27.516667 0.000000

Список элементов:

1 3 2
2 3 4

4 6 7
4 3 6
7 6 8
3 5 6
6 5 9
6 9 8

Стоит отметить, что есть несколько способов задать один и тот же элемент. Индексы узлов можно перечислять по часовой стрелке или против часовой стрелки. Обычно используется перечисление против часовой стрелки, поэтому мы будем предполагать, что все элементы заданы именно таким образом.

Давайте создадим некий входной файл с полным описанием задачи. Чтобы избежать путаницы и не усложнять, лучше использовать индексацию, которая начинается с нуля, а не с единицы, так как в C/C ++ массивы индексируются с нуля. Для первого тестового входного файла, мы создадим самую простую сетку из возможных:



Пусть первая строка будет описание свойств материала. Например, для стали, коэффициент Пуассона ν = 0,3 и модуль Юнга Е = 2000 МПа:

0.3 2000

Затем следует строка с количеством узлов и затем сам список узлов:

4
0.0 0.0

1.0 0.0
0.0 1.0
1.0 1.0

Затем следует строка с количеством элементов, далее список элементов:

2
0 1 2

1 3 2

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

• 
  1 — запрещено перемещение в направлении оси x
  
• 
  2 — запрещено перемещение в направлении оси y
  
• 
  3 — запрещено перемещение в обоих, х и у направлениях
  

Первая строка определяет количество ограничений:

2
0 3

1 2

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

2
2 0.0 1.0

3 0.0 1.0

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



Входной файл целиком:

0.3 2000
4

0.0 0.0
1.0 0.0
0.0 1.0
1.0 1.0
2
0 1 2
1 3 2
2
0 3
1 2
2
2 0.0 1.0
3 0.0 1.0