Stiffness Matrix for Linear Tetrahedral Element

[Hassan2]

Dear all,

I am writing my finite element code for linear elasticity analysis with linear tetrahedral elements. I am not new to the finite element method and my code already works well with element types of triangles, quadrangles in 2D and hexahedrons in 3D but I don't get the correct result with tetrahedrons. I doubt the calculation of the element stiffness matrix and need a valid stiffness matrix of a single tetrahedron for comparison. Does anyone has hands on such software ?

Attachment contains the stiffness matrix and its eigenvalues from my code for a unit tetrahedron whose base are at (1,0,0), (0,1,0), and (0,0,1) and whose apex at (0,0,0), with Young's modulus of 200 MPa and Poisson's ratio of 0.3.

Your help is highly appreciated.

Hassan

 

[AlephZero]

This has an example on page 17: http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf

 

[Hassan2]

Thank you very much AlephZero.

My code calculates the matrix correctly. When I run the code for a beam whose mesh is fine enough, the displacement vectors look fine but their magnitudes are scaled down with a factor usually between 1.0 and 2.0 and the factor depends on the mesh and loading. I don't think the error is due to element type.

 

[AlephZero]

A linear tetrahedron with exact integration is not a good element for continuum mechanics. You can probably find examples where the errors are much worse than 50%, if you search the literature. It is a better element for scalar field problems like heat transfer, especially for nonlinear problems where the gradient discontinuities between the elements may be physically realistic (e.g. they correspond to changes of state in the material).

On the other hand the quadratic (10 node) tetrahedron is very well behaved for structural analysis, if you want to use an automatic tetrahedral mesh generator.

 

[Hassan2]

That could be the reason. I have been spending several days on this problem. Before trying higher order tetrahedrons, I should explain that in 2D, I didn't find much difference between triangular and quad elements so I didn't expect much differences between tetrahedrons and beam elements too.

My test model is a cantilever beam with rectangular cross-section. The original mesh is constructed with beam elements and I believe the mesh is fine enough. By dividing each element into 5 tetrahedrons, I expect the results to be nearly the same but the difference is more than 50 %.

 

[AlephZero]

I should explain that in 2D, I didn't find much difference between triangular and quad elements.

That seems rather surprising, unless you were lucky in your choice of test problems, or your linear quad element was formulated to avoid problems like shear locking etc.

By dividing each element into 5 tetrahedrons, I expect the results to be nearly the same but the difference is more than 50 %.

If you think how five constant-strain tetrahedra could approximate the strain variation in one linear brick element, (or in the actual cantilever beam) you might see why there is a big difference.