Numerical Integration in Finite Element Method
Numerical integration is a method in which approximate solution of definite integral is evaluated using numerical techniques.
In this method the polynomial function is evaluated at some known base points (sampling points) with weighting co-efficient.
There are a wide range of methods available for numerical integration. Most common methods are Newton Cotes and Gauss quadrature.
Gauss quadrature is widely used in finite element method since it produces the most accurate numerical approximations possible.
Gauss Quadrature Method
The basic assumption in Gauss numerical integration is that the polynomial function is evaluated in the interval -1 to +1 and the weight coefficients depend upon the number of base points known as gauss integration points.
w = weight coefficients and are base points.
The order of Gauss integration depends upon the number of gauss points selected. By increasing the number of gauss points, the accuracy of the solution increases but the cost of the solution also increases.
Hence in Gauss quadrature method we select n Gauss points and n weights such that the integration formula provides accurate solution for the polynomial function
Gauss Points and Weights for Gaussian Quadrature
Multiple integration can be done by extending the one-dimensional integration formula.
In finite element formulation, we encounter integrals which are to be integrated for each element individually. These integrals are evaluated numerically using Gauss quadrature formula.
For an example, let’s consider a four node quadrilateral element
Stiffness matrix integral of the element in local coordinate system is given as :
where
B: Strain displacement matrix
D : Material Matrix
te : thickness of element
For numerically integrating the stiffness matrix, we transform the integral in to Gauss integration formula
Now, using one-point, two-point and three-point Gauss integration rule we integrate the stiffness matrix integral.
The order of numerical integration plays a vital role in the accuracy of the solution in finite element analysis. Lower order integration causes the stiffness matrix rank deficient i.e. it will have more zero eigenvalues than the no.of physical rigid body modes. Such modes are termed as spurious zero energy modes.
Higher order integration will increase solution accuracy but will increase the computational cost which will have a significant impact on solution time for 3D analysis.
Case Study
Considering the rectangular element as shown in Fig. Test quadrilateral element , it is a rectangle of base 50 and height 25. The plate has unit thickness and isotropic material with E = 96 and ν = 1/3
Test quadrilateral element
Now we compute the stiffness matrix of the test element in SCIFESOL solver using the above three integration schemes and calculate eigenvalues of the stiffness matrices to observe the effect of order of numerical integration.
In SCIFESOL we can select stiffness matrix and mass matrix integration rules using element properties dialog box. In the script file we can directly modify the global variables ITGRUL and ITGMRUL to specify integration rules for stiffness and mass matrix.
Element Properties dialog box
We plot mode shapes of the element calculated at different integration rules.
One-Point Integration (1x1 rule)
Mode Shapes 1, 2, 3 for 1x1 integration rule
Mode Shapes 4, 5, 6 for 1x1 integration rule
Mode Shapes 7 and 8 for 1x1 integration rule
It is observed that with one-point integration, the rank of the stiffness matrix is only 3 instead of 5 (8 degrees of freedom minus 3 rigid body modes). The element has 2 spurious zero energy modes which leads to ill-conditioning of the stiffness matrix.
TWO-Point Integration (2x2 rule)
Mode Shapes 1, 2, 3 for 2x2 integration rule
Mode Shapes 4, 5, 6 for 2x2 integration rule
Mode Shapes 7 and 8 for 2x2 integration rule
Using two-point integration, the stiffness matrix has correct rank of 5. In three-point integration, the eigenvalues obtained are same as two-point integration due to negligible difference in stiffness matrix for 4 node quadrilateral element.
Eigenvalues of test element
Integration order of the finite element matrices should be chosen with care so that good numerical accuracy is achieved. Selecting too low integration order will give spurious zero energy modes and too high integration order will increase the computation time of the matrices, eventually leading to higher solution time.
note:Script file for above analysis is available in the SCIFESOL package.