FEM Stiffness Matrix is always close to Singular or Badly Scaled

My comments apply specifically to a 3D structural finite element but the principles are applicable to more general elements.

A test you will want to execute as part of element development is to calculate the eigenvalues and eigenvectors for the element stiffness matrix. For a 3D structural element, you should get exactly six eigenvalues that are zero (or very close to zero). The eigenvectors for these should be rigid body displacements in the three directions and rigid body rotations about the three axes (or linear combinations of these motions). If you get more than six, your stiffness matrix doesn't have the necessary rank and is therefore defective. It sounds like this might be your problem. (However, in general, if you get less than six zero-eigenvalues, your formulation is also defective.)

So what can cause this problem?

The most likely cause is that the shape functions are invalid in some way. A second possibility is that you are not accurately integrating the terms in the stiffness matrix. In my previous post on a similar topic (3D Solid 8 Node FEM Matlab Code), I alluded to the problems with using too few integration points. Specifically, the problem is that this can cause zero-frequency (zero eigenvalue) displacement modes, often referred to as "hourglass" modes due to the shape of the eigenvector.

The simple solution is to use enough integration points to integrate the stiffness matrix terms exactly. There are also clever numerical tricks that can be employed to remove these hourglass modes but still allow the under-integration. Also, interestingly, sometimes a mesh of defective elements with appropriate boundary conditions will stabilize each other so that a solution of the complete model is possible. But, in general, this is not something you want to rely on.