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.
