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.

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 [

Various 3D strain smoothing approach variants have been introduced. Firstly, Nguyen-Thoi et al. [

The application of main interest of this work, i.e., soft materials, has been also considered by Jiang et al. [

Instead, in the present work, we employ a cubic bubble function to enrich S-FEM [

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 bubble-enhanced 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.

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 sub-domains 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 [

The salient abilities of S-FEM are its insensitive to locking and to element distortion [

Four different types of S-FEM have been introduced, depending on how the smoothing domains are built [

In the cell-based S-FEM four sub-cells are constructed by splitting the tetrahedron about its centroidal point.

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

In the same manner to CS-FEM, the tetrahedral cell is split into four tetrahedral sub-domains by the centroidal point of the cell as shown in

Assembling the edge-based smoothing domain in three-dimensions is rather more complex than in its two-dimensional counterpart.

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

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 [

Note that red circles are Gauss points located at the centroid of each surface of the element and

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

where

where a point

The smoothed strain can be obtained as follows, applying the divergence theorem:

where _{k}

The discrete trial and test functions are expanded in terms of linear Lagrange shape function

When the strain in each sub-domain is constant, the smoothed strain can be written as:

where _{SD}

To solve the non-linear equilibrium equations, we adopt the Newton–Raphson iterative method. The detailed information could be found in [

where

where _{k}

and

where

To solve

therefore, the energy functional

where

where the smoothed tangent stiffness matrix

where _{t}

The load vector

where

The smoothed global system of equations for

thus the displacement

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 [

For this reason, Nguyen-Xuan et al. [

Since

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.

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 [

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},

In this section we investigate the performance of the proposed method for quasi-incompressibility, with Poisson’s ratio ^{1}

Firslty, the non-homogeneous Dirichlet boundary conditions (BCs) problem in quasi-incompressible limit is considered. Simple torsion of a rectangular block is used with human brain-like material [

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 [

where the torsion per unit length

Then, cylindrical coordinates are formulated to the following Cartesian coordinates:

where

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

Method | Strain energy |
Relative error (%) |
---|---|---|

Analytical solution | 2.0613 | – |

FEM T4 | 3.4950 | 69.55 |

CS-FEM T4 | 3.4950 | 69.55 |

FS-FEM T4 | 2.7071 | 31.33 |

ES-FEM T4 | 1.8473 | −10.38 |

NS-FEM T4 | 1.2833 | −37.74 |

FS-FEM with the bubble | 2.2481 | 9.06 |

where strain energy relative error is obtained as follows:

Method | Computational time ( |
---|---|

FEM T4 | 1.005538 |

CS-FEM T4 | 1.418809 |

FS-FEM T4 | 1.870802 |

ES-FEM T4 | 2.757793 |

NS-FEM T4 | 2.215367 |

FS-FEM with the bubble | 147.109921 |

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

In this section, the well-known Cook’s membrane is investigated. The geometry is shown in

The structure is a cantilever, therefore fixed boundary conditions are implemented as:

Left surface

and a bending force (

The deformed configurations of Cook’s membrane obtained by reference solution (MINI element) and the proposed strain smoothing approaches are shown in

Method | Strain energy |
Relative error (%) |
---|---|---|

Analytical solution | 3.9652 | – |

FEM T4 | 3.5883 | −9.5046 |

CS-FEM T4 | 3.5883 | −9.5046 |

FS-FEM T4 | 3.7302 | −5.9251 |

ES-FEM T4 | 3.9895 | 0.6139 |

NS-FEM T4 | 4.0481 | 2.0914 |

FS-FEM with the bubble | 3.9329 | −0.8116 |

ES-FEM with the bubble | 4.1875 | 5.6075 |

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 details of the device are given in

To model polymer polyurethane material, neo-Hookean hyperelastic limit is again used. Shear and bulk moduli are

Bottom surface

Top surface

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.

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

Method | Strain energy ( |
Relative error (%) |
---|---|---|

Analytical solution | 10.5202 | – |

FEM T4 | 11.1317 | 5.81 |

CS-FEM T4 | 11.1317 | 5.81 |

FS-FEM T4 | 10.3227 | 3.68 |

ES-FEM T4 | 10.9084 | −1.88 |

FS-FEM with the bubble | 10.4747 | −1.39 |

Lastly, in this section, the mesh distortion sensitivity for the proposed approach for quasi-incompressible problem is investigated. For this test, regularly distributed finite element cell meshes are artificially distorted by the following terms [

where _{c}

A cube model under compression, with edge lengths 0.1 m (see in

Bottom surface

and a pressure of

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

Method | Displacement | Relative error (%) | Strain energy ( |
Relative error (%) |
---|---|---|---|---|

Analytical solution | −0.0143 | – | 3.09 | – |

MINI | −0.0134 | −6.38 | 3.02 | −2.54 |

CS-FEM T4 | −0.0132 | −7.87 | 2.71 | −12.52 |

FS-FEM T4 | −0.0134 | −6.14 | 2.82 | −8.89 |

ES-FEM T4 | −0.0137 | −4.11 | 3.01 | −2.69 |

FS-FEM with the bubble | −0.0136 | −4.97 | 3.05 | −1.29 |

ES-FEM with the bubble | −0.0152 | 6.54 | 3.27 | 5.63 |

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,

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 cell-based 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 [

ESCO RTS Co., Ltd., (