Abstract
Elastomeric components are widely used in the engineering field since their mechanical properties can vary according to a specific condition, enabling their applications under large deformations and multiaxial loading. In this context, the present study seeks to investigate the main challenges involved in the finite element hyperelasticity simulation of rubber-like material components under different cases of multiaxial loading and precompression. The complex geometry of a conical rubber spring was chosen to deal with several deformation modes; this component is in the suspension system placed between the frame and the axle for railway vehicles. The framework of this study provides the correlation between axial and radial stiffness under precompression obtained by experimental tests in prototypes and virtual modeling obtained through a curve fitting procedure. Since the material approaches incompressibility, different shape functions were adopted to describe the fields of pressure and displacements according to the finite element hybrid formulation. The material parameters were accurately adjusted through an optimization algorithm implemented in Python program language which calibrates the finite element model according to the prototype test data. However, as an initial guess, the proper constitutive model and its parameters were first defined based only on the uniaxial tensile test data, since this test is easy to perform and well understood. The validation of the simulation results in comparison with the experimental data demonstrated that care should be given when the same component is subjected to different multiaxial loading cases.
1. Introduction
Since new products development has always been increasing, several efforts have been focused on reducing time and costs. As the usual prototype construction is generally very expensive and time consuming, numerical simulation tools like the finite element method should become more reliable.
Even though a number of materials are available, elastomers play an important role in the product properties definition since it fulfills a wide range of functional tasks. Their typical applications are related to vibration damping, acoustic insulation, sealing, coupling, tires, consumer goods, and so on. The high elasticity and the ability to support large deformations under external forces are their basic characteristics [1, 2].
In railway industry the rubber springs are very useful, since they accommodate oscillatory motions and allow axes misalignments. They are located in the primary suspension, between the bogie frame and the axle box. Although a wide range of geometries can be found in order to fulfill user requirements, the conical rubber spring is one of the most used designs, since it acts as a universal damping and guiding element [3, 4]. It consists of rubber pads bonded on its inner and outer metal surfaces layers conically arranged. Thus, it is worth mentioning that this component needs to combine high stiffness for axle-guiding and an optimally vertical rigidity due to damping purposes.
Besides the conical rubber spring, the most common rubber components are subjected to different multiaxial loading cases in the real working conditions. For this reason, choosing a proper material model to be implemented in the numerical simulations is necessary to better describe the component behavior under each loading case [5].
Since rubbers are materials that exhibit nearly incompressible behavior and often experience high strains in service, their strain states are usually very complex. Thus, their mechanical behavior should be defined based on strain energy potentials formulated according to the hyperelasticity theory. They are a mixture of tension, compression, and shear with a very small amount of volume change.
The strain energy density function () defines the strain energy stored in the material per unit of reference volume. This function depends on the principal stretches or invariants of the strain tensor and it is directly linked to the material’s stress-strain relationship which depends on a series of parameters (material constants) [6]. Therefore, in order to assess the rubber-like material constitutive law, some experiments should be performed for input data in curve fitting procedures for the Finite Element Analysis (FEA).
Several published works related to rubber bushing have been extensively researched since last decades. The first studies have been performed by Adkins and Gent [7], where the authors found a formulation to predict the stiffness variation according to the bushing length for four principal modes of deflection termed: torsional, axial, radial, and tilting. Hill [8–10] performed studies related to the radial deformations of bonded cylindrical rubber bush mountings and its precompression effect. Petek and Kicher [11] experimentally investigated the nonlinear behavior of a shear bushing with conical ends focusing on the axial quasi-static load/deflection properties.
In addition to the empirical bushing models, constitutive modeling through FEA has also been studied. Despite the main challenges about computational efforts and time-consuming process, Morman and Pan [12] performed studies comparing the closed-form analytical equations developed for the application in the design of elastomeric components and simulation response through FEA. They argued that the closed-form equations should not provide accurate results in case of more complicated geometries and complex boundary conditions.
Besides the studies presented above, Horton et al. [13–15] derived more accurate expressions for annular rubber bush mountings subjected to radial loading and tilting deflection based on the classical theory of elasticity. In the meantime, Luo et al. [3] obtained the stress state of a simplified conical rubber spring by FEA according to the von Mises criterion. They compared the prototype experimental data with simulations, but only focusing on strength and durability of metal parts of the system. Kadlowec et al. [16] performed studies in which annular bushings were subjected to radial, torsional, and coupled radial-torsional modes of deformation. They compared the elastic bushing response obtained experimentally with finite element results. As far as the torsional mode is concerned, Horton and Tupholme [17] derived new closed-form expressions for the torsional stiffness of spherical rubber bush mountings in the two principal modes of angular deformation. Despite the reasonably agreement of some analytical relationships, there have been some limitations in the use of closed-form equations, being not possible to cover all the real situations.
Thus, for the last ten years the modeling of rubber components through the finite element method has been increasingly used. Gil-Negrete et al. [18] predicted the dynamic stiffness of filled rubber isolators using a finite element (FE) code. Olsson [19] presented a method to analyze the dynamic behavior of rubber components under radial loading by considering an overlay of viscoelastic and elastoplastic finite element models. Gracia et al. [20] studied two types of filled rubber industrial components subjected to several loads using FE analysis with the overlay model. Finally, Lee et al. [21] proposed a hybrid method based on FEA and empirical modeling to obtain the hysteresis of suspension rubber bushings and predict the dynamic stiffness without performing iterative experiments and avoiding high computational costs.
Since the rubber properties for FEA should be defined from experimental data, Seibert and his coworkers [22] studied an optimized specimen’s shape through biaxial extension to calibrate the hyperelastic material parameters for the simulation of an engine mount placed in a car under multiaxial loading. Furthermore, Kaya et al. [23] simulated the behavior of a vehicular rubber bushing using FEA and performed the shape optimization to redesign its geometry and to meet the target static stiffness curve. Lalo and Greco [6] compared the rubber bushing behavior extracted experimentally with hyperelastic models based on shore hardness and uniaxial extension.
Based on the apparent extent of earlier works in modeling bushing response through FEA, it is possible to state that the hyperelastic constitutive models predefined in finite element codes are able to simulate geometry-independent bushings. However, some challenges need to be overcome when dealing with different multiaxial loading cases. In this situation, the present paper seeks to implement an optimization algorithm to characterize the hyperelastic constitutive model according to the applied multiaxial loading. Although the loading directions are axial and radial, the complex geometry of the component results in multiaxial modes of deformation on the rubber pads. It will be shown that care should be taken when the load direction is changed and the previous characterization cannot be valid anymore.
2. Hyperelasticity Theory
Hyperelastic material models are very usual to define the mechanical behavior of elastomers, foams, and many biological tissues. Since elastomers are materials that exhibit nearly incompressible behavior and often experience high strains in service, its states of strain are usually very complex [6].
The constitutive laws for hyperelastic materials are basically defined in two main theories: phenomenological or micromechanical. The phenomenological models capture the overall behavior of some polymers and fit their experimental data with reasonable accuracy aiming to minimize computational effort. The first representative formulations were developed in the works of Mooney [24] and Rivlin [25]. After that, Rivlin and Saunders [26] have proposed a new constitutive model, currently known as Mooney-Rivlin. In the 90s, Yeoh [27, 28] proposed a model by truncating the polynomial series to the first invariant and added more terms in order to increase accuracy. The simplest form of Rivlin’s strain energy function is the neo-Hookean, based only on the first term [28]. The phenomenological models can also be formulated in terms of principal stretches as in Valanis and Landel [29] and Ogden [30].
On the other hand, the micromechanical models are based on statistical mechanics to capture the network evolution of cross-linked polymer chains and to predict the three-dimensional material response. This model can be applied to unusual polymers types, being first studied in the works of Flory and Rehner [31], James and Guth [32], Treloar [33–35], and Wang and Guth [36]. Later, the classical models such as van der Waals proposed by Kilian et al. [37], Arruda-Boyce [38], and Gent [39] were developed and nowadays are implemented in several FEA commercial software. In 2002, Pucci and Saccomandi [40] proposed some modifications in the Gent model seeking to obtain a very good fit to uniaxial data over the full range of deformations. Apart from these models many others have been developed over the years.
The strain energy functions (called “”) used in hyperelastic models are derived according to the invariants of Green strain tensor which are expressed in terms of principal stretch ratios:where the three invariants are given in terms of the principle stretch ratios , , and by
Aiming to minimize the numerical instabilities, when the material compressibility is a concern, it is necessary to split out the deviatoric (represented by “superscript bar”) and volumetric terms of the strain energy function, represented by and , respectively. As a result, the volumetric term is a function of the volume ratio only and the deviatoric part depends on the modified stretches ():
For large deformations, it is usual to work with the Cauchy stress tensor [41], where its representation can be given through the principal stretches as in where the first term is relative to the deviatoric (the deviatoric stress is related to shape change and it is what is left after subtracting out the hydrostatic stress) part of stress, while refers to the hydrostatic (the hydrostatic stress is related to volume change) part of stress which is associated with the incompressibility condition.
2.1. Constitutive Models
The strain energy function particular models are defined according to a number of material parameters. Thus, the nominal stress versus nominal strain data obtained from experimental tests is required to fit these parameters in most models’ theoretical behavior available. Although several constitutive models have been developed over the years, each one has its peculiarities of application [42, 43]. For this reason, the present paper will be focused on some classical constitutive models which are already implemented in the commercial software of FEA Abaqus®.
2.1.1. Neo-Hooke
The neo-Hookean material model was formulated by truncating the power series of a polynomial form [44]. This is the simplest phenomenological model since it is derived only in terms of the first invariant ().where is a material parameter for data fitting and relates to the compressibility ratio.
Treloar [45] also formulated the neo-Hookean model according to the statistical theory of rubber elasticity. Even though the phenomenological and statistical theories were formulated from quite different premises, they are considered equivalent.
This model can be a good starting point. However, the strains should be limited up to 30-40% in uniaxial extension and up to 80-90% in shear deformation [46].
2.1.2. Mooney-Rivlin
The Mooney-Rivlin model is very usual for representing the large strain nonlinear behavior of elastomers. It is also a phenomenological model and gives a good response for moderately strains in uniaxial extension and shear deformation [47]. The first order Mooney-Rivlin model can be defined by
This model is based on a polynomial expansion. Thus, to obtain a order, more terms can be added to (9). However, according to Sasso et al. [48] the results do not show major improvements.
2.1.3. Ogden
Proposed by Ogden [30], this is a phenomenological model based on the principal stretches rather than strain invariants. This model has been effective for very large strains and it is written aswhere and are real material parameters, being positive or negative.
If a two-term expansion () is applied with and , it coincides with the 2 parameter Mooney-Rivlin form for typical values of and . If is set to 1 with and , it degenerates to the neo-Hookean form [41]. The Ogden model is very useful for simple extension up to 700% [46].
2.1.4. Yeoh
Yeoh [27, 28] introduced a phenomenological model only based on the . Since a polynomial form of is used but the other deviatoric strain invariants are not, the Yeoh model can be also called reduced polynomial model and it has the general form:
According to Yeoh [28], the terms containing were eliminated from the equation since its variation in the sensitivity of function is negligible if compared to . Thus, it is possible to improve the model’s ability when predicting the behavior of deformation states over a large strain range even when limited test data is available.
The Yeoh model is commonly considered with , and the upturn of the stress-strain curve can be captured. This model is generally applied in the characterization of carbon-black filled rubbers [27] but is not useful when low strains are involved.
2.1.5. Arruda-Boyce
Proposed by Arruda and Boyce [38], this model was developed based on a statistical treatment of the non-Gaussian chains, and for this reason it is also known as the eight-chain model. As this model represents the physics of network deformation, they are micromechanical models and can be described as follows:where the material constants are predefined functions of the limiting network stretch “” as follows:
From (12), it is possible to note that the Arruda-Boyce model is similar to the Yeoh form, but with .
According to Arruda and Boyce [38], this model is unique since it is able to provide great accuracy for multiple modes of deformation based only on the standard uniaxial tensile test.
2.2. Finite Element Hybrid Formulation
In this approach, pressure is treated as an uncoupled variable and it must receive a suitable formulation through the weak form of the Finite Element Method (FEM). Since in almost incompressible problems the bulk modulus is relatively high, its relation to the pressure is given through the Cauchy stress tensor:
Linearizing (15) in relation to the displacement field [41], there is
Therefore, it is worth mentioning that pressure has less variation than the displacements field and, for this reason, it must be treated as an uncoupled term. That is, it must be treated separately as an independent variable in the FEM formulation [49]. If the problem is effectively incompressible, a relatively high value of the variable is assigned, forcing it to an incompressibility condition solved by a penalty procedure, in which the pressure variables act as Lagrange multipliers to force the condition of incompressibility.
In view of the above, the hybrid formulation can be used by FEM, where different shape functions are adopted to describe the fields of pressure and displacement. In the case of pressure, a lower order variation must be attributed:where refers to the nodal pressure variables constituting the vector and the nodal displacements that constitute the vector .
Thus, (17) can be rewritten as
All in all, the boundary value problem can be expressed in the weighted form, integrating the solid domain through a weight function. Therefore, through the Principle of Virtual Work (PVW) the equilibrium condition of the internal forces and external in the volume can be obtained [41].
The forces and pressures variables coupling vector is called and can be derived aswhere is the pressure variation in terms of nodal variables and . The subscript shows that the tensor refers to the 2nd Piola-Kirchoff stress tensor. Thus, the differentiation of the matrix is given by
In addition, the pressure-displacement relationship of (15) should still be considered. In this case, the Galerkin Method must be applied to obtain the weak form of (15) by multiplying it by and integrating it over the element.
This relationship should hold for any , where represents the lack of pressure compatibility:
The equations above represent the governing relationships between the displacement and nodal pressure variables. For a purely incompressible problem, the term from (23) would disappear and this equation could be used in the weak form of the incompressibility constraint.
It is worth mentioning that the governing equation (23) involves the pressure and not its derivative, so the pressure to be continuous between the elements is not necessary. Hence, they can be treated as internal variables of the system. As a way to reduce the number of additional unknowns, in general one of these conditions (force or displacement) is admitted in the strong form, while the second is the weak form obtained by weighting. Since the formulations were adopted in the almost incompressible materials, the weak form has been attributed to the pressure.
3. Data Fitting Optimization Algorithm
Based on optimization techniques the simulation models can be calibrated by minimizing different error measures related to force-displacement experimental response.
The problem arises from a set of arguments that give a minimal function value under a set of constraints. The general mathematical description of an optimization problem can be expressed as indicated in (24), where by convention the problem is often represented as the minimization of an objective function.where the vector refers to the material parameters assigned as design variables and is the objective function which represents the error function between experimental data and simulation model. The constraints mean the set of requirements the project must meet to be admissible as a solution, being in this case, represented by and , respectively. The parameters and are vectors that respectively correspond to the values of the equality and inequality constraints for each position and .
The implemented pattern search optimization method was proposed by Hooke and Jeeves [50] and does not require a gradient calculation. This algorithm examines points near the current point by perturbing design variables until an improved point is found. It then follows the favorable direction until no more design improvement is possible. The convergence is detected according to a Step Size Reduction Factor until the step-length is sufficiently small or when a maximum number of runs are reached. The implemented algorithm can be found in [51].
Since the nature of the objective function is unknown a priori, this is a well-suited method for the curve fittings obtained by finite element simulation. Performing some adjustments of the tuning parameters like step sizes and number of iterations, it finds the optimal solution even if multiple local minima are present. However, if a quadratic function is used, it is obviously less efficient than gradient methods. The optimization process coupled with finite element simulation is depicted in Figure 1. The design variables change automatically throughout the process and the force-displacement curve depends on the material parameters for accounting the error.
This process was implemented through Python and Abaqus® scripting interface. The files containing the Python commands run with the extension “.py”. This option became interesting for the proposed problem since it can be used to perform the following tasks:(i)Automation of repetitive tasks, such as the sequence of structural analyses via FEM;(ii)Parametric studies which modify the model, such as the attribution of different physical properties;(iii)Access to an output database, such as reading the postprocessing results.
4. Prototype Experimental Tests
Conical rubber springs are used as vibration isolators combining the function of primary damping and axle linkage in a single component. Its design is compact and light, making it easy for retrofitting. It consists of one or several conically arranged elastomeric layers bordered by adequately molded metallic components, the conical metallic molded parts should provide a separation between the various elastomeric layers as shown in Figure 2.
The conical spring application and its performance depend on the number of elastomeric layers, the cone tilt angle, and the selection of an optimum elastomeric compound. As the rubber is in both compression and shear when loaded, the conical springs allow for considerable excursion and have a progressive characteristic curve. It allows the transmission of longitudinal traction/braking loads and transverse guidance loads, as well as the vertical flexibility. It must fulfill with the requirements of the standard NF EN 13913 [52] for rubber-based parts.
The vulcanization characteristic curve should be in accordance with standard NFT 43-015 [53] or ISO 6502 [54] and the mechanical characteristics on test pieces in conformity with standard NFT 46-002 [55] or ISO 37 [56].
The maximum load of 40kN must be supported by the spring without any damage.
4.1. Static Axial Stiffness
The static axial stiffness is measured in the vertical position. First two preload cycles are applied from 0 to 40kN to stabilize the rubber properties due to stress softening. Then, the load/deflection curve is recorded during the third cycle.
The testing machine is illustrated according to Figure 3, where two power screws are responsible for the translational motion to compress the spring through a metal block controlled by a load cell.
4.2. Static Radial Stiffness
To account for static radial stiffness, an axial preload equivalent to the tare load of 20kN is first applied (Figure 4). Following, two preload cycles ranging from equivalent to -6 to 6mm are applied in the radial direction (Figure 5). Then, the load/deflection curve is recorded during the third cycle.
5. Finite Element Method for Analysis Correlation
The finite element method is a popular tool for designing elastomeric components. This numerical method is able to approximate the stress-strain behavior of rubber components based on a theoretical material model.
The least squares method is used to determine the best curve fitting for each constitutive model since the number of stress equations will be greater than the number of unknown constants. Thus, for each “” stress-strain pair which makes up the test data, the following error measure “” is minimized [57]:where is the stress value from the test data and is the theoretical stress from the fitted curve.
In many cases, the material curve is only fitted according to the uniaxial test data. This happens due to test limitations for the other states of pure deformation. Therefore, the present study carried out uniaxial tensile tests for 5 samples from the same batch according to ASTM D412 [58], in order to obtain the mean uniaxial curve as shown in Figure 6.
The correlation between uniaxial test data and theoretical curves fitted according to the classical constitutive models in study is shown in Figure 7. The biaxial and planar shear behaviors could also be estimated by each constitutive model strain energy, and they are described according to Figures 8 and 9, respectively. The material was assumed to be almost incompressible with a Poisson’s ratio equal to 0.4995 (. The relative errors should also be accounted in order to infer the constitutive models performance as in [59, 60], and they are depicted according to Figure 10 with respect to the uniaxial experimental data over the entire strain range.
It is possible to note that the uniaxial data could be well fitted by higher order phenomenological constitutive models and micromechanics. The lower order models (Mooney-Rivlin, Neo-Hooke, and Ogden-N1) tried to fit the data linearly and for this reason the relative error had only two minimum points over the strain. The accumulated relative error for each model is presented according to Table 1.
5.1. Axial Stiffness FEA Modeling
For the case of axial stiffness simulation, the conical rubber spring under investigation has been treated as an axisymmetric model due to its rotational symmetry and axial load applied in the center point.
The FEA was performed using the commercial software Abaqus®. To represent the rubber material, a 4-node bilinear quadrilateral axisymmetric element with hybrid formulation, constant pressure, and reduced integration with hourglass control (CAX4RH) has been used, and the metal parts were modeled by using a 3-node linear triangle axisymmetric element (CAX3). The hybrid formulation is suitable for strictly incompressible materials, such as rubber, and it is represented by the letter “H” at the end of element name. The capabilities of this implementation with its boundary conditions can be shown in Figure 11, where a simplified 3D model is discretized into 2978 2D axisymmetric elements.
5.2. Radial Stiffness FEA Modeling
For the case of radial stiffness simulation, the conical rubber spring has been treated as a 3D solid model. Here, it is not possible to use the axisymmetric approach since the load is not applied in the axial direction.
In the same way, as in axial stiffness, FEA was performed using the commercial software Abaqus®. The rubber was modeled by an 8-node linear hex element with hybrid formulation and constant pressure (C3D8RH); the reduced integration and hourglass control was also inputted to reduce the volumetric constraints, avoiding the overly stiff behavior due to volumetric locking. The metal parts were modeled by using a 10-node tetrahedron with an improved surface stress formulation (C3D10I).
The implementation must be performed in two analysis steps. The first step refers to the vertical preload to put the component in a prestressed condition. Afterward, in the second step the radial force is applied through an imposed displacement of 6mm.
The 3D model was discretized into 213814 solid elements. The implementation with its boundary conditions for each step is shown in Figure 12.
(a)
(b)
6. Finite Element Analysis Validation according to Prototype Experimental Data
The performance of some classical constitutive models implemented in the commercial software Abaqus (Yeoh, Mooney-Rivlin, Ogden, Neo-Hooke, and Arruda-Boyce) had been first investigated through uniaxial test data only. Next, the constitutive models with their respective material parameters were implemented in the conical rubber spring model to evaluate its behavior in the axial deflection mode. The comparison between FEA and prototype test data is shown in Figure 13. The relative errors were also accounted over the displacement according to Figure 14 in order to analyze the performance of each constitutive model, but considering now that the same rubber is subjected to multiaxial loading. Table 2 shows the accumulated relative error for each model.
During the simulations, Mooney-Rivlin and Ogden-N1 models were not able to achieve the high axial deflection performed by the prototype test. The Mooney-Rivlin model became unstable under deflections higher than approximately 37mm; on the other hand, Ogden-N1 could overcome 37mm but was only stable for deflections up to 55mm. In both cases, the analyses were aborted due to an artificial increasing in load-deflection curve at a structural level and an overprediction of the collapse load tending to volumetric locking, despite the hybrid formulation and the reduced integration adopted in the modeling process.
In spite of the good fitting of experimental simple tension data, Ogden-N3 provided very high relative errors. The load-deflection curve presented in Figure 13 showed a much higher stiffness than the experimental response. This may happen, since the planar shear and biaxial curves had been defined with a much higher stiffness for this model. Thus, if test data for some of the other deformation modes is missing, care should be taken with this model when trying to predict the rubber component behavior under multiaxial loading.
According to Figure 15 it is possible to compare the undeformed (a) and deformed (b) shapes. Maximum displacement value occurs around the metal axle, more exactly at the point where the concentrated load was applied. Although the rubber section undergoes large strains, it still keeps a reasonable shape, without excessive distortion.
(a)
(b)
Thus, by evaluating the relative error curves obtained in Figure 14, it can be observed that the Arruda-Boyce model is the one which best correlated the prototype load-deflection experimental curve. From this model, an optimization study was performed based directly on the prototype test data in order to recalibrate the material parameters, aiming to increase the accuracy degree in the component behavior that depends on the multiaxial deformation mode in which the rubber is submitted.
Table 3 shows the constants of the Arruda-Boyce model obtained by uniaxial stretching test and which was assigned as initial guesses of the design variables in the pattern search optimization process, as well as their final values after the simulation process is completed.
From the results in Figure 16, we note that the implemented optimization algorithm allows a much more accurate adjustment. However, the prototype test data will always be necessary for the material characterization, aiming as an objective function minimizing the error between the curves computed by the area difference. Figure 17 shows the relative errors over the displacement to prove the optimization algorithm efficiency in the adopted model. The accumulated relative error dropped from 76.8 to 18.1.
Through the material model obtained by the implemented optimization algorithm for axial deflection, the component mechanical behavior under precompression is verified by applying the loading in the radial deflection direction.
Figure 18 shows the force-displacement curve for the second analysis step, taking into account the precompression condition when applying the radial imposed displacement. It can be seen that the material model used in this second study is no longer capable of reproducing the results with the same accuracy as the previous loading mode. This is because the rubber has different properties and consequently different behaviors for each mode of deformation.
Figure 19 shows the component in the undeformed condition (a), in the deformed condition after precompression (b), and under the action of the radial loading (c). It is worth noting that according to Figure 16, the axial displacement should be approximately 40mm for an axial load of 20kN. Therefore, the precompression behavior was obtained with accuracy, but the radial deflection in the second step lost this accuracy, evidencing the change of rubber properties for a different multiaxial loading mode.
(a)
(b)
(c)
7. Conclusions
In order to illustrate the load-deflection relations for a rubber spring component, the influence of multiaxial loading conditions of a typical railway industry component: the conical rubber spring was inquired.
Although several techniques and assumptions have been developed to simplify the rubber spring behavior, the formulated equations became inapplicable as the geometry becomes more complex. For this reason, the using of FEA was considered the most accurate, versatile, and comprehensive method for solving the presented complex design problem.
The finite element simulations introduced in this work dealt with nearly incompressible hyperelastic material under quasistatic condition. Since Poisson’s ratio for elastomers is between 0.49 and 0.5, the elements used in FEA need to be reformulated to accommodate this high value of Poisson’s ratio. For this kind of situation, the bulk modulus presents a relatively high value, much larger than the shear modulus (μ); consequently, a small change in displacement produces high changes in pressure, so the solution based purely on the displacement field becomes very numerically sensitive. Due to this instability, the present work applied the hybrid element formulation aiming to avoid the called volumetric locking. This problem may occur when the mesh can no longer represent the strain field due to an excessively stiff behavior.
Concerning the constitutive parameters obtained by the optimization process, it is possible to understand why curve fittings on elastomers are based on the main pure deformation modes (uniaxial, biaxial, and planar shear). Each deformation mode can influence the calibration constants and modify the model fitted curve. This behavior also accounts the compressibility degree since it is related to the hydrostatic part of the strain energy potential.
When the output of the simulation is compared to the experimentally obtained load-deflection curves, the huge influence of the constitutive constants under different multiaxial loading cases can be observed. That is the reason why very complex rubber component geometries show a great challenge for numerical simulations. Even with experimental test data for the main pure deformation modes, it is not always possible to obtain the component response with a high accuracy degree.
In the future, the investigations regarding cyclic multiaxial deformations taking into account the Mullins effect and viscoelastic behavior will be executed. Other deformation modes and constitutive models including the GG-model [40], which states to be very effective during the entire range of uniaxial tension data [60] and simpler than Arruda-Boyce [40], will also be evaluated in order to compare the simulation response with prototypes test and better clarify the rubber complex behavior.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors are grateful for the financial support granted by Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Programa de Pós-Graduação em Engenharia de Estruturas (PROPEEs-UFMG), and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), under grant numbers TEC-PPM-00409-16, TEC-PPM-00444-18, and 3023762/2016-0. We also acknowledge Vibtech Industrial Ltda. for the gently provided materials.