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

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

• What does $\zeta_i$ stand for? You have $\zeta$ as a scalar integration variable, so it can't be referring to vector components. Feb 9, 2014 at 23:58
• oh, sorry, $$\zeta_i$$ are the z-coordinate of standard triangular prism, i.e., $$\zeta_i=-1, i=1,2,3, \\ \zeta_i=1, i=4,5,6.$$ Thanks Feb 10, 2014 at 0:07
• Ok, so $\zeta_i$ is the $z$-coordinate of the $i$th vertex of the prism? Are the $x$ and $y$-coordinates of the vertices not used anywhere, only the $z$-components are used? And what are the $x_i,y_i,z_i$? Are they arbitrary constants? And what does $J^{-T}$ represent? Is it shorthand for $(J^{-1})^T$? Feb 10, 2014 at 0:55
• I have a feeling it won't be possible to get an analytic form, since $J^{-1}$ is going to have entries which are fractions that have polynomials in their denominators, which might be problematic. But I might be wrong. Feb 10, 2014 at 1:01
• Hi DumpsterDoofus, thanks for your kindly reply. The $$x_i, y_i, z_i$$ are coordinates of vertices of original triangular prism. The x and y-coordinates of the vertices are used in the integration. Yes, $$J^{-T}$$ is the shorthand for $$(J^{-1})^{T}$$Thanks Feb 10, 2014 at 1:02

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.

• Thanks for your so beautiful code:) Feb 10, 2014 at 2:28
• Could you please also give me the code to output the whole matrix K by integrating? Thanks Feb 10, 2014 at 2:34
• Yeah, I wish it gave a useful answer but unfortunately it doesn't. Kind of surprising; did you say that the tetrahedron matrix had a closed-form? I'd be curious to see how they calculated it, and why a tetrahedron is so much less complicated than a triangular prism. Feb 10, 2014 at 2:41
• Thanks for your kindly reply. There are many papers about tetrahedron, for example: dspace.uta.edu/bitstream/handle/10106/933/… Feb 10, 2014 at 2:48
• Well, if the tetrahedron case took an entire PhD thesis to find, I imagine the triangular prism case will be similarly hard. Feb 10, 2014 at 4:39