pymecht.MatModel module

class pymecht.MatModel.MatModel(*modelsList)[source]

Bases: object

The main material model class that allows adding different models together. This will usually be the starting point in using pyMechT.

The models can be provided as a list of model objects, or as a list of strings.

The passed strings can have lower/upper case. The following models are available:

‘NH’: Neo-Hookean model

‘MR’: Mooney-Rivlin model

‘YEOH’: Yeoh model

‘LS’: Lee-Sacks model

‘MN’: May-Newman model

‘GOH’: Gasser-Ogden-Holzapfel model

‘HGO’: Holzapfel-Gasser-Ogden model

‘expI1’: A model with an exponential of I1

‘polyI4’: A model with a polynomial of I4

‘HY’: Humphrey-Yin model

‘Holzapfel’: Holzapfel model

‘volPenalty’: A penalty model for volumetric change

‘ArrudaBoyce’: Arruda-Boyce model

‘Gent’: Gent model

‘splineI1’: A spline model of I1

‘splineI1I4’: A spline model of I1 and I4

‘StructModel’: A structural model with fiber distribution

‘ARB’: Arbitrary model with user-defined strain energy density function

Multiple models can be added together to create a composite model, with each model having its own fiber(s) and parameters.

Example

>>> from MatModel import *
>>> mat1 = NH() #A neo-Hookean model
>>> mat2 = GOH([np.array([1,0,0]),np.array([0.5,0.5,0])]) #A GOH model with two fiber families
>>> mat_model = MatModel(mat1,mat2) #can provide as many models as needed as long as their parameters are uniquely named
>>> mat_model = MatModel(mat1) + MatModel(mat2) #can also add them after creating the MatModel object
>>> mat_model.models #provides a list of the models included
>>> mat_model.parameters #provides a dictionary of parameters and their default values
>>> mat_model = MatModel('GOH','NH') #can also provide a list of model names, however one has to assign fiber directions afterwards
>>> mats = mat_model.models
>>> mats[0].fiber_dirs = [np.array([1,0,0]),np.array([0.5,0.5,0])]

property parameters: Parameters of the model

property models: List of component models

energy(F=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), theta=None)[source]

Returns the energy density of the material model

Parameters:

F (np.array([3,3])) – the deformation gradient, default: identity matrix (no deformation)
theta (ParamDict or dict) – the parameters of the model, if None, then the default values are used

Returns:

the energy density

Return type:

float

stress(F=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), theta=None, stresstype='cauchy', incomp=False, Fdiag=False)[source]

Returns the stress tensor of the material model

Parameters:

F (np.array([3,3])) – the deformation gradient, default: identity matrix (no deformation)
theta (ParamDict or dict) – the parameters of the model, if None, then the default values are used
stresstype (str) –
the type of stress tensor to return with the following options (case insensitive)
- ’cauchy’: Cauchy stress,
- ’1pk’ or ‘1stpk’ or ‘firstpk’: 1st Piola-Kirchoff stress
- ’2pk’ or ‘2ndpk’ or ‘secondpk’: 2nd Piola-Kirchoff stress
- default: Cauchy stress
incomp (bool) – if True, then the material is assumed to be incompressible, default: False
Fdiag (bool) – if True, then it is assumed that F is diagonal (for faster computation), default: False

Returns:

the stress 3X3 tensor

Return type:

np.array([3,3])

energy_stress(F=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), theta=None, stresstype='cauchy', incomp=False, Fdiag=False)[source]

Returns the energy density and stress tensor of the material model

Parameters:

F (np.array([3,3])) – the deformation gradient, default: identity matrix (no deformation)
theta (ParamDict or dict) – the parameters of the model, if None, then the default values are used
stresstype (str) –
the type of stress tensor to return with the following options (case insensitive)
- ’cauchy’: Cauchy stress,
- ’1pk’ or ‘1stpk’ or ‘firstpk’: 1st Piola-Kirchoff stress
- ’2pk’ or ‘2ndpk’ or ‘secondpk’: 2nd Piola-Kirchoff stress
- default: Cauchy stress
incomp (bool) – if True, then the material is assumed to be incompressible, default: False
Fdiag (bool) – if True, then it is assumed that F is diagonal (for faster computation), default: False

Returns:

first element is the energy density second element is the stress 3X3 tensor

Return type:

tuple of (float,np.array([3,3]))

class pymecht.MatModel.InvariantHyperelastic[source]

Bases: object

An abstract class from which all the invariant-based hyperelastic models should be derived. Currently, it allows for models that depend on I1, I2, J, and I4 with multiple fiber families It should not be used directly, but rather derived from. Use MatModel class to use the derived classes.

property fiber_dirs: Fiber direction vector(s)

class pymecht.MatModel.NH[source]

Bases: InvariantHyperelastic

Neo-Hookean model

\[\Psi = \frac{\mu}{2}(I_1-3)\]

class pymecht.MatModel.MR[source]

Bases: InvariantHyperelastic

Mooney-Rivlin model

\[\Psi = c_1(I_1-3) + c_2(I_2-3)\]

class pymecht.MatModel.YEOH[source]

Bases: InvariantHyperelastic

Yeoh model

\[\Psi = \sum_{i=1}^4 c_i(I_1-3)^i\]

class pymecht.MatModel.LS(M=[])[source]

Bases: InvariantHyperelastic

Lee–Sacks model

\[\Psi = \frac{k_1}{2(k_4k_2+(1-k_4)k_3)}(k_4\exp(k_2(I_1-3)^2)+(1-k_4)\exp(k_3(I_4-1)^2)-1)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> k1,k2,k3,k4,I1bar,I4bar = symbols('k1 k2 k3 k4 I1bar I4bar')
>>> Psi = k1/2/(k4*k2+(1-k4)*k3)*(k4*exp(k2*(I1bar-3)**2) + (1-k4)*exp(k3*(I4bar-1)**2)-1)
>>> diff(Psi,I1bar)
>>> diff(Psi,I4bar)

class pymecht.MatModel.MN(M=[])[source]

Bases: InvariantHyperelastic

MayNewman model

\[\Psi = \sum_{i=1}^N \frac{k_1}{k_2+k_3}(\exp(k_2(I_1-3)^2+k_3(\sqrt{I_{4i}}-1)^4)-1)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> k1,k2,k3,I1bar,I4bar = symbols('k1 k2 k3 I1bar I4bar')
>>> Psi = k1/(k2+k3)*(exp(k2*(I1bar-3)**2+k3*(sqrt(I4bar)-1)**4)-1)
>>> diff(Psi,I1bar)
>>> diff(Psi,I4bar)

class pymecht.MatModel.GOH(M=[])[source]

Bases: InvariantHyperelastic

Gasser-Ogden-Holzapfel model (without the Neo-Hookean term)

\[\Psi = \sum_{i=1}^N \frac{k_1}{2k_2}(\exp(k_2(k_3I_1+(1-3k_3)I_{4i}-1)^2)-1)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> k1,k2,k3,I1bar,I4bar = symbols('k1 k2 k3 I1bar I4bar')
>>> Psi = k1/2/k2*(exp(k2*(k3*I1bar+(1-3*k3)*I4bar-1)**2)-1)
>>> diff(Psi,I1bar)
>>> diff(Psi,I4bar)

class pymecht.MatModel.Holzapfel(M=[])[source]

Bases: InvariantHyperelastic

Holzapfel (2005) model without the neo-Hookean term

\[\Psi = \sum_{i=1}^N \frac{k_1}{2k_2}(\exp(k_2(k_3(I_1-3)^2+(1-k_3)(I_{4i}-1)^2)-1)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> k1,k2,k3,I1bar,I4bar = symbols('k1 k2 k3 I1bar I4bar')
>>> Psi = k1/2/k2*(exp(k2*(k3*(I1bar-3)**2+(1-k3)*(I4bar-1)**2))-1)
>>> diff(Psi,I1bar)
>>> diff(Psi,I4bar)

class pymecht.MatModel.expI1[source]

Bases: InvariantHyperelastic

Model with exponential of I1

\[\Psi = \frac{k_1}{k_2}(\exp(k_2(I_1-3))-1)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> k1,k2,I1bar = symbols('k1 k2 I1bar')
>>> Psi = k1/k2*(exp(k2*(I1bar-3))-1)
>>> diff(Psi,I1bar)

class pymecht.MatModel.HGO(M=[])[source]

Bases: InvariantHyperelastic

HGO 2000 model’s fiber part

\[\Psi = \sum_{i=1}^N \frac{k_3}{2k_4}(\exp(k_4(I_{4i}-1)^2)-1)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> k1,k2,k3,k4,I1bar,I4bar = symbols('k1 k2 k3 k4 I1bar I4bar')
>>> Psi = k3/2./k4*(exp(k4*(I4bar-1)**2)-1)
>>> diff(Psi,I1bar)
>>> diff(Psi,I4bar)

class pymecht.MatModel.HY(M=[])[source]

Bases: InvariantHyperelastic

Humphrey-Yin model’s fiber part

\[\Psi = \sum_{i=1}^N \frac{k_3}{k_4}(\exp(k_4(\sqrt{I_{4i}}-1)^2)-1)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> k1,k2,k3,k4,I1bar,I4bar = symbols('k1 k2 k3 k4 I1bar I4bar')
>>> Psi = k3/k4*(exp(k4*(sqrt(I4bar)-1)**2)-1)
>>> diff(Psi,I1bar)
>>> diff(Psi,I4bar)

class pymecht.MatModel.volPenalty[source]

Bases: InvariantHyperelastic

Volumetric penality

\[\Psi = \frac{\kappa}{2}(J-1)^2\]

class pymecht.MatModel.polyI4[source]

Bases: InvariantHyperelastic

Polynomial in I4 model

\[\Psi = \sum_{i=1}^3 d_i(I_{4i}-1)^i\]

class pymecht.MatModel.ArrudaBoyce[source]

Bases: InvariantHyperelastic

Arruda-Boyce model

\[\Psi = \sum_{i=1}^5 C\frac{\alpha_i}{N^i}(I_1^i-3^i)\]

where \(\alpha_i\) are [1,1/20,11/1050,19/7000,519/673750]

class pymecht.MatModel.Gent[source]

Bases: InvariantHyperelastic

Gent model

\[\Psi = -\frac{\mu J_m}{2}\ln\left(1-\frac{I_1-3}{J_m}\right)\]

Derivation

The partial derivative expressions are obtained using the following code

>>> from sympy import *
>>> I1, mu, Jm = symbols('I1,mu,Jm')
>>> Psi = -mu*Jm/2.*ln(1-(I1-3)/Jm)
>>> diff(Psi,I1)

class pymecht.MatModel.splineI1[source]

Bases: InvariantHyperelastic

Spline-based model in I1 for data driven models (without an analytical expression)

Psi is loaded from a scipy spline object

set(W, alpha=1)[source]

Set the spline function and the weight

Parameters:

W (scipy.interpolate._bsplines.BSpline) – The spline function
alpha (float) – The weight of W in the energy function (i.e., \(\Psi = \alpha W(I1)\))

Return type:

None

class pymecht.MatModel.splineI1I4[source]

Bases: InvariantHyperelastic

Spline-based model in I1 and I4 for data driven models (without an analytical expression)

Psi is loaded from a scipy spline object

set(W, alpha=1)[source]

Set the spline function and the weight

Parameters:

W (scipy.interpolate.fitpack2.RectBivariateSpline) – The spline function
alpha (float) – The weight of W in the energy function (i.e., \(\Psi = \alpha W(I1,I4)\))

Return type:

None

Notes

The spline function should be defined on a rectangular grid

class pymecht.MatModel.StructModel[source]

Bases: object

A class for simplified structural model with integration over fiber directions

\[\Psi = \int_{-\pi/2}^{\pi/2} \Gamma(\theta) \psi(\theta) d\theta\]

where \(\Gamma(\theta)\) is the fiber distribution function and \(\psi(\theta)\) is the fiber energy function.

The fiber distribution function is assumed to be a truncated Gaussian distribution

\[\Gamma(\theta) = d_e\frac{1}{P}\exp\left(-\frac{(\theta-\theta_0)}{2\sigma^2}\right) + \frac{1-d_e}{\pi}\]

where \(P\) is the normalization factor and \(d_e\) is the anisotropy fraction.

The fiber energy function is assumed to be an exponential function

\[\psi(\epsilon) = \frac{A}{B}\left(\frac{\exp(B \epsilon)-1}{B}-\epsilon\right)\]

class pymecht.MatModel.ARB(_W='', _init_guess='', _low_bound='', _up_bound='', dparams=False)[source]

Bases: InvariantHyperelastic

A material model class that allows the arbitrary definition of any strain energy density function (SEDF). Sympy’s symbolic differentiation is used to calculate the partial derivatives of the SEDF with respect to the invariants I1, I2, J, and I4. The SEDF is provided as a string as a function of the invariants I1, I2, J, and I4, alongside initial guesses and upper/lower bounds for the parameters. The SEDF an initial guess for the parameters must be provided to identify parameters of the material model. If the strings are not provided, the user will be prompted to provide them.

The parsed strings have the following functions replaced by their respective numpy equivalents:

exp -> np.exp
sqrt -> np.sqrt
log -> np.log
log10 -> np.log10
log2 -> np.log2
sin -> np.sin
cos -> np.cos
tan -> np.tan
arcsin -> np.arcsin
asin -> np.asin
arccos -> np.arccos
acos -> np.acos
arctan -> np.arctan
atan -> np.atan
hypot -> np.hypot
arctan2 -> np.arctan2
sinh -> np.sinh
cosh -> np.cosh
tanh -> np.tanh
arcsinh -> np.arcsinh
asinh -> np.asinh
arccosh -> np.arccosh
acosh -> np.acosh
arctanh -> np.arctanh
atanh -> np.atanh

Example

>>> from MatModel import *
>>> model = pmt.ARB('mu/2.*(I1-3)','mu=1.','mu=0.01','mu=10.')
>>> mat = pmt.MatModel(model)
>>> F = np.random.rand
>>> mat.stress(F)

partial_deriv(**extra_args)[source]

Returns the partial derivatives of the SEDF. In the case of dPsidIi(I4), for i < 4, the partial derivative is summed over all fibre directions

Returns:: the partial derivative of the SEDF with respect to I1, I2, J, and I4.
Return type:: list of floats