# Trouble with creating correct element matrices in Finite Element Analysis for a cantilever beam

I'm trying to solve for displacements of a cantilever beam numerically with FEA. The beam is modeled as a 3D-solid made up of a set of 8-noded hexahedral elements, which are in their undeformed state, constructed as cubes. Every element has the same volume before solving for the displacements, where each side is 1 cm. See picture below:

I apply a uniformly distributed load at the 4 outermost nodes on the top of the beam in the positive z-direction (global coordinates). And the result I get is this horribly deformed structure!

Obviously there is something faulty with my solving, the beam should simply bend without these awful artifacts. I will now try to explain the things I have considered so far, bare with me:

Node numbering

For each element, I have numbered the 8 nodes (global node numbers) counter-clockwise starting from the bottom left doing a full rotation of the "bottom face" and then the same for the "top face", (see attached image below). With this numbering schema I've managed to get a positive Jacobian for every element stiffness matrix. It does become negative, as expected if the numbering is instead done clock-wise. However the Jacobian is very small, in fact it is almost constant for every element stiffness matrix, being for example 1.250000000000003e-07, only differing slightly with respect to the last decimal.

Eigenvalues of element matrices

When testing for the eigenvalues of a stiffness matrix many of the element stiffness matrices have a total of 6 eigenvalues that are zero, corresponding to the rigid body displacements stated in an answer to a similar post here. In that post the answerer states that:

If you get more than six, your stiffness matrix doesn't have the necessary rank and is therefore defective.

My element matrices are 24x24, as each node has 3 DOF's. Should there be a total of 6 eigenvalues that are zero if I have modeled everything only with 3 degrees of freedom? I reckon that the additional 3 values are in respect to the rotational degrees of freedom as well, even though I have only modeled my problem with the 3 translational DOF's.

I'm creating the element matrices with a total 2 Gauss points in each of the 3 directions: since the shapes functions are linear the number of Gauss points should be enough.

for i = 1:numgp
for j = 1:numgp
for k = 1:numgp
xi = gp(i);
eta = gp(j);
zeta = gp(k);
%In here I compute Jacobian, strain matrix
mat = mat + w(i)*w(j)*w(k)*B'*c_mat*B*J;
end
end
end


and compute the derivatives in natural coordinates for each of the 8 shape functions: for example the derivative of the first shape function:

dN1dXI = -(1/8)*(eta-1)*(zeta-1);
dN1dETA = -(1/8)*(xi-1)*(zeta-1);
dN1dZETA = -(1/8)*(xi-1)*(eta-1);


And construct the gradient matrix as follows:

%Gradient matrix of shape functions N1-N8
GN = [
dN1dXI, dN2dXI, dN3dXI, dN4dXI, dN5dXI, dN6dXI, dN7dXI, dN8dXI;
dN1dETA, dN2dETA, dN3dETA, dN4dETA, dN5dETA, dN6dETA, dN7dETA, dN8dETA;
dN1dZETA, dN2dZETA, dN3dZETA, dN4dZETA, dN5dZETA, dN6dZETA, dN7dZETA, dN8dZETA;
];

%Jacobian matrix
%C is a 8x3 matrix with the global (cartesian) coordinates
%the columns x y and z where the i:th row is the i:th node in the element.
Jmat = GN*C;
%Jacobian
J = det(Jmat)

Bmat = Jmat\GN;

B1x = Bmat(1, 1); B2x = Bmat(1, 2); B3x = Bmat(1, 3); B4x = Bmat(1, 4);
B5x = Bmat(1, 5); B6x = Bmat(1, 6); B7x = Bmat(1, 7); B8x = Bmat(1, 8);

B1y = Bmat(2, 1); B2y = Bmat(2,2); B3y = Bmat(2,3); B4y = Bmat(2, 4);
B5y = Bmat(2, 5); B6y = Bmat(2, 6); B7y = Bmat(2,7); B8y = Bmat(2,8);

B1z = Bmat(3, 1); B2z = Bmat(3, 2); B3z = Bmat(3, 3); B4z = Bmat(3, 4);
B5z = Bmat(3, 5); B6z = Bmat(3, 6); B7z = Bmat(3, 7); B8z = Bmat(3,8);

B1 = [B1x 0 0; 0 B1y 0; 0 0 B1z; 0 B1z B1y; B1z 0 B1x; B1y B1x 0];
B2 = [B2x 0 0; 0 B2y 0; 0 0 B2z; 0 B2z B2y; B2z 0 B2x; B2y B2x 0];
%..and so on for all the 8 nods

B = [B1 B2 B3 B4 B5 B6 B7 B8];


This will however not result in numerically symmetric element stiffness matrices. They are almost symmetric but differ slightly in the values causing them not to be exactly symmetric.

Boundary conditions

Before solving for the displacements I simply remove the first 48 rows and columns, in my assembled (global stiffness matrix), that correspond to the Dirichlet boundary conditions of the cantilever beam, that is, there are a total of 16 nodes (each having 3 translational DOF's) that should have zero displacements.

I appreciate any help I can get, I could include more code snippets but as you might understand there is a whole lot of code so I'll keep it rather stripped down as for now.

• From your plots it looks like all elements are the same. So each element matrix should be the same There should be 6 zero eigenvalues and the non-zero eigenvalues should be the same for all elements. It appears you have a basic programming error. Your impression is correct that, on this site, uploading code and expecting others to debug it for you is strongly frowned upon. Jun 7, 2019 at 23:21
• I suggest that you first solve for a single element under uniaxial stress. You can compare your matrix with the values presented in this post: scicomp.stackexchange.com/a/28475/9667 Jun 8, 2019 at 3:41
• Also, this scicomp.stackexchange.com/questions/27093/… may be relevant. Jun 8, 2019 at 7:52
• The stiffness matrices should be identical and symmetric to a very close tolerance-- around 1e-16 assuming double precision computations. If this is not the case, you can limit your debugging to a single element. Jun 8, 2019 at 13:14
• I suggest you solve the 1-element problem suggested by @nicoguaro above. Jun 10, 2019 at 10:19