Constitutive Models

Principal Cauchy stress constitutive models

Hyperelastic material models often rely upon the strain energy density function W formulated in terms of strain invariants. Under the assumption of an isotropic incompressible behavior of the material under uniaxial loading, the principal Cauchy stress constitutive models as a function of invariants yields:

$\sigma_{\text{true uniax}} = 2 \left(\lambda^{2} -\dfrac{1}{\lambda}\right) \cdot \left(\dfrac{\partial W}{\partial I_{1}} + \dfrac{1}{\lambda}\dfrac{\partial W}{\partial I_{2}}\right)$

$\sigma_{\text{eng uniax}} = 2 \left(\lambda -\dfrac{1}{\lambda^{2}}\right) \cdot \left(\dfrac{\partial W}{\partial I_{1}} + \dfrac{1}{\lambda}\dfrac{\partial W}{\partial I_{2}}\right)$

where, the Cauchy-Green invariants I₁ and I₂, are in this particular case a function of the stretch lambda:

Ogden Model

Strain energy density function

Principal Cauchy stress

Python implementation

import numpy as np

class Hyperelastic:

    def __init__(self, model, parameters, order, data_type):
        self.model = model
        self.order = order
        self.parameters = parameters
        self.param_names = []
        self.data_type = data_type

        if model == 'Ogden':
            initialGuessMu = np.array([0.1]*self.order)
            initialGuessAlpha = np.array([0.2]*self.order)
            self.initialGuessParam = np.append(initialGuessMu,initialGuessAlpha)
            self.nbparam = self.order*2
            muVec_names = ["µ1","µ2","µ3"][0:self.order]
            alphaVec_names = ["α1","α2","α3"][0:self.order]
            self.param_names = np.append(muVec_names,alphaVec_names)

def OgdenModel(self, parameters, Strain):
    """Ogden hyperelastic model (incompressible material under uniaxial tension)"""

    # parameter is a 1D array : [mu0,mu1,...,mun,alpha0,alpha1,...,alphan] 
    muVec = parameters.reshape(2, self.order)[0]    # muVec is a 1D array : [mu0,mu1,...,mun]
    alphaVec = parameters.reshape(2, self.order)[1] # muVec is a 1D array : [alpha0,alpha1,...,alphan]

    if self.data_type == 'True':
        lambd = np.exp(Strain)   # lambd i.e lambda
    elif self.data_type == 'Engineering':
        lambd = 1 + Strain
    else:
        print("Data type error. Data is neither 'True' or 'Engineering'. ") 

    # broadcasting method to speed up computation
    lambd = lambd[np.newaxis, :]
    muVec = muVec[:self.order, np.newaxis]
    alphaVec = alphaVec[:self.order, np.newaxis]

    if self.data_type == 'True':
        Stress = np.sum(2*muVec*(lambd**(alphaVec - 1) - lambd**(-((1/2)*alphaVec + 1))), axis=0)
    elif self.data_type == 'Engineering':
        Stress = np.sum((2*muVec*(lambd**(alphaVec - 1) - lambd**(-((1/2)*alphaVec + 1)))/lambd), axis=0)

    return Stress

Example

import numpy as np
import matplotlib as plt
from Hyperelastic import Hyperelastic

# Hyperelastic object
hyperelastic = Hyperelastic('Ogden', np.array([0]), 3, 'True')  # Instantiate an Hyperalastic object : Ogden model of order 3 using True data type

model_param = np.array([0.354, -0.129, -0.226, 3.316,  3.278, 3.278])   # [mu1, mu2, mu3, alpha1, alpha2, alpha3]
trueStrain = linspace(0, 1.5, 50)
trueStress = hyperelastic.ConsitutiveModel(model_param, trueStrain)

plt.plot(trueStrain, trueStress, 'r')
plt.xlabel('True Strain')
plt.ylabel('True Stress (MPa)')
plt.grid('on')
plt.show() 

Yeoh Model

Strain energy density function

Principal Cauchy stress

Python implementation

def YeohModel(self, cVec, Strain):
    """Yeoh hyperelastic model (incompressible material under uniaxial tension)"""
    
    if self.data_type == 'True':
        lambd = np.exp(Strain)
    elif self.data_type == 'Engineering':
        lambd = 1 + Strain
    else:
        print("Data type error. Data is neither 'True' or 'Engineering'. ")

    I1 = lambd**2 + 2/lambd

    Stress = np.zeros((self.order,len(Strain)))

    for i in range (0,self.order):
        if self.data_type == 'True':
            Stress[i,:] = 2*(lambd**2 - 1/lambd)*(i+1)*cVec[i]*((I1-3)**(i)) # true
        elif self.data_type == 'Engineering':
            Stress[i,:] = (2*(lambd - 1/(lambd**2))*(i+1)*cVec[i]*(I1-3)**(i))/lambd # eng
        else:
            print("Data type error. Data is neither 'True' or 'Engineering'. ")

        Stress_sum = np.sum(Stress, axis=0)
    return Stress_sum

Example

import numpy as np
import matplotlib as plt
from Hyperelastic import Hyperelastic

# Hyperelastic object
hyperelastic = Hyperelastic('Yeoh', np.array([0]), 2, 'Engineering')  # Instantiate an Hyperalastic object : Yeoh model of order 2 using Engineering data type

model_param = np.array([3.278, 0.384])   # [C1, C2]
engStrain = linspace(0, 1.5, 50)
engStress = hyperelastic.ConsitutiveModel(model_param, engStrain)

plt.plot(engStrain, engStress, 'r')
plt.xlabel('Engineering Strain')
plt.ylabel('Engineering Stress (MPa)')
plt.grid('on')
plt.show() 

Relevant references

Martins PALS, Natal Jorge RM, Ferreira AJM. A Comparative Study of Several Material Models for Prediction of Hyperelastic Properties: Application to Silicone-Rubber and Soft Tissues. Strain. 2006;42(3):135–147.

Rackl M. Material testing and hyperelastic material model curve fitting for Ogden, Polynomial and Yeoh models. In: ScilabTEC (7th International Scilab Users Confer ence). Paris; 2015.

Ogden RW, Saccomandi G, Sgura I. Fitting hyperelastic models to experimental data. Computational Mechanics. 2004 nov;34(6):484–502.

Bergström J. Elasticity/Hyperelasticity. In: Mechanics of Solid Polymers: Theory and Computational Modeling. William Andrew Publishing; 2015. p.209–307.