Перейти к содержанию

Численное интегрирование

Векторы внутренних сил элемента \(\boldsymbol{q}^e, \boldsymbol{Q}^e\) и матрица жёсткости элемента \(\boldsymbol{K}^e\), полученные в разделах Дискретизация внутренней виртуальной работы и Внешняя виртуальная работа и сборка глобального уравнения, имеют вид интегралов по области элемента \(\Omega^e\) или \(\Omega^e_0\). FrontISTR вычисляет их численно с помощью квадратур Гаусса.

Квадратура Гаусса и замена переменных

Квадратура Гаусса аппроксимирует интеграл по эталонной области \(\Xi\) линейной комбинацией значений подынтегральной функции в точках интегрирования \(\boldsymbol{\xi}_i \in \Xi\) с весами \(w_i\). При применении к области элемента \(\Omega^e\) выполняется замена переменных через отображение \(\boldsymbol{x}: \Xi \to \Omega^e\), что даёт

\[ \int_{\Omega^e} f(\boldsymbol{x})\, dv \approx \sum_{i=1}^{n_q} w_i\, f(\boldsymbol{x}(\boldsymbol{\xi}_i))\, J_{\xi_i}, \qquad J_{\xi_i} = \left.\det\!\left(\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{\xi}}\right)\right|_{\boldsymbol{\xi}_i} \]

где \(n_q\) — число точек интегрирования, а \(J_{\xi_i}\) — определитель якобиана преобразования. Эталонная область \(\Xi\) определяется для каждого типа элемента (для гексаэдра — \([-1,1]^3\); для треугольников, тетраэдров и клиньев — соответствующие эталонные формы), а точки интегрирования \(\boldsymbol{\xi}_i\) и веса \(w_i\) задаются числовыми таблицами. Поверхностные интегралы обрабатываются в той же форме посредством отображения грани элемента из двумерной эталонной области.

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

Тип элемента Квадратурная формула Число точек интегрирования
4-узловой тетраэдр (tet4n) 1-точечная формула 1
10-узловой тетраэдр (tet10n) 4-точечная формула 4
6-узловая треугольная призма (prism6n) 2-точечная формула 2
15-узловая треугольная призма (prism15n) 9-точечная формула 9
8-узловой гексаэдр (hex8n) Гаусс—Лежандр 2×2×2 8
20-узловой гексаэдр (hex20n) Гаусс—Лежандр 3×3×3 27
4-узловой четырёхугольник (quad4n) Гаусс—Лежандр 2×2 4
8-узловой четырёхугольник (quad8n) Гаусс—Лежандр 3×3 9
3-узловой треугольник (tri3n) 1-точечная формула 1
6-узловой треугольник (tri6n) 3-точечная формула 3

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

Применение к интегрированию по элементу

При интегрировании по элементу естественные координаты \(\boldsymbol{r}\) используются как координаты эталонной области (\(\boldsymbol{r} = \boldsymbol{\xi}\)), а отображение в физические координаты задаётся интерполяцией узловых координат функциями формы. В зависимости от выбора отсчётной конфигурации (Схема инкрементного анализа) используются следующие формулировки.

Полная лагранжева формулировка (интегрирование по исходной конфигурации \(\Omega^e_0\)): отображение и якобиан имеют вид

\[ \boldsymbol{X}^e(\boldsymbol{r}) = \sum_{\alpha} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{X}^e_\alpha, \qquad J_{r_i} = \left.\det\!\left(\frac{\partial \boldsymbol{X}}{\partial \boldsymbol{r}}\right)\right|_{\boldsymbol{r}_i} \]

а внутренние силы элемента и матрица жёсткости аппроксимируются как

\[ \boldsymbol{Q}^e \approx \sum_{i=1}^{n_q} w_i\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \boldsymbol{S}\, J_{r_i}, \qquad \boldsymbol{K}^e \approx \sum_{i=1}^{n_q} w_i\, \boldsymbol{K}^e_{x}\, J_{r_i} \]

Все величины \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}, \boldsymbol{S}, \boldsymbol{K}^e_{x}\) вычисляются в точке интегрирования \(\boldsymbol{r}_i\).

Обновлённая лагранжева формулировка (интегрирование по текущей конфигурации \(\Omega^e\)): отображение и якобиан имеют вид

\[ \boldsymbol{x}^e(\boldsymbol{r}) = \sum_{\alpha} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{x}^e_\alpha = \sum_{\alpha} N_\alpha^e(\boldsymbol{r})\, (\boldsymbol{X}^e_\alpha + \boldsymbol{u}^e_\alpha), \qquad J_{r_i} = \left.\det\!\left(\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{r}}\right)\right|_{\boldsymbol{r}_i} \]

и используются следующие аппроксимации

\[ \boldsymbol{q}^e \approx \sum_{i=1}^{n_q} w_i\, \boldsymbol{B}_L^T\, \boldsymbol{\sigma}\, J_{r_i}, \qquad \boldsymbol{K}^e \approx \sum_{i=1}^{n_q} w_i\, \boldsymbol{K}^e_{x}\, J_{r_i} \]

Единственное различие между двумя формулировками состоит в том, какие узловые координаты подаются на отображение — \(\boldsymbol{X}^e_\alpha\) или \(\boldsymbol{x}^e_\alpha\); точки интегрирования, веса и структура цикла по точкам интегрирования являются общими.

Полное и пониженное интегрирование

Интегрирование с количеством точек, достаточным для точного интегрирования полиномиальной степени подынтегральной функции, называется полным интегрированием, а интегрирование с числом точек на один уровень меньше — пониженным интегрированием. Пониженное интегрирование используется для уменьшения сдвигового и объёмного запирания, но требует обработки паразитных мод деформации, например мод «песочных часов». Число точек интегрирования для каждого типа элемента и выбор между полным и пониженным интегрированием рассматриваются в разделе Схема нумерации элементов и библиотека функций формы и последующих разделах, а также в разделе Расширенные формулировки элементов.

Связанные разделы