Finite element analysis

Igor Kokcharov

FINITE ELEMENT METHOD

The theories of elasticity, plasticity, plates and other analytical theories can be used to solve many engineering problems. Frequently, practical engineering problems cannot be solved analytically due to complexity of the structure's geometry and boundary conditions. The simple examples given in A, B and C can be solved to obtain inner stresses and displacements with analytical methods. More complicated geometries such as the propeller in example D is usually treated with a numerical method such as finite element method (FEM).

FEM is applied in the following manner:

1. Identify the problem, sketch the structure and loads.

2. Create the geometry with the FE package solid modeler or a CAD system.

3. Mesh the model.

4. Apply boundary conditions (constraints and loads) on the model.

5. Solve numerical equations.

6. Evaluate the results.

Steps 1, 2, 3, 4 are known as preprocessing, the solution of equations in step 5 is the processor and step 6 is considered postprocessing.
The FE model is normally subdivided into finite elements of a specific and simple shape. A typical 3D finite element may be a brick or a wedge with nodes representing the vertices. The displacement of the element is determined by nodal displacements and simple polynomial shape functions that describe the assigned shape of the element. The strains and stresses are calculated by the unknown nodal displacements. Once the nodal displacements are known, element stresses and strains can be calculated.

The most difficult and lengthy step of FEM is the preprocessing, or creating the finite element model. This step includes defining and generating the mesh and applying the correct loading and displacement boundary conditions. Automatic meshing is not always simple, especially in very small features or at the edges and corners. It can be difficult to apply boundary conditions that correspond to the real situation. However, FEM solvers that process the equations in step 5 work automatically and can be rather fast depending on the number of nodes. Powerful and robust visualization tools can allow for a very thorough analysis in step 6.

Degrees of freedoms

Degrees of freedoms (DOF) are associated with each unknown nodal displacement. Each node of a 3D tetrahedral element has 3 DOF representing 3 translational motions. The equations of equilibrium are assembled in a matrix form. Problems with well over 100,000 DOF can be solved with a notebook computer.
Equilibrium equations refer to the equilibrium of each node in each direction:

The sum of all forces at an axis is equal to zero.

The sum of inner forces is equal to the sum of external forces.

The number of nodes is usually bigger than the number of elements for structured 3D models. The number of degrees of freedom is 3 times the number of nodes less the number of kinematic boundary conditions.

The stiffness matrix [K] is the relationship between the vectors of nodal displacements {D} and forces {F}. The stiffness matrix is a diagonal-dominant matrix and is symmetric. It solves for nodal displacement given the loading scheme.

FINITE ELEMENTS

The mathematical model is subdivided by finite elements which are connected to each other with nodes. Forces act at the nodes. The finite element is not a rigid body, the model assumes stresses and strains exist inside the finite elements.

There are a few commonly-used finite elements: beam A, truss B, thin shell C, 2D or 3D continuum (D). It is possible to use several types of elements in one model.

Reliability of the FEA predictions depends on the number of finite elements. If the inner stresses do not vary greatly then the number of elements does not have a significant effect on the accuracy.

Solid elements may be linear (first-order elements) or parabolic (second-order elements). Linear elements have corner nodes only and their edges are straight. The minimum number of nodes for 3D elements is 4. Parabolic elements can have a node placed centrally along each edge and therefore the edges are parabolic. Given the same number of elements, the higher order elements are more accurate because they have a more complex mathematical formulation to describe the element shape (shape function). Also, they represent curved geometry more accurately. An analysis involving higher order elements requires more computational resources.

There are 3 degrees of freedom in a node and 8 nodes for a brick (hexahedral) element. This means that 24 nodal displacements and 24 nodal forces must be considered. The size of the stiffness matrix that relates the nodal displacement vector with the nodal forces vector is [24*24].

The stiffness matrix components are inversely proportion to the modulus of elasticity. Zero modulus of elasticity means that there is no finite elements. Division by zero modulus of elasticity leads to numerical errors in the FEM procedures. Infinity modulus of elasticity means that a part of structure is absolutely rigid.

Although in the theory of elasticity the tensile stress in the crack tip is equal to infinity, all stresses are finite in FEM.

Long elements can be used if there is not a large gradient of displacements, strains and stresses. Such elements can be used in fields of uniform stress but not for stress concentrations. The stress gradient is high but still finite in zones of stress concentration.

If a structure and its loads are symmetrical about the axis the problem can be solved using axisymmetric 2D finite elements.

Five tetrahedral elements are enough to form a cube. Parabolic pyramid elements provide results that are at least as accurate as linear brick elements.

MESHING

A key step in finite element analysis procedure is to mesh the model. Meshing is process of breaking the model into small pieces (finite elements). The network of nodes and elements is called a mesh.

There are two broad types of mesh-generation methods: structured and unstructured meshes.

A structured mesh B is formed by grid-based subdividing of the geometry.

Unstructured mesh are formed automatically. The size of neighboring elements can be significantly different for unstructured mesh. There are no "rows and columns" for such mesh, A. There are more nodes than elements for unstructured mesh with a large number of elements. The ratio between elements and nodes is approximately 2:1 for 2D unstructured mesh and 6:1 for 3D unstructured mesh with tetrahedral elements.

Smaller mesh size h corresponds to a larger number of finite elements in the model. The calculation time increases exponentially as size decreases. The errors decrease for finer mesh but never fall to zero since FEM is always an numerical approximation.

A linear element requires a finer mesh than a parabolic (quadratic) or a cubic element.
Structured mesh B is preferable over an unstructured mesh, A.
Rectangular 4 node elements, C are more preferable than triangular elements, B.
Quadratic (second order) triangular elements, D have at least the same accuracy as first order 4 node elements, C.
Rectangular 8 node elements, E are preferable over triangular second order elements, D despite their larger size.
Cubic displacement approximation F does not need fine meshing.

The FEM is an approximate method. The accuracy of the predictions depends on the assumptions made within the element types and the mesh. A fine mesh is required where there are stress and strain gradients (rates of change). A coarse mesh can be used in areas of reasonably constant stress or regions that are not of user's interest. Users must be able to identify regions of stress concentration. Points of interest may consist of fracture points of previously tested structure, holes, fillets, corners, contact zones, complex details, and high stress areas.

The accuracy decreases if the sizes of the neighboring elements near stress concentrators are significantly different.

The shape of finite elements affects the accuracy. It is preferable not to have sharp corners in finite elements. Elements with similar sides produce smaller errors.

The FE mesh is built without gaps between elements. Both triangular and rectangular elements can be used in the same FE model.

The nodes are numbered sequentially for manual meshing. It is forbidden to build four node elements with an obtuse (> 180^o) inner corner.

BOUNDARY CONDITIONS

Application of boundary conditions is the most critical stage of the finite element analysis.

To simulate the constraints imposed on the physical motion of the structure, displacement boundary conditions A, B must be defined. Prescribed displacements can have zero A or non-zero B values. There are also load boundary conditions, C.

The boundary conditions (fixation or force for a direction) are applied at the nodes only. Maximum number of boundary conditions for a node is equal to the degrees of freedom: 3 restraints or forces for a node in this example.

Great care must be taken so that the finite element model is neither under constrained nor over-constrained.

It is not possible to fix all degrees of freedom (all nodal displacements) for an element. It is better to remove the element from the model.

It is not possible to fix a node and to apply force in the same direction.

Absence of restraints along an axis can lead to a shift along the axis due to errors in numerical calculations. The correct boundary conditions must have at least one restraint for each axis.

Different sets of load and displacement boundary conditions can be implemented to represent model tension, pure bending or shear.

There are three planes of symmetry in this example. If the cube is more flexible than the platens then there is no need to model the entire structure. All points of the upper cube surface will have virtually the same vertical displacements.

St Venants principle:

Two sets of statically equivalent forces produce the same stress field at distance that is large compared to linear dimensions of cross section: b > a.

This principle is often utilized to replace complex boundary conditions with statically equivalent loads.

The figure shows two equivalent loading schemes.

Tensile stresses are frequently the reason for failure in a structure. If the region of maximum tensile stress extends beyond the region of the applied force it is not necessary to have a very fine mesh in that region. There are compressive stresses in region of applied force.

This is an example how the distributed load are spread over the nodes. The sum of the force is equal to 18 kg for a half of the plate. The force can be distributed as the weight over 6 nodes.

DEFORMATION

Most finite element procedures are based on the "displacement method". From the law of equilibrium, the sum of the forces (internal and external) on a node must equal zero. The unknown variables are the displacements. The following is the matrix form of the equilibrium equations:

[K] {D} = {F}

[K] = global stiffness matrix;
{D}= displacement vector;
{F} = load vector.

The stiffness matrix [K] is symmetrical about the diagonal. There are two main types of solvers: direct and iterative. Direct solvers are usually based on Gaussian elimination technique. The direct solvers are more robust but can be slow and require large amounts of disk space for very large problems. Iterative solvers can be extremely fast and require small disk space.

The matrix equation is solved for the displacement vector {D}. Strains can be computed from the displacement results.

The figure shows how the strain depends on the nodal displacements. There is a shift along axis x by 5mm : e_x = 0. Displacement v increases along axis y: e_y is positive.

For first order elements the tensile strain is defined by the difference between corresponding displacements.

There is a linear function for displacements inside the 4-node element. This means that the maximum displacement is only in nodes. The strain is constant for the element. The strain function is quadratic for the second order element shown.

The post-processor visualizes the static or animated description of the deformed shape. The deformed shape of the structure is obtained by summing the nodal coordinates and the nodal displacements multiplied by factor k. The deformed shape helps us to understand many things about the structure such as the position of a region under deflection, where maximum distortion takes place, the accuracy of the applied restraints, and other features of the structural deformation.

The figure shows the initial and deformed shape of a thin plate. The magnification factor is k=1. It is possible to obtain large deformation even if the material properties remains linear (elastic). The problem was solved by a (geometrically) nonlinear structural procedure. The step-by-step loading solution is implemented in this case. The low bending rigidity assists in large linear deformation. The bending rigidity is low if the elastic modulus or thickness t are small. Theoretically, the size of the finite elements do not effect rigidity.

The deformed shape can help the user to decide if the boundary conditions are properly prescribed. No rotations is allowed at edges 1 and 2, only 1 rotation componet at edges 3 and 4.

All displacements are linear for these four-nodes plane elements. There are no gaps between the finite elements. The edges can be polygonal lines if the number of finite elements in the model is large.

This figure shows the deformed shapes for different loading schemes. The uniform stress and deformation fields can be obtained for the second loading scheme. It is better if the loads at the edges are half of those in the center.

The maximum deformation is in the point where the force is applied. The maximum shear stress in the body occurs in the surface element under the force. The rigid plate redistributes force in the body. The contact stress in the body is smaller for the second example. The calculated value of maximum shear stress depends on the size of the elements.

The left end of the shaft is fixed in all nodes. Two nodes on the opposite end are shifted by 2 mm. This loading situation corresponds to torsion. There is a stress irregularity at the right end. If the ends of the shaft are made of relatively rigid material, the local deformation will be distributed over the entire shaft. The stress field is more homogeneous for the second instance.

ACCURACY

Numerical modeling requires approximations and simplifications. The results of FEA are not error-free. Using such a powerful tool as a "black box" without proper understanding of the its principal features can lead users into serious mistakes. Unfortunately human error is inevitable.

An design engineer must understand:

- Which FE analysis is appropriate for an engineering problem;
- What part of the structure must be studied in detail;
missing value

Formulation errors take place if the finite elements do not precisely describe the behavior of the physical problem. Selecting the proper element type and mesh size will reduce formulation errors. Formulation of boundary conditions is critical for the analysis. The success of FEA depends on the how closely the boundary conditions, geometry, and material behavior of the model match the actual situation.

FEM's approximation of a real engineering structure with a finite number of finite elements, as well as the size and shape of the elements can cause discretization error.

Numerical errors are usually rare in comparison with the discretization and formulation errors.

Displacements are the primary unknowns. The FEM solution is usually obtained as a vector of nodal displacement {v}. The solution at other locations throughout the element is generally calculated by interpolation. After approximating the displacement field with shape functions the strain and stress can be calculated. This means that accuracy is at a maximum for nodal displacements.

Strain is calculated by determining the difference in displacement between corresponding points. This is why the accuracy for strain and stress is at a maximum in the central part of the finite elements.

The figure shows a stress pattern for bending. The theoretical and numerical solutions coincide at the center point of the elements.

The type and number of finite elements affect the accuracy of the modeling.

For the nonlinear analysis where force is calculated the number of elements in the model affects the value of the force for small numbers only. The force becomes more stable with an increase in the number of elements.

The are two methods to increase the accuracy for the solution shown:

H-method increases the number of elements, h is the length of one side of an element.
P-method increase the polynomial order of the element.

Second-order elements, B demonstrate the higher accuracy. For bending of the cantilever beam one element over the height of the beam is not the best choice. The more parabolic elements in the length, the better.

Mesh A is preferred for the problem of computational fluid dynamics (CFD).

A. Structured rectangular mesh. Numerical approximations are centered about the centroid of the rectangle element.
B. Structured deformed mesh. There are many elements where numerical approximations are not centered (or symmetrical).
C. Block-structured mesh. Elements are condensed at one point along the straight border.

Nonlinear structural analysis predicts the behavior of the mild steel specimen. The coarse mesh at the left end does not have a significant effect on results for the region with the finer mesh such as the value of maximum stress, the strain pattern in the net section Q-Q, or deformed shape of the central hole. The condition of force applied on the left side of the model roughly agrees with the real situation. The errors in the numerical diagram "Force - Displacement F-v" are possible due to the coarse mesh at the left end of the model.

Errors due to numerical analysis increases if small values of volumetric strain e_V is multiplied by large values of the bulk modulus K. Such analysis requires additional tests. The situation occurs when Poisson's ratio approaches a value of 0.5.