Abstract
We propose a fast stiffness matrix calculation technique for nonlinear finite element method (FEM). Nonlinear stiffness matrices are constructed using Green-Lagrange strains, which are derived from infinitesimal strains by adding the nonlinear terms discarded from small deformations. We implemented a linear and a nonlinear finite element method with the same material properties to examine the differences between them. We verified our nonlinear formulation with different applications and achieved considerable speedups in solving the system of equations using our nonlinear FEM compared to a state-of-the-art nonlinear FEM.
1. Introduction
Mesh deformations have widespread usage areas, such as computer games, computer animations, fluid flow, heat transfer, surgical simulations, cloth simulations, and crash test simulations. The major goal in mesh deformations is to establish a good balance between the accuracy of the simulation and the computational cost; achieving this balance depends on the application. The speed of the simulation is far more important than the accuracy in computer games. The simulation needs to be real-time in order to be used in games so free-form deformation or fast linear FEM solvers can be used. However, high computation cost gives much more accurate results when we are working with life-critical applications such as car crash tests, surgical simulators, and concrete analysis of buildings; even linear FEM solvers are not adequate enough for these types of applications in terms of the accuracy.
For realistic and highly accurate deformations, one can use the finite element method (FEM), a numerical technique to find approximate solutions to engineering and mathematical physics problems. FEM could be used to solve problems in areas such as structural analysis, heat transfer, fluid flow, mass transport, and electromagnetics [1, 2].
We propose a fast stiffness matrix calculation technique for nonlinear FEM. We derive nonlinear stiffness matrices using Green-Lagrange strains, themselves derived from infinitesimal strains by adding the nonlinear terms discarded from infinitesimal strain theory.
We mainly focus on the construction of the stiffness matrices because change in material parameters and change in boundary conditions can be directly represented and applied without choosing a proper FEM [3]. Joldes et al. [4] and Taylor et al. [5] achieve real-time computations of soft tissue deformations for nonlinear FEM using GPUs; however, they do not describe how they compute stiffness matrices; thus, we cannot implement their methods and compare them with our proposed method. Cerrolaza and Osorio describe a simple and efficient method to reduce the integration time of nonlinear FEM for dynamic problems using hexahedral 8-noded finite elements [6]. We compare our stiffness matrix calculations with Pedersen’s method [3] to measure performance and verify correctness. We achieve a speedup in calculating the stiffness matrices and a speedup in solving the whole system on average, compared to Pedersen’s method.
2. The Nonlinear FEM with Green-Lagrange Strains
We use tetrahedral elements for modeling meshes in the experiments. Overall, there are 12 unknown nodal displacements in a tetrahedral element. They are given by [2] In global coordinates, we represent displacements by linear function by
For all 4 vertices, (2) is extended as Constants can be found as where is given by is , where is the volume of the tetrahedron. If we substitute (4) into (2), we obtain , , , , and the volume are calculated by Because of the differentials in strain calculation, is not used in the following stages. If we expand (6), we obtain For tetrahedral elements, to express displacements in simpler form, shape functions are introduced (). They are given by
In our method, nonlinear stiffness matrices are derived using Green-Lagrange strains (large deformations), which themselves are derived directly from infinitesimal strains (small deformations), by adding the nonlinear terms discarded in infinitesimal strain theory. The proposed nonlinear FEM uses the linear FEM framework but it does not require the explicit use of weight functions and differential equations. Hence, numerical integration is not needed for the solution of the proposed nonlinear FEM. Instead of using weight functions and integrals, we use displacement gradients and strains to make the elemental stiffness matrices space-independent in order to discard the integral. We extend the linear FEM to the nonlinear FEM by extending the linear strains to the Green-Lagrange strains.
We constructed our linear FEM by extending Logan’s 2D linear FEM to 3D [2]. To understand Green-Lagrange strains, we must first see how they differ from the infinitesimal strains used to calculate the global stiffness matrices in a linear FEM. Figure 1 shows a 2D element before and after deformation, where the element edge with initial length becomes . The engineering normal strain is calculated as the change in the length of the line The final length of the elemental edge can be calculated using By neglecting the higher-order terms in (11), 2D infinitesimal strains are defined by
By the definition, the nonlinear FEM differs from the linear FEM because of the nonlinearity that arises from the higher-order term neglected in calculation of strains. The strain vector used in the linear FEM relies on the assumption that the displacements at the -axis, -axis, and -axis are very small. The initial and final positions of a given particle are practically the same; thus, the higher-order terms are neglected [7]. When the displacements are large, however, this is no longer the case and one must distinguish between the initial and final coordinates of the particles; thus the higher-order terms are added into the strain equations. By adding these high-order terms, 3D strains are defined as [8] which leads to The Green-Lagrange strain tensor is represented in matrix notation as where is the nodal displacement, [] is the linear, and is the nonlinear part of the matrix [3]. For a specific element, [] and are constant, as with the matrix in the linear FEM. With the variation of [9], (15) becomes
Gathering the strain components together, we can rewrite (15) and (16) as
The linear part of the matrix () is the same as the matrix in the linear FEM. Calculating becomes more complex with the introduction of the nonlinear terms. After finding the nonlinear strains, these equations are combined with the shape functions to find matrix : The most frequently used terms, which are the nine displacement gradients for calculating the nonlinear strains, are , , , , , , , , and . Using (8) for displacements, we can construct the displacement gradients using the partial derivatives of the shape functions. They are represented by where represents .
We can evaluate the partial derivatives of the shape functions as follows (for the 1st node of ): Using (17) and (20), we obtain for the 1st node (21). Similarly, using (17) and (20), we obtain for the 1st node (22).
Consider
The FEM is derived from conservation of the potential energy, which is defined by where is the strain energy of the linear element and is the work potential. They are given by where the engineering strain vector is From (24), the engineering stress vector is related to the strain vector by The secant relations are described by the matrix . We substitute (15) and (26) into (24), obtaining the element stiffness matrix We can discard the integrals as we did for the linear FEM. , , and are constant for the four-node tetrahedral element, so (27) is rewritten as The secant stiffness matrix which is is nonsymmetric because of the fact that .
Introducing nodal forces, we obtain With the equilibrium equation and cancelling , the whole system for one element reduces to By substituting with , we obtain Finally, only nonlinear displacement functions remain, which are solved with Newton-Raphson to find the unknown displacements [10].
Element residuals are necessary for the iterative Newton-Raphson method. The element residual is a vector for a specific element. The residual for a specific element is defined as Having determined , we can now express (32) in expanded vector form as
The tangent stiffness matrix is also necessary for the iterative Newton-Raphson method. The tangent stiffness matrix is also matrix, like the elemental stiffness matrix. However, the tangent stiffness matrix depends on residuals, unlike the elemental stiffness matrix. Elemental stiffness matrices are used to construct residuals and the derivatives of the residuals are used to construct the elemental tangent stiffness matrices. We can express the elemental tangent stiffness matrix for a specific element as
Newton-Raphson method is a fast and popular numerical method for solving nonlinear equations [10], as compared to the other methods, such as direct iteration. In principle, the method works by applying the following two steps (cf. Algorithm 1): (i) check if the equilibrium is reached within the desired accuracy; (ii) if not, make a suitable adjustment to the state of the deformation [11]. An initial guess for displacements is needed to start the iterations. The displacements are updated according to
| ||||||||||||
In the proposed nonlinear FEM, is the vector that keeps the information of the nodal displacements. Instead of making only one assumption, we make whole vector initial guess in order to start the iteration.
Consider where is residual of the global stiffness matrix calculated in (33) and is the tangent stiffness matrix calculated in (34).
At every step, the vector and the matrix are updated for every element with the new values. Then, and are assembled as we did with for the global stiffness matrix and the global force vector in linear FEM. Boundary conditions are applied to the global vector and the global matrix. Using the global vector and the global matrix, we have is updated with the solution of (37). Consider Then, we check if the equilibrium is reached within the desired accuracy defined by as After the desired accuracy is reached, the unknown nodal displacements are found.
3. Experimental Results
The proposed nonlinear FEM and Pedersen’s nonlinear FEM were implemented using MATLAB programming language. The visualizer was implemented with C++ language and connected to the solver using the MATLAB engine [12], which allows users to call the MATLAB solver from C/C++ or Fortran programs. The simulation results, interaction with the 3D model, and the 3D models themselves were visualized using OpenGL, and the nose experiment was visualized using 3ds Max [13]. To speed up the nonlinear FEM, we used MAPLE’s symbolic solver [14], which is integrated into MATLAB. We conducted all the experiments on a desktop computer with a Core i7 3930 K processor overclocked at 4.2 GHz with 32 GB of RAM. We used linear material properties for the models in the experiments. We used 1 GPa for Young’s modulus because polypropylene has Young’s modulus between and 2 GPa and polyethylene HDPE has Young’s modulus 0.8 GPa, which shows plastic properties and they are close to 1 GPa. Because most steels and plastic materials undergo plastic deformation near the value of , we used for Poisson’s ratio (). Our simulation is static so we used a single load step in all experiments. Multiple load steps are used when the load forces are time-dependent or the simulation is dynamic [15].
We conducted four experiments, each having different number of elements to observe the speedup for both stiffness matrix calculation and for the solution of the system. As expected, the proposed method and Pedersen’s method produced same amount of nodal displacements in all experiments.
The first experiment was conducted for a cube mesh with eight nodes and six tetrahedral elements. Figure 2 shows that the cube is constrained at the upper four nodes and pulled downwards with a small amount of force (one unit force for each of the upper four nodes). This experiment was conducted with a small mesh in order to carefully examine the nodal displacements and strains for each element. Figure 3 shows the initial and final positions of the nodes for the linear and nonlinear FEMs, respectively. As seen in Figure 3, the linear and nonlinear methods produce similar displacements when the force magnitude is small. Tables 1 and 2 show the displacements at force applied nodes (first, second, third, and fourth) using the linear and nonlinear FEMs, respectively. Figure 4 shows that displacement increases linearly with force magnitude. However, as expected, the nonlinear FEM behaves quadratically due to the nonlinear strain definitions. Figure 5 depicts the convergence of the Newton-Raphson method for the nonlinear FEM.
(a)
(b)
The second experiment was conducted on a beam with 90 nodes and 216 tetrahedral elements. Figures 6(a) and 6(b) show that the beam is constrained at the blue nodes and twisted at both ends. Figure 7 shows the final shape of the beam mesh for both the proposed and Pedersen’s methods. Tables 3 and 4 show the displacements at force applied nodes (green nodes) for the second experiment using the linear and nonlinear FEMs, respectively.
(a)
(b)
(a)
(b)
(c)
(d)
The third experiment was conducted with a cross mesh of 159 nodes and 244 tetrahedral elements. We aimed to observe if there is a root jump occurring when solving the system for a high amount of force (50 N units), and its effect on the computation times for both methods. Figure 8 shows that the cross shape is constrained at the blue nodes and pushed towards the green nodes. Figure 9 shows the final shape of the beam mesh for both the proposed and Pedersen’s methods. Tables 5 and 6 show the displacements at force applied nodes (green nodes) using the linear and nonlinear FEMs, respectively.
(a)
(b)
We conducted fourth experiment with a liver mesh of 465 nodes and 1560 tetrahedral elements. Figure 10 shows that the mesh is constrained at the blue nodes and pulled from the green node (30 N units) in the direction of the arrow. We aimed to observe the similar amount of speedup like previous experiments for a high density mesh. Figure 11 shows the final shape of the beam mesh for both the proposed and Pedersen’s methods. Tables 7 and 8 show the displacements at force applied node (node number 271) using the linear and nonlinear FEMs, respectively.
(a)
(b)
(a)
(b)
Computation times of the finite element experiments are required to compare how much faster our proposed method is than Pedersen’s. When comparing nonlinear FEMs, we calculated the computation times to construct the stiffness matrices as well as the computation times of the nonlinear FEM solutions to determine how different calculations affect them. Table 9 depicts the computation times for the stiffness matrix calculation and Table 10 depicts the computation times for the system solution. Table 11 shows the iteration counts to solve the system using the Newton-Raphson procedure.
The speed-up columns of Tables 9 and 10 depict the speedups of the proposed method compared to Pedersen’s method for the stiffness matrix calculation and the system solution using a single thread, respectively. The speedup is calculated as follows:
Our proposed method outperforms Pedersen’s method. On the average, it is faster at computing stiffness matrices because Pedersen’s method uses more symbolic terms. However, both methods use Newton-Raphson to solve nonlinear equations, which takes approximately of the computation time. Thus, the overall speedup decreases to on average. In experiment 3, because of more iterations due to root jumps, there was a performance loss against Pedersen’s method.
We also applied our method for corrective operation on the misshapen nose of a head mesh. The head mesh is composed of 6709 nodes and 25722 tetrahedral elements (see Figure 12(a)); all the operations were performed in the nose area of only 1458 tetrahedral elements. Figure 12(b) shows that the head mesh is constrained at the blue nodes and pushed upwards at the green nodes and Figure 12(c) shows the result of the nonlinear FEM.
(a)
(b)
(c)
4. Conclusions and Future Work
We propose a new stiffness matrix calculation method for nonlinear FEM that is easier to analyze in terms of constructing elemental stiffness matrices and is faster than Pedersen’s method. The proposed method is approximately times faster, on average, at computing stiffness matrices and faster at computing the whole system than Pedersen’s method.
Although the proposed nonlinear FEM has significant advantages over Pedersen’s nonlinear FEM, there is still room for the following development.(1)Heuristics could be applied to avoid root jumps.(2)Although we decreased system memory usage by simplifying the solution process for the nonlinear FEM, a significant amount of system memory is still used. The solution process could be further optimized to decrease memory usage.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.