1 INTRODUCTION

Finite element (FE) analysis is widely used to predict deformations and stress distributions, including for hyper-elastic materials. There are many material models that describe these materials, for which Hooke’s law does not apply. Moreover, proper analysis of rubber-like components requires nonlinear FE analysis tools, which differ greatly from those used for metallic components. The unique properties of rubber-like materials include the following (MSC 2010):

an ability to undergo large deformations under load, sustaining strains of up to 500% in engineering applications;

markedly non-linear load-extension behaviour; and

near incompressibility, with no appreciable change in volume in response to stress.

The successful modelling and design of rubber components depends on the selection of an appropriate strain energy function and accurate determination of the material constants in this function. However, the testing of hyper-elastic materials for the purpose of defining material models is often very difficult, with experimental requirements that differ from those for most standardized test methods. Moreover, national and international standards organizations have yet to clearly define appropriate experiments (MSC 2010).

The aim of this study was to select an appropriate material model for a polyurethane rail pad. The main role of the elastic rail pads applied in fastening systems is to reduce the dynamic effect of a rolling railway vehicle on concrete sleepers. The pads also provide electrical insulation between the rail foot and the sleeper. The stiffness of a rail pad depends on the track class and the maximum track capacity (^{Szurgott et al. 2012}). The rail pad considered in this study was designed using computer-aided engineering (CAE) systems as a result of the modernization of existing solutions and was optimized mainly in terms of ensuring the appropriate stiffness. The basic pads were made of polyurethane without any specified material parameters, and unspecified material was applied for the modernized pad as well. The conducted studies were limited to static analysis, and therefore cyclical deformation modes were omitted.

When selecting a constitutive model, the degree of complexity is an important factor, and compromise is often required. More complex models will better approximate the experimental data but may result in instability of the analysis and increased central processing unit (CPU) time. Among constitutive models, hyper-elastic Mooney-Rivlin (M-R) based material models have been applied in numerical simulations of various problems.

For example, ^{Martins et al. (2006}) considered seven material models, including some described in the current paper, i.e. Neo-Hookean, M-R and Yeoh. Material constants were obtained via the inverse method, and the fitting accuracy of the experimental stress-strain relation was evaluated. The correlation coefficient exceeded 0.995 for the M-R material model and was close to 1 for the Yeoh model. Uniaxial stretching and equibiaxial tests were proposed by ^{Sasso et al. (2008}) to fit the model and extract material parameters. Numerical simulations of the test procedures gave a useful comparison between the numerical and experimental data. Four material models from the M-R group of models were taken into consideration, and the second order M-R model was found to be the most accurate. ^{Berselli et al. (2011}) presented accurate, hyper-elastic constitutive models describing photopolymer behaviour. The authors fit the numerical data to the experimental data by applying the Neo-Hookean and 2-term M-R models. Uniaxial compression and tension tests were taken into account. The Neo-Hookean model could not adequately capture the material behaviour within the stretch range under consideration, whereas the M-R model provided acceptable fitting accuracy within 10% relative error.

With respect to rubbery materials, ^{Hamzah and Razaq (2013}) took into account the complex behaviour of these materials in theoretical and experimental works in order to attain an appropriate constitutive model of rubber components for a specific application. Uniaxial tensile tests were carried out to obtain stress-stretch curves. The 2-term M-R model adequately described the hyper-elasticity of the material under consideration, and the numerical results exhibited good agreement with the experimental data. ^{Jahani and Mahmoodzade (2014}) identified dynamic material constants of the 2-term M-R constitutive model for elastomeric components. The material constants were obtained by curve fitting of a uniaxial stress-strain curve obtained from an experimental tension test. The selection of hyper-elastic material models for a natural rubber reinforced with carbon black was proposed by ^{Shahzad et al. (2015}). Four experimental tests were performed: uniaxial, biaxial, planar shear and volumetric. The results showed that the Yeoh model was the most suitable choice, taking into account all four tests. ^{Kumar and Rao (2016}) performed experimental curve fitting using four M-R material models with 2, 3, 5 or 9 *C*
_{
ij
} . The 5-term model yielded the most accurate approximation while maintaining the correct behaviour of the. The 9-parameter M-R model was used to describe high-temperature silicone rubber in (^{Lobdell and Croop, 2016}). Stress-strain curves were obtained from the uniaxial tension/compression and planar tension tests, and a biaxial curve was selected on the basis of the uniaxial compression data. ^{Abdelsalam et al. (2017}) selected material models to appropriately describe styrene butadiene rubber based on a uniaxial tensile test only and various elastomeric material models available in MSC.Marc software. The relative error was calculated to determine the accuracy of data fitting. The third order deformation model proved to be the best because it generated the smallest relative percentage error in experimental data fitting. Finally, ^{Łagan and Liber-Kneć (2017}) applied an M-R material model for swine skin tissue modelling and obtained good matching to the experimental data from uniaxial tensile tests for the 2-parameter M-R and Yeoh models.

To obtain reliable results, the constitutive model must be selected carefully, as demonstrated by ^{Tobajas et al. (2018}). The authors considered six material models widely cited in the literature, including the Neo-Hookean, 2-term M-R and Yeoh models. To determine the quality of each of the models, the *R*
^{2} correlation coefficient between the experimental and calculated values was employed. The 2-parameter M-R model was most capable of reproducing the real behaviour of the material under consideration. Frequently, when obtaining the material constants for an M-R model, only a uniaxial data set is used. ^{Keerthiwansa et al. (2018}) described the risk of using such an approach for fitting the material model, considering the 2- and 3-parameter M-R models and the Yeoh model. The authors discussed the variations related to three basic deformation modes (uniaxial, equibiaxial, pure-shear) and stated that using only the uniaxial set of data may be insufficient.

Based on this analysis of the literature, it is clear that relatively simple M-R models-those with two or three *C*
_{
ij
} constants-dominate numerical analyses of hyper-elastic materials. Applying more complex models with a large number of *C*
_{
ij
} constants does not guarantee better correlation and higher accuracy. The selection of material constants is usually based on the correct approximation of experimental results, and other evaluation criteria are rarely considered. Moreover, uniaxial tests are often the only available experimental data. To address these shortcomings, the current paper presents a broader view of the selection process for hyper-elastic material constants.

2 MODELLING OF HYPER-ELASTIC MATERIALS

Numerical simulations were carried out using MSC.Marc/Mentat computer code. Therefore, the terminology used here corresponds to the MSC software nomenclature. Selection of the material model was limited only to the M-R material models available in the mentioned software. Other models such as Ogden (principal stretch model) or Arruda-Boyce and Gent (micro-mechanical models) were not taken into account at this stage.

The Mooney-Rivlin constitutive models are based on the strain energy function *W*. This function, disregarding terms associated with volumetric deformation and temperature changes, can be written in the following polynomial form (MSC 2010):

where *I*
_{1} and *I*
_{2} are deformation invariants and *C*
_{
ij
} are constants defined for a particular material. The deformation invariants are given as follows:

where *λ*
_{1}, *λ*
_{2}, *λ*
_{3} are the elongations of the considered element along the respective directions.

M-R constitutive models differ from each other in the degree of complexity of the strain energy function resulting from the number of *C*
_{
ij
} constants. Older versions of the MSC.Marc/Mentat software gave the possibility of using one of seven different models due to the combination of five *C*
_{
ij
} constants. In the most recent versions, the family of polynomial strain energy functions can be generalized to a complete 5th order as a combination of constants *C*
_{
ij
} , where *i*= 1, 2, 3, 4, 5 and *j*= 1, 2, 3, 4, 5. However, the current paper will show that a higher-order fit does not necessarily mean a more correct model.

The *C*
_{
ij
} constants expressed in Eq. (1) can be determined on the basis of the stress-strain curves obtained from experimental tension/compression tests. In the current study, the *C*
_{
ij
} constants were determined using the Experimental Data Fitting algorithm included in the MSC.Mentat preprocessor, which is based on the method of least squares (MSC 2010). Since uniaxial tension/ compression test results were available, it was necessary to verify the model behaviour in other states, e.g. biaxial, planar shear etc. (Table 1). Moreover, the predicted response of the model beyond the range of strains experienced during the experimental test was checked.

Assuming incompressibility of the material:

and particular elongations *λ*
_{1}, *λ*
_{2}, *λ*
_{3} as a function of elongation *λ*= 1_{ }+_{ }
*ε*, the deformation invariants can be expressed as follows:

• for the uniaxial deformation mode

• for the biaxial deformation mode

• for the planar shear deformation mode

Having defined the strain energy function W, one can determine the stress components along the direction of interest as a derivative σ(ε) = ∂W/∂ε, where σ is the engineering stress and ε is the engineering strain. For each of the M-R models, the strain energy function W was converted to a form including the Cij constants and elongation λ. The σ(ε) relation was not derived explicitly due to the complexity of the mathematical transformations. The σ(λ) relations were limited to the uniaxial deformation mode only (MSC 2010):

1. Neo-Hookean model

2. Mooney(2) model - 2-term M-R model

3. Mooney(3) model - 3-term M-R model

4. Signiorini model

5. Second order invariant model

6. Third order deformation model

7. Yeoh model

The number of components (terms) in the σ(λ) relation increases with the number of constants Cij. The polynomial form of the stress-strain relations (7)-(13) would be even more complex after expansion of the (1 + ε)n expression.

3 EXPERIMENTAL TESTS

Experimental tests were carried out to determine the engineering stress-strain characteristic in the uniaxial tension/compression state for the polyurethane material under consideration. These tests included tensile tests of 5 dumbbell specimens and compression tests of 10 cylindrical specimens. Both tests were carried out using an INSTRON 8802 universal testing machine as described in detail previously (^{Szurgott et al. 2012}), (MUT 2010). The tension tests were conducted according to the ASTM (2002) standard. The tension force and vertical displacement of the machine crosshead were registered during the test. Moreover, the longitudinal strain was measured using an extensometer attached to the specimen. Each test was conducted once the longitudinal strain reached 200%. The compression test was carried out according to the Polish standard (PCS 1954). The diameter-to-height ratio for the specimen corresponded to the values recommended in (ASTM 2001).

The results of the tension and compression tests were averaged and combined into the final curve depicted in Figure 1. This curve was used to determine the material constants in the constitutive models. The number of σeng - εeng characteristic points was limited to 53 with intervals of 0.05 for the strain. This set of points was intended for approximation using the EXPERIMENTAL DATA FIT module in the MSC.Mentat pre-processor.

4 SELECTION OF THE MATERIAL MODEL

4.1 Experimental data fitting

The Cij constants were determined using the EXPERIMENTAL DATA FIT algorithm included in the MSC.Mentat pre-processor, which is based on least squares fitting (MSC 2010). The least squares error is given by either

or

where ΔR and ΔA are the relative and absolute least squares errors, respectively (the relative error is the default in MSC.Mentat), n is the total number of data points, σc is the calculated stress, and σm is the measured engineering stress (taken from the experimental data).

Figures 2a to
8a (^{Jarzębski 2015}) show the σeng - εeng curves obtained as a result of data fitting for each material model. The results relate to the uniaxial mode only since these data were available from the experimental test. The graphs are presented for a deliberately increased range of strain to permit an assessment of the behaviour of the model beyond the range of strain received during the experimental test.

The material Cij constants in the constitutive models under consideration determined using the EXPERIMENTAL DATA FIT algorithm are provided in Tables 2 and 3 for the relative and absolute error criteria, respectively. As the material Cij constants and, in turn, the complexity of the constitutive model increase, the relative and absolute errors decrease. The approximation of the experimental data is more accurate, but the behaviour of the model is not appropriate in some cases. Several material models are unstable since they do not meet ^{Drucker’s stability postulate (Drucker 1957}):

where dσ is the stress increment tensor associated with the strain increment tensor dε through the constitutive relation. Relation (16) states that the incremental internal energy of a material can only increase. Stability was assessed based on the waveforms from Figures 2a- 8a.

Material model | C10 | C01 | C11 | C20 | C30 | Error |
---|---|---|---|---|---|---|

Neo-Hookean | 2.38916 | - | - | - | - | 5.66059 |

Mooney(2) | 1.52181 | 1.51389 | - | - | - | 2.92835 |

Mooney(3) | 2.21179 | 1.40132 | -0.204704 | - | - | 1.97711 |

Signiorini | 2.26824 | 1.05553 | - | -0.0785725 | - | 2.34187 |

2nd order invariant | 1.70769 | 2.21650 | -0.537369 | 0.171096 | - | 1.70834 |

3rd order deformation | 2.03154 | 1.90897 | -0.420575 | 0.0371954 | 0.00786259 | 1.68343 |

Yeoh | 3.76949 | - | - | -0.415479 | 0.0289463 | 2.49475 |

Material model | C10 | C01 | C11 | C20 | C30 | Error |
---|---|---|---|---|---|---|

Neo-Hookean | 2.72904 | - | - | - | - | 703.676 |

Mooney(2) | 1.69702 | 0.947423 | - | - | - | 81.0011 |

Mooney(3) | 2.02289 | 1.18392 | -0.124802 | - | - | 33.3222 |

Signiorini | 2.30745 | 0.768876 | - | -0.0676634 | - | 44.4381 |

2nd order invariant | 1.98840 | 1.21917 | -0.134267 | 0.0065626 | - | 33.2525 |

3rd order deformation | 2.60906 | 0.856871 | -0.023915 | -0.206314 | 0.0145353 | 25.3836 |

Yeoh | 4.17015 | - | - | -0.320778 | 0.0136151 | 286.421 |

The approximations of the experimental values by the different material models according to the model stability criterion are summarized in Table 4. Only two material models-the Neo-Hookean and 2-term Mooney-are stable in the uniaxial tension/compression state regardless of the accepted error criterion. The next two models-the Yeoh and 3rd order deformation-are stable during compression and tension only when one error criterion is used. The remaining models are generally unstable in or beyond the range of strain received during the experimental test.

Material model | Relative error criterion | Absolute error criterion | ||
---|---|---|---|---|

Compression | Tension | Compression | Tension | |

Neo-Hookean | stable | stable | stable | stable |

Mooney(2) | stable | stable | stable | stable |

Mooney(3) | unstable | unstable | unstable | unstable |

Signiorini | stable | unstable | stable | unstable |

2nd order invariant | unstable | stable | unstable | stable |

3rd order deformation | unstable | stable | stable | stable |

Yeoh | stable | stable | stable | unstable |

4.2 Approximation error analysis

Error analysis is performed for each model to more clearly evaluate the match between the calculated and experimental values, as depicted in Figures 2b- 8b. The normalized error (%) is defined as follows:

where σm(ε) and σc(ε) are the measured and calculated engineering stress, respectively. The normalized error for the strain ε = 0 is omitted in the graphs. Each graph includes curves for the relative and absolute error criteria. The values correspond to the engineering stress-engineering strain (σeng - εeng)curves from Figures 2a- 8a.

The values of the normalized errors are averaged according to the following relation:

for compression and tension independently (Figure 9a and 9b, respectively). A better approximation is reached for tension. The average of the normalized errors for the whole range of strain (Figure 9c) is proposed to qualitatively compare the approximations obtained using the different material models and assess the best compatibility between the experimental and calculated results.

Based on the comparison depicted in Figure 9, a better approximation is obtained using the default relative error criterion. The weighted average does not exceed 3.7% except for the Neo-Hookean model (9%). Once the absolute error criterion is applied, the average of the normalized errors increases and varies between -12% and 11.4%. Therefore, the 3-term, Signiorini, 2nd and 3rd order material models approximate the experimental values most accurately. However, some models have been eliminated due to instability.

4.3 Correlation between measured and calculated values

Another determinant used in this matching process was the Pearson correlation coefficient between the measured and calculated values of engineering stress. It can be defined as follows (^{Cohen 1998}):

where σm, σc are the measured and calculated engineering stress, respectively, and
*σ*
_{
m
} , *σ*
_{
c
} . As mentioned before, the number of *σ*(ε) characteristic points was limited to 53.

It can be assumed that the correlation between the two considered variables is large if the correlation coefficient ranges between 0.5 and 1.0 (^{Cohen 1998}). In the current study, the correlation is very large since the obtained values of the correlation coefficient (Table 5) are close to unity. However, the correlation for the Neo-Hookean model clearly stands out from the rest of the models as the lowest. The highest correlation is observed for the 3rd order deformation model when the absolute error criterion is taken into account.

4.4 Verification of material models in states other than uniaxial tension/compression

Only the uniaxial tension/compression test results were taken into consideration above, and the model behaviour in other states should be assessed. The MSC.Mentat pre-processor permits the verification of the material model behaviour under biaxial and planar shear (Figure 1) for the assumed *C*
_{
ij
} constants. The waveforms shown in Figures 10-
16 (^{Jarzębski 2015}) should be considered only in qualitative terms because neither biaxial nor planar shear tests were conducted. Nevertheless, the waveforms permit an assessment of the model behaviour in such states. The non-physical behaviour of the higher-order models with a larger number of *C*
_{
ij
} constants confirms the previous analyses. Table 6 provides information on the computational stability of the material model-understood as a lack of non-physical behaviour-in all three deformation modes.

The results presented in Table 6 clearly show that the simplest models-the Neo-Hookean and the 2-term M-R-most accurately describe the polyurethane under consideration. Furthermore, there was no non-physical behaviour of the 3rd order deformation model with *C*
_{
ij
} constants determined on the basis of the smallest absolute error criterion. The above-mentioned material models can thus be considered stable in all three deformation modes during both tension and compression.

5 RESULTS OF FINITE ELEMENT ANALYSES

The numerical analyses in the current study were limited only to the static compression test of the pad. An FE model of a quarter of the pad was developed due to the symmetry of the problem. The pad was modelled as a defor mable body, whereas the top surface of the sleeper and the bottom face of the rail were treated as rigid bodies (Figure 17). The main aim of the FE analyses was to estimate the CPU time depending on the applied material model and to verify the stability of the FE model and convergence of the process. The use of the CPU time criterion for such simple analyses allows the qualitative assessment of the effectiveness of a given material model. It can be assumed that the relationship between the CPU time and the applied material model will be similar in more complex analyses.

Material model | Uniaxial | Biaxial | Planar shear | Overall rating | ||||
---|---|---|---|---|---|---|---|---|

-σ(ε) |
+σ(ε) |
-σ(ε) |
+σ(ε) |
-σ(ε) |
+σ(ε) |
|||

Relative error | Neo-Hookean | positive | positive | positive | positive | positive | positive | positive |

Mooney(2) | positive | positive | positive | positive | positive | positive | positive | |

Mooney(3) | negative | negative | negative | negative | negative | negative | negative | |

Signiorini | positive | negative | negative | negative | positive | positive | negative | |

2nd order invariant | negative | positive | positive | negative | negative | negative | negative | |

3rd order deformation | negative | positive | positive | negative | negative | negative | negative | |

Yeoh | positive | positive | positive | negative | positive | negative | negative | |

Absolute error | Neo-Hookean | positive | positive | positive | positive | positive | positive | positive |

Mooney(2) | positive | positive | positive | positive | positive | positive | positive | |

Mooney(3) | negative | negative | negative | negative | negative | negative | negative | |

Signiorini | positive | negative | negative | negative | negative | negative | negative | |

2nd order invariant | negative | positive | negative | negative | negative | negative | negative | |

3rd order deformation | positive | positive | positive | positive | positive | positive | positive | |

Yeoh | positive | negative | positive | negative | positive | negative | negative |

The FE model consisted of 10,368 solid TYPE 84 eight-node isoparametric elements with an additional ninth node for the pressure (MSC 2011), as recommended for the analysis of hyper-elastic materials. Compression of the pad was carried out by vertical displacement of the rail surface (25% of the pad thickness). The displacement was increased from 0 to 2 mm with a constant step. Because the main aim of the current paper was to identify the most accurate material model, the results of the FE analyses are limited to the contours of equivalent stress depicted in Figures 18-
24 (^{Jarzębski 2015}). Some cases led to non-convergence in the last calculation step.

The results of the maximum equivalent stress in the last three calculation steps are summarized in Figure 25. In general, the stress values increase in subsequent steps except for the Mooney(3) and 2nd order invariant material models with the relative error criterion.

The results obtained for the relative error criterion (Figure 25a) vary more greatly than those for the absolute error criterion (Figure 25b); the latter yields similar values regardless of the applied material model. Based on the results presented in Figure 25, the absolute error criterion seems to be a better criterion for selecting *C*
_{
ij
} constants in material models due to the full convergence in each of the seven cases and the similar values of equivalent stress obtained.

The CPU time was the last criterion for material model selection. All FE analyses were carried out on the same average-class computer. The shortest calculation time, 397 s, was obtained with the 2nd order invariant material model, whereas the longest analysis-excluding non-convergent cases-was obtained with the Neo-Hookean material model, with a time of 626 s. In both cases, the absolute error criterion for selecting *C*
_{
ij
} constants was assumed.

6. CONCLUSIONS

The selection of a material model on the basis of the most accurate approximation of the experimental results is an important and very difficult problem. Ultimately, more complex material models are not necessarily preferable. Increasing the number of criteria for model selection allows researchers to make the right choice.

The current study takes into account the following selection criteria:

the most accurate approximation of the experimental results in the uniaxial tension/compression state, including approximation error analysis and the correlation between the calculated and experimental values;

the complexity of the material model, including the number of

*C*_{ ij }constants and the complexity of the stress-strain relation;verification of the material models in states other than uniaxial tension/compression;

the results obtained from preliminary FE analyses of a simple component made of the hyper-elastic material, including the convergence of these analyses and CPU time.

Taking into account all of these factors allows appropriate selection of the material model. Based on the proposed comprehensive selection of a material model, the following conclusions can be drawn. Material models with a larger number of *C*
_{
ij
} constants approximate the stress-strain relation more accurately. However, they are less computationally stable and may lead to non-convergence in analysis. Application of the absolute error criterion during the determination of the *C*
_{
ij
} constants seems to be more reasonable because similar values are obtained regardless of the applied material model. A large number of *C*
_{
ij
} constants does not increase the CPU time, which remains at a similar level for most material models.

Observations during analyses should not be generalized to the whole group of hyper-elastic materials since they relate only to the specific material under consideration. However such observations may provide a foundation for researchers involved in the analysis of these types of engineering materials. In conclusion, the material model describing the polyurethane under consideration in the most appropriate manner-satisfying the above conditions and guaranteeing the stability of the model and convergence of the analysis-is the Mooney(2) model with only two constants, *C*
_{10} and *C*
_{01}.