A. Comparison of the rubber models

The characteristics of rubber can be decomposed into time-independent and time-dependent behavior. Time-independent behavior considering finite deformation is described by isotropic hyperelastic models such as the neo-Hookean, the Mooney–Rivlin, the Ogden, the Yeoh, the Arruda–Boyce, and the Van der Waals model.Reference Yeoh^4, Reference Mooney^6–14 The models of rubber elasticity such as the statistical mechanics and phenomenological continuum mechanics theories have been developed by some different bases.Reference Treloar^8, Reference Boyce and Arruda¹³ The statistical theory assumes randomly oriented long molecular chains of rubber to be idealized structures, whereas the phenomenological theory has been developed to conveniently describe the properties based on the continuum mechanics. As described by Treloar,Reference Treloar⁸ the phenomenological theories such as the Mooney–Rivlin, the Ogden, and the Yeoh models are based not on molecular concepts but on purely mathematical reasoning. In this study, we only deal with the phenomenological theory because it is convenient to describe the behavior of rubber materials in the continuum mechanics, and most studies including experimental works to evaluate the stress–strain properties of rubber use the phenomenological properties.

The invariant-based constitutive equations of rubber material are expressed by the relationship between strain energy density and strain invariants as follows:

(1)$$W\, = f(I_1 ,I_2 ,I_3 ),$$

where W is the strain energy density (or strain energy potential), which means the stored strain energy per unit volume, and I ₁, I ₂, and I ₃ are the three invariants of the Green deformation tensor given in terms of the principal stretches (or extension ratios) λ₁, λ₂, and λ₃as follows:

(2)$$I_1 \, = \,{\rm{\lambda }}_1^2 \, + \,{\rm{\lambda }}_2^2 \, + \,{\rm{\lambda }}_3^2 ,$$

(3)$$I_2 \, = \,{\rm{\lambda }}_1^2 {\rm{\lambda }}_2^2 \, + \,{\rm{\lambda }}_2^2 {\rm{\lambda }}_3^2 \, + \,{\rm{\lambda }}_3^2 {\rm{\lambda }}_1^2 ,$$

(4)$$I_3 \, = \,{\rm{\lambda }}_1^2 {\rm{\lambda }}_2^2 {\rm{\lambda }}_3^2 .$$

The strain energy density function in terms of the invariants is written as

(5)$$W = \mathop\sum\limits_{i + j + k = 0}^\infty {C_{ijk} (I_1 - 3)^i (I_2 - 3)^j (I_3 - 3)^k } .$$

Since rubber material is almost incompressible (I ₃ = 1), Eq. (5) can be reduced to

(6)$$W = \mathop\sum\limits_{i + j = 0}^\infty {C_{ij} (I_1 - 3)^i (I_2 - 3)^{\,j} } .$$

Several reduced models of the polynomial form Eq. (6) have been suggested, such as the neo-Hookean model, the Mooney–Rivlin model, and the Yeoh model. The Mooney–Rivlin form is a linear function of the invariants I ₁ and I ₂, i.e.,

(7)$$W\, = \,C_{10} (I_1 - 3)\, + \,C_{01} (I_2 - 3).$$

The Mooney–Rivlin model is widely used because of its simplicity, but it is effective only in relatively small strain. For larger deformation, the Ogden modelReference Ogden¹¹ is more appropriate than the Mooney–Rivlin form and can be expressed as

(8)$$W\, = \,\mathop\sum\limits_{i = 1}^n {{{{\rm{\mu }}_i } \over {{\rm{\alpha }}_i }}({\rm{\lambda }}_1^{{\rm{\alpha }}_i } + {\rm{\lambda }}_2^{{\rm{\alpha }}_i } + {\rm{\lambda }}_3^{{\rm{\alpha }}_i } - 3)} ,$$

where n, μ, and α are material parameters. It is noteworthy that the neo-Hookean and the Mooney–Rivlin models are special cases of the Ogden form.Reference Ali, Hosseini and Sahari¹⁵ The Ogden model, however, requires lots of fitting parameters, so that this method is effective only if sufficient experimental data are can be obtained. The Yeoh modelReference Yeoh^3, Reference Yeoh⁴ and its extended modelReference Yeoh and Fleming^5, Reference Beda¹⁶ have been thus suggested to be applicable for larger deformation with simple experimental data. The constitutive equation is represented by

(9)$$W\, = \,C_1 (I_1 - 3)\, + C_2 (I_1 - 3)^2 \, + \,C_3 (I_1 - 3)^3 .$$

To compare the constitutive models, a uniaxial tensile test and a shear test have been performed for a vulcanized rubber. The Ogden, the Mooney–Rivlin, and the Yeoh coefficients are obtained from the fitting of the tensile test data. By using the coefficients, the nominal tensile stress–strain curves are generated and compared with the test data as shown in Fig. 1. Compared to the Mooney–Rivlin form, the Yeoh model is more suitable for describing the behavior of rubber materials. Moreover, although the Yeoh model is simpler than the Ogden model, it can follow well the experimental data including large deformation. Most of all, whereas the Mooney–Rivlin and the Ogden models often inaccurately predict the behavior for different modes by a single experiment, the Yeoh model can predict well the behavior in different deformation modes from data gained in just one deformation mode.Reference Ali, Hosseini and Sahari¹⁵ The shear stress–shear strain curves are then generated by the Yeoh and the Mooney–Rivlin coefficients, which are obtained from tensile test data, to compare them with experimental shear test data as shown in Fig. 2. The Mooney–Rivlin model shows large deviation from experimental data, but the Yeoh model can predict well the curve up to the shear strain is around 1.

FIG. 1.

Nominal stress–strain data predicted by the Ogden, the Mooney–Rivlin, and the Yeoh models of which coefficients are obtained from a uniaxial tensile test.

FIG. 2.

Nominal shear stress–shear strain data predicted by the Yeoh and the Mooney–Rivlin models of which coefficients are obtained from a uniaxial tensile test.

FIG. 9.

(a) W versus r/D and (b) I ₁ versus r/D for various values of l/D where l is the distance from specimen surface to data acquisition depth before deformation and r is the projected distance from axisymmetric center after deformation.

As we assume incompressible materials, i.e., I ₃ = 1, the strain energy density for an isotopic incompressible elastomer can be expressed as a function of I ₁ and I ₂ [W = f (I ₁, I ₂)]. In a strict sense, ∂W/∂I ₂ has finite values.Reference McKenna and Zapas^17, Reference McKenna and Zapas¹⁸ However, prior studies showed that ∂W/∂I ₂ is sufficiently smaller than ∂W/∂I ₁,Reference Yeoh and Fleming^5, Reference Gregory¹⁹ and therefore, ∂W/∂I ₂ is often assumed to be 0.Reference Yeoh and Fleming^5, Reference Beda^16, Reference Qi and Boyce²⁰ If ∂W/∂I ₂ is assumed to be 0, W can be represented as a function of only I ₁. Coefficients for rubber materials can then be obtained from one single test, e.g., uniaxial tension test.Reference Ali, Hosseini and Sahari¹⁵ Therefore, a neo-Hookean model, a Yeoh model, and an Arruda–Boyce model are developed by neglecting I ₂. We also use a strain energy density function based on I ₁ [Eq. (9)]. In this study, the Yeoh model is thus selected as the constitutive equation.

B. The Yeoh model

The Yeoh model describes the behavior of rubber materials by W, I ₁, and three coefficients C ₁, C ₂, and C ₃. Hence, it is very important to investigate the effect of each coefficient on the stress–strain curve.

For a uniaxial tensile/compression test of incompressible rubber material, I ₃ is equal to 1 (${\rm{\lambda }}_1^2 {\rm{\lambda }}_2^2 {\rm{\lambda }}_3^2 $ = 1). In addition, the principal stretches perpendicular to the tensile/compression direction are identical (λ₂ = λ₃), the following equations can be easily obtained:

(10)$${\rm{\lambda }}_2^2 = {\rm{\lambda }}_3^2 = 1/{\rm{\lambda }}_1 ,$$

(11)$$I_1 = {\rm{\lambda }}_1^2 + 2/{\rm{\lambda }}_1 ;\,\,\,{\rm{\lambda }}_1 \, = 1 + {\rm{\varepsilon }}_1 .$$

For uniaxial tensile/compression tests, the nominal stress, σ₁, can be written as follows:

(12)$$W = \int {{\rm{\sigma }}_1 } \,d{\rm{\varepsilon }}_1 \, \leftrightarrow \,{\rm{\sigma }}_1 = {{\partial W} \over {\partial {\rm{\varepsilon }}_1 }} = {{\partial W} \over {\partial I_1 }}{{\partial I_1 } \over {\partial {\rm{\lambda }}_1 }}{{\partial {\rm{\lambda }}_1 } \over {\partial {\rm{\varepsilon }}_1 }}.$$

The first derivative of Eq. (9) with respect to I ₁is

(13)$${{\partial W} \over {\partial I_1 }} = C_1 + 2C_2 (I_1 - 3) + 3C_3 (I_1 - 3)^2 .$$

Substituting this equation [Eq. 13] and Eqs. (11) and (12), we obtain

(14)$${\rm{\sigma }}_1 = {{\partial W} \over {\partial I_1 }}{{\partial I_1 } \over {\partial {\rm{\lambda }}}} = [C_1 + 2C_2 (I_1 - 3) + 3C_3 (I_1 - 3)^2 ]\left[ {2\left( {{\rm{\lambda }} - {1 \over {{\rm{\lambda }}^2 }}} \right)} \right] \,\,\,\, = \left[ {C_1 + 2C_2 \left( {{\rm{\lambda }}^2 + {2 \over {\rm{\lambda }}} - 3} \right) + 3C_3 \left( {{\rm{\lambda }}^2 + {2 \over {\rm{\lambda }}} - 3} \right)^2 } \right]\left[ {2\left( {{\rm{\lambda }} - {1 \over {{\rm{\lambda }}^2 }}} \right)} \right], } $$

where we simply set ${\rm{\sigma }}_1 = {\rm{\sigma }},\,\,{\rm{\varepsilon }}_1 = {\rm{\varepsilon }},$ and ${\rm{\lambda }}_1 = {\rm{\lambda }}$. Since λ = 1 + ε, Eq. (14) can be rewritten as follows:

(15)$$\eqalignb{ {\rm{\sigma }} =\, [ {C_1 + 2C_2 ( {(1 + {\rm{\varepsilon }})^2 + {2 \over {(1 + {\rm{\varepsilon }})}} - 3} ) \cr+ 3C_3 ( {(1 + {\rm{\varepsilon }})^2 + {2 \over {(1 + {\rm{\varepsilon }})}} - 3} )^{2} } ] \cr \!\times [ {2 ( {(1 + {\rm{\varepsilon }}) - {1 \over {(1 + {\rm{\varepsilon }})^2 }}} )} ]. \cr}$$

The Yeoh model is defined by the relationship between the strain energy density and the first invariant, i.e., W–I ₁. It can be expressed in terms of nominal stress (σ) and nominal strain (ε) as shown in Eq. (15).

Stress–strain curves obtained by Eq. (15) are shown in Fig. 3 for three Yeoh coefficients, C ₁, C ₂, and C ₃. The bold line in each figure indicates the stress–strain curve obtained by tentatively selected reference properties (C ₁ = 1.5 MPa, C ₂ = −0.6 MPa, and C ₃ = 0.3 MPa). One of the three coefficients is varied, while the other two values are fixed as the reference values. In the initial stage of deformation, C ₁ dominates the stress–strain curve. As the strain is increased, C ₂ becomes a second influential variable, followed by C ₃. Consequently, to accurately acquire C ₁ and C ₂, the maximum strain should be at least 0.4 (Fig. 3), whereas it should be at least 1.0 to accurately acquire C ₃.

FIG. 3.

Stress–strain curves for various values of three coefficients of the Yeoh model.

The behavior of rubber, which is usually explained by a phenomenological constitutive equation, is expressed as the relationship between the strain energy density and the strain invariants (or the principal stretches), and the constitutive equation should express well various deformation modes. If the tensile and compression tests are performed ideally, the material constants obtained from tensile and compression tests with an equivalently softened specimen must be almost identical. Figure 4(a) shows the stress–strain curves of a rubber material, SMR-CV60, obtained from the tensile and compression tests in a reference.Reference Kim, Kim, Kim, Woo and Lee²¹ The figure indicates that the tensile and compression data can be fitted with a set of the Yeoh coefficients, which are obtained from the compression data only. As described in Sec. II. A, the behavior in different deformation modes can be accurately predicted by the Yeoh model from data with one deformation mode. Figure 4(b) shows the true stress–strain curves of the material. The true stress–strain curve is not an odd function, whereas the true stress–strain curves of many metallic materials usually assume to be point symmetry when the Bauschinger effect is negligible. It means the tensile and compression behaviors of rubber material are not equivalent, but the Yeoh model can express these nonsymmetric behaviors with a single set of coefficients. Consequently, considering the mixed tensile, compression, and shear stress fields in the deformed indentation specimen, the Yeoh model is adequate to compare indentation properties with tensile properties.

FIG. 4.

(a) Nominal stress–strain curves of a natural rubber, SMR-CV60, obtained from tensile and compression tests. The stress–strain curve of the Yeoh model is predicted by the Yeoh coefficients of the compression data. (b) The true stress–strain curves of the tensile and compression test data. The point symmetry of the compression curve is plotted to compare the tensile and compression curves.

C. The ranges of the Yeoh coefficients

To set up the range of the three coefficients of the Yeoh model, C ₁, C ₂, and C ₂, we carry out tensile tests for three kinds of rubbers, nitrile butadiene rubber (NBR), chloroprene rubber (CR), and silicon rubber (SIU) according to ISO37:2005. Around 7–10 cycles of loading–unloading are performed before the actual test to preclude Mullins softening effectReference Mullins²² (Fig. 5). We study then the dependence of the load–displacement curve on the strain rate using four different strain rates, 0.0067, 0.013, 0.13, and 0.8 s⁻¹. From Fig. 6, it is confirmed that the load–displacement curves converge when the tensile strain rate is lower than 0.013/s. Hence, we define 0.013/s as the quasistatic strain rate. The corresponding material properties are regarded as quasistatic material properties. By using Eqs. (11) and (12), the relationship between W and I ₁ is obtained from the tensile test data, and the Yeoh coefficients are acquired from fitting the W–I ₁ curve. The ranges of the three coefficients are confined in C ₁ = 0.6–2.2 MPa, C ₂ = −1.1–0.0 MPa, and C ₃ = 0.1–0.6 MPa. Table I shows the Yeoh coefficients of the three rubber materials. In Sec. V, we perform instrumented indentation tests with these three rubber materials. These material properties obtained from tensile tests are used as reference values.

FIG. 5.

The effect of Mullins softening on uniaxial tensile tests for three kinds of rubber materials, NBR, CR, and SIU.

FIG. 6.

The effect of strain rates on tensile tests for three kinds of rubber materials, NBR, CR, and SIU.

TABLE I.

Yeoh coefficients of three kinds of rubber materials, NBR, CR, and SIU, obtained from tensile tests.

A. FE modeling for spherical indentation

Figure 7 shows an axisymmetric finite element model. A commercial code ABAQUS¹⁴ is used to simulate spherical indentation tests for rubber materials. Since load and geometry are axisymmetric and the rubber shows incompressible behavior, four-node axisymmetric hybrid elements (CAX4H) can be used in analysis.¹⁴ To minimize the constraint caused by the bilinear shape function and to detect sufficiently small deformation, fine elements with size 1% of indenter diameter are placed around the contact surface of specimen. As the elastic modulus of indenter is relatively high compared to that of rubber, the indenter is assumed to be rigid; therefore, rigid surface¹⁴ is used. The finite element (FE) model consists of about 12,300 elements and 12,500 nodes. Contact surfaces¹⁴ are placed on both contact surfaces of specimen and indenter. Roller boundary conditions are imposed on the bottom and the axisymmetric axis of the specimen. The bottom of the indenter is a hemispherical rigid body with diameter D = 1 mm, and the upper part of the indenter is a cylindrical rigid body with radius R = 0.5 mm (Fig. 7). It should be noted that our FE model can be applied to any indenter diameter because all numerical formulas normalized by the indenter diameter make the method universally applicable irrespective of its size.Reference Lee, Kim, Lee and Nahm^2, Reference Lee, Lee and Pharr²³ To analyze and quantify the effects of indentation variables on the load–displacement curve, we set up the ranges for the main variables, i.e., the material constants of the Yeoh model, as described in Sec. II. C. To precisely extract the property C ₃ of the Yeoh form, the maximum indentation depth, h _max, is set as h _max = D, which is deeper than that of our previous work.Reference Lee, Kim, Lee and Nahm² The friction coefficient used in this study is explained in Sec. III. B.

FIG. 7.

The axisymmetric finite element model for the spherical indentation.

To investigate the size effects of the indentation specimen on the load–displacement curve for a given indenter diameter, we perform FE analyses with selected material properties (C ₁ = 1.0 MPa, C ₂ = −0.4 MPa, and C ₃ = 0.4 MPa) in the defined property range. The height (H) and radius (W) of the cylindrical specimen are varied. In total, nine cases are applied in finite element analysis (FEA) with three ratios for H/D and W/D (24, 48, and 96). From this, we conclude that H/D = 48 and W/D = 48 are the optimal specimen size because the effect of specimen size is saturated here. The z-directional displacement at the end of the top surface can be regarded as 0. The effects of the specimen size on the indentation results with some other combinations of the coefficients are checked, and we conclude that H/D = 48 and W/D = 48 are sufficient to use as saturated dimensions in this study. It should be noted that in indentation study, 10% rule (the ratio of the indentation depth to the specimen size) is conventionally used, but it is about 2% in this study because we set h _max = D, which is much smaller than the conventional value.

B. Friction coefficient

Friction forces occur at an interface where tangential and normal forces are applied. To define friction, the Coulomb friction model is generally used. The friction coefficient, f, can change depending on the state of lubrication, contamination, and materials. Kim et al.Reference Kim, Kim, Na and Lee²⁴ have shown that the friction coefficient between structural steel and vulcanized natural rubber is affected by normal load, deformation rate, state of lubrication, and the thickness of the rubber plate. In their study, the friction coefficient was estimated at 1.6–3.9 in the case of unlubricated surfaces and 0.1–0.2 for grease lubrication. It is therefore required to quantitatively evaluate the effect of the friction coefficient on the indentation analysis. Thus, we perform FEA for various friction coefficients between indenter and rubber specimen. Figure 8(a) shows the indentation load–displacement curves for a wide variety of friction coefficients when the Yeoh coefficients are assumed to be C ₁ = 1.0 MPa, C ₂ = −0.4 MPa, and C ₃ = 0.4 MPa. Here, h is the indentation depth, and P is the indentation load. Lee et al.Reference Lee, Kim, Lee and Nahm² used the values of 0.09–0.13 for the friction coefficient assuming well-lubricated conditions; however, the assumption can be impractical. Figure 8(b) shows an example where two sets of material properties and friction coefficients can show almost the same P–h curve; this means that the previous methodReference Lee, Kim, Lee and Nahm² could produce wrong material properties for unlubricated conditions. Hence, we concluded that there is no reason to consider lubrication between an indenter and a specimen, which makes this indentation testing simpler. It is therefore reasonable to suggest an updated indentation approach with properly assumed values of the normalized indentation depth h/D and the friction coefficient. As shown in Fig. 8(a), the indentation load–displacement curves converge as the friction coefficient increases, and when it is greater than 0.5, they are not affected by the friction coefficient any more. The effects of the friction on the indentation responses with other combinations of the Yeoh coefficients are tested, and 0.5 is suitable to use as a sufficiently converged friction coefficient. Hence, we select 0.5 for the friction coefficient between indenter and rubber specimen and perform static FEA by using this value. Therefore, the static indentation approach to evaluate rubber material properties can be assumed to be independent of the frictional coefficient.

FIG. 8.

(a) Indentation load–displacement curves of a rubber material (C ₁ = 1.0 MPa, C ₂ = −0.4 MPa, and C ₃ = 0.4 MPa) with various friction coefficients. (b) Two set of material properties and friction coefficients can generate almost the same P–h curve.

A. The optimal data acquisition point

The distributions of nominal stress–strain (and W–I ₁) data under the indenter are definitely complicated. To estimate rubber properties from indentation tests, it is therefore necessary to define the representative stress–strain or the representative strain energy density–first invariant corresponding to the indentation load–displacement. The detailed basis of the selection is well described in our preceding work.Reference Lee, Lee and Pharr^23, Reference Lee, Kim and Lee²⁵ Applying these indentation methodologies to the present study, we obtain the relation between (representative) W–I ₁ data at a selected data acquisition point. To find the optimal data acquisition point, we investigate the distributions of W [Fig. 9(a)] and I ₁ [Fig. 9(b)] along radial direction (r/D) for four l/D values; here, l is the distance from specimen surface to data acquisition depth before deformation, and r is the projected distance from axisymmetric center after deformation. When h/D = 1, the strain energy density and the first principal invariant I ₁ reach a maximum at l/D = 35%. In this study, we select r/D = 0 and l/D = 35% as the optimal data acquisition point.

B. Numerical formulas for spherical indentation

The constitutive equation for the Yeoh modelReference Yeoh^3, Reference Yeoh⁴ is expressed in terms of the strain energy density W and the first invariant I ₁. Lee et al.Reference Lee, Kim, Lee and Nahm² adopted nondimensional functions such as Eqs. (1) and (2) to map load–displacement onto W–I ₁:

(16)$${\rm{\psi }}\, \equiv \,WhD/P\quad,$$

(17)$$\phi \, \equiv \,I_1\quad .$$

They indented materials up to 40% of indenter diameter (h/D = 0.4), which is the maximum depth where the effects of friction and C ₃ on load–displacement curve can be neglected. For commercial vulcanized rubber, the next equations for ψ and $\phi $ as a function of C ₁, C ₂, and indentation depth were obtained through FEA when C ₃= 0.005 MPa as a representative constant value:

(18)$${\rm{\psi }} = \mathop\sum\limits_{i = 0}^3 {f_i ^{\rm{\psi }} (C_1 ,\,C_2 )\left( {{h \over D}} \right)^{i} \quad,} $$

(19)$$\phi = \mathop\sum\limits_{i = 0}^3 {f_i ^\phi (C_1 ,\,C_2 )\left( {{h \over D}} \right)^{i} }\quad .$$

Additional indentation up to h/D = 0.8 was then performed to evaluate C ₃. Equation (20) represents the relationship between C ₃ and P.

(20)$$C_3 = f^C (P) = \mathop\sum\limits_{i = 0}^2 {\left[ {\mathop\sum\limits_{j = 0}^2 {\left( {\mathop\sum\limits_{k = 0}^2 {\eta _{ijk} C_1^k } } \right)C_2^j } } \right]} \left( {{P \over {D^2 }}} \right)^{i}\quad ,$$

where η_ijk are the coefficients of the polynomial function. This method is however likely to have a large error for very large deformation severely affected by the coefficient C ₃. This is because to estimate large deformation well, the simultaneous variation of the three coefficients, C ₁, C ₂, and C ₃, should be considered. In the present work, a new method to overcome this limitation is suggested.

From Fig. 3, we can see that the effects of the coefficients C ₂ and C ₃ of Yeoh model in the initial stage of stress–strain curves are negligible. Hence, we fix C ₁ and vary C ₂ and C ₃ one by one to identify the range in which C ₁ independently affects the load–depth curve. We first fix C ₁ and C ₃ and vary C ₂, as shown in Fig. 10(a), and then fix C ₁ and C ₂ and vary C ₃, as shown in Fig. 10(b). Figure 10 makes it obvious that indentation load–displacement curves at h/D < 0.1 are determined by C ₁ only irrespective of C ₂ and C ₃ values. We can thus evaluate C ₁ first for h/D < 0.1 and then obtain C ₂ and C ₃ from triple fitting functions where h > 0.1. The relationship between indentation load $P_{h/D\,\, = \,\,0.1} $ and C ₁ is shown in Fig. 11. Based on the fitting of data in the figure, we can write

(21)$$\eqalignb{ C_1 = f(P_{h/D\,\, = \,\,0.1} ) = H_0 + H_1 P_{h/D\,\, = \,\,0.1} /D\,^2 \, \cr (H_0 ,\,\,H_1 ) = (0.01663\,\,{\rm{MPa}},\,\,4.3758)\quad. \cr} $$

FIG. 10.

Indentation P–h curves for the various Yeoh coefficients of (a) C ₂ and (b) C ₃. In each case, other coefficients are fixed as reference values.

FIG. 11.

Comparison of C ₁–P_{h/D = 0.1}/D ² curve and its regression curve (h/D = 0.1).

To obtain C ₂ and C ₃, by the same approach as Lee et al.,Reference Lee, Kim, Lee and Nahm² we newly define normalized functions and match each load–displacement data to its equivalent W–I ₁ data. In the present work, contrary to Eqs. (16) and (17), we define new normalized mapping functions as follows:

(22)$${\rm{\psi }} = {{Wh^3 } \over {PD(I_1 - 3)}}\quad,$$

(23)$$\phi \,\, = \,{{2I_1 D} \over {2D - h}}\quad.$$

The functions are obtained from fitting FE data for a total of 2244 cases: 17 cases for C ₁ (from 0.6 to 2.2 MPa with an increment of 0.1), 12 cases for C ₂ (from −1.1 to 0.0 MPa with an increment of 0.1), and 11 cases for C ₃ (from 0.1 to 0.6 MPa with an increment of 0.05). The mapping functions ψ and $\phi $can be expressed as

(24)$$\eqalignb{ {\rm{\psi }} \, = \,f_i^{\rm{\psi }} (C_1 ,\,C_2 ,\,C_3 )\left( {{h \over D}} \right)^i ;\;\,\,\crf_i^{\rm{\psi }} (C_1 ,\,C_2 ,\,C_3 ) = {\rm{\alpha }}_{ij}^{\rm{\psi }} (C_1 ,\,C_2 )\,C_3^j \quad, \cr \phi \, = \,3\, + \,f_i^\phi (C_1 ,\,C_2 ,\,C_3 )\left( {{h \over D}} \right)^i ;\,\,\,\crf_i^\phi (C_1 ,\,C_2 ,\,C_3 ) = {\rm{\alpha }}_{ij}^\phi (C_1 ,\,C_2 )\,C_3^j \quad, \cr i\, = \,1,\,2,\,3;\,\,\,j\, = \,0,\,1,\,2,\,3\quad, \cr} $$

where $f_i^{\rm{\psi }} $ and $f_i^{\phi} $ are the fitting coefficients with functions of C ₁, C ₂, and C ₃. For each C ₁ and C ₂, the coefficients $\alpha _{ij}^{\rm{\psi }} $ and $\alpha _{ij}^\phi $ in Eq. (24) are induced from the fitting with respect to C ₃. For each C ₁, the next coefficients $\beta _{ijk}^{\rm{\psi }} $ and $\beta _{ijk}^\phi $ in Eq. (25) are then obtained for the variation of C ₂ with the above $\alpha _{ij}^{\rm{\psi }} $ and $\alpha _{ij}^\phi $coefficients as follows:

(25)$$\eqalign{ {\rm{\alpha }}_{ij}^{\rm{\psi }} (C_1 ,\,C_2 ) = {\rm{\beta }}_{ijk}^{\rm{\psi }} (C_1 )\,C_2^k ;\,\,\,{\rm{\alpha }}_{ij}^\phi (C_1 ,\,C_2 ) = {\rm{\beta }}_{ijk}^\phi (C_1 )\,C_2^k\quad , \cr i = 1,\,2,\,3\,;\,\,j,k = 0,\,1,\,2,\,3\quad. \cr} $$

Finally, ${\rm{\gamma }}_{ijkl}^{\rm{\psi }} $ and ${\rm{\gamma }}_{ijkl}^\phi $ are produced for the variation of C ₁ with the above ${\rm{\beta }}_{ijk}^{\rm{\psi }} $ and ${\rm{\beta }}_{ijk}^\phi $coefficients as follows:

(26)$$\eqalignb{ {\rm{\beta }}_{ijk}^{\rm{\psi }} (C_1 ) = {\rm{\gamma }}_{ijkl}^{\rm{\psi }} C_1^l ;\,\,{\rm{\beta }}_{ijk}^\phi (C_1 ) = {\rm{\gamma }}_{ijkl}^\phi C_1^l \quad, \cr i = 1,\,2,\,3;\,\,\,j,k,l = 0,\,1,\,2,\,3\quad. \cr} $$

In Fig. 12 and 13, FE data for $\psi - h/D$ and $\phi - h/D$ (symbols) are compared with curves generated from Eq. (24) for C ₁ = 1.0 MPa and C ₂= −0.3, −0.4, and −0.5 MPa. The developed fitting functions closely follow the FE results.

FIG. 12.

Comparison of regressed ψ curves with FEA data.

FIG. 13.

Comparison of regressed $\phi $curves with FEA data.

C. Evaluation of the Yeoh properties

In the preceding section, the normalized functions, ψ and $\phi $, which map the indentation load–displacement curve onto the W–I ₁ curve, were established from FE data for a property range of the generally used rubber materials. The flowchart of the property evaluation program by using the normalized mapping functions is shown in Fig. 14. The detailed procedure can be described as follows. First, indentation tests are simulated and a load–displacement curve up to h/D = 1.0 is generated from FE analysis. C ₂ and C ₃ are initially assumed to be −0.4 and 0.4 MPa, respectively; second, C ₁ is determined by substituting the load at h/D = 0.1 into Eq. (21). Using C ₁, C ₂, and C ₃, the ψ and $\phi $ functions in Eq. (24) are determined; the W–I ₁ curve corresponding to the load–displacement curve is obtained from Eqs. (22) and (23). Regression of this W–I ₁ curve produces updated coefficients, C ₂ and C ₃. These updated properties are compared with the previous ones. The process is repeated until the updated C ₂ and C ₃ satisfy a convergence limit. The predicted properties are compared with given material properties for 45 cases (Table II). The average errors of C ₁, C ₂, and C ₃are 0.6%, 7.0%, and 7.6%, and the standard deviations are 0.5, 5.8, and 5.6%, respectively.

FIG. 14.

Flowchart for the determination of the Yeoh coefficients.

TABLE II.

Comparison of computed material property values to those given.

To investigate the effects of prediction errors on W–I ₁ and nominal stress–strain curves, we compare W–I ₁, stress–strain, and load–displacement curves for given properties with those of predicted properties, especially for the two kinds of materials having relatively large prediction errors in 45 cases (Fig. 15). The prediction errors are about 20.8% and 27.3% for C ₂ and 11.0% and 17.6% for C ₃. However, as we cannot distinguish the differences in the curves in Fig. 15, we may state that the quantitative errors of the properties are of secondary importance. The accuracy of the developed method should be therefore confirmed by comparisons of these three curves. From these comparisons, our new technique is highly accurate for the selected ranges of material properties of hyperelastic rubber.

FIG. 15.

Comparison of the predicted data to those given for (a) W–I ₁, (b) nominal stress–nominal strain, and (c) P–h curves.