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.