I am trying to model a basic solid cube that rests on a ground (z = 0) with gravity pointing in the negative z direction. To do this, I have a basic cube that is divided into 6 tetrahedrons. The entire solid is made up of these basic cubes so that the tetrahedrons are just repeated.
So far, I managed to find the global element stiffness matrix Ke and global element load vector Fe but I have a question about that. However, my main problem is assembling them into the global stiffness matrix K and global load vector F for the entire solid, and then doing some modifications to take into account the floor before using a solver to solve for KU = F.
Let me summarise the steps I've done.
1. Element stiffness matrix and load vector for each tetrahedral element
Reference tetrahedron Consider firstly a reference tetrahedron that is made up of nodes (0,0,0), (1,0,0), (0,1,0) and (0,0,1). The general formula for the volume of a tetrahedron made of nodes (x_i, y_i, z_i) , i = 1,2,3,4 is given by \begin{equation} \mathcal{V} = \frac{1}{6} \det \begin{pmatrix} 1&1&1&1\\x_1&x_2&x_3&x_4\\y_1&y_2&y_3&y_4\\ z_1&z_2&z_3&z_4 \end{pmatrix} \end{equation} so that in the case of the reference tetrahedron, \begin{equation} \mathcal{V} = \frac{1}{6} \det \begin{pmatrix} 1&1&1&1\\x_1&x_2&x_3&x_4\\y_1&y_2&y_3&y_4\\ z_1&z_2&z_3&z_4 \end{pmatrix} \end{equation} so that in the case of the reference tetrahedron, \begin{equation} \mathcal{V}_{ref} = \frac{1}{6} \end{equation}
The displacement of any point in the tetrahedron can be expressed as a combination of the displacements of each node of the tetrahedron as such:
\begin{equation} U_h(x,y,z) = N(x,y,z) d_e \end{equation}
where $U_h \in {R}^3, N(x,y,z) \in \mathcal{M}_{3 \times 12}({R}), d_e \in {R}^{12}$, \begin{equation} N(x,y,z) = \begin{pmatrix} N_1 & 0 & 0 & & N_4 & 0 & 0\\ 0 & N_1 & 0 &\dots& 0 & N_4 & 0\\ 0 & 0 & N_1 & & 0 & 0 & N_4\\ \end{pmatrix} \end{equation} $N_i$ is called a shape function and be found easily by geometry for a reference tetrahedron.
\begin{equation} N_i = \frac{1}{6 \mathcal{V}} (a_i + b_i x + C_i y+ d_i z) \end{equation} The coefficients are given as \begin{align} a_i & =\quad \det \begin{pmatrix} x_j & y_j & z_j \\ x_k & y_k & z_k \\x_l & y_l & z_l\\ \end{pmatrix}\\\notag b_i & = -\det \begin{pmatrix} 1 & y_j & z_j \\1 & y_k & z_k \\1& y_l & z_l\\ \end{pmatrix}\\\notag c_i & = -\det \begin{pmatrix} y_j & 1& z_j \\y_k & 1& z_k \\ y_l & 1&z_l\\ \end{pmatrix}\\\notag d_i & = -\det \begin{pmatrix} y_j & z_j & 1\\ y_k & z_k &1\\ y_l & z_l &1\\ \end{pmatrix}\\\notag \end{align} With $N_i$ all known, we can calculate the strain matrix $B = LN $ which will help us in calculating the strain in three dimensions, $\epsilon = LU = LN d_e = B d_e$.
The stress relationship with the displacement $U$ is already known. \begin{align} \epsilon_{xx} &= \frac{\partial u}{\partial x}\\\notag \epsilon_{yy} &= \frac{\partial v}{\partial y}\\\notag \epsilon_{zz} &= \frac{\partial w}{\partial w}\\\notag \epsilon_{xy} &= \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\\\notag \epsilon_{xz} &= \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}\\\notag \epsilon_{yz} &= \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}\\\notag \end{align} By inspection, the matrix $L$ is \begin{equation} L = \begin{pmatrix} \partial x & 0 & 0\\ 0 & \partial y & 0\\ 0 & 0 & \partial z \\ 0 & \partial z & \partial y\\ \partial z & 0 & \partial x\\ \partial y & \partial x & 0\\ \end{pmatrix} \end{equation} Hence we have the strain matrix $B = LN$
\begin{equation} B = \frac{1}{6\mathcal{V}} \begin{pmatrix} b_1 & 0 & 0 & b_2 & 0 & 0 & b_3 & 0 & 0 & b_4 & 0 & 0\\ 0 & c_1 & 0 & 0 & c_2 & 0 & 0 & c_3 & 0 & 0 & c_4 & 0\\ 0 & 0 & d_1 & 0 & 0 & d_2 & 0 & 0 & d_3 & 0 & 0 & d_4\\ c_1 & b_1 & 0 & c_2 & b_2 & 0 & c_3 & b_3 & 0 & c_4 & b_4 & 0\\ 0 & d_1 & c_1 & 0 & d_2 & c_2 & 0 & d_3 & c_3 & 0 & d_4 & c_4\\ d_1 & 0 & b_1 & d_2 & 0 & b_2 & d_3 & 0 & b_3 & d_4 & 0 & b_4\\ \end{pmatrix} \end{equation} Plugging in the coordinates of the reference tetrahedron to find $a_i, b_i, c_i d_i, i = 1,2,3,4$, B for the reference tetrahedron is \begin{equation} B = \frac{1}{6\mathcal{V}_{e,ref}} \begin{pmatrix} -1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ -1 & -1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\ -1 & 0 & -1 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ \end{pmatrix} \end{equation} One remark is that the strain matrix is constant within the tetrahedron, hence the stress as well. Finally, we can obtain the element stiffness matrix for a reference tetrahedron. \begin{align} k_{e,ref} &= \int_{\mathcal{V}_{e,ref}} B^T C B dV\\ & = V_{e,ref} B^T C B \end{align} where $C$ is the elasticity matrix \begin{equation} C = \frac{E}{(1+ \nu)(1-2 \nu)} \begin{pmatrix} 1-\nu &\nu & \nu & 0 & 0 & 0\\ \nu & 1-\nu & \nu 0 & 0 & 0 \\ \nu & \nu & 1-\nu 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0\\ 0 & 0 & 0 &0 & 0 &\frac{1-2\nu}{2} \end{pmatrix} \end{equation} The element load vector (global right?) for a body force $(b_x, b_y, b_z)$ is simply \begin{equation} f_e = \frac{1}{4}\rho g \mathcal{V} \begin{pmatrix} 0\\ 0\\ 1\\ \vdots \\0\\0\\1 \end{pmatrix} \end{equation}
So the $f_e$ is actually the global element load vector, but $k_{e, ref}$ is only the local element stiffness matrix. To get the global element stiffness matrix, I decided that I can either calculate $B$ separately for each tetrahedron, or I could find a transformation matrix $T$ such that \begin{align} d_e &= TD_e \\ f_e &= TF_e \\ \end{align} and the element equations would be \begin{align} k_{e,ref} d_e &= f_e \\ K_e D_e &= F_e \\ \end{align} where \begin{align} K_e &= T^T k_{e,ref} T \\ F_e &= T^T f_{e} \\ \end{align}
My first question is, how do I find this transformation matrix $T$? As the matrix $k_{e,ref}$ is a 12 by 12 matrix, it is hard to figure it out. Here is what I've tried.
Let $X_0 = (0,0,0) , X_1 = (1,0,0), X_2 = (0,1,0), X_3 = (0,0,1)$ be the coordinates of the reference tetrahedron. Let $X_0', X_1',X_2', X_3'$ be the coordinates of the arbitrary tetrahedron. Then \begin{equation} T^{-1} \begin{pmatrix} X_0^T\\ X_1^T\\ X_2^T\\ X_3^T\\ \end{pmatrix} = \begin{pmatrix} X_0'^T\\ X_1'^T\\ X_2'^T\\ X_3'^T\\ \end{pmatrix} \end{equation} Here, I'm stuck. To move on, I just calculated B for every single tetrahedron to find the global element sitffness matrix $K_e$ for each tetrahedral element.
2. Assembling the global stiffness matrix $K$ and global load vector $F$.
To assemble the global $K$ and $F$, it suffices to add up each contribution of $K_e$ and $F_e$ by a node in an element.
Pseudo algorithm:
-for node in the globalnodelist
-----for elem in globalelementlist
---------if node is found in elem
----------------run through all other nodes of this element and add up the contribution from Ke_{i,j} where i,j are the local node positions
I do a similar algorithm to add up the load vector. This step should not be much of a problem as it is well documented in textbooks.
3. Taking into account the floor as a Dirichlet boundary condition##
Now that I have $KU = F$, I still have to modify it to take into account the floor (z= 0). My cube is in the positive z space and is resting on z= 0, with gravity pointing in the negative z direction. To do this, I decided that for all nodes $U_i$ with z coordinate = 0, I will set $U_{i,z} = 0$ and this translates to converting the entire row $K(i,:) = 0$ and $F(i) = 0$, and set $K(i,i) = 1$. Is this correct?
Once this is done, I then try to solve for $U$ using a QR solver. My results look terrible. Is there something I am missing? The results look like the nodes at $z=0$ have definitely stayed at $z=0$, but the other nodes with $z>0$ are passing through the floor. Is there a possible problem of reaction forces lacking in my setup? If so, how do I fix that?
Sample of K modified if $U_1$ is at $z=0$, so that $U_{1,z} = 0$
\begin{equation} K = \begin{pmatrix} 0 & 0 & 1 & \dots & 0 \\ \vdots & \dots & \dots &\dots & \vdots \\ \end{pmatrix} \end{equation}
(note: in black are the nodes at the initial positions. In red are the ending positions of each node once solved) Result for a 5 by 5 cube
Thank you in advance if anyone can help point me in the right direction!
