Using MatModel
[1]:
import pymecht as pmt
import numpy as np
Single material model
The core module in pyMechT is MatModel which allows us to add models together to simulate the stress-strain behavior. A basic usage is when there is only one model, for example we can create a material with Yeoh model, which is isotropic and only a function of the first invariant:
[2]:
mat1 = pmt.MatModel('yeoh')
print(mat1.parameters) # Returns parameters as a dictionary
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
c1_0 1.00 No 1.00e-04 1.00e+02
c2_0 1.00 No 0.00 1.00e+02
c3_0 1.00 No 0.00 1.00e+02
c4_0 0.00 No 0.00 1.00e+02
------------------------------------------------------------------
In principle, one can directly use this material to calculate stresses given a deformation gradient using stress method
[3]:
mat1.stress?
Signature:
mat1.stress(
F=array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]]),
theta=None,
stresstype='cauchy',
incomp=False,
Fdiag=False,
)
Docstring:
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
-------
np.array([3,3])
the stress 3X3 tensor
File: /usr/local/lib/python3.9/site-packages/pymecht/MatModel.py
Type: method
Thus, we can calculate the Cauchy stress at an identity deformation gradient. However, since this model does not have a volumetric part, we must set incompressibility to be true to get a zero stress. This is done using a Lagrange multiplier, which is calculated internally in pyMechT by setting the normal stress along the third axis equal to zero.
[4]:
mat1.stress(np.eye(3),incomp=True)
[4]:
array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
Or, we can calculate other stresses at non-identity deformation gradient
[5]:
F = np.array([[1.2,0,0],[0.1,1.2,0],[0,0,1]])
print(mat1.stress(F,stresstype='1pk',incomp=True))
print(mat1.stress(F,stresstype='2pk',incomp=True))
[[3.78128667 0.71615278 0. ]
[1.03126 3.78128667 0. ]
[0. 0. 0. ]]
[[3.15107222 0.59679398 0. ]
[0.59679398 3.10133939 0. ]
[0. 0. 0. ]]
One can also change the parameters and pass them while calculating the stress.
[6]:
params = mat1.parameters
params['c1_0'].set(2)
print(mat1.stress(F,params,stresstype='1pk',incomp=True))
[[4.51462 0.85504167 0. ]
[1.23126 4.51462 0. ]
[0. 0. 0. ]]
Alternatively, instead of passing the parameters, one can set the values and call the stress function without the parameters
[7]:
mat1.parameters = params
print(mat1.stress(F,stresstype='1pk',incomp=True))
[[4.51462 0.85504167 0. ]
[1.23126 4.51462 0. ]
[0. 0. 0. ]]
Adding materials
A convenient feature of pyMechT is that we can easily add different models together (potentially, the same model repeated). For example
[8]:
mat1 = pmt.MatModel('yeoh','nh')
mat2 = pmt.MatModel('goh','goh','nh')
print(mat1, mat1.parameters)
print(mat2, mat2.parameters)
Material model with 2 components:
Component1: YEOH
Component2: NH
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
c1_0 1.00 No 1.00e-04 1.00e+02
c2_0 1.00 No 0.00 1.00e+02
c3_0 1.00 No 0.00 1.00e+02
c4_0 0.00 No 0.00 1.00e+02
mu_1 1.00 No 1.00e-04 1.00e+02
------------------------------------------------------------------
Material model with 3 components:
Component1: GOH
Component2: GOH
Component3: NH
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
k1_0 10.00 No 0.10 30.00
k2_0 10.00 No 0.10 30.00
k3_0 0.10 No 0.00 0.33
k1_1 10.00 No 0.10 30.00
k2_1 10.00 No 0.10 30.00
k3_1 0.10 No 0.00 0.33
mu_2 1.00 No 1.00e-04 1.00e+02
------------------------------------------------------------------
We can even add them afterwards, although this is rarely needed
[9]:
mat3 = mat1 + mat2
print(mat3,mat3.parameters)
Material model with 5 components:
Component1: YEOH
Component2: NH
Component3: GOH
Component4: GOH
Component5: NH
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
c1_0 1.00 No 1.00e-04 1.00e+02
c2_0 1.00 No 0.00 1.00e+02
c3_0 1.00 No 0.00 1.00e+02
c4_0 0.00 No 0.00 1.00e+02
mu_1 1.00 No 1.00e-04 1.00e+02
k1_2 10.00 No 0.10 30.00
k2_2 10.00 No 0.10 30.00
k3_2 0.10 No 0.00 0.33
k1_3 10.00 No 0.10 30.00
k2_3 10.00 No 0.10 30.00
k3_3 0.10 No 0.00 0.33
mu_4 1.00 No 1.00e-04 1.00e+02
------------------------------------------------------------------
To get the model components back, we can use the .models
[10]:
mat2.models
[10]:
(<pymecht.MatModel.GOH at 0x14d4ca730>,
<pymecht.MatModel.GOH at 0x14d4ca9a0>,
<pymecht.MatModel.NH at 0x14d4caac0>)
Setting fiber directions
Note that, for anisotropic models (such as goh), fiber directions need to be specified before they can be called. That is, mat1 above can be used but not mat2
[11]:
print(mat1.stress(np.eye(3),incomp=True))
try:
mat2.stress(np.eye(3),incomp=True)
except ValueError as e:
print("Value error occured", e)
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
Value error occured GOH model class uses I4 but no fiber directions have been defined. Did you forget to set the fiber directions?
The fiber directions are set for each component model, and each component can have multiple fiber directions (the response is summed over them). This means, that the individual components can have different fiber directions (which is why one might want to add two GOH models together). Below we set two fiber directions for each of the GOH component.
[12]:
model_comps = mat2.models
model_comps[0].fiber_dirs = [ np.array([1,1,0]), np.array([0,1,0])]
model_comps[1].fiber_dirs = [ np.array([1,0,0]), np.array([0.5,1,0])]
print(mat2)
print(mat2.stress(np.eye(3),incomp=True))
print(mat2.stress(F,incomp=True))
Material model with 3 components:
Component1: GOH with fiber direction(s):[array([0.70710678, 0.70710678, 0. ]), array([0., 1., 0.])]
Component2: GOH with fiber direction(s):[array([1., 0., 0.]), array([0.4472136 , 0.89442719, 0. ])]
Component3: NH
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
[[ 84.11760455 66.84260966 0. ]
[ 66.84260966 126.67097888 0. ]
[ 0. 0. 0. ]]
Note that the fiber direction vectors we supplied were not unit vectors. Internally they are made unit vectors. Thus, when we print them, we see some differences. Lastly, if we also set fiber directions to isotropic materials, it does not affect their response.
[13]:
model_comps[2].fiber_dirs = [ np.array([1,0,0]), np.array([0.5,1,0])]
print(mat2)
print(mat2.stress(np.eye(3),incomp=True))
print(mat2.stress(F,incomp=True))
Material model with 3 components:
Component1: GOH with fiber direction(s):[array([0.70710678, 0.70710678, 0. ]), array([0., 1., 0.])]
Component2: GOH with fiber direction(s):[array([1., 0., 0.]), array([0.4472136 , 0.89442719, 0. ])]
Component3: NH with fiber direction(s):[array([1., 0., 0.]), array([0.4472136 , 0.89442719, 0. ])]
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
[[ 84.11760455 66.84260966 0. ]
[ 66.84260966 126.67097888 0. ]
[ 0. 0. 0. ]]
However, for using pyMechT, one would seldom use the MatModel directly. Instead, it would be simply created and then used to create a SampleExperiment instance, which internally sets the deformation gradient based on the mode of deformation, incompressibility condition, and even a helper function to set the fiber directions more easily. These aspects are covered in the next examples.
Arbitrary ARB model
If the desired hyperelastic model is not implemented in pymecht, then one can use the arbitrary ARB model. It uses symbolic toolbox sympy to convert a string in terms of I1, I2, J and I4 into a material model. An example is follows.
[14]:
pmt.ARB?
Init signature: pmt.ARB(_W='', _init_guess='', _low_bound='', _up_bound='', dparams=False)
Docstring:
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)
File: /usr/local/lib/python3.9/site-packages/pymecht/MatModel.py
Type: type
Subclasses:
[15]:
model_I1I2 = pmt.ARB('mu/2.*(I1-3.)+nu*(I1-3)*(I2-3)','mu=1.,nu=1.','mu=0.01, nu=0.01','mu=10., nu=10.')
mat = pmt.MatModel(model_I1I2)
print(mat, mat.parameters)
print(mat.stress(np.eye(3),incomp=True))
print(mat.stress(F,incomp=True))
Material model with 1 component:
Component1: ARB
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
mu_0 1.00 No 1.00e-02 10.00
nu_0 1.00 No 1.00e-02 10.00
------------------------------------------------------------------
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
[[2.27637222 0.55893333 0. ]
[0.55893333 2.32295 0. ]
[0. 0. 0. ]]
Spline based data-driven model
Instead of an analytical expression, if one has data of the strain energy density function, it can be interpolated with a spline and used as a material model. pymecht has splineI1 (a function of I1 alone) and splineI1I4 (a function of I1 and I4) implemented to achieve this. These can be used as follows:
[16]:
from scipy.interpolate import make_interp_spline, RectBivariateSpline
I1 = np.linspace(3,3.1,10)
Psi = np.tan(I1-3)
sp = make_interp_spline(I1,Psi) #create an interpolating spline
mat = pmt.MatModel('splineI1')
mat.models[0].set(sp) #set the spline for the model
print(mat,mat.parameters)
print(mat.stress(F,incomp=True))
Material model with 1 component:
Component1: splineI1
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
alpha_0 1.00 No -1.00e+01 10.00
------------------------------------------------------------------
[[1.10778312 0.30212267 0. ]
[0.30212267 1.13296001 0. ]
[0. 0. 0. ]]
As we can see, there is a parameters alpha_0 which is the coefficient of the spline. This can be used if the function is expressed as a linear combination of more than one splines. For example:
[17]:
sp2 = make_interp_spline(I1,np.sin(I1-3)) #create another interpolating spline
mat = pmt.MatModel('splineI1','splineI1')
models = mat.models
models[0].set(sp)
models[1].set(sp2)
params = mat.parameters
params.set('alpha_0',0.4)
params.set('alpha_1',0.6)
mat.parameters = params
print(mat, mat.parameters)
print(mat.stress(F,incomp=True))
Material model with 2 components:
Component1: splineI1
Component2: splineI1
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
alpha_0 0.40 No -1.00e+01 10.00
alpha_1 0.60 No -1.00e+01 10.00
------------------------------------------------------------------
[[0.6650294 0.18137165 0. ]
[0.18137165 0.6801437 0. ]
[0. 0. 0. ]]
With a spline, an interpolation is used within the range of defined values. Therefore, to track if we are extrapolating during the stress calculation, we can set the _warn flag to be True.
[18]:
models[0]._warn = True
models[1]._warn = True
print(mat.stress(F,incomp=True))
[[0.6650294 0.18137165 0. ]
[0.18137165 0.6801437 0. ]
[0. 0. 0. ]]
/usr/local/lib/python3.9/site-packages/pymecht/MatModel.py:946: UserWarning: Outside the training range; be careful interpreting the results 3.8899999999999997
3.0 3.1
warnings.warn(w)
The same procedure can be used for a bivariate spline that is a function of I1 and I4 defined on a rectangular grid as shown below.
[19]:
mat = pmt.MatModel('splineI1I4')
I4 = np.linspace(0.9,1.1,15)
I1grid, I4grid = np.meshgrid(I1,I4)
Psi = (I1grid-3) + (I4grid-1)**2 + (I4grid-1)*(I1grid-3)
sp = RectBivariateSpline(I1, I4, Psi.T, s=0)
mat.models[0].set(sp)
print(mat,mat.parameters)
mat.models[0].fiber_dirs = np.array([1,0,0])
print(mat.stress(F,incomp=True))
Material model with 1 component:
Component1: splineI1I4
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
alpha_0 1.00 No -1.00e+01 10.00
------------------------------------------------------------------
[[1.27222222 0.23333333 0. ]
[0.23333333 0.69166667 0. ]
[0. 0. 0. ]]