Strojniški vestnik - Journal of Mechanical Engineering 64(2018)5, 303-309 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5081 Original Scientific Paper Received for review: 2017-11-13 Received revised form: 2018-02-25 Accepted for publication: 2018-03-19 Consideration of Body Forces within Finite Element Analysis Christian Glenk* - Florian Hüter - Daniel Billenstein - Frank Rieg University of Bayreuth, Faculty of Engineering Science, Germany The finite element analysis (FEA) is an established numerical method for the mechanical analysis of structural components, which cannot be analytically described sufficiently well. Despite ever-increasing computing power, the efficiency enhancement of the calculation methods and their computational implementation is one of the highest development goals in this field. This article shows an efficient method for the calculation of weight forces that is based on a factorization of the volume integral so that the integration is independent of the global node coordinates of the finite element mesh and solely dependent on the shape functions of the element type. The numerical integration can thus be done in advance, and the resultant numerical values can be tabulated in the finite element software. No quadrature needs to be performed during the program runtime. This approach leads to a reduction of the computing effort and time. Keywords: finite element analysis, body force, weight force, numerical integration, Gaussian quadrature Highlights • Efficient method for the calculation of body forces within finite element analysis was implemented. • Factorization of the volume integral makes the integration independent of the global node coordinates. • Most of the computational effort, i.e. numerical integration, does not need to be performed during the FE-program runtime. • Computing effort and computation time can be reduced. • The introduced method enables exact numerical integration without affecting computation time. 0 INTRODUCTION The swift development in recent decades towards more powerful computers has allowed the widespread use of complex simulation programs based on numerical methods. Thus, finite element analysis (FEA) is now an established standard procedure in many areas of numerical strength analysis. Despite the increasing computing power, the endeavour is to develop more efficient algorithms and thus to minimize the computation time during the FEA. The main idea of the FEA is to partition a continuum into a finite number of finitely large subcontinua [1] to [3]. Such discretized elements are referred to as finite elements (FE). Their nodes act as joints between neighbouring elements and as target points for the specification of loads and displacement boundary conditions. The solid deforms under mechanical load, resulting in a deformation of each finite element. The displacement of any point x within an element is specified by the vector field u . The location and displacement of any point inside a finite element are interpolated from the position p and the displacement di of the element nodes. These nodal displacements are computed via a discretized formulation of the mechanical momentum equation. Starting from the geometrical discretization of the continuum, one possible approach for the derivation of the discretized momentum equation is the Galerkin method [1], [2] and [4]. In the case of a static equilibrium problem, this ultimately leads to a system of linear equations for determining the nodal displacements: K • d = Rs + Rb. (1) Matrix K is the stiffness matrix and the vectors RS and Rb represent the load vectors of surface forces and body forces. Body forces, such as the weight force, affect almost all technical components and cannot be neglected in stress analysis, e.g. in areas such as structural engineering. They affect the whole element from a physical point of view, but they are proportionately distributed to the element nodes within the finite element method. For each element, the equivalent nodal forces can be calculated as follows [5]: R =J]J N'pgdV, (2) Ve where N is the matrix of the shape functions and pg describes the mass density and the gravitational constant. If the mass density and the gravitational factor are assumed to be constant for each finite element, both quantities can be put outside the integral, leading to the following expression: Rb =j]J NT dV -ps. (3) Ve The vector Rb in Eq. (1) can be calculated by assembling each elemental force vector into one global force vector [2]: *Corr. Author's Address: University of Bayreuth, Faculty of Engineering Science, Department of Engineering Design and CAD, Universitätsstrasse 30, 95444 Bayreuth, Germany, christian.glenk@uni-bayreuth.de e In this paper, we present a mathematical approach to a time-efficient calculation of the volume integral in Eq. (3). 1 ELEMENT DESCRIPTION AND NUMERICAL INTEGRATION Interpolation polynomials are used to describe the displacement field u within a finite element, which are composed of the shape functions N and the displacement vectors d of each element node. Typically, the following formula is used in the FEM literature [1], [5] and [6]: u(r,s,t) N(r,s,t) • di. (5) The location dependency of the displacement field is included in the shape functions. If these functions were dependent on the global Cartesian coordinates (x,y, z), the shape functions N would have to be determined individually for each finite element of the structure. To avoid this effort, a coordinate transformation into a generally curvilinear element coordinate system (r, 5, t) is typically performed within the finite element method. In that way, the finite elements are transformed into their undistorted geometry, which is why this coordinate system is also called the natural coordinate system. An interpolation approach similar to Eq. (5) is typically used to establish a relation between the global coordinates (x, y, z) and the curvilinear, natural coordinates (r, 5, t): x(r, s, t) M (r,s, t)P , (6) where Ni denotes the shape functions for the interpolation of the geometry and p the node coordinates in the global Cartesian coordinate system. Generally, the supporting points for the interpolation polynomials of u and do not have to be coincidental [6]. Other or fewer nodes can be used for the calculation of x . If the supporting points are chosen to be identical, the shape functions NNi and Ni coincide. In this case, the element formulation is called isoparametric [2]. According to Eq. (3), a volume integral must be solved for the calculation of the gravitational nodal forces. If the element formulation is assumed to be isoparametric, the integration can be performed in the natural coordinate system [3]: R = Hi NTdetJ drdsdt ■ pg. (7) (r ,s,t) The term det J denotes the determinant of the Jacobian matrix, also called functional determinant, and is defined as: dx dx dx dr ds ~dt dy dy dy dr ds dt dz dz dz dr ds dt The evaluation of the volume integral, shown in Eq. (7), is usually carried out by numerical integration [3]. Typically, Gaussian quadrature is used in this context [1], which is based on the idea of evaluating the integrand at certain grid points, so-called Gauss points, and calculating the weighted sum [3]: R Z Z a«aa &etj}{ra, sb, tc) -pi. (9) a=1 b=1 c=1 The parameters (aa, ab, ac) denote the Gaussian weights, and (ra, sb, tc) are the natural coordinates of the Gauss points [3]. These parameters might be different for various element geometries. For example, there are Gaussian weights and Gauss points specially developed for triangles and tetrahedrons [1], [7] and [8]. If the polynomial degree is denoted as p, at least m = (p+1) / 2 supporting points are required to obtain an exact integration [9]. 2 FACTORIZATION OF THE JACOBIAN DETERMINANT During the calculation of the Jacobian determinant, the partial derivatives of the Cartesian coordinates (x, y, z) with respect to the natural coordinates (r, 5, t) must be evaluated. If an isoparametric element formulation is assumed, the following equation holds [3]: y — m dr ' y — Tf ds ' y - Tf dt ' det J = yy-N-y- jl dr yd— jf ds yd— jf dt . (10) t=i dr k ±Nzk ds k ydN^z h dtZk This determinant can be expanded to the following expression: det J = ¿ ¿ ¿ N N ^ i=i j=i k=1 dr ds dt x¡yjzk +ttt i=1 j=1 k=1 +t t t i=1 j=1 k=1 -t It It i=1 j=1 k=1 -t t t DN, dN, dNj r s t N j dNk dN, r s t Nk dNj dN, r s t dN 3Nt dNj r s t dNj dN dNk -ttt i=1 j=1 ,=1 dr ds dt xiy,zk x¡yjzk x¡yjzk xJjzk x¡yjzk (il) By bracketing out the node coordinates, this equation can be rewritten as: det J = XZZ i=i j=i k=i dNi dN, dN dr ds dt dNj dNj dNj dr ds ~dt~ N N dNk dr ds dt • xiyjzk. (12) According to this equation, the Jacobian determinant can be "factorized" in such a way that the node coordinates (xi,y, zk) and the partial derivatives of the shape functions are separated: det J = ^rjrj^Hijk(r,s,t) ■ Xijk. (13) i=1 j=1 k=1 Based on the general properties of determinants [10], it can be shown that the determinant of the partial differential derivatives Hijk possess the following properties: 1) Hijk=0 for i = j and/or j = k and/or k = i, 2) Hijk = -Hjik = Hjki = -Hkji = Hkij = -Hikj■ These properties result from the fact that the rows of Hijk are identical and that interchanging rows leads to a change of the sign of the determinant. According to these properties, the following statements can be concluded: 1. To determine the value of all Hjk, it suffices to calculate solely those j with: i = 1, ..., n-2, j = i + 1, ..., n - 1, k= j + 1, ..., n■ 2. The sign of all Hijk can be determined via property 2). 3. The number of the determinants Hijk to be calculated decreases from n3 to To proof these statements, Hijk is considered as a tensor of third order with n elements per row, i.e. it describes a cube of "edge length" n. According to property 1), all components of this tensor with i = j and/or j = k and/or k = i are equal to zero. All these zero-values Hiik, Hj and Hkjjk are arranged in such a manner that they form three planes within the cube, which intersect along the body diagonal of the cube, where i = j=k holds. The cube is thereby divided into six tetrahedrons, Fig. 1. These tetrahedrons can be described by the intervals of indices a = 1, ..., n - 2, b = a + 1, ..., n - 1 and c = b + 1, ..., n, where the parameters (a, b, c) are placeholders for each of the six possible permutations of (i,j, k). One possible tetrahedron is defined as i = 1, ..., n - 2, j = i + 1, ..., n - 1, and k =j + 1, ..., n, which is shown in Fig. 1. plane diagonal ►k body diagonal relevant area of the tensor Fig. 1. Composition of the tensor Hyk According to property 2) all other five tetrahedrons have the same entries regarding their absolute value. They only differ in sign. Therefore, it suffices to calculate only one of these six tetrahedrons. All the others can be derived from the first tetrahedron based on property 2), which is why the second statement is valid. The third statement is verified with an example. If n = 3, i.e. i,j,k= 1, ..., 3, then there are n3 = 27 different "combinations" for the index triple: (i,j, k) = (1, 1, 1), (1, 1, 2), (1, 2, 1), (2, 1, 1) etc. According to property 1), all combinations comprised of at least two identical indices, such as (1, 1, 2) or (3, 3, 3), are of no interest. Consequently, the number of relevant combinations decreases to n! = 6. Furthermore, the sequence of the indices within each index triple is irrelevant according to property 2), i.e. | Hijk\ = | Hjik\ = | Hkij \, etc. The number of relevant (n - 3)! 3!' Consideration of Body Forces within Finite Element Analysis index combinations is thereby reduced to only n n! , -.x, = I- F°r this reason, (1,2, 3) is the only (n - 3)! 3! interesting combination in this example. If the value of H123 is calculated, all other 26 determinants Hijk can be derived via properties 1) and 2). Therefore, the n! calculation of n3 - --„. = 26 integrals can be (n - 3)! 3! omitted, resulting in a reduction of computing effort n! 1 to -= —. n3(n - 3)! 3! 27 In this example, the combination (1,2,3) was arbitrarily chosen as the "only" interesting combination. Likewise, any other permutation could be used. The relevant point is that calculating one of them is sufficient. This example shows that the computing effort for the evaluation of the Jacobian determinant (Eq. (13)) can be significantly reduced by exploiting the properties of Hijk. 3 FACTORIZATION OF THE VOLUME INTEGRAL The factorization of the Jacobian determinant according to Eq. (13) can be utilized for an efficient evaluation of the volume integral, Eq. (7). The basic idea is to make the integral independent of the global coordinates of the element nodes by exploiting the factorization of the Jacobian determinant. Inserting Eq. (13) into Eq. (7) leads to the following expression: R = JIf NTZZ 1^HijkXijkdrdsdt • pg. (14) (r,s,t) =1 j=1 k=1 This equation can be rewritten as follows: R=ZEE HI nth^ i=1 j=1 k=1 (r ,s,t) drdsdt ■ pgXi ijk =ZZZ I i =1 j =1 k=1 X ijk ijk (15) As a result, the volume integral in Eq. (14) has been split up into a sum of two factors. The first factor \ijk consists of an integral, which is independent of the node coordinates (xi,yi, z). The node coordinates are included in the second factor Xjk. This factorization of the original volume integral leads to a significant increase of the total computing effort because several integrals need to be calculated instead of one. Furthermore, since several integrals must be calculated and summed up, a full integration order should be used to avoid summing up numerical errors, which makes the numerical integration more expensive. This effect is particularly noticeable for elements with a high number of nodes, such as quadratic hexahedrons, for which many integrals must be calculated. However, the advantage of this approach is that these integrals are independent of the node coordinates, the mass density, and the gravitational factor. They only depend on the shape functions and their partial derivatives. Therefore, these integrals are specific for an element type just like the shape functions. They only need to be calculated and tabulated once for a specific element type. Afterwards, they can be used to calculate the nodal forces for all finite elements of the same element type. In practice, the integrals \ijk are calculated only once and their numerical values are stored hard coded within the finite element program. During the program runtime these numerical values of \ijk only need to be multiplied with the vectors Xjk and summed up for each element according to Eq. (15). No Gaussian quadrature must be performed. Therefore, both the computing effort and the computing time can be reduced in practice. 4 COMPARISON OF METHODS In this section, the efficiency of the introduced factorization-based method is considered. For this purpose, the calculation of the weight force of a tetrahedron element with quadratic shape functions is performed, see Fig. 2. The required shape functions and their partial derivatives can be taken from Rieg et al. [3]. [1] (0.0|0.0|0.0) [3] (0.0|1.0|0.0) [5] (0.5|0.0|0.0) [7] (0.0|0.5|0.0) [9] (0.0|0.5|0.5) [2] (1.0|0.0|0.0) [4] (0.0|0.0|1.0) [6] (0.5|0.5|0.0) [8] (0.5|0.0|0.5) [10] (0.0|0.0|0.5) Fig. 2. Quadratic tetrahedron element In this purely academic example, the mass density is assumed to be p = 1 and the gravitational factor is set to g = (1,0,0). First, the volume integral (Eq. (3)) is solved via the standard approach using the Gaussian quadrature according to Eq. (9). The total order of the polynomial inside the integrand of the volume integral is p = 5 because quadratic shape functions are used. An integration order of at least m = 3 must be chosen to ensure an exact numerical integration. In this case, the numerical integration leads to the following value for the weight force on the tetrahedron rounded to six decimal places: ||Re|| = 0.166667. The efficiency of this standard approach is compared to the factorization-based method according to Eq. (15). As mentioned in Section 3, the integrals \ijk are not calculated during the finite element analysis, which is why no numerical integration must be performed. These integrals are solved in advance, and their numerical values are stored hard coded within the program. Only the summation needs to be performed during runtime, which finally leads to the same result as above: ^Reb || = 0.166667. This result can be easily verified by an analytical calculation of the weight force on the tetrahedron: 1 F = pgv = - 0.166667. (16) The results of all three methods match exactly on the full mantissa length. A test program has been written to compare the performance of both integration methods. In Table 1, a comparison of the computation times of the standard integration approach tGQ using Gaussian quadrature and the factorization-based method tF is given. The computation time for the calculation of the integrals Ijk is not taken in to account, because the numerical integration is not performed during runtime. Different integration orders have been applied during Gaussian quadrature. The computing time of the factorization-based method is lower in each case for the test setting considered here. Especially for higher integration orders, the advantage is significant according to Table 1. Table 1. Comparison of computation times for different integration orders Integration order m Ratio of computation times tGQ / tF 2 1.7 3 5.8 4 14.1 5 27.2 This observation can be explained by the fact that a higher integration order implies a significantly higher number of Gaussian integration points, which must be considered during numerical integration. This affects the computation time tGQ. In contrast, the computation time of the factorization-based method tF does not change, because no numerical integration is performed during runtime. The factorization-based method can be applied to every continuum element. Fig. 3 shows a cantilever beam under dead load. lg z y —>y axa L v. s g = 9. 81 • 103 N/t p = 7.85 • 10"9 t/mm3 E = 206 • 103 N/mm2 u = 0.3 L — 2 • 103 mm a — 100 mm Fig. 3. Cantilever beam under dead load The aim of this example is to demonstrate the applicability of this method to different finite element types. Both linear and quadratic hexahedrons and tetrahedrons are considered. The results of the finite element simulations are summarized in Table 2. Each finite element model considered the weight force correctly compared to the analytical solution: F = pgLa2 = 1540.17 N. (17) The displacement field of the quadratic hexahedron model is shown in Fig. 4. The values of the maximum total displacement are summarized in Table 2. Table 2. Results of the maximum displacement de-pending on the element type and number of elements Element type Number of elements 2Â [N] Maximum displacement [mm] Linear hexahedron 20000 1540.17 0.891 Quadratic hexahedron 20000 1540.17 0.896 Linear tetrahedron 86011 1540.17 0.871 Quadratic tetrahedron 86011 1540.17 0.896 The quality of the calculated maximum total displacements depends on the element type and on the number of finite elements. An analytical calculation based on the Bernoulli beam theory [11] leads to a value of 0.897 mm for the displacement of the end of the cantilever beam. The comparison illustrates that the factorization-based method is applicable to large finite element structures of different finite element types. m + 0 00E+000 + 8 14E-002 +8 14E-002 + 1 ¿3E-001 + 1 Ó3E-001 - +2 44E-001 +2 44E-001 - +3 26E-001 + 3 2dE-001 - + 4 07E-001 + 4 07E-001 - + 4 89E-001 +4 89E-001 - + 5 70E-001 +5 70E-001 +6 51E-001 +6 51E-001 - + 7 33E-001 ■ +7 33E-001 + 8 14E-001 1 +8 14E-001 + 8 96E-001 Fig. 4. Total displacement of the cantilever beam 5 DISCUSSION Using the factorization-based method for the calculation of weight forces within the finite element analysis can lead to a noticeable reduction of the computation time, as shown in Section 3. However, this approach is not limited to the calculation of weight forces but can be generalized as follows: jjj f (r,s,t)dV = £