Is it possible to give the closed-form of the stiffness matrix of triangular prism element?

Question

Suppose the vertices of triangular prism are 1, 2, 3 (bottom triangle), 4, 5, 6 (top triangle), the coordinates of the vertices are $$ (x_i, y_i, z_i), i=1,2,3,4,5,6. $$

We can write the stiffness matrix of FEM of the triangular prism as follows:

$$ K_{i,j} = \int_0^1\int_0^{1-\xi}\int_{-1}^1 \begin{pmatrix}\frac{\eth N_i}{\eth\xi} & \frac{\eth N_i}{\eth\eta} & \frac{\eth N_i}{\eth\zeta}\end{pmatrix} (\textbf{J}^{-1})^T\textbf{J}^{-1}\begin{pmatrix}\frac{\eth N_j}{\eth\xi} \\ \frac{\eth N_j}{\eth\eta} \\ \frac{\eth N_j}{\eth\zeta}\end{pmatrix}det(\mathbf{J})d\xi d\eta d\zeta\,\\ i,j=1,2,...,6. $$

where

$$ \textbf{J}=\begin{pmatrix}\frac{\eth x}{\eth\xi} & \frac{\eth y}{\eth\xi} & \frac{\eth z}{\eth\xi} \\ \frac{\eth x}{\eth\eta} & \frac{\eth y}{\eth\eta} & \frac{\eth z}{\eth\eta} \\ \frac{\eth x}{\eth\zeta} & \frac{\eth y}{\eth\zeta} & \frac{\eth z}{\eth\zeta} \end{pmatrix} $$

and $$ \textbf{N}=\begin{pmatrix}N_1 & N_2 & N_3 & N_4 & N_5 & N_6\end{pmatrix}^T\\ \left\{\begin{array}{} N_1=\frac{1}{2}(1-\zeta)(1-\xi-\eta)\\ N_2=\frac{1}{2}\xi(1-\zeta)\\ N_3=\frac{1}{‌​2}\eta(1-\zeta)\\ N_4=\frac{1}{2}(1+\zeta)(1-\xi-‌​‌​\eta)\\ N_5=\frac{1}{2}\xi(1+\zeta)\\ N_6=\frac{1}{2}\eta(1+\zeta) \end{array} \right.\ \\ \left\{\begin{array}{} x=\mathbf{N^T}\mathbf{x_e}\\ y=\mathbf{N^T}\mathbf{y_e}\\ z=\mathbf{N^T}\mathbf{z_e}\\ \end{array} \right.\ \\ \left\{\begin{array}{} \mathbf{x_e}=\begin{pmatrix}x_1 & x_2 & x_3 & x_4 & x_5 & x_6\end{pmatrix}^T\\ \mathbf{y_e}=\begin{pmatrix}y_1 & y_2 & y_3 & y_4 & y_5 & y_6\end{pmatrix}^T\\ \mathbf{z_e}=\begin{pmatrix}z_1 & z_2 & z_3 & z_4 & z_5 & z_6\end{pmatrix}^T \end{array} \right.\ $$

Is it possible to give the closed-form of the integration by the powerful mathematica?

Thanks, Tang Laoya

Answer 1

So with the following code:

n = 1/2 {(1 - \[Zeta]) (1 - \[Xi] - \[Eta]), \[Xi] (1 - \[Zeta]), \
\[Eta] (1 - \[Zeta]), (1 + \[Zeta]) (1 - \[Xi] - \[Eta]), \[Xi] (1 + \
\[Zeta]), \[Eta] (1 + \[Zeta])};
xe = {x1, x2, x3, x4, x5, x6};
ye = {y1, y2, y3, y4, y5, y6};
ze = {z1, z2, z3, z4, z5, z6};
x = n.xe;
y = n.ye;
z = n.ze;
J = D[{x, y, z}, {{\[Xi], \[Eta], \[Zeta]}, 1}]\[Transpose];
s[i_] := {D[n[[i]], \[Xi]], D[n[[i]], \[Eta]], D[n[[i]], \[Zeta]]};
k[i_, j_] := s[i].Inverse[J.J\[Transpose]].s[j] Det[J];

Executing k[1,1] yields an enormous expression which is 38 pages long when printed. Integrating it would be extremely difficult, unless there is some hidden way to simplify it. So I'm not sure it's possible.