# Stiffness matrix in finite element method

• Hassan2
In summary: B matrix (e.g. the Euler-Lagrange or the Hamilton-Jacobi), and use a different rule for the D matrix (e.g. the Neumann or the Crank-Nicolson). You can also use a "finite difference" rule for the B and D matrices, but I'm not sure what that would be. I suspect that your integration rule is incorrect, and that is why the stiffness matrix is behaving strangely.
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.

Thanks

#### Attachments

• stiffness.JPG
68.4 KB · Views: 2,759
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?

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?

-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

Last edited:
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.

## What is a stiffness matrix?

A stiffness matrix is a square matrix that represents the relationship between the applied forces and resulting displacements in a finite element model. It is a key component in the finite element method, which is a numerical technique for solving complex engineering problems.

## How is a stiffness matrix calculated?

A stiffness matrix is calculated by combining the stiffness values of individual elements in a finite element model. These stiffness values are obtained from the material properties and geometry of each element, and are then assembled into a global stiffness matrix.

## What is the significance of the stiffness matrix in the finite element method?

The stiffness matrix is a crucial part of the finite element method as it allows for the prediction of the behavior of a structure under various loading conditions. It also enables engineers to make design decisions and optimize the performance of a structure.

## What are the properties of a stiffness matrix?

A stiffness matrix is a symmetric, positive definite matrix, meaning it is equal to its transpose and has all positive eigenvalues. It is also a sparse matrix, meaning it has a large number of zero elements due to the localized nature of finite element calculations.

## What are the limitations of using a stiffness matrix in the finite element method?

One limitation of using a stiffness matrix is that it assumes linear behavior of the structure, which may not always be accurate for more complex structures. Additionally, the accuracy of the results depends on the mesh size and element types used in the finite element model.

• Mechanical Engineering
Replies
1
Views
866
• Mechanical Engineering
Replies
1
Views
1K
• Mechanical Engineering
Replies
0
Views
288
• Mechanical Engineering
Replies
3
Views
2K
• Mechanical Engineering
Replies
5
Views
6K
• Mechanical Engineering
Replies
9
Views
2K
• Mechanical Engineering
Replies
3
Views
2K
• Mechanical Engineering
Replies
2
Views
1K
• Mechanical Engineering
Replies
4
Views
1K
• Mechanical Engineering
Replies
1
Views
996