Вектор нагрузки 3D

Теперь вычислим вектор нагрузки для трёхмерного случая. Вектор нагрузки связан с интегралом произведения функции источника на пробную функцию. Рассмотрим интеграл для одного тетраэдра с вершинами

Пробная функция на тетраэдре имеет вид $\upsilon_{(i)(i+3)}(x, y, z)$ $= q_i \cdot \phi_i(x, y, z)$ $+ q_{i+1} \cdot \phi_{i+1}(x, y, z)$ $+ q_{i+2} \cdot \phi_{i+2}(x, y, z)$ $+ q_{i+3} \cdot \phi_{i+3}(x, y, z)$. Для вычисления интеграла необходимо интерполировать функцию источника $f(x, y, z)$ на тетраэдре линейной функцией. Представим $f(x, y, z)$ в виде $\widetilde{f}(x, y, z)$ $= e_1$ $+ e_2 \cdot x$ $+ e_3 \cdot y$ $+ e_4 \cdot z$. Коэффициенты $e_1, e_2, e_3, e_4$ определяются из условий интерполяции в узлах

Поскольку функции «крышки» $\phi_n$ равны единице в своём узле и нулю в остальных, линейный интерполянт функции источника, выраженный через них, принимает вид

где $f_i, f_{i+1}, f_{i+2}, f_{i+3}$ — значения функции источника в вершинах тетраэдра.

Подставим интерполированную функцию источника и пробную функцию в (6.45$$\int_{\text{тет}} f(x, y, z) \cdot \upsilon \,dV.$$). Примем во внимание (6.30$\displaystyle \phi_i(x, y, z)$ $\displaystyle = a_i$ $\displaystyle + b_i \cdot x$ $\displaystyle + c_i \cdot y$ $\displaystyle + d_i \cdot z$$\displaystyle \phi_{i+1}(x, y, z)$ $\displaystyle = a_{i+1}$ $\displaystyle + b_{i+1} \cdot x$ $\displaystyle + c_{i+1} \cdot y$ $\displaystyle + d_{i+1} \cdot z$$\displaystyle \phi_{i+2}(x, y, z)$ $\displaystyle = a_{i+2}$ $\displaystyle + b_{i+2} \cdot x$ $\displaystyle + c_{i+2} \cdot y$ $\displaystyle + d_{i+2} \cdot z$$\displaystyle \phi_{i+3}(x, y, z)$ $\displaystyle = a_{i+3}$ $\displaystyle + b_{i+3} \cdot x$ $\displaystyle + c_{i+3} \cdot y$ $\displaystyle + d_{i+3} \cdot z$) из раздела о функциях «крышек»

где индексы $m$ и $n$ пробегают вершины тетраэдра $\{i, i$ $+1, i$ $+2, i$ $+3\}$.

Воспользуемся соотношениями (6.31$\displaystyle \int_{\text{тет}} \phi_m \cdot \phi_n \,dV$ $\displaystyle = \begin{cases} \frac{\displaystyle V_{\text{тет}}}{\displaystyle 10}, & m = n\\ \frac{\displaystyle V_{\text{тет}}}{\displaystyle 20}, & m \neq n \end{cases}$) для интегралов произведений функций «крышек», выведенными при вычислении матрицы демпфирования, где $V_{\text{тет}}$ — объём тетраэдра, который вычисляется по формуле (6.15$\displaystyle V_{\text{тет}}$ $\displaystyle = \frac{\displaystyle |\Delta|}{\displaystyle 6},$). После приведения подобных получаем

Введём обозначения для элементов локального вектора нагрузки тетраэдра

Таким образом, локальный вектор нагрузки для трёхмерного тетраэдрального элемента имеет вид

Заметим, что полученный вектор нагрузки совпадает с произведением локальной матрицы демпфирования (6.33$$\begin{aligned}\mathbf{C}_{\text{тет}} = \begin{bmatrix} c_{(i)(i)} & c_{(i)(i+1)} & c_{(i)(i+2)} & c_{(i)(i+3)}\\ c_{(i+1)(i)} & c_{(i+1)(i+1)} & c_{(i+1)(i+2)} & c_{(i+1)(i+3)}\\ c_{(i+2)(i)} & c_{(i+2)(i+1)} & c_{(i+2)(i+2)} & c_{(i+2)(i+3)}\\ c_{(i+3)(i)} & c_{(i+3)(i+1)} & c_{(i+3)(i+2)} & c_{(i+3)(i+3)} \end{bmatrix} = \frac{\displaystyle V_{\text{тет}}}{\displaystyle 20} \begin{bmatrix} 2 & 1 & 1 & 1\\ 1 & 2 & 1 & 1\\ 1 & 1 & 2 & 1\\ 1 & 1 & 1 & 2 \end{bmatrix}.\end{aligned}$$) на вектор узловых значений источника, то есть $\mathbf{R}_{\text{тет}}$ $= \mathbf{C}_{\text{тет}} \cdot \mathbf{f}$.

Глобальный вектор нагрузки $\mathbf{R}$ получается путём суммирования вкладов от всех тетраэдральных элементов сетки методом сборки: элементы локальных векторов добавляются к соответствующим элементам глобального вектора согласно глобальной нумерации узлов. Размерность глобального вектора нагрузки равна $N$, где $N$ — общее количество узлов сетки. Важно отметить, что для внутренних узлов происходит суммирование вкладов от всех смежных тетраэдральных элементов, содержащих данный узел. Для граничных узлов элементные вклады собираются так же, но при последующем учёте граничных условий соответствующие компоненты системы могут изменяться.