0
$\begingroup$

I am making a code for an 18-node (3x3x2) 3D element FEM. However, even though I am (pretty) sure that all the shape functions are correct and whatnot, whenever I try and invert the stiffness matrix to solve for displacement, I get the warning message from Matlab telling me that my matrix is close to singular.

Is there a checklist of things I should test for to fix this problem? What kind of error commonly causes these kind of errors? Or is it just too broad to even answer meaningfully?

EDIT: Thanks for all your input. Sorry for responding so late, but here is the link to my code. I originally didn't want to post it, it seemed a bit much for me to ask strangers to subject themselves to such annoying debugging.

Improve this question
$\endgroup$
6
  • $\begingroup$ @Additional_Pylons when you try to solve for the displacement Do you invert explicitly the stiffness matrix? From a numerical point of view you avoid always to invert matrix, but is preferred solve the linear system. $\endgroup$
    – Mauro Vanzetto
    Nov 30, 2016 at 13:21
  • $\begingroup$ Which PDE are you solving? Laplace/Poisson? Helmholtz? $\endgroup$
    – coolguy1000000
    Nov 30, 2016 at 13:50
  • 1
    $\begingroup$ Do you have boundary conditions? $\endgroup$
    – nicoguaro ♦
    Nov 30, 2016 at 18:55
  • $\begingroup$ What kind of finite element is this? A hexahedron? To me this sounds like a problem with reduced integration, errors with the integration or not applying the boundary conditions in the right way. $\endgroup$
    – P. G.
    Dec 1, 2016 at 23:58
  • 1
    $\begingroup$ The derivatives of shape functions should always be 'hardcoded'. Deriving the shape functions while running the code is not efficient, especially when calculating bigger structures. Maybe you should post the code to get some more help. $\endgroup$
    – P. G.
    Dec 3, 2016 at 13:49

2 Answers 2

Reset to default
7
$\begingroup$

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.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thank you! I did not know about the eigenvalue test for stiffness matrices. I knew that something was wrong with the code, but I just didn't know what. Turns out, the problem was with Matlab. Because the element was a 18-node 3D element, I simply input the shape functions with symbols and had Matlab diff() and then later subs() the Gauss Point values. From there, I assume many various slight calculation errors caused it to be singular. I have thus fixed it and it gives reasonable values, but I still get 14 zero eigenvalues. Perhaps I need more Gauss points? Would that help? $\endgroup$
    – Additional Pylons
    Dec 3, 2016 at 8:22
  • $\begingroup$ What integration orders are you using in the three directions? $\endgroup$
    – Bill Greene
    Dec 3, 2016 at 12:46
  • $\begingroup$ I am using just 2-point Gauss point integration in all three directions. I understand that 2-point integration is enough for 3 node elements, and since my element is 3x3x2, it should be fine. $\endgroup$
    – Additional Pylons
    Dec 4, 2016 at 14:05
  • $\begingroup$ Okay, I will definitely test for that. And while I did check the Jacobian at all points and noticed nothing awry, I would like to ask what entails "something wrong" with the Jacobian? Since I am using a 3x3x2 element, I don't think the Jacobian should be symmetric, like for the 8-node case. Is there some checklist I should know about? $\endgroup$
    – Additional Pylons
    Dec 5, 2016 at 3:20
  • $\begingroup$ The Jacobian simply scales the parametric to real-space coordinates. So, for a simple "rectangular" element, it should be diagonal with the terms easily-calculated from the side lengths. It should be identical to your 8-node element for the same size element. $\endgroup$
    – Bill Greene
    Dec 5, 2016 at 13:08
2
$\begingroup$

The elemental stiffness matrix must be always singular because while deriving it we do not impose any constraints or boundary conditions. Thus inverting the stiffness matrix to solve for displacements/position-vectors/degrees-of-freedom should yield indeterminate results.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.