Blog

Modeling Elastomers using FEM: Do’s and Don’ts

Elastomers are increasingly used in industrial applications and understanding their mechanical behavior is even more important than earlier. Tires being one of the largest application of rubbery materials, some of the other very common industrial applications include seals, bearings, dampeners, O-rings etc.

O-rings were invented in the 1890s and used to seal fluids. O-rings are almost present in all industrial and commercial applications and work on the concept of balancing the fluid pressure with the contact stresses. As the fluid pressure increases the O-ring is compressed in the direction normal to the fluid surface as shown in Fig. 01. Elastomers being incompressible, they expand in the transverse direction to increase the contact stresses and effectively seal the interface. The sealing does not fail as long as there is no mechanical failure of the O-ring. This is one of the primary reasons they are still widely used.

Fig 01: O-ring subjected to fluid pressure. The fluid is stopped by the O-ring in (a). As the fluid pressure increases in (b) and (c), the O-ring is pushed in to increase the contact stress and seals the interface.

In this article, we discuss a series of tips for modeling elastomeric materials using the Finite Element Method (FEM) to improve the accuracy and reliability of obtained solutions.

Choosing the right material model for the elastomers

Elastomeric and rubbery materials typically demonstrate:

Isotropic and elastic behavior without any permanent set (like plasticity)

Nonlinear stress-strain behavior as shown in Fig 02.

Different tensile and compressive behavior and comparatively softer in tension

Most of these characteristics are captured well through Hyperelastic models and depending on the material behavior a suitable model should be selected. The material parameters are obtained through experiments like uniaxial & biaxial tensile test, uniaxial compression test, planar & simple shear and volumetric test. A more detailed discussion on the selection of hyperelastic material model can be found on the SimScale blog.

Tip 01: Using the parameters highlighted in this blog article, choosing the right hyperelastic material model is half the job done well.

Meshing the structure

At present, SimScale allows for the meshing of the structures with Tetrahedrons. Hence, in this article, we will limit the discussion to meshing 3D structures with tetrahedrons. A general overview on meshing & meshing tutorial is available in the SimScale documentation. Presently, SimScale allows tetrahedral meshing with following options:

Automatic

Parametric

Tetrahedral with local refinement

As shown in Fig. 03, all these cases allow for the specification of mesh order: first or second. Using a first-order generates a mesh with four-noded tetrahedral elements and second-order a mesh with ten-noded tetrahedral elements. For more detailed tips on meshing for structural mechanics problems, watch out for upcoming articles.

Fig 03: Order of interpolation

Tip 02: For most elastomeric structures second order would be the safe choice. However, in some cases, a first-order mesh could also be sufficient and save significant computational time.

Considering geometry & loading

The dimensions of the structure also play a pertinent role in the accuracy of solutions obtained. Many of the applications of elastomers, like conveyor belts, have one or more dimensions (like length & breadth) much larger than the others (thickness).

For example, as shown in Fig. 04, if the thickness and breadth of the structure are much larger than its length, then it can be considered as a beam. Similarly, if the thickness is much smaller than both length & breadth, then the structure can be considered as a plate or shell. If all the dimensions are comparable, then the structure can be considered to be a full solid. Irrespective of the structure type, solid elements are being used to mesh the structure. Thus, is becomes important to understand the form of a structure being considered so that effective meshes can be generated.

Fig 04: Illustration of beam, plate, shell & solid structures

It is important to first identify the category of the structure. For solid structures, the automatic tetrahedral meshing can be directly used and both the first and second order could yield good results. However, the complexity arises for structures that resemble beams / plates / shells.

Secondly, in addition to geometry, it is also important to identify the type of loading these structures are subjected to. For example, the O-ring illustration shows that the loading on the elastomer structure is of bending-type in nature.

When the structures are of the form of beams-plates-shells, they can be considered as thin structures. These thin-structures demonstrate a locking effect (well known as shear locking) especially when a first-order mesh is used and structure is subjected to bending loads. The easiest way to circumvent the shear locking is through mesh refinement or using higher-order elements. A more detailed discussion on the shear locking is dealt in detail in our upcoming blog article on meshing tips for structural mechanics problems.

Tip 03: If you are using a solid structure, first-order might still work. If the structure can be classified as one of the thin structures and/or subjected to bending loads, then second-order mesh is strongly recommended.

Dealing with incompressibility in elastomers

Most elastomers and rubbers have a Poisson ratio of almost 0.5. Physically, this means that when the material undergoes elongation in the x-direction, then there is enough compression in the other two directions such that the initial and final volumes remain same. Note that a Poisson ratio of 0.5 is equivalent to saying that the Bulk modulus is infinity.

Or otherwise, this means that it would take an infinite amount of energy to compress the material. To get a feel of the idea, one could consider water. Water is a perfect example of an incompressible material. If one would fill a can with water and try to compress, it can be seen that the volume of water remains unchanged. So one could say that the bulk modulus of water is almost infinity. Since infinity cannot be used numerically, a very large bulk modulus is used to signify Poisson ratio tending to 0.5.

When the bulk modulus is considered such that the Poisson ratio is larger than 0.4, most normal elements show a locking effect. This is commonly known as volumetric locking. Unlike the earlier discussed shear locking, volumetric locking cannot be circumvented through mesh refinement. Volumetric locking is not only common to elastomers but also pertinent to any elastoplastic material.

There are several approaches available in the literature:

Selectively reduced integration (B-Bar method)

Uniformly reduced integration (or more commonly just known as reduced integration)

Enhanced strain formulation

Mixed formulations

First, it would be important to understand the cause of this volumetric locking before finding a solution. Considering a simple example of a loaded structure made of elastomer in 2-D, as shown in Fig 05, the cause of volumetric deformation is observed. The mesh is made of two triangles (or could be one quadrilateral either) with linear interpolation (first-order mesh). A small force (or alternatively displacement) is applied on the free node (node number 01) along the direction of 45-deg.

Fig 05: Illustration of volumetric locking

Since the material is volume preserving and the elements use linear interpolation (i.e. lines 1-2 and 1-4 need to be straight), the only direction the free node can move is along the circle shown in red color. However, this would mean moving perpendicular to the effective reaction force and such a motion cannot be allowed. Eventually, the only possible solution is the trivial solution of nodal displacements being zero (or very small numerically). This is clearly disastrous!

One of the simplest ways to solve the problem is using second (or higher) order elements. When second-order elements are used, the edges do not need to be straight anymore. This would facilitate motions that could lead to non-zero nodal displacements and still conserve the volume. However, as the Poisson ratio approaches 0.499, 0.4995 etc, even the second-order elements do not behave well. To understand this behavior of second-order elements, it is important to understand the mathematics behind the interface.

In problems involving elastomers (or nonlinear FEM), matrix equations [K]{u} – {f} = {R} are solved such that {R} tends to zero. The stiffness matrix [K] is first computed locally for each element and then assembled to a bigger global system. Locally, in each element, [K] is computed at fictitious points commonly known as Gauss (or integration) points. The above problem with second-order elements is caused since the volumetric strain is being enforced to be zero at all of these integration points. The total energy from volumetric strain can be approximately written as a product of Bulk modulus x (volumetric strain)^2. Since the bulk modulus is very large, even a small volumetric strain can lead to large extra addition to the energy. This renders the total equation system to behave very stiff and the solvers more often converge to the trivial solution (of zero displacements).

An alternative and very commonly used approach, also incorporated in SimScale, is using reduced integration with second-order elements. This means that the deviatoric strains are computed as earlier but the zero volumetric strain is enforced only at certain points. Though this might sound to be inaccurate, reduced integration has been observed to be extremely robust and yield quite accurate solutions across a wide range of Poisson ratio.

Note that first-order tetrahedrons use only one integration point and thus a reduced integration with these elements is not possible.

Tip 04: If the Poisson ratio is less than or equal to 0.4, first order elements can be used. For Poisson ratio’s even up to 0.45, second-order tetrahedrons could work well. For much higher Poisson ratios, it is strongly recommended to use second-order elements with reduced integration.

Force vs. displacement boundary condition for modeling elastomers

One of the common causes of non-convergence is unbalanced residuals. When modeling elastomers, the system of equation being solved is nonlinear, i.e. [K(u)] {u} – {f} = {R}. Here the stiffness matrix [K(u)] is a function of the unknowns {u}. Thus, small displacements / forces are applied in small increments to reach the desired value.

Using smaller time steps (i.e. smaller force increments) can significantly improve the convergence properties. Alternatively, if possible, displacement controlled boundary conditions are most likely a safe option. As shown in Fig 05, where a very highly nonlinear problem is considered, for every displacement there is a unique solution of the force. However, the same is not true vice-versa.

Fig 05: Nonlinear force-displacement curve showing that for a particular force there can be multiple displacements (or non-unique solution). But for each displacement, there is only one force (or a unique solution).

This means that when a force boundary condition is applied, it need not lead to a unique displacement. Such non-uniqueness can easily lead to a lack of convergence. Alternatively, using a displacement boundary condition can likely lead to better convergence. In such cases, advanced techniques like line search etc are employed to improve the overall convergence of the system. Watch out our SimScale blog for the upcoming article related to solvers for structural mechanics problems that will delve more into areas like line search etc.

Tip 05: If possible, use displacement boundary conditions.

Convergence issues

There are several reasons why a simulation can fail and show a lack of convergence. Some of the very common ones pertaining to elastomers include, but not limited to,

Material instability

Time step size

Inelastic effects

Material instability

The material parameters are generally obtained by fitting the experimental data. The fitting is limited by the maximum stretch to which the data is available. If the model is used beyond this limit, the model might not be unconditionally stable.

For example: Say if the data was only available for a stretch to 30% and the model parameters were fitted with this. When a simulation is done using this model where the strains are much larger (> 30%), it is possible that the nonlinear behavior beyond 30% is not accurately captured.

Secondly, depending on the number of tests used for fitting, it is possible that the material demonstrates an instability when subjected to a different type of loading. For example: Say material parameters for Mooney-Rivlin model are fitted using only uniaxial tests. When a biaxial simulation is done, it is possible that the results show significant instabilities.

A possible solution to this would be to test the material model with uniaxial, biaxial, shear, volumetric tests for up to failure limits. Using these simulations, possible points of material instabilities can be identified a priori.

Time step size

The second possible reason for the lack of convergence can be time step size. If a manual time step is used and if the time step is too large then it can lead to a lack of convergence. SimScale allows for auto time-stepping and it is recommended to use auto-stepping.

Fig 06: Screenshot for auto time stepping

As shown in Fig. 06, a retiming event is used to kick-in to reduce the time step. When the time step is too large and leads to non-convergence, this kicks in a retiming event. In such a case, the time step is reduced and the simulation restarts from the last converged time step. This significantly reduces the failure of simulations. There are several options available for time step calculation type and “Newton Iterations Target” option is a good place to start.

Inelastic effect

In this article, we have considered the elastomers to be purely hyperelastic and we will discuss the effects of inelasticity in an upcoming blog article

Tip 06: Check the material parameters thoroughly for instability before using it under general loading conditions. Auto time stepping is a great option to prevent failure of simulations!

Overall, summarizing the type of problem or procedure, types of element recommended and to avoid are:

SimScale is the world's first cloud-based simulation platform, enabling you to perform CFD, FEA, or thermal analyses. Sign up for the 14-day free trial and join the community of 70 000 engineers and designers. No payment data required.

Ajay Harish

Ajay Harish is a Doctoral candidate with Prof. habil. Dr.-Ing. Dr. h.c.mult. Dr.-Ing. E.h. Peter Wriggers at the Institute of Continuum Mechanics (Leibniz Universität Hannover, Germany). His primary focus of research includes constitutive modeling of soft materials - granular materials - composites, nonlinear fracture mechanics, contact mechanics and automation of computational modeling. His official profile can be found on IKM Website (https://goo.gl/hKWoMa). More information can be found on his personal website (http://www.ajaybharish.com).