Mathematical Problems in Engineering

Volume 2017, Article ID 1467356, 14 pages

https://doi.org/10.1155/2017/1467356

## A Smoothed Finite Element-Based Elasticity Model for Soft Bodies

^{1}College of Computer Science and Technology, Beijing Normal University, Beijing 100875, China^{2}Engineering Research Center of Virtual Reality Applications, Ministry of Education of the People’s Republic of China, Beijing 100875, China

Correspondence should be addressed to Mingquan Zhou; nc.ude.unb@uohzqm

Received 8 December 2016; Revised 9 February 2017; Accepted 14 February 2017; Published 29 March 2017

Academic Editor: Giovanni Garcea

Copyright © 2017 Juan Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

One of the major challenges in mesh-based deformation simulation in computer graphics is to deal with mesh distortion. In this paper, we present a novel mesh-insensitive and softer method for simulating deformable solid bodies under the assumptions of linear elastic mechanics. A face-based strain smoothing method is adopted to alleviate mesh distortion instead of the traditional spatial adaptive smoothing method. Then, we propose a way to combine the strain smoothing method and the corotational method. With this approach, the amplitude and frequency of transient displacements are slightly affected by the distorted mesh. Realistic simulation results are generated under large rotation using a linear elasticity model without adding significant complexity or computational cost to the standard corotational FEM. Meanwhile, softening effect is a by-product of our method.

#### 1. Introduction

Physically based dynamic simulation of deformation has been an active research topic in computer graphics community for more than 30 years. Since Terzopoulos and his colleagues [1] introduced this field into graphics in 1980s, a large body of works have been published in pursuit of visually and physically realistic animation of deformable objects. According to the method for discretization of partial differential equations (PDEs) governing dynamic elasticity deformation, physically based methods are divided into mesh-based and mesh-free methods. Our method falls mainly within the first category by using the finite element method (FEM).

The FEM is one of the most widely used mesh-based methods in solving PDEs. It discretizes the objects into a mesh of finite connected nodes and approximates the field within an element by interpolation of the values at the nodes of the element. As a result, the quality of the mesh plays an important role in the FEMs. Adopting either explicit integration or implicit integration methods, badly shaped elements directly lower down the accuracy and speed of numerical solutions of governing PDEs [2]. To make it worse, the mesh is not static and is changing with the object deformation. The traditional and most promising way to improve mesh quality is the adaptive remeshing. However, remeshing methods are daunting because they involve tedious mesh topology and geometry operations and projections of field variables from the previous mesh [3, 4]. Even excellent field variable projection schemes might lead to significant errors in displacements and velocities [5]. Remeshing methods include refinement and coarsening. Both of them modify the edge length of the mesh, which disturbs the Courant-Friedrichs-Lewy (CFL) condition and introduces stability problems in turn. Our approach departs from this traditional viewpoint by smoothing the strain field on the mesh instead of smoothing the mesh. This idea is adopted by Liu and his colleagues [6] in their smoothed finite element methods (S-FEM). This new point of view avoids the problems caused by remeshing such as instability. To make the simulation fast, we combine the strain smoothing technique with the stiffness warping approach [7]. Our results show that this method is capable of producing correct simulation results using either well-shaped or ill-shaped meshes under large rotations.

In summary, our main contributions are as follows:(i)A novel smoothed pseudolinear elasticity model is presented for deformable solid simulation. It adopts a linear elasticity model for nonlinear simulation and compensates for the error of linear calculation in the nonlinear simulation using stiffness warping. Different from the previous approaches, the pseudolinear elasticity model is built on the smoothed finite element method. It is the first time that the smoothed finite element method has been used for deformable solid simulation in computer graphics.(ii)A novel smoothing domain-based stiffness warping approach is proposed to accommodate the change of integration domains in our smoothed pseudolinear elasticity model.(iii)The strain smoothing technique adopted provides an alternative way to avoid problems related to mesh distortion met in deformable solid simulation in computer graphics.

The paper is structured as follows. First, closely related researches are presented in Section 2. Then the S-FEM proposed by Liu et al. [6] is briefly reviewed in Section 3. Our smoothed and corotational elasticity model (CS-FEM for short) is presented in Section 4. Next, the experiment results and analysis are delivered in Section 5. Finally, conclusions and future studies are sketched in Section 6.

#### 2. Related Work

Several researchers have developed physically based models for deformable solid simulation since Terzopoulos and his colleagues introduced methods for simulating elastic [1] and inelastic [8] materials into the graphics community. We refer the reader to the survey papers that focus on deformation modeling in computer graphics for in-depth review [9–12]. Here, we only focus on works on the FEM.

The FEM, one of the most popular methods, has been widely used in both elastic material (summarized in [10] and later reviews) and inelastic material simulations [13, 14]. Linear elasticity FEM models are applicable to interactive application because of their stability and efficiency. However, they are not suitable for large rotational deformations because of their well-known geometric distortion. Capell et al. proposed dividing an object into small parts based on its skeleton manually [15]. Müller and his colleagues took a different approach: stiffness warping [16]. Then, they further improved the vertex-based stiffness warping to the element-based stiffness warping and extended their approach to inelastic material simulation [7]. Their method produces fast and robust simulation results and has been widely used in interactive simulation [7, 14, 16–18]. Later, Chao et al. [19], McAdams et al. [20], and Barbic [21] gave an exact corotational FEM stiffness matrix for a linear tetrahedral element by adding the higher-order terms of element rotation. Recently, Civit-Flores and Susín proposed a method to handle degenerate elements in isotropic elastic materials [22]. Our method extends the element-based stiffness warping method to the face-based stiffness warping method.

Because mesh quality is an important performance factor in mesh-based simulations, there is a huge body of literature on spatial or geometric adaptive methods. Remeshing is a widely used approach in maintaining mesh quality. Basis function refinement [23, 24], mesh embedding [25, 26], and other mixed models also work well in general. Manteaux and his colleagues gave a thorough review on them [4] recently in computer graphics. To avoid element distortion encountered in the FEM, a mesh-free method (MFM) [27] was developed which took point-based representations for both the simulation volume and the boundary surface of elastically and plastically deforming solids. Later, many useful techniques have been developed in mesh-free methods (summarized in [10]). In computational science, Liu et al. combined the strain smoothing technique [28] used in mesh-free methods [29] into the FEM and formulated a cell/element-based smoothed finite element method (S-FEM) [6]. Their method inherits the mesh distortion insensitivity of mesh-free methods and the accuracy of the FEM. After the theoretical aspects of S-FEM were clarified and its properties were confirmed by numerical experiments [30], the concept of smoothing was extended to formulate a series of smoothed FEM models such as the node-based S-FEM (NS-FEM) [31, 32], the alpha-FEM (-FEM) [33], the edge-based S-FEM (ES-FEM) [34, 35], and the face-based S-FEM (FS-FEM) [36]. As the earliest S-FEM model, the cell-based FEM converges in higher rate than the standard FEM in both displacement and energy norms. The NS-FEM can alleviate the volumetric locking problem effectively. But it is usually less computationally efficient and temporally unstable. The -FEM was proposed by Liu et al. to avoid the spurious nonzero energy modes in NS-FEM for dynamic problems. It proves to be stable and convergent, but its variational consistency depends on how it is formulated [37]. The ES-FEM-T3 often offers superconvergent and better solution than the standard FEM. It also excels in handling the volumetric and/or bending locking problem using bubble functions [38]. It is immune from the “overly soft” problem met in NS-FEM and the cell-based S-FEM. The ES-FEM was further extended for 3D problems to form the FS-FEM. Similar to ES-FEM, the FS-FEM is more accurate than the standard FEM using the same T4 mesh for dynamic problems [39] and both linear and nonlinear problems [36]. Because of its excellent properties, S-FEM has been applied to a wide range of practical mechanics problems such as fracture mechanics and fatigue behavior [40–43], nonlinear material behavior analysis [35, 38, 44–48], plates and shells [49–52], piezoelectric structures [43, 53–55], heat transfer and thermomechanical problems [56–59], vibration analysis and acoustics problems [39, 58, 60–62], and fluid and structure interaction problems [63–66]. We refer the reader to [67, 68] for recent in-depth reviews of S-FEM. As an alternative, Leonetti and Aristodemo proposed a composite mixed finite element model (CM-FEM) to generate new smoothed operators [69]. A composite triangular mesh is assumed over the domain. An element is subdivided into three triangular regions and linear stress/quadratic displacement interpolations are used to approximate the exact stress/displacement fields. Their method is also insensitive to mesh distortion and applicable to the elastic and plastic analysis in plane problems [70]. Bilotta and his colleagues adopted the same idea and proposed a composite mixed finite element model for 3D structural problems with small strain fields [71, 72].

Our method is based on FS-FEM because it is stable, accurate, and insensitive to mesh distortion in linear elasticity simulation using tetrahedral meshes. Essentially, our method uses both FS-FEM in linear elasticity modeling and stiffness warping in rotational distortion elimination. Like FS-FEM, it makes use of strain smoothing technique to avoid the problem related to element distortion instead of remeshing. Different from the original FS-FEM, stiffness warping is combined to produce 3D nonlinear deformable model using the corotational linear elasticity model. Moreover, the stiffness warping is converted from element-based stiffness warping to smoothing domain-based stiffness warping.

#### 3. The Smoothed Finite Element Method

In this section, we will present the main idea of S-FEM, which is the foundation of our model. The common concepts and key equations are also listed here.

##### 3.1. The Idea of S-FEM Models

The S-FEM, a numerical and computational method, was proposed by Liu et al. [6] in 2007 and based on FEM and some mesh-free techniques. In the standard FEM, once the displacement is properly assumed, its strain field is available using the strain-displacement relation, which is called fully compatible strain field. Then the standard FEM is formulated using the standard Galerkin formulation. The fully compatible FEM model leads to locking behavior for many problems. The assumed piecewise continuous displacement field induces a discontinuous strain field on all the element interfaces and inaccuracy in stress solutions in turn. In addition, the Jacobian matrix related to domain mapping becomes badly conditioned on distorted elements, leading to deterioration in solution accuracy. It also prefers quadrilateral elements and hexahedral elements and gives poor accuracy, especially for stresses, for easily obtained triangular and tetrahedral elements.

To avoid the upper problems of the FEM, the S-FEM uses the smoothing technique to modify the fully compatible strain field or construct a strain field without computing the compatible strain field over all the smoothing domains. Then the smoothed Galerkin weak form, instead of the Galerkin weak form used in the standard FEM, is used to establish the discrete linear algebraic system of equations. Except the strain field construction and discrete linear algebraic system establishment, both the S-FEM and the FEM follow the same steps. S-FEM models building on different smoothing domains have different features and properties. Next, we will introduce the common concepts and steps in S-FEM models. In Section 4, we will take the FS-FEM as an example to illustrate the smoothing domain building, the strain filed construction, and discrete linear algebraic system establishment of S-FEM.

##### 3.2. Smoothing Operation

We first introduce the integral representation of the approximation of function :where is a prescribed smoothing function defined in the smoothing domain of point . Smoothing domain can be moved and overlapped for different . must follow the conditions of the partition of unity, positivity, and decay. Similarly, the integral form of the gradient of function can be represented asNote that (1) and (2) are the standard forms of smoothing operation and were widely used in the smoothed particle hydrodynamics and mesh-free methods for field approximation. Applying Green’s divergence theorem on (2), we getwhere is the boundary of and is the outward normal on . Equation (3) holds if is continuous and at least piecewise differentiable.

For simplicity, a local constant smoothing function is used.. Substituting (4) into (3), we get the smoothed gradientThe integrals involving gradient of a function are recast into the boundary integrals involving only the function and the boundary normals by the gradient smoothing technique.

##### 3.3. Strain Smoothing

The S-FEM models apply the smoothing operation on a strain field and get where is the smoothed strain-displacement matrix of smoothing domain , is the nodal displacement matrix, and is the number of nodes in a smoothing domain. The smoothed strain is assumed to be constant in each smoothed domain in S-FEM.

For a linear elasticity model, the standard compatible strain-displacement matrix iswhere is the symmetric gradient operator, is the finite element shape function matrix, and is the number of nodes in an element. Note that, applying smoothing operation on , we get the smoothed strain-displacement matrix:where and where is the number of nodes in a smoothing domain . Compared to the standard FEM models, the smoothed strain of the S-FEM models only depends on the nodal displacements, the normal of the element boundary, and the shape function matrix. Usually one Gaussian point is used for line integration along each segment of boundary . From (8), we can see that the smoothed strain-displacement matrix is an average of the standard strain-displacement matrices over the smoothing domain .

##### 3.4. Smoothed Stiffness Matrix

Then the S-FEM models use the smoothed strain to compute the potential energy and define the smoothed Galerkin weak form:where is the material constant matrix, is a trial function, is a test function, is the external body force applied over the problem domain, and is the external force applied on the natural boundary.

The S-FEM models formulated using (10) are variationally consistent and spatially stable if the number of linear independent smoothing domains is sufficient and converge to the exact solution of a physically well-posed linear elasticity problem. Substituting the assumed approximations and ,into (10), we have the standard discretized algebraic system of equations:where is the nodal displacement, is the nodal displacements for all the nodes in the S-FEM model, and is the load vector:and is the global smoothed stiffness matrix defined by is a sparse matrix with nonzero entries if node and node share the same smoothing domain. For a static analysis problem, the nodal displacements can be obtained by solving (12).

#### 4. Smoothed Pseudolinear Elasticity Model

In this section, our smoothed corotational elasticity model is formulated. In S-FEM models, the FS-FEM is both stable and computationally efficient for 3D problems, compared to the NS-FEM, the cell-based S-FEM, and the ES-FEM. So our method is based on the linear FS-FEM (LFS-FEM). Under large rotations, the linear elasticity model produces geometric distortion. Stiffness warping is good at solving this kind of problem. In the following, LFS-FEM is outlined first. Then stiffness warping is abstracted and extended from an element-based method to a smoothing domain-based method. Next, the dynamics of the system are introduced to simulate the dynamic behavior of an object. In the end, the procedure of simulation is summed up.

##### 4.1. LFS-FEM

This section formulates the LFS-FEM presented by Nguyen-Thoi et al. [36] for 3D problems using tetrahedral elements. We first give the construction of face-based smoothing domain on a tetrahedral mesh and then show the equations of smoothed strain corresponding stiffness matrix definition.

*Smoothing Domain Construction*. In the FS-FEM, both the strain smoothing and the stiffness matrix integration are based on local smoothing domains. These local smoothing domains are constructed based on faces of the neighboring elements such that and , , in which is the total number of smoothing domains (i.e., faces here) in . For tetrahedral elements, a smoothing domain is created by connecting three nodes of its associated face to the central point of two adjacent tetrahedral elements as shown in Figure 1.