Damping matrix 3D

Let us proceed to computing the damping matrix for the three-dimensional case. The damping matrix is related to the integral of the square of the trial function. Consider the integral for a single tetrahedron with vertices

(xi,yi,zi),\displaystyle (x_i, y_i, z_i),​(xi+1,yi+1,zi+1),\displaystyle (x_{i+1}, y_{i+1}, z_{i+1}),​(xi+2,yi+2,zi+2),\displaystyle (x_{i+2}, y_{i+2}, z_{i+2}),​(xi+3,yi+3,zi+3)\displaystyle (x_{i+3}, y_{i+3}, z_{i+3})
∫tetυ2 dV.\int_{\text{tet}} \upsilon^2 \,dV.
(6.29)

The trial function on the tetrahedron has the form υ(i)(i+3)(x,y,z)\upsilon_{(i)(i+3)}(x, y, z)​=qi⋅ϕi{} = q_i \cdot \phi_i​+qi+1⋅ϕi+1{} + q_{i+1} \cdot \phi_{i+1}​+qi+2⋅ϕi+2{} + q_{i+2} \cdot \phi_{i+2}​+qi+3⋅ϕi+3{} + q_{i+3} \cdot \phi_{i+3}. Analogously to the two-dimensional case, the hat functions for the tetrahedron are linear according to (6.13)

ϕi(x,y,z)\displaystyle \phi_i(x, y, z) =ai\displaystyle {} = a_i​+bi⋅x\displaystyle {} + b_i \cdot x​+ci⋅y\displaystyle {} + c_i \cdot y​+di⋅z\displaystyle {} + d_i \cdot zϕi+1(x,y,z)\displaystyle \phi_{i+1}(x, y, z) =ai+1\displaystyle {} = a_{i+1}​+bi+1⋅x\displaystyle {} + b_{i+1} \cdot x​+ci+1⋅y\displaystyle {} + c_{i+1} \cdot y​+di+1⋅z\displaystyle {} + d_{i+1} \cdot zϕi+2(x,y,z)\displaystyle \phi_{i+2}(x, y, z) =ai+2\displaystyle {} = a_{i+2}​+bi+2⋅x\displaystyle {} + b_{i+2} \cdot x​+ci+2⋅y\displaystyle {} + c_{i+2} \cdot y​+di+2⋅z\displaystyle {} + d_{i+2} \cdot zϕi+3(x,y,z)\displaystyle \phi_{i+3}(x, y, z) =ai+3\displaystyle {} = a_{i+3}​+bi+3⋅x\displaystyle {} + b_{i+3} \cdot x​+ci+3⋅y\displaystyle {} + c_{i+3} \cdot y​+di+3⋅z\displaystyle {} + d_{i+3} \cdot z
(6.30)

Substitute the trial function into (6.29)

∫tetυ(i)(i+3)2 dV\displaystyle \int_{\text{tet}} \upsilon_{(i)(i+3)}^2 \,dV =∫tet[qi⋅ϕi(x,y,z)+qi+1⋅ϕi+1(x,y,z)+qi+2⋅ϕi+2(x,y,z)+qi+3⋅ϕi+3(x,y,z)]2 dV\displaystyle {} = \int_{\text{tet}} \Big[ 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) \Big]^2 \,dV

For linear hat functions on the tetrahedron the following relations hold

∫tetϕm⋅ϕn dV\displaystyle \int_{\text{tet}} \phi_m \cdot \phi_n \,dV ={Vtet10,m=nVtet20,m≠n\displaystyle {} = \begin{cases} \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}, & m = n\\ \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}, & m \neq n \end{cases}
(6.31)

where VtetV_{\text{tet}} is the volume of the tetrahedron, which is computed by formula (6.15).

Expand the square of the trial function

∫tetυ(i)(i+3)2 dV\displaystyle \int_{\text{tet}} \upsilon_{(i)(i+3)}^2 \,dV =∫tet[qi2⋅ϕi2+qi+12⋅ϕi+12+qi+22⋅ϕi+22+qi+32⋅ϕi+32+2⋅qi⋅qi+1⋅ϕi⋅ϕi+1+2⋅qi⋅qi+2⋅ϕi⋅ϕi+2+2⋅qi⋅qi+3⋅ϕi⋅ϕi+3+2⋅qi+1⋅qi+2⋅ϕi+1⋅ϕi+2+2⋅qi+1⋅qi+3⋅ϕi+1⋅ϕi+3+2⋅qi+2⋅qi+3⋅ϕi+2⋅ϕi+3] dV\displaystyle {} = \int_{\text{tet}} \Big[ q_i^2 \cdot \phi_i^2 + q_{i+1}^2 \cdot \phi_{i+1}^2 + q_{i+2}^2 \cdot \phi_{i+2}^2 + q_{i+3}^2 \cdot \phi_{i+3}^2 + 2 \cdot q_i \cdot q_{i+1} \cdot \phi_i \cdot \phi_{i+1} + 2 \cdot q_i \cdot q_{i+2} \cdot \phi_i \cdot \phi_{i+2} + 2 \cdot q_i \cdot q_{i+3} \cdot \phi_i \cdot \phi_{i+3} + 2 \cdot q_{i+1} \cdot q_{i+2} \cdot \phi_{i+1} \cdot \phi_{i+2} + 2 \cdot q_{i+1} \cdot q_{i+3} \cdot \phi_{i+1} \cdot \phi_{i+3} + 2 \cdot q_{i+2} \cdot q_{i+3} \cdot \phi_{i+2} \cdot \phi_{i+3} \Big] \,dV

Apply the formulas (6.31)

∫tetυ(i)(i+3)2 dV\displaystyle \int_{\text{tet}} \upsilon_{(i)(i+3)}^2 \,dV =qi2⋅Vtet10\displaystyle {} = q_i^2 \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}​+qi+12⋅Vtet10\displaystyle {} + q_{i+1}^2 \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}​+qi+22⋅Vtet10\displaystyle {} + q_{i+2}^2 \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}​+qi+32⋅Vtet10\displaystyle {} + q_{i+3}^2 \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}​+2⋅qi⋅qi+1⋅Vtet20\displaystyle {} + 2 \cdot q_i \cdot q_{i+1} \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}​+2⋅qi⋅qi+2⋅Vtet20\displaystyle {} + 2 \cdot q_i \cdot q_{i+2} \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}​+2⋅qi⋅qi+3⋅Vtet20\displaystyle {} + 2 \cdot q_i \cdot q_{i+3} \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}​+2⋅qi+1⋅qi+2⋅Vtet20\displaystyle {} + 2 \cdot q_{i+1} \cdot q_{i+2} \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}​+2⋅qi+1⋅qi+3⋅Vtet20\displaystyle {} + 2 \cdot q_{i+1} \cdot q_{i+3} \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}​+2⋅qi+2⋅qi+3⋅Vtet20\displaystyle {} + 2 \cdot q_{i+2} \cdot q_{i+3} \cdot \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}

Simplify the expression

∫tetυ(i)(i+3)2 dV\displaystyle \int_{\text{tet}} \upsilon_{(i)(i+3)}^2 \,dV =Vtet20\displaystyle {} = \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}​⋅[2⋅qi2+2⋅qi+12+2⋅qi+22+2⋅qi+32+2⋅qi⋅qi+1+2⋅qi⋅qi+2+2⋅qi⋅qi+3+2⋅qi+1⋅qi+2+2⋅qi+1⋅qi+3+2⋅qi+2⋅qi+3]\displaystyle {} \cdot \Big[ 2 \cdot q_i^2 + 2 \cdot q_{i+1}^2 + 2 \cdot q_{i+2}^2 + 2 \cdot q_{i+3}^2 + 2 \cdot q_i \cdot q_{i+1} + 2 \cdot q_i \cdot q_{i+2} + 2 \cdot q_i \cdot q_{i+3} + 2 \cdot q_{i+1} \cdot q_{i+2} + 2 \cdot q_{i+1} \cdot q_{i+3} + 2 \cdot q_{i+2} \cdot q_{i+3} \Big]

We introduce the notation for the elements of the local damping matrix of the tetrahedron

c(i)(i)\displaystyle c_{(i)(i)} =Vtet10\displaystyle {} = \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}c(i+1)(i+1)\displaystyle c_{(i+1)(i+1)} =Vtet10\displaystyle {} = \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}c(i+2)(i+2)\displaystyle c_{(i+2)(i+2)} =Vtet10\displaystyle {} = \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}c(i+3)(i+3)\displaystyle c_{(i+3)(i+3)} =Vtet10\displaystyle {} = \frac{\displaystyle V_{\text{tet}}}{\displaystyle 10}c(m)(n)\displaystyle c_{(m)(n)} =c(n)(m)\displaystyle {} = c_{(n)(m)}​=Vtet20,\displaystyle {} = \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20},​m\displaystyle m​≠n,\displaystyle {} \neq n,​m,n\displaystyle m, n​∈{i,i\displaystyle {} \in \{i, i​+1,i\displaystyle {} +1, i​+2,i\displaystyle {} +2, i​+3}\displaystyle {} +3\}
(6.32)

Thus, the local damping matrix for the tetrahedral element has the form

Ctet\displaystyle \mathbf{C}_{\text{tet}} =[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)]\displaystyle {} = \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}​=Vtet20\displaystyle {} = \frac{\displaystyle V_{\text{tet}}}{\displaystyle 20}​[2111121111211112].\displaystyle \begin{bmatrix} 2 & 1 & 1 & 1\\ 1 & 2 & 1 & 1\\ 1 & 1 & 2 & 1\\ 1 & 1 & 1 & 2 \end{bmatrix}.
(6.33)

The local damping matrix is symmetric. The global damping matrix C\mathbf{C} is obtained by summing the contributions from all tetrahedral elements of the mesh using the assembly method: the elements of the local matrices are added to the corresponding elements of the global matrix according to the global node numbering. The dimension of the global damping matrix is N×NN \times N, where NN is the total number of mesh nodes.