Bubble-Enriched Smoothed Finite Element Methods for Nearly-Incompressible Solids

This work presents a locking-free smoothed finite element method (S-FEM) for the simulation of soft matter modelled by the equations of quasi-incompressible hyperelasticity. The proposed method overcomes well-known issues of standard finite element methods (FEM) in the incompressible limit: the over-estimation of stiffness and sensitivity to severely distorted meshes. The concepts of cell-based, edge-based and node-based S-FEMs are extended in this paper to three-dimensions. Additionally, a cubic bubble function is utilized to improve accuracy and stability. For the bubble function, an additional displacement degree of freedom is added at the centroid of the element. Several numerical studies are performed demonstrating the stability and validity of the proposed approach. The obtained results are compared with standard FEM and with analytical solutions to show the effectiveness of the method. three times and it also shows around 26% less error in strain energy compared to ES-FEM.


Introduction
The deformation of soft matter, which can be modelled using the equations of hyperelasticity, is an important problem with applications throughout the physical sciences and engineering. Such soft matter, e.g., rubbers, polymers and human soft tissues, behave in a geometrically and materially non-linear manner characterized by large deformations and nearly incompressible behaviour. It is important to design numerical approaches that provide robust and accurate simulation for scientists and engineers [1]. The first step to numerically solve such problem is to obtain a subdivision of the problem domain as automatically (without user intervention) and rapidly as possible. In the context of finite element methods, this operation is called mesh generation. Meshes consisting of low-order simplex tetrahedral cells (triangulations) are best suited to automatic mesh generation and there exist schemes that allow the generation of millions of tetrahedral elements per second [2]. However, finite elements defined on simplex tetrahedral cells are also notorious for the excessive stiffness and the sensitivity to volumetric locking. Rabczuk et al. [3] discussed a possibility of non-mesh based method to handle incompressibility undergoing large deformations. However, this work focuses on development of an alternative approach which leads to use simplex tetrahedral cells. Therefore, the strain-smoothed approach in the framework of finite element approximation [4][5][6] (so-called smoothed finite element methods-S-FEM) is an another possible option to address the following issues in the context of simulating soft matter: 1) volumetric locking associated with nearly-and fully-incompressible materials and 2) the difficulty of generating quality subject-specific meshes. The method relies on the smoothed deformation gradient, which therein is assumed homogeneous over smoothing domains.
Various 3D strain smoothing approach variants have been introduced. Firstly, Nguyen-Thoi et al. [7] used the face-based smoothed finite element method (FS-FEM) for linear and geometrically non-linear problems. Viscoelastoplastic problems are successfully addressed by the face-based strain smoothing [8] and 3D extension of node-based smoothed finite element method (NS-FEM) [9], both of which show higher efficiency than FEM and the edge-based smoothing method. Cazes et al. [10] presented 3D-extended edge-based smoothed finite element method (ES-FEM) for solid problems. Not only solid problems but also coupled structural-acoustic problems are tackled by smoothed finite element methods. He et al. [11] used a coupled 3D ES-FEM/FEM formulation, with ES-FEM and FEM applied to structure and fluid parts, respectively, while Li et al. [12] used a coupled ES/FS-FEM approach for solid and fluid. Moreover, to handle heat transfer and thermal mechanical problems, Feng et al. [13] used face-based S-FEM (FS-FEM) and Cui et al. [14] used cell-based S-FEM (CS-FEM). In particular, Cui and colleagues employed cell-based strain smoothing approach with the meshless radial point interpolation method to obtain the radial basis function on smoothing domains. Other hybrids of node-and edge-based smoothed finite element methods with meshless point interpolation methods for heat transfer problems are presented by Wu et al. [15,16]. Natarajan et al. [17] constructed cell-based polyhedral S-FEM with one sub-cell. To stabilize the polyhedral S-FEM, they adapted the stabilizing term of the virtual element method to S-FEM without increasing the number of smoothing domains. Additionally, further applications for linear fracture problems were investigated by 3D ES-FEM by Zeng et al. [18], employing singular 7 node tetrahedral elements for the crack front.
The application of main interest of this work, i.e., soft materials, has been also considered by Jiang et al. [19,20]. They employed the selective 3D-ES/NS-FEM and FS/NS-FEM to solve nearly-incompressible hyperelastic models. To use the selective S-FEM, the strain energy density is split into two parts: the deviatoric and volumetric. Since it is known that ES-FEM is immune to the distorted meshes and NS-FEM can avoid the over-stiffness phenomenon [21], the deviatoric part of the strain energy density is solved by 3D-ES-and FS-FEM while NS-FEM tackles the volumetric part.
Instead, in the present work, we employ a cubic bubble function to enrich S-FEM [22,23]. The study of Ong et al. [23] showed that the inf-sup condition for the bubble-enriched S-FEM was satisfied in small deformation cases and Lee et al. [22] used the concept of the bubbleenhanced S-FEM to analyze large deformation problems in the quasi-incompressible hyperelastic limits. In this paper, the suggested bubble function is integrated into two strain-smoothed approaches, 3D-edge-based and face-based smoothing methods to deal with the distorted meshes in near-incompressibility.
An outline of this paper as follows: in the following section, the basic properties of cell-based, edge-based, face-based and node-based gradient smoothing approaches extended to three-dimension in the framework of finite elasticity are recalled. In the next section the bubbleenhanced smoothed FEM is formulated for the hyperelastic neo-Hookean material. In the penultimate section numerical tests including large deformation and insensitivity to mesh distortion in incompressible models are studied with several cases of soft materials parameters, e.g., polymers and human organs. In the final section, we draw a conclusion to remarks of the proposed methods and give some possible ideas for future work.

Review of Gradient-Smoothed Finite Element Approximation
Many existing studies of strain-smoothed finite element approaches describe their basic formulation and properties. We briefly revisit the background of the strain-smoothed method in this section. The basic idea behind S-FEM is to divide the computational domains into subdomains wherein gradients are smoothed. Gradients are constant over the smoothing domains but are discontinuous across the boundaries of smoothing domains. Sub-domains are usually constructed using the topology provided by the mesh. Depending on this construction, various solution behaviors are observed, offering the use of a spectrum of methods with particular properties [5,24].
The salient abilities of S-FEM are its insensitive to locking and to element distortion [22]. The integration for S-FEM is performed in the global space, hence the Jacobian transformation from the local to global space is obviated. Moreover, since the bandwidth of the stiffness in the smoothed-strain method is wider than in the existing FEM, there is a potential to design methods without over-estimated stiffness.
Four different types of S-FEM have been introduced, depending on how the smoothing domains are built [21]: cell-based (CS-FEM), edge-based (ES-FEM), node-based (NS-FEM) and face-based (FS-FEM). In this section, we extend the idea of two-dimensional cell-based, edgebased and node-based smoothing to three-dimension and Figs. 1-4 illustrate how to construct the smoothing domains for each strain-smoothed technique.

Cell-Based Smoothing Domain
In the cell-based S-FEM four sub-cells are constructed by splitting the tetrahedron about its centroidal point. Fig. 1 shows sub-cells and smoothing domains for CS-FEM. In Fig. 1a, it is shown that how a tetrahedral cell with vertices labelled 1234 can be divided into four tetrahedral sub-cells by point 5, the centroid of tetrahedron. These sub-cells are new domains where the socalled smoothed strain can be constructed. Fig. 1b shows an example of the selection of the cellbased smoothing domain when two cells, tetrahedron 1234 and tetrahedron 2345, are neighbored. The cells 1234 and 2345 share a triangular face 234, and points 6 and 7 are the centroids of the respective tetrahedra. When the target cell 2357 is chosen, its smoothing domain is constructed as itself, a sub-cell 2357. The concept of the construction of the cell-based smoothing domain has strong similarities with the standard finite element method.

Face-Based Smoothing Domain
The basic idea of FS-FEM is the extension of the 2D edge-based smoothing method to 3D, meaning that the target edge of a triangular cell is transferred to the faces of the tetrahedral cell as shown in Fig. 2.

Edge-Based Smoothing Domain
Assembling the edge-based smoothing domain in three-dimensions is rather more complex than in its two-dimensional counterpart. Fig. 3 illustrates the division of a cell into sub-cells geometrically and assembly of sub-cells to form a smoothing domain around a target edge. Firstly, to create sub-cells for ES-FEM, five new points must be generated at the centroid of the cell and the centroids of the cell's four faces. Then, as shown in Fig. 3a, sub-cells associated with each edge are created from the nodes of those edges, the cell centroid, and the centroids of the two respective adjoining faces. The number of sub-cells for ES-FEM is six since each tetrahedral cell has six edges. Fig. 3b depicts the construction of the smoothing domain associated with the target edge.

Node-Based Smoothing Domain
For NS-FEM, the tetrahedral cell is divided into four hexahedral sub-domains based around the four nodes. To create the hexahedral sub-domains, we construct 11 additional points: the centroid of the cell (1), the centroid of each face (4), and centroid of each edge (6). Then the sub-cells are constructed from the four original vertices and 11 additional points as shown in Fig. 4a. Assembly of sub-cells from adjoining elements to form the smoothing domain around a target node is illustrated (for the case of node 2) in Fig. 4b.

Numerical Integration
Integrals are evaluated over smoothing domains using a Gauss integration scheme. Sample points, at which shape functions and outward normals are computed, are located in the middle of each face of the smoothing domain [25]. Note that red circles are Gauss points located at the centroid of each surface of the element and n is the outward normal vectors.

Strain Smoothing Approximation
To formulate the non-linear strain smoothing approximation in the following section, we recall the fundamentals of smoothed finite element approximation. The infinitesimal strain tensor ε h ij is given as: where u is the displacement field and X is the initial configuration. The infinitesimal strain tensor is smoothed over the smoothing domain s k as: where a point x k is located in the smoothing domain s k and (x) is the weight function of gradient-smoothed approximation. The properties of the weight function (x) are given as: The smoothed strain can be obtained as follows, applying the divergence theorem: where V k is the volume of smoothing domain s k , s k is the boundary of the smoothing domain s k and n is the outward normal on boundary s k . The 3D form of the outward normal vector matrix is given as: The discrete trial and test functions are expanded in terms of linear Lagrange shape function I with degrees of freedom associated with vertex I of a simplex mesh: When the strain in each sub-domain is constant, the smoothed strain can be written as: where ε (x) is the constant compatible strain of the qth sub-domain of the smoothing domain s k , V k, q is the smoothing domain volume of the qth sub-domain and N SD is a number of subdomains of the smoothing domain s k .

Smoothed Finite Element Approximation in Finite Elasticity, Enhanced by the Bubble Function
To solve the non-linear equilibrium equations, we adopt the Newton-Raphson iterative method. The detailed information could be found in [22,26], but we herein briefly explain the smoothed Galerkin weak form and its linearization. ∂W where v is the set of admissible test functions andF is the smoothed deformation gradient computed on each smoothing domain. The smoothed deformation gradientF is defined as: where G k is a set of nodes and ∂ /∂X is the shape function derivatives. Note that, 1 (10) can be replaced by a set of the outward normal vectors n and the shape functions as shown in Eq. (5).
Eq. (9) can be defined to the energy functionR (u) and its directional derivatives DR (u) · u which are given as: and where i, j, k, l ∈ {1, 2, 3} for three dimensional problem.
To solve Eq. (9) approximately, the Newton-Raphson method is utilized in this work: therefore, the energy functional Eq. (11) and its directional derivatives Eq. (12) can be rewritten: whereC is the smoothed right Cauchy-Green deformation tensor defined asC =F TF and i, j, k, l, p, s ∈ {1, 2, 3}.
Eq. (13) can be assembled into matrix form at every iteration given as: ⎡ where the smoothed tangent stiffness matrixK tan =K mat +K geo can be defined by the smoothed discretized gradientsB 0 andB NL : where N t is the number of target cells, edges, vertices and faces,S sym is a symmetric matrix form of the smoothed second Piola-Kirchhoff stress (PK2)S, andC is the smoothed fourth-order elasticity tensor of the stored strain energy which are defined as: The load vectorb can be expressed as: where S is the vector form of the PK2 stress. Note that the matrix form of the smoothed discretized gradientsB 0 andB NL are the same as the conventional FEM replacing the smoothed deformation gradient.
The smoothed global system of equations for Eq. (13) at each iteration can be written as: thus the displacement u is obtained by the Newton-Raphson iteration:

The Bubble Enrichment
The domain-based selective strain smoothing finite element approximation was suggested due to the fact that the edge-based strain-smoothed method does not fully overcome locking [11]. Recent study of the domain-based S-FEM showed that the selective S-FEM can alleviated locking [19,20]; however, an additional treatment is required the discretized equations for the each strain smoothing approach.
For this reason, Nguyen-Xuan et al. [27] proposed the mixed displacement-pressure formulation in the framework of the edge-based smoothing method using the cubic bubble function and its mathematical properties about the stability of the inf-sup condition were provided by Ong et al. [23]. The bubble function requires the supplement of an additional degree of freedom for the displacement field approximation at the centroid of the element. At the centroid of the element, only the linear displacement field exists as the unknown field in the bubble-enriched gradient smoothing because the strain smoothing approach is based on displacement-based formulation. The bubble function is given with four linear basis functions and the cubic bubble function for the tetrahedral element as follows: Since Eq. (21) does not satisfy the partition of unity, the basis function necessarily requires the transformed form can be obtained as: Hence the transformed bubble basis function can be expressed as: where the properties of the bubble is given as: meaning that the bubble function has a value one at the center of the element and the value is zero on the boundaries of the element.

Numerical Tests
A series of test problems is used to assess the effectiveness of the proposed method in dealing with the shortcoming explained in the previous section: simple torsion of a rectangular block, Cook's membrane, a hollow cylinder and a cube.
To test near-incompressibility we use the following soft materials: polymers, e.g., polyurethane, and human tissues. Both soft media exhibit the nonlinear strain-stress behaviour commonly associated with highly deformable materials, and their Poisson's ratios are close to 0.5. Polyurethane is mainly used in energy absorption applications, such as packaging or the cushion of padded chairs [28]. In this work, we model a Polyurethane spring which is currently used in a base isolation device to regulate mass energy from earthquakes using its reversible deformation property. Human soft tissues can sustain large deformations and, due to their high water content, are nearly incompressible. Correspondingly, they are frequently modelled as hyperelastic materials. In our test case, we employ the quasi-incompressible neo-Hookean model used by Taylor et al. [29] to model brain material.
Another aspect that we consider is the mesh distortion sensitivity. To study the influence of the mesh distortion, finite element meshes are artificially distorted, changing the position of a vertex randomly in each direction.
For all numerical examples the following parameters for the Newton iteration method are used: tolerance 10 −9 , 15∼300 load steps and 4∼6 iterations at each load step, including quadratic convergence of the Newton-Raphson algorithm.

Near-Incompressibility
In this section we investigate the performance of the proposed method for quasiincompressibility, with Poisson's ratio ν → 0.5. The DOLFIN finite element software [30,31] is used as a source of pseudo-analytical solutions for the described simple tension of a rectangular block, Cook's membrane, hollow cylinder and a cube. Python scripts implementing the models are derived from the DOLFIN python demos 1 , and the MINI element [32], which is a finite element with linear basis functions enriched with cubic bubbles for displacements, and linear basis function for pressure, is used.

Simple Torsion of a Rectangular Block
Firslty, the non-homogeneous Dirichlet boundary conditions (BCs) problem in quasiincompressible limit is considered. Simple torsion of a rectangular block is used with human brain-like material [29]: the shear modulus μ = 1000 Pa and the bulk modulus κ = 50000 Pa (equivalent to Poisson's ratio ν ≈ 0.49). The geometry of the rectangular block is given in Fig. 6  For this test, cylindrical coordinates and Dirichlet BCs based on a solid cylinder are needed; however cylindrical coordinates and Dirichlet BCs are transformed to Cartesian coordinates and imposed to a rectangular solid (further details can be found in [26]). The deformation in cylindrical polar coordinates are given as: where the torsion per unit length τ = π/3 is used in this test. The deformation gradient F and the right Cauchy-Green tensor C are defined as: Then, cylindrical coordinates are formulated to the following Cartesian coordinates: where R = X 2 + Z 2 and tan θ = Z X . Therefore, the following non-homogeneous Dirichlet BCs can be obtained: Fig. 7 shows deformed configurations of the rectangular blocks obtained by each method and Fig. 8 illustrates the convergence of strain energy errors for the block, respectively.
The number of elements along to x-and z-axes is 2 and the numbers of elements in y-axis are 4, 8, 16, 24 and 32. As shown in Fig. 7 and Tab. 1, the bubble-enhancement leads to improve the accuracy" FEM and CS-FEM provide identical results. 3D edge-based and node-based smoothing strain methods show the minus values of the relative error in strain energy, meaning their strain energies are smaller than the analytical solution. FS-FEM shows around 31% of the relative error, while its bubble-enhanced approach gives only about 9% of the error. It is clearly shown that the bubble-enhancement is able to reduce almost 71% of the error and it helps to increase.
where strain energy relative error is obtained as follows: Tab. 2 provides the computational time for the standard FEM, the standard S-FEM and the proposed enriched S-FEM. To obtain the computational time, the following computer system is used: Matlab R2019b running on Windows 10 with 3.40 GHz Intel Core i7-6800K and 32 GB memory. It is known that the bandwidth of the stiffness matrix for S-FEM is larger than the standard FEM; therefore, in general, S-FEM needs more computational time than the standard FEM. As given in Tab. 2, FS-FEM and 3D-NS-FEM are around two times slower than the standard FEM and 3D-ES-FEM needs almost three times more computational effort to converge the problem. On the other hand, the computational cost of the bubble-enriched FS-FEM is the most expensive among the present approaches. However, as shown in Fig. 8 and Tab. 1, the bubble-enhancement leads to improve the quality of result meaning that with coarser meshes the bubble-enriched S-FEM provides the relatively accurate results.

Cook's Membrane
In this section, the well-known Cook's membrane is investigated. The geometry is shown in Fig. 9 with L x = 48 m, L y1 = 44 m, L y2 = 16 m and L z = 10 m. Human brain-like materials are also used given as: μ = 1000 Pa and κ = 50000 Pa equivalent to Poisson's ratio ν ≈ 0.49.
The structure is a cantilever, therefore fixed boundary conditions are implemented as: 0, 0, 0), and a bending force (F = 1/16 N) is imposed on the right-hand side surface where X = 48. For this test, four different numbers of DOFs are used: 903, 1320, 2589 and 6678. DOLFIN finite element software is also used to obtain the analytical solution with finer MINI element mesh (40941 DOFs). The deformed configurations of Cook's membrane obtained by reference solution (MINI element) and the proposed strain smoothing approaches are shown in Fig. 10. Fig. 11 illustrates the convergence rate of the relative error in strain energy of Cook's membrane and the detailed values are provided in Tab. 3. As shown in Fig. 11 and Tab. 3, FEM and CS-FEM show identical results and this phenomenon is also shown in the aforementioned non-homogeneous Dirichlet BCs problem. The most accurate result is obtained by ES-FEM with about 0.6% of error and the bubble-enhanced FS-FEM provides around 0.81% of error in strain energy. As shown in Tab. 3, the bubble-enhanced FS-FEM gives relatively accurate results among others, while FS-FEM shows almost 6% of error. It is clearly seen that the bubble-enhancement leads to an increase in the accuracy of results by about 86% compared to the standard FS-FEM. The notable result obtained in the present test is the bubble-ES-FEM which produces a worse result that is around eight times higher error than ES-FEM and as well as NS-FEM shows three times higher error.

Pressure of a Hollow Cylinder
In this section, Polyurethane MER (Mass Energy Regulator) spring used for the seismic base isolation unit (EQS, Eradi Quake System), developed by ESCO RTS 2 . The basic idea of the EQS unit is that the seismic force can be absorbed by the disk bearing and the horizontal force is dissipated by the MER spring.   Considering the instability of NS-FEM and bES-FEM observed in the previous section, only CS-FEM, ES-FEM, FS-FEM and the bubble-FS-FEM are examined in this section. Fig. 13 illustrates deformed configurations of the polyurethane spring for respective S-FEM approaches and for the MINI element obtained by DOLFIN software.
The deformation predicted by bFS-FEM is similar to that of the reference solution, while other methods produced somewhat different patterns. Distinct improvements in the strain energy results were also observed, as shown in Tab. 4. CS-FEM and FEM, again, produce the same results, and cannot avoid locking in quasi-incompressibility. The edge-based smoothing is more accurate than the face-based smoothing; however, the bubble-enrichment notably increases the accuracy retaining the stability. The bubble-face-based smoothed finite element method is able to improve the accuracy about three times and it also shows around 26% less error in strain energy compared to ES-FEM.

Mesh Distortion Sensitivity
Lastly, in this section, the mesh distortion sensitivity for the proposed approach for quasiincompressible problem is investigated. For this test, regularly distributed finite element cell meshes are artificially distorted by the following terms [24]: where X , Y , Z and X , Y , Z are the initial and updated coordinates respectively, r c is a random number −1 ≤ r c ≤ 1, α distortion is the distortion magnitude and x, y, z are the initial mesh sizes.   The pseudo-analytical solution obtained with DOLFIN software using MINI element with regular fine meshes. The convergence of displacement and strain energy errors and the detailed relative errors of displacement and strain energy for the mesh distortion problem can be found in Fig. 16 and Tab. 5, respectively. MINI element gives lower strain relative error but higher displacement error than cell-based, face-based and edge-based smoothing methods. This result shows that the MINI element effectively handles volumetric locking in nearly-incompressible media, compared to normal smoothed-strain approaches. However, MINI element is still suffered from highly distorted meshes. The bubble-enriched-FS-FEM shows the improvement of accuracy in the both displacement and strain energy rather than the standard FS-FEM: around 20% of improvement in displacement and 85% of improvement in strain energy. Moreover, bFS-FEM provides more accurate result in strain energy although displacement relative error is slightly higher than ES-FEM, reducing almost half of strain energy error. A remarkable result we obtain from this problem is the performance of bubble-enhanced face-based smoothing approach is rather better than ES-FEM for the combined mesh distortion and nearly-incompressible problem.

Conclusions
In this work, the strain-smoothed finite element approaches were extended to the problem of three-dimensional hyperelasticity. On each smoothing domain, the smoothed deformation gradients are constructed, allowing expressions for the smoothed strain-displacement matrix, the smoothed tangent stiffness matrix and the smoothed internal force vector to be obtained. The cubic bubble enrichment function for the 3D face-based and the edge-based strain smoothing methods is also introduced. The proposed bubble function is straightforwardly implemented into the existing smoothed-strain framework.
The proposed S-FEMs and their bubble-enhancement are examined by numerical benchmark tests to evaluate their immunity to volumetric locking and mesh distortion. Comparisons are made with solutions obtained by DOLFIN using the well-known MINI element. Quasi-incompressible media, ν > 0.49, under large deformations, and with severely distorted finite element cells are the main considerations of the proposed methods.
From the results carried out in numerical tests, the following conclusions are obtained: • the behaviour of CS-FEM is similar to FEM in quasi-incompressible media; that, the cellbased strain smoothing method is not able to avoid the over-estimation of stiffness; • the accuracy and stability of the edge-based and node-based smoothed-strain approaches are problem-dependent, which means ES-FEM and NS-FEM occasionally fail to solve the problem; and • the bubble-enrichment in the framework of the strain smoothing approach shows more accurate and reliable results, especially in the quasi-incompressible case, than do to the standard S-FEM.
A notable issue of the instability of 3D-extended NS-FEM and ES-FEM is found through certain numerical tests, for which those methods either failed to converge or produced substantially higher error. FS-FEM, especially the bubble-enhancement version, however, provides a very stable solution and relatively accurate results compared with FEM or even other S-FEMs. Because the concept behind the face-based smoothing method has been geometrically derived from the 2D edge-based S-FEM, bFS-FEM is more reliable than the direct transfer of ES-FEM and NS-FEM into three dimensions.
Overall, the bubble-enriched face-based strain-smoothed approach with low order tetrahedral cells provides a computationally efficient, yet stable and accurate solution for large deformation quasi-incompressible problems.
In future work, we plan to integrate the proposed method with the efficient total Lagrangian Explicit Dynamics (TLED) algorithm [33] which should further reduce solution times, especially when deployed on massively parallel GPU hardware [33].

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.