FEM: Rank deficiency and hourglassing

[bda23]

Hello, I am having somewhat difficulty understanding the concepts of rank deficiency and hourglassing in finite element methods. Essentially, I have been reading a book outlining this very briefly on half a page and I need a bit more information. As an example:

If we have a 2D elasticity problem with a bilinear element, the element stiffness matrix "k" will be an 8x8 matrix, composed of the integral of B(T) D B, where B is a 3x8 matrix and D is a 3x3 matrix. So far so good. In 2D, this matrix "k" is said to be of rank 5 because you subtract 3 rigid body motions from this 8x8 matrix. While I see that you can have these 3 motions, I don't quite understand why this reduces the rank of the matrix. Could anyone give me some insight into how this works? I know, this probably sounds ridiculously simple, but somehow I am having problems picturing it.

Secondly, when it comes to hourglassing: If you use one Gauss point only to integrate this "k"-matrix numerically, you have: k = 4J B(T)DB | xi, eta = 0, which makes sense. However, this k-matrix has rank 3, the same as the D-matrix, and I don't really know why it isn't 5, like before. This rank deficiency obviously then leads to 2 hourglassing modes. Does anyone have any advice of how to determine, not just here, but also in other cases, the rank of a stiffness matrix? That would be very helpful, thank you.

 

[AlephZero]

Consider the statics equation Kx = F for a single element. If the displacement vector x corresponds to a rigid body motion, then F must be 0. There are 3 linearly indepedent non-zero vectors x corresponding to the 3 rigid body motions, so the rank of K should be 3 less than the order of K (i.e rank 5 if K is of order 8).

One way to think about hourglassing modes is to consider the strain energy in the element which is $x^T K x / 2 = x^T F / 2$. The energy can also be written in terms of the stresses and strains as $(1/2)\iint \sigma^T \epsilon\, \text{dA}$.

If you use numerical intergation to compute K, you not considering the stresses and strains everywhere in the element, but only at the integration ponits. If you use a one-point integration rule, you are only using the strain at the element centroid. For the hourglassing shapes, the element shape functions give all the strain components at the centroid as zero. So if you calculate Kx = F where x is an hourglass shape, you will get F = 0.

A general way to find the hourglass modes of an element is to find the eignenvalues and vectors of K. The number of zero eigenvalues equals the number of rigid body modes + the number of hourglass modes, and the eigenvectors give the corresponding sets of displacements.

 

[bda23]

Thanks a lot for the explanation! This reminds me of some linear algebra I did a while back: The rank of a linear transformation = dimension of transformation - dim(kernel), is that essentially what you're saying about ranks?

As far as hourglassing goes, I can follow your explanations, makes sense. However, when I turn back to the example I mentioned in my first post, I still don't quite see why the element stiffness matrix "k", when integrated numerically with one GP, is reduced in rank by two. Basically, I understand what you're saying, but I still can't make sense of what the author of the book (where I got this example from) meant, I can't follow his logic. Any help is greatly appreciated. Thanks again for your reply, it was very useful already.

 

[AlephZero]

Thanks a lot for the explanation! This reminds me of some linear algebra I did a while back: The rank of a linear transformation = dimension of transformation - dim(kernel), is that essentially what you're saying about ranks?

Yes.

As far as hourglassing goes, I can follow your explanations, makes sense. However, when I turn back to the example I mentioned in my first post, I still don't quite see why the element stiffness matrix "k", when integrated numerically with one GP, is reduced in rank by two.

When you form K using numerical integration, you only "sample" the stresses and strains at the integration points. So for 1-point integration, you only use the straiin at the element centroid.

Draw a rectangular element and join the mid-points of the opposite sides with two lines that intersect at the center of the element. If the element has linear shape functions, these lines remain straight when the element deforms. The change in length of the lines measures the direct strain components at the center, and the change in angle between the lines measures the shear strain.

You can deform the element in two ways ( the two hourglass modes) so the midpoints of all four sides don't move. and therefore the strains at the integration point in the center are zero. The strains everywhere else in the elements are non-zero, but one-point numerical integration doesn't "see" that. So if you find the strain energy for an hourglass mode using the one-point stiffness matrix, you get zero. Of course for the hourglass modes this is the "wrong" answer, unlike the rigid body modes where the energy should be zero.

If you try to fix this problem by using a different integrration rule, for example 2x2 point gaussian integration, you over-compensate and get a different problem called "shear locking" where the calculated energy is muck bigger than the "right" answer, or in other words the element is much too stiff for some types of deformation. You didn't mention shear locking in your OP, but your book will probably deal with it, and it might help to think about shear locking and hourglassing together, since they are two different symptoms of the same basic problem, which is that linear shape functions don't correctly represent the physical behavior that you want to model.

 

[bda23]

Yes.



When you form K using numerical integration, you only "sample" the stresses and strains at the integration points. So for 1-point integration, you only use the straiin at the element centroid.

Draw a rectangular element and join the mid-points of the opposite sides with two lines that intersect at the center of the element. If the element has linear shape functions, these lines remain straight when the element deforms. The change in length of the lines measures the direct strain components at the center, and the change in angle between the lines measures the shear strain.

You can deform the element in two ways ( the two hourglass modes) so the midpoints of all four sides don't move. and therefore the strains at the integration point in the center are zero. The strains everywhere else in the elements are non-zero, but one-point numerical integration doesn't "see" that. So if you find the strain energy for an hourglass mode using the one-point stiffness matrix, you get zero. Of course for the hourglass modes this is the "wrong" answer, unlike the rigid body modes where the energy should be zero.

If you try to fix this problem by using a different integrration rule, for example 2x2 point gaussian integration, you over-compensate and get a different problem called "shear locking" where the calculated energy is muck bigger than the "right" answer, or in other words the element is much too stiff for some types of deformation. You didn't mention shear locking in your OP, but your book will probably deal with it, and it might help to think about shear locking and hourglassing together, since they are two different symptoms of the same basic problem, which is that linear shape functions don't correctly represent the physical behavior that you want to model.

Thanks a lot, that was very helpful!