3. Hyperelastic material models for elastomers

Elastomers used in seals, bushings, tires, and vibration isolators can sustain engineering strains well beyond 100% and recover almost completely. Linear elasticity cannot represent this behavior, so large-deformation FEA relies on hyperelastic constitutive models defined by a strain energy density function W. The stress is obtained by differentiating W with respect to the deformation, which guarantees path-independence and thermodynamic consistency. This section reviews the five most widely used formulations — neo-Hookean, Mooney-Rivlin, Yeoh, Ogden, and Arruda-Boyce — and shows how to identify their parameters from test data.

Kinematics and the Strain Energy Framework

For an isotropic hyperelastic solid, the strain energy is expressed through the invariants of the right Cauchy–Green tensor C = F^TF:

I₁ = tr(C),   I₂ = ½[(tr C)² − tr(C²)],   J = det(F)

Nearly incompressible elastomers are usually split into deviatoric and volumetric parts, W = W_dev(Ī₁, Ī₂) + W_vol(J), where Ī₁ = J^−2/3I₁ and Ī₂ = J^−4/3I₂. The volumetric term is typically W_vol = (1/D₁)(J − 1)², with D₁ controlling bulk stiffness. In practice, D₁ is set from a target bulk modulus (e.g., K ≈ 1000·C₁₀) rather than fitted to mechanical data.

Model Formulations

Neo-Hookean. The simplest form, W = C₁₀(Ī₁ − 3), has a single deviatoric constant. It captures the initial shear modulus μ = 2C₁₀ but cannot reproduce the S-shaped upturn at large strain, so it is best limited to moderate deformations (roughly <50% strain) or used as a baseline.

Mooney-Rivlin. Adding a second invariant term gives W = C₁₀(Ī₁ − 3) + C₀₁(Ī₂ − 3). It improves the fit for uniaxial and equibiaxial data simultaneously and remains popular for rubber components with moderate strain. However, it can exhibit unstable (non-monotonic) stress predictions in some deformation modes if the constants are poorly conditioned.

Yeoh. This model retains only the first invariant but adds higher-order terms: W = Σ_i=1..3 C_i0(Ī₁ − 3)ⁱ. Because it depends solely on I₁, it is robust and captures the strong stiffening at large strain. It is the default choice for filled rubbers and for cases where only uniaxial data are available.

Ogden. Formulated directly in principal stretches λ_i, W = Σ_p (2μ_p/α_p²)(λ₁^α_p + λ₂^α_p + λ₃^α_p − 3). The fractional exponents α_p allow excellent fits across multiple loading modes, but the parameters are not directly tied to physical quantities and require careful regularization to avoid instability.

Arruda-Boyce. Based on an eight-chain network representation, W = μ Σ_i (C_i/λ_m²ⁱ⁻²)(Ī₁ⁱ − 3ⁱ), where λ_m is the locking stretch. It predicts the characteristic stress upturn as chains approach full extension and needs only two physical parameters (μ, λ_m), making it attractive when data are limited.

_m μ = 0.8 MPa, λ_m = 3.2 Limited data, physical insight

Parameter Identification from Test Data

Reliable parameters come from fitting the model to experimental stress–strain curves. The standard procedure uses three homogeneous deformation modes: uniaxial tension, equibiaxial tension, and planar (pure) shear. For each mode, the nominal stress P is expressed analytically in terms of the model constants and the stretch λ. For example, for uniaxial tension of an incompressible neo-Hookean solid:

P_uni = 2·C10·(λ − λ⁻²)

For Mooney-Rivlin, the uniaxial nominal stress becomes:

P_uni = 2·(λ − λ⁻²)·(C10 + C01/λ)

Similar expressions exist for equibiaxial and planar shear. A least-squares minimization over all three datasets simultaneously yields the constants. In commercial codes such as Abaqus, this is done with the *HYPERELASTIC material option and the built-in material evaluation tool; in open-source workflows, a Python script using scipy.optimize.least_squares can perform the same fit.

Worked Example: Fitting a Yeoh Model

Suppose a filled natural rubber specimen is tested in uniaxial tension up to λ = 4.0. The measured nominal stress at selected stretches is given below (example data, not from a real test):

λ	P (MPa)
1.5	0.62
2.0	1.05
2.5	1.48
3.0	2.10
3.5	3.25
4.0	5.40

Using the Yeoh uniaxial expression P = 2(λ − λ⁻²)·(C₁₀ + 2C₂₀(Ī₁ − 3) + 3C₃₀(Ī₁ − 3)²) with Ī₁ = λ² + 2/λ, a least-squares fit yields approximate values C₁₀ ≈ 0.38 MPa, C₂₀ ≈ −0.02 MPa, C₃₀ ≈ 0.003 MPa. The negative C₂₀ is typical for filled rubbers and reflects the initial softening region. Plotting the fitted curve against the data should show residuals below 5% of the peak stress; if not, add equibiaxial data or switch to Ogden.

Model Selection Guidance

There is no universal best model. The choice depends on the available test data, the strain range, and the deformation modes present in the FE simulation. For small strains (<50%) and limited data, neo-Hookean or Arruda-Boyce suffice. For moderate strains with multiaxial data, Mooney-Rivlin is adequate. For large strains (>200%) and filled rubbers, Yeoh or Ogden are preferred. Always validate the fitted model against an independent test mode before committing to a full simulation, and monitor the FE solution for convergence issues that may indicate an unstable parameter set.