Ansys, Abaqus

Eigen — математическая библиотека

Перед тем, как приступить к написанию кода, я хотел бы рассказать о математической библиотеке — Eigen. Это зрелая, надежная и эффективная библиотека, которая обладает очень чистым и выразительным API. В ней есть множество compile-time оптимизаций, библиотека способна выполнять явную векторизацию и поддерживает SSE 2/3/4, ARM и NEON инструкции. Как по мне, эта библиотека замечательна хотя бы тем, что позволяет реализовывать сложные математические вычисления кратким и выразительным образом.
Я хотел бы сделать краткое описание части API Eigen, которое мы будем использовать далее. Те читатели, которые знакомы с этой библиотекой, могут пропустить этот раздел.

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

Плотные матрицы

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

Матрицы фиксированного размера могут быть определены следующим образом:

Eigen::Matrix m;

Где m — float матрица c фиксированный размеров 3 × 3. Также можно использовать следующие предопределенные типы:

Eigen::Matrix2i

Eigen::Matrix3i
Eigen::Matrix4i
Eigen::Matrix2f
Eigen::Matrix3f
Eigen::Matrix4f
Eigen::Matrix2d
Eigen::Matrix3d
Eigen::Matrix4d


Есть немного больше предопределенных типов, эти являются основными. Цифра обозначает размерность квадратной матрицы и буква обозначает тип элемента матрицы: i — целое число, f — число с плавающей точкой одинарной точности, d — число с плавающей точкой двойной точности.

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

Eigen::Matrix v;

Или можно использовать один из предопределенных типов (обозначения те же, что и для матриц):

Eigen::Vector2i

Eigen::Vector3i
Eigen::Vector4i
Eigen::Vector2f
Eigen::Vector3f
Eigen::Vector4f
Eigen::Vector2d
Eigen::Vector3d
Eigen::Vector4d


Для быстрого ознакомления, я написал следующий пример:

#include

int main()

{
//Fixed size 3x3 matrix of floats
Eigen::Matrix3f A;
A << 1, 0, 1,
2, 5, 0,
1, 1, 2;

std::cout << "Matrix A:"

<< std::endl
<< A
<< std::endl;

//Access to matrix element [1, 1]:

std::cout << "A[1, 1]:"
<< std::endl
<< A(1, 1)
<< std::endl;

//Fixed size 3 vector of floats

Eigen::Vector3f B;
B << 1, 0, 1;

std::cout << "Vector B:"

<< std::endl
<< B
<< std::endl;
//Access to vector element [1]:
std::cout << "B[1]:"
<< std::endl
<< B(1)
<< std::endl;

//Multiplication of vector and matrix

Eigen::Vector3f C = A * B;
std::cout << "C = A * B:"
<< std::endl
<< C
<< std::endl;
//Dynamic size matrix of floats
Eigen::MatrixXf D;
D.resize(3, 3);

//Let matrix D be an identity matrix:

D.setIdentity();

std::cout << "Matrix D:"

<< std::endl
<< D
<< std::endl;

//Multiplication of matrix and matrix

Eigen::MatrixXf E = A * D;
std::cout << "E = A * D:"
<< std::endl
<< E
<< std::endl;

return 0;

}


Вывод:

Matrix A:

1 0 1
2 5 0
1 1 2
A[1, 1]:
5
Vector B:
1
0
1
B[1]:
0
C = A * B:
2
2
3
Matrix D:
1 0 0
0 1 0
0 0 1
E = A * D:
1 0 1
2 5 0
1 1 2


Для более подробной информации лучше обратиться к документации.

Разреженные матрицы

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

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



Размер памяти, необходимый для хранения матрицы огромный! Разреженная матрица потребуется в 10000 раз меньшее количество памяти.

Разреженные матрицы используют память более эффективно, так как хранят только ненулевые элементы. Разреженные матрицы могут быть представлены ​​двумя способами: сжатым и несжатым. В несжатом виде легко добавить или удалить элемент из матрицы, но такой вид представления не является оптимальным с точки зрения использования памяти. Eigen, при работе с разреженной матрицей, использует своеобразное сжатое представление, несколько более оптимизированное для динамического добавления/удаления элементов. Eigen также умеет конвертировать такую матрицу в сжатый вид, более того он это делает прозрачно, делать это в явном виде не нужно. Большинство алгоритмов требуют сжатый вид матрицы, и применение любого из этих алгоритмов автоматически преобразует матрицу в сжатый вид. И наоборот, добавление/удаление элемента автоматически конвертирует матрицу в оптимизированное представление.

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

Например, мы имеем следующую разреженную матрицу:

0 3 0 0 0

22 0 0 0 0
0 5 0 1 0
0 0 0 0 0
0 0 14 0 8


Первое, что мы должны сделать, это включить соответствующий заголовок библиотеки Eigen: . Затем, мы объявляем пустую разреженную матрицу размерности 5х5. Далее, мы определяем вектор triplets и заполняем его значениями. Конструктор triplet'а принимает следующие аргументы: i-индекс, j-индекс, значение.

#include

#include

int main()

{
Eigen::SparseMatrix A(5, 5);

std::vector< Eigen::Triplet > triplets;

triplets.push_back(Eigen::Triplet(0, 1, 3));

triplets.push_back(Eigen::Triplet(1, 0, 22));
triplets.push_back(Eigen::Triplet(2, 1, 5));
triplets.push_back(Eigen::Triplet(2, 3, 1));
triplets.push_back(Eigen::Triplet(4, 2, 14));
triplets.push_back(Eigen::Triplet(4, 4, 8));

A.setFromTriplets(triplets.begin(), triplets.end());

A.insert(0, 0);

std::cout << A;

A.makeCompressed();

std::cout << std::endl << A;

}


Вывод будет следующим:

Nonzero entries:

(0,0) (22,1) (_,_) (3,0) (5,2) (_,_) (_,_) (14,4) (_,_) (_,_) (1,2) (_,_) (_,_) (8,4) (_,_) (_,_)

Outer pointers:

0 3 7 10 13 $
Inner non zeros:
2 2 1 1 1 $

0 3 0 0 0

22 0 0 0 0
0 5 0 1 0
0 0 0 0 0
0 0 14 0 8

Nonzero entries:

(0,0) (22,1) (3,0) (5,2) (14,4) (1,2) (8,4)

Outer pointers:

0 2 4 5 6 $

0 3 0 0 0

22 0 0 0 0
0 5 0 1 0
0 0 0 0 0
0 0 14 0 8

Структуры данных

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

struct Element

{
void CalculateStiffnessMatrix(const Eigen::Matrix3f& D, std::vector >& triplets);

Eigen::Matrix B;

int nodesIds[3];
};

Еще одна простая структура которая нам понадобится — это структура для описания закреплений. Это простая структура состоящая из перечисляемого типа, который определяет тип ограничения, и индекса узла на который накладывается ограничение.

struct Constraint

{
enum Type
{
UX = 1 << 0,
UY = 1 << 1,
UXY = UX | UY
};
int node;
Type type;
};


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

• 
  Количество узлов
  
• 
  Вектор с х-координатой узлов
  
• 
  Вектор с y-координатой узлов
  
• 
  Вектор элементов
  
• 
  Вектор закреплений
  
• 
  Вектор нагрузок
  

В коде, мы определим их следующим образом:

int nodesCount;

int noadsCount;
Eigen::VectorXf nodesX;
Eigen::VectorXf nodesY;
Eigen::VectorXf loads;
std::vector< Element > elements;
std::vector< Constraint > constraints;

Чтение входного файла

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

Первое, что мы делаем — создаем потоки для ввода/вывода:

int main(int argc, char *argv[])

{
if ( argc != 3 )
{
std::cout << "usage: " << argv[0] << "

\n";
return 1;
}
std::ifstream infile(argv[1]);
std::ofstream outfile(argv[2]);


В первой строке входного файла находятся свойства материала, прочитаем их:

float poissonRatio, youngModulus;

infile >> poissonRatio >> youngModulus;


Этого достаточно, чтобы построить матрицу упругости изотропного материала для плоско-напряженного состояния, которая определена следующим образом:





Откуда взялись эти выражения? Они легко получаются из закона Гука в общем виде, действительно, мы можем найти выражение для матрицы D из следующих соотношений:



Стоит отметить, что плоско-напряженное состояние означает, что σ_Z равно нулю, но не ε_Z. Плоско-напряженная модель хорошо подходит для решения ряда инженерных задач, где рассматриваются плоские структуры и все силы действуют в плоскости. Во времена, когда расчеты объемных задач были слишком дорогими, множество задач удавалось свести к плоским, жертвуя при этом точностью.

При плоско-напряженном состоянии, ничто не мешает телу деформироваться в направлении нормали к его плоскости, поэтому деформация ε_Z в общем случае не равна нулю. Она не появляется в уравнениях выше, но может быть легко получена из следующего соотношения, учитывая что напряжение σ_Z равно нулю:


У нас есть все исходные данные чтобы задать матрицу упругости, давайте сделаем это:

Eigen::Matrix3f D;

D <<
1.0f, poissonRatio, 0.0f,
poissonRatio, 1.0, 0.0f,
0.0f, 0.0f, (1.0f - poissonRatio) / 2.0f;

D *= youngModulus / (1.0f - pow(poissonRatio, 2.0f));

Далее, мы читаем список с координатами узлов. Сначала читаем количество узлов, затем задаем размер динамических векторов х и у. Далее, мы просто читать координаты узлов в цикле, строка за строкой.

infile >> nodesCount;

nodesX.resize(nodesCount);
nodesY.resize(nodesCount);

for (int i = 0; i < nodesCount; ++i)

{
infile >> nodesX[i] >> nodesY[i];
}

Затем, мы читаем список элементов. Все то-же самое, читаем количество элементов, а затем индексы узлов для каждого элемента:

int elementCount;

infile >> elementCount;

for (int i = 0; i < elementCount; ++i)

{
Element element;
infile >> element.nodesIds[0] >> element.nodesIds[1] >> element.nodesIds[2];
elements.push_back(element);
}

Далее, читаем список закреплений. Все то же самое:

int constraintCount;

infile >> constraintCount;

for (int i = 0; i < constraintCount; ++i)

{
Constraint constraint;
int type;
infile >> constraint.node >> type;
constraint.type = static_cast(type);
constraints.push_back(constraint);
}


Вы могли заметить static_cast, он нужен чтобы преобразовать int тип, в объявленный нами ранее перечисляемый тип для закреплений.

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



Таким образом, для каждого узла, у нас есть два элемента в векторе нагрузки, которые представляют х и у компоненты усилия. Таким образом, х-компонента усилия в определенном узле будет храниться в элементе с индексом id = 2 * node_id + 0, а у-составляющая усилия для этого узла будет храниться в элементе с индексом id = 2 * node_id + 1.

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

int loadsCount;

loads.resize(2 * nodesCount);
loads.setZero();

infile >> loadsCount;

for (int i = 0; i < loadsCount; ++i)

{
int node;
float x, y;
infile >> node >> x >> y;
loads[2 * node + 0] = x;
loads[2 * node + 1] = y;
}
Расчет глобальной матрицы жесткости

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



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

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





Где: u_i и v_i — u-компонента и v-компонента перемещения i-го узла.

Вектор внешних сил:





Где: R_xi и R_yi — х-компонента и y-компонента внешнего усилия приложенного к i-му узлу.

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



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

Как построить матрицу жесткости? Суть глобальной матрицы жесткости есть суперпозиция матриц жесткости каждого элемента. Если мы рассмотрим один элемент отдельно от остальных, то мы можем определить такую-же матрицу жесткости, но только для данного элемента. Эта матрица определяет соотношения между перемещениями узлов и узловыми силами. Например для 3-узлового элемента:





Где: [k]^е — матрица жесткости e-го элемента; [δ]^е — вектор перемещений узлов е-го элемента; [F]^е — вектор узловых сил е-го элемента; i, j, m — индексы узлов элемента.

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

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

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

Если просуммировать все уравнения для каждого узла, меж-элементные узловые силы уйдут. С правой стороны уравнения останутся только внешние силы — нагрузки. Таким образом, учитывая этот факт, мы можем написать:



Рисунок ниже иллюстрирует вышеприведенное выражение:



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

std::vector > triplets;

for (std::vector::iterator it = elements.begin(); it != elements.end(); ++it)
{
it->CalculateStiffnessMatrix(D, triplets);
}


Где std::vector > — вектор triplets; CalculateStiffnessMatrix — метод класса элемента, который заполняет вектор triplets.

Как уже упоминалось ранее, мы можем построить глобальную разреженную матрицу прямо из вектора triplets:

Eigen::SparseMatrix globalK(2 * nodesCount, 2 * nodesCount);

globalK.setFromTriplets(triplets.begin(), triplets.end());

Расчет матрицы жесткости отдельных элементов

Мы выяснили, как собрать глобальную матрицу жесткости (ансамбль) из матриц элементов, Осталось выяснить как получить матрицу жесткости элемента.

Функции формы

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



Где [N] представляет собой матрицу функций координат (х, у). Эти функции называются функциями формы. Каждую компоненту перемещения, u и v можно интерполировать независимо, тогда уравнение выше, можно переписать в следующем виде:



Или можно записать в раздельной форме:




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



Если значения f в узлах известны, то мы можем задать систему из трех уравнений:



Или в матричной форме:



Из этой системы уравнений, мы можем найти неизвестный вектор коэффициентов для интерполирующего выражения:



Введем обозначение



Наконец, получим интерполяционное выражение:



Возвращаясь к перемещениям, мы можем сказать, что:





Тогда, функции формы будут иметь вид:

Деформации и градиентная матрица

Зная поле перемещений, мы можем найти поле деформаций, по соотношениям из теории упругости:



Теперь мы можем заменить u и v, на функции формы:






Или мы можем написать то же самое в комбинированной форме:



Матрица с правой стороны, обычно обозначается как матрица