Stiffness matrix in finite element method

[Hassan2]

Dear all,

I have written a code for structural analysis using the finite element method. For some reason, I directly started with 3D elements ( hexahedron). I used to believe that the code was fine but recently i realized that the results ( deformation, natural frequency,..) strongly depend on the mesh size even for very small meshes. I have been trying for several days to find the problem in the code in vain. The problem might be in matrix assembling but before that, I need to make sure that my elemental stiffness matrix is correct. Does anyone have a reference with a numerical example of 8-node 3D elemental stiffness matrix?

The figure shows the matrix construed by the code. Perhaps you can judge the correctness by of the matrix by looking at the entries.

Your help would be appreciated.

Thanks

 

[AlephZero]

Some basic checks that you should do:

1. Find the eigenvalues of an element stiffness matrix. There should be 6 zeros corresponding to the rigid body motions. Do this check for a rectangular brick shape, and also for an arbitrary distorted shape.

2. For a rectangular brick shaped elements, simulate some "constant stress field" sitations. For example if there is a constant stress in the x direction, the extension in x should correspond to an axial bar with axial stiffness = EA/L, and there should be a contraction in the other two directions according to Poisson's ratio (Run the checks for Posson's ratio = 0, and also for a typical nonzero value like 0.3). Work out the correct displacement vector for the 24 degrees of freedom, multiply by the stiffness, and check you get the correct nodal loads (i.e. equal and opposite loads in the x direction and zero in y and z). Repeat in for constant stress in the Y and Z directions. (Note, there is a reason for doing the test this way rather than just restraining the structure and solving the equations: doing it this way removes the possibiltiy that you restrained the structure incorrectly, and therefore you had both errors in the element formulation and errors in the FE model!).

Then do the same checks for a constant shear stress.

You need to do these sort of checks anyway, but I'll make a guess at the real cause of your problem (ignoring the possibility of bugs in your code). I'm guessing the element has linear shape functions (and is isoparametric, for distorted element shapes). What integration rule are you using to form the stiffness?

 

[Hassan2]

I use Gauss integration rule and with two points I get the accurate result.

The eigenvalues obtained for the element is as following. Here I used the diagonal mass matrix. It means a 24x24 matrix with all diagonal elements equal to the element mass. I hope the mass matrix is fine. I see some odd things in the eigenvalues, like lack of symmetry. Any thoughts?

Added: The element is cubic.

-0.00000
-0.00000
0.00000
0.00000
0.00000
0.00000
0.02654
0.02654
0.03539
0.03539
0.03539
0.07962
0.07962
0.07962
0.07962
0.07962
0.07962
0.10616
0.15924
0.15924
0.15924
0.15924
0.15924
0.15924

 

[AlephZero]

Hm... as you say it looks a bit strange that for a cube you have a pair of eigenvalues 0.026 and a single one 0.106, but it's not necessarily wrong.

Check the eigenvector for the 0.106 mode. You might find the displacements have have a lot of "internal symmetry" and you can't create a different shape that is "symmetrical" with it. An example would be if every node was moving radially in and out from the center of the element. Another possibility is when 4 nodes (forming the vertices of an equilateral tetrahedron) move radially out, and the other 4 nodes move radially inwards.

I use Gauss integration rule and with two points I get the accurate result.

That seems like a good idea, but it isn't. Elememts with linear shape functions work fine for applications like heat transfer with only one degree of freedom per node, but they don't work very well for structural mechanics without a bit of tweaking. The "obvious" formulation that you tried suffers from what is called "shear locking" which makes is much to stiff (possibly orders of magnitude too stiff!) if the structure is trying to bend. The simplest cure for that is to use a 1-point integration rule, but that introduces a different problem called "hourglassing" where the FE model can have spurious zero-energy deformation modes.

One way to fix this is to split the stiffness matrix into two parts and use different integration rules for each part. If you write the stiffness as
$$\iiint B^T D B\,dV$$
you can split the 6x6 material properties matrix $D$ into 3x3 sub-matrices. The two diagonal sub-matrices represent the relationship between the direct stresses and strains (call that $D_\epsilon$), and the shear stresses and strains $D_\gamma$. The two off-diagonal blocks are conveniently zero, for isotropic matierals. So, you can write the element stiffness matrix as
$$\iiint B^T (D_\epsilon + D_\gamma) B\,dV = \iiint B^T D_\epsilon B\,dV + \iiint B^T D_\gamma B\,dV = K_\epsilon + K_\gamma$$
and then integrate $K_\epsilon$ using 2-point integration and $K_\gamma$ using 1-point integration. (You also partition the 24x6 B matrix into two 24x3 matrices).

The mass matrix doesn't suffer from this problem, and a 2x2 itegration rule should work fine.

It means a 24x24 matrix with all diagonal elements equal to the element mass.

That's not quite right. The diagonal elements should be 1/8 of the element mass, for an 8 node element. (But that doesn't affect the pattern of eigenvalues of the stiffness, so long as all the mass matrix diagonals are equal).

Keywords to Google are "shear locking", "hourglasing", and "selective reduced integration" for more on this. The same problem affects 4-node plane elements, and you will probably find more explanations and pictures of those than of 3D elements.