Skip to content

Constitutive Models

Jump to bottom
Luc Marechal edited this page Jun 15, 2023 · 110 revisions

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 I1 and I2, are in this particular case a function of the stretch lambda:

$I_{1} = \lambda^{2} + \dfrac{2}{\lambda}$, $I_{2} = 2\lambda + \dfrac{1}{\lambda^{2}}$, with $\lambda=exp(\epsilon_{\text{true uniax}})$

For the case of an incompressible Mooney–Rivlin material under uniaxial elongation, $\lambda_{1} = \lambda$, $\lambda_{2}=\lambda_{3}=1/{\sqrt{\lambda}}$

Neo-Hookean Model

Strain energy density function

$W = C_{10}(I_{1}−3)$

Where, $I_{1}$ is the first invariant of the right Cauchy-Green Tensor

With the initial shear modulus: $\mu = 2 C_{10} $

Principal Cauchy stress for uniaxial tension

$\sigma_{\text{true uniax}} = \mu \left(\lambda^{2} -\dfrac{1}{\lambda}\right)$

$\sigma_{\text{eng uniax}} = \mu \left(\lambda - \dfrac{1}{\lambda^{2}}\right)$

The Neo-Hookean model is a simple model that is typically only accurate for strains less than 20%.

Mooney-Rivlin Model

Strain energy density function

$W = C_{10}(I_{1}-3)+C_{01}(I_{2}-3)+C_{20}(I_{1}-3)^{2}$

Where, $I_{1}$ and $I_{2}$ are respectively the first and second invariants of the right Cauchy-Green Tensor

Principal Cauchy stress for uniaxial tension [1]:

$\sigma_{\text{true uniax}} = 2\left(\lambda^{2} -\dfrac{1}{\lambda}\right)\left[C_{10} + \dfrac{C_{01}}{\lambda} + 2C_{20}\left(\lambda^{2} + \dfrac{2}{\lambda}-3\right)\right]$

$\sigma_{\text{eng uniax}} = 2\left(\lambda -\dfrac{1}{\lambda^{2}}\right)\left[C_{10} + \dfrac{C_{01}}{\lambda} + 2C_{20}\left(\lambda^{2} + \dfrac{2}{\lambda}-3\right)\right]$

With the initial shear modulus: $\mu = 2\left(C_{01}+C_{10}\right)> 0$

The Mooney-Rivlin material law is typically accurate for strains up to 100%

Python implementation

def MooneyRivlinModel(cVec, Strain, order, data_type):
    """Mooney Rivlin hyperelastic model (incompressible material under uniaxial tension)"""
    
    # cVec = [C10,C01,C20]
    cVec = np.append(cVec, np.zeros(3-order) ) #To ensure CXX is zero if unsed
    C10 = cVec[0]
    C01 = cVec[1]
    C20 = cVec[2]

    if data_type == 'True':
        lambd = np.exp(Strain)       # lambd i.e lambda
        Stress = 2*(lambd**2 - 1/lambd)*(C10 + C01/lambd + 2*C20*(lambd**2 + 2/lambd -3))
    elif data_type == 'Engineering':
        lambd = 1 + Strain
        Stress = 2*(lambd - 1/lambd**2)*(C10 + C01/lambd + 2*C20*(lambd**2 + 2/lambd -3))
    else:
        print("Data type error. Data is neither 'True' or 'Engineering'. ")

    return Stress

Ogden Model

Strain energy density function

$W=\sum\limits_{p=1}^{n}{\frac{\mu_{p}}{\alpha_{p}}}\left[\lambda_{1}^{\alpha_{p}}+\lambda_{2}^{\alpha_{p}}+\lambda_{3}^{\alpha_{p}}-3\right]$

Principal Cauchy stress for uniaxial tension

Under the assumption of an incompressible material under uniaxial tension [2]:

$\sigma_{\text{true uniax}}=\sum\limits_{p=1}^{N}{\mu_{p}}\left(\lambda^{\alpha_{p}} - \lambda^{-\frac{1}{2}\alpha_p}\right)$

$\sigma_{\text{eng uniax}}=\sum\limits_{p=1}^{N}{\mu_{p}}\left(\lambda^{(\alpha_{p}-1)} - \lambda^{(-\frac{1}{2}\alpha_{p}-1)}\right)$

For material stability, it is required that each material constant pair : $\alpha_{p}.\mu_{p} > 0$

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        
        self.fitting_method = 'lm'

        if model == 'Ogden':
            initialGuessMu = np.array([1.0]*self.order)
            initialGuessAlpha = np.array([1.0]*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)
            self.fitting_method = 'trust-constr'
        elif model == 'Yeoh':
            self.initialGuessParam = np.array([0.1]*self.order)
            self.nbparam = self.order
            self.param_names = ["C1","C2","C3"][0:self.order]
            self.fitting_method = 'lm'
        else:
            print("Error")


    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[0:self.order]  # [mu0,mu1,...,mun]
        alphaVec = parameters[self.order:]  # [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(muVec*(lambd**alphaVec - 1/(lambd**(alphaVec/2))), axis=0)
        elif self.data_type == 'Engineering':
            Stress = np.sum((muVec*(lambd**alphaVec - 1/(lambd**(alphaVec/2)))/lambd), axis=0)

        return Stress

Defining Linear Constraints:

The Mooney Rivlin model requires the following linear constraints for initial shear modulus: $\mu = 2\left(C_{01}+C_{10}\right)> 0$

The constraints $C_{01}+C_{10} > 0$ can be written in the linear constraint standard format:

$$\left[\begin{array}{c} 0 \\ 0 \end{array}\right] \leq \left[\begin{array}{cc} 1 & 1 \\ 0 & 0 \end{array}\right] \left[\begin{array}{c} C_{01} \\ C_{10} \end{array}\right] \leq\left[\begin{array}{c} +\infty \\ +\infty \end{array}\right]$$

and defined using a LinearConstraint object.

from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint([[1, 1], [0, 0]], [0, 0], [np.inf, np.inf])

Defining Non Linear Constraints:

The Ogden model requires the following non linear constraints for material stability : $\alpha_{p}.\mu_{p} > 0$

We use the trust-region constrained method to deal with constrained minimization problems. The minimize from _scipy.optimize.minimize_function provides algorithms for constrained minimization, namely 'trust-constr'. The method 'trust-constr' requires the constraints to be defined as a sequence of objects LinearConstraint and NonlinearConstraint.

$$x=\left[\begin{array}{c}\mu_1\\ \mu_2\\ \mu_3\\ \alpha_1\\ \alpha_2\\ \alpha_3 \end{array}\right]$$

The nonlinear constraint:

$$c(x)=\left[\begin{array}{c} \mu_1 . \alpha_1 \\ \mu_2 . \alpha_2 \\ \mu_3 . \alpha_3 \\ \end{array}\right] \geq\left[\begin{array}{l} 0 \\ 0 \\ 0 \end{array}\right]$$

with Jacobian matrix:

$$J(x)=\left[\begin{array}{cccccc} \alpha_1 & 0 & 0 & \mu_1 & 0 & 0\\ 0 & \alpha_2 & 0 & 0 & \mu_2 & 0 \\ 0 & 0 & \alpha_3 & 0 & 0 & \mu_3 \end{array}\right]$$

is defined using a NonlinearConstraint object.

# x = [mu1, mu2, mu3, alpha1, alpha2, alpha3]
def cons_f(x):
    return [x[0]*x[3], x[1]*x[4], x[2]*x[5]]

def cons_J(x):
    return [[x[3], 0, 0, x[0], 0, 0], [0, x[4], 0, 0, x[1], 0], [0, 0, x[5], 0, 0, x[2]]]

from scipy.optimize import NonlinearConstraint
nonlinear_constraint = NonlinearConstraint(cons_f, 0, np.inf, jac=cons_J)

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

equation

Principal Cauchy stress for uniaxial tension [3]:

equation

equation

With the initial shear modulus: $\mu = 2 C_{1}$

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

[1] Rivlin, R. S., Large elastic deformations of isotropic materials. IV. Further developments of the general theory, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 241(835), 1948. pp. 379–397.

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

[3] Yeoh, O. H., Some forms of the strain energy function for rubber, Rubber Chemistry and technology, Volume 66, Issue 5, November 1993, Pages 754-771.

[3] 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.

[4] 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.

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

Clone this wiki locally