DAMASK_EICMD/python/damask/mechanics.py

308 lines
8.0 KiB
Python
Raw Normal View History

import numpy as np
def Cauchy(P,F):
2019-10-30 22:35:44 +05:30
"""
Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
2019-10-30 22:35:44 +05:30
Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
F : numpy.ndarray of shape (:,3,3) or (3,3)
Deformation gradient.
2020-03-13 05:00:49 +05:30
P : numpy.ndarray of shape (:,3,3) or (3,3)
First Piola-Kirchhoff stress.
2019-10-30 22:35:44 +05:30
"""
if np.shape(F) == np.shape(P) == (3,3):
sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T)
else:
sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F)
return symmetric(sigma)
2019-11-27 16:49:37 +05:30
2020-03-03 03:41:05 +05:30
def deviatoric_part(T):
2019-11-27 16:49:37 +05:30
"""
2020-02-15 18:40:16 +05:30
Return deviatoric part of a tensor.
2019-11-27 16:49:37 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the deviatoric part is computed.
2019-11-27 16:49:37 +05:30
"""
2020-03-03 03:41:05 +05:30
return T - np.eye(3)*spherical_part(T) if np.shape(T) == (3,3) else \
T - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[T.shape[0],3,3]),spherical_part(T))
2019-11-27 16:49:37 +05:30
2020-03-03 03:41:05 +05:30
def eigenvalues(T_sym):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return the eigenvalues, i.e. principal components, of a symmetric tensor.
2020-02-15 18:40:16 +05:30
The eigenvalues are sorted in ascending order, each repeated according to
its multiplicity.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvalues are computed.
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
return np.linalg.eigvalsh(symmetric(T_sym))
2020-03-03 03:41:05 +05:30
def eigenvectors(T_sym,RHS=False):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return eigenvectors of a symmetric tensor.
The eigenvalues are sorted in ascending order of their associated eigenvalues.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvectors are computed.
RHS: bool, optional
Enforce right-handed coordinate system. Default is False.
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
(u,v) = np.linalg.eigh(symmetric(T_sym))
if RHS:
2020-03-03 03:41:05 +05:30
if np.shape(T_sym) == (3,3):
if np.linalg.det(v) < 0.0: v[:,2] *= -1.0
else:
v[np.linalg.det(v) < 0.0,:,2] *= -1.0
2020-02-15 18:40:16 +05:30
return v
2020-03-03 03:41:05 +05:30
def left_stretch(T):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return the left stretch of a tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the left stretch is computed.
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
return __polar_decomposition(T,'V')[0]
2020-03-03 03:41:05 +05:30
def maximum_shear(T_sym):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return the maximum shear component of a symmetric tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the maximum shear is computed.
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
w = eigenvalues(T_sym)
return (w[0] - w[2])*0.5 if np.shape(T_sym) == (3,3) else \
2020-02-15 18:40:16 +05:30
(w[:,0] - w[:,2])*0.5
2019-10-19 12:21:51 +05:30
def Mises_strain(epsilon):
2019-10-30 22:35:44 +05:30
"""
Return the Mises equivalent of a strain tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
epsilon : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric strain tensor of which the von Mises equivalent is computed.
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
return __Mises(epsilon,2.0/3.0)
2020-02-15 18:40:16 +05:30
def Mises_stress(sigma):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return the Mises equivalent of a stress tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
sigma : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric stress tensor of which the von Mises equivalent is computed.
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
return __Mises(sigma,3.0/2.0)
def PK2(P,F):
2019-10-30 22:35:44 +05:30
"""
Calculate second Piola-Kirchhoff stress from first Piola-Kirchhoff stress and deformation gradient.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
P : numpy.ndarray of shape (:,3,3) or (3,3)
First Piola-Kirchhoff stress.
2020-03-13 05:00:49 +05:30
F : numpy.ndarray of shape (:,3,3) or (3,3)
Deformation gradient.
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
if np.shape(F) == np.shape(P) == (3,3):
S = np.dot(np.linalg.inv(F),P)
else:
S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P)
return symmetric(S)
2020-03-03 03:41:05 +05:30
def right_stretch(T):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return the right stretch of a tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the right stretch is computed.
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
return __polar_decomposition(T,'U')[0]
2020-03-03 03:41:05 +05:30
def rotational_part(T):
"""
2020-02-15 18:40:16 +05:30
Return the rotational part of a tensor.
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the rotational part is computed.
"""
2020-03-03 03:41:05 +05:30
return __polar_decomposition(T,'R')[0]
2020-03-03 03:41:05 +05:30
def spherical_part(T,tensor=False):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return spherical (hydrostatic) part of a tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the hydrostatic part is computed.
2020-02-15 18:40:16 +05:30
tensor : bool, optional
Map spherical part onto identity tensor. Default is false
2019-10-19 16:40:46 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
if T.shape == (3,3):
sph = np.trace(T)/3.0
2020-02-15 18:40:16 +05:30
return sph if not tensor else np.eye(3)*sph
else:
2020-03-03 03:41:05 +05:30
sph = np.trace(T,axis1=1,axis2=2)/3.0
2020-02-15 18:40:16 +05:30
if not tensor:
return sph
else:
2020-03-03 03:41:05 +05:30
return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(T.shape[0],3,3)),sph)
2019-10-25 17:00:20 +05:30
2020-02-15 18:40:16 +05:30
def strain_tensor(F,t,m):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return strain tensor calculated from deformation gradient.
For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
https://de.wikipedia.org/wiki/Verzerrungstensor
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
F : numpy.ndarray of shape (:,3,3) or (3,3)
Deformation gradient.
2020-02-15 18:40:16 +05:30
t : {V, U}
Type of the polar decomposition, V for left stretch tensor and U for right stretch tensor.
2020-02-15 18:40:16 +05:30
m : float
Order of the strain.
2019-10-25 17:00:20 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-03-17 16:52:48 +05:30
F_ = F.reshape(1,3,3) if F.shape == (3,3) else F
2020-02-15 18:40:16 +05:30
if t == 'V':
B = np.matmul(F_,transpose(F_))
w,n = np.linalg.eigh(B)
elif t == 'U':
C = np.matmul(transpose(F_),F_)
w,n = np.linalg.eigh(C)
2019-10-25 17:00:20 +05:30
2020-02-15 18:40:16 +05:30
if m > 0.0:
eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
- np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
elif m < 0.0:
eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
+ np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
else:
eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n))
2019-10-25 17:00:20 +05:30
2020-03-17 16:52:48 +05:30
return eps.reshape(3,3) if np.shape(F) == (3,3) else \
2020-02-15 18:40:16 +05:30
eps
2020-03-03 03:41:05 +05:30
def symmetric(T):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return the symmetrized tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the symmetrized values are computed.
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
return (T+transpose(T))*0.5
2020-03-03 03:41:05 +05:30
def transpose(T):
2019-10-30 22:35:44 +05:30
"""
2020-02-15 18:40:16 +05:30
Return the transpose of a tensor.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the transpose is computed.
2019-10-25 17:00:20 +05:30
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
return T.T if np.shape(T) == (3,3) else \
np.transpose(T,(0,2,1))
2019-10-25 17:00:20 +05:30
2020-03-03 03:41:05 +05:30
def __polar_decomposition(T,requested):
2019-10-30 22:35:44 +05:30
"""
Singular value decomposition.
2019-10-30 22:35:44 +05:30
Parameters
----------
2020-03-13 05:00:49 +05:30
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the singular values are computed.
requested : iterable of str
Requested outputs: R for the rotation tensor,
V for left stretch tensor and U for right stretch tensor.
2019-10-30 22:35:44 +05:30
"""
2020-03-03 03:41:05 +05:30
u, s, vh = np.linalg.svd(T)
R = np.dot(u,vh) if np.shape(T) == (3,3) else \
2019-10-30 22:35:44 +05:30
np.einsum('ijk,ikl->ijl',u,vh)
2019-10-30 22:35:44 +05:30
output = []
if 'R' in requested:
output.append(R)
if 'V' in requested:
2020-03-03 03:41:05 +05:30
output.append(np.dot(T,R.T) if np.shape(T) == (3,3) else np.einsum('ijk,ilk->ijl',T,R))
2019-10-30 22:35:44 +05:30
if 'U' in requested:
2020-03-03 03:41:05 +05:30
output.append(np.dot(R.T,T) if np.shape(T) == (3,3) else np.einsum('ikj,ikl->ijl',R,T))
2019-10-30 22:35:44 +05:30
return tuple(output)
2020-02-15 18:40:16 +05:30
2020-03-03 03:41:05 +05:30
def __Mises(T_sym,s):
2020-02-15 18:40:16 +05:30
"""
Base equation for Mises equivalent of a stres or strain tensor.
Parameters
----------
2020-03-13 05:00:49 +05:30
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the von Mises equivalent is computed.
2020-02-15 18:40:16 +05:30
s : float
Scaling factor (2/3 for strain, 3/2 for stress).
2020-03-03 03:41:05 +05:30
2020-02-15 18:40:16 +05:30
"""
2020-03-03 03:41:05 +05:30
d = deviatoric_part(T_sym)
return np.sqrt(s*(np.sum(d**2.0))) if np.shape(T_sym) == (3,3) else \
2020-02-15 18:40:16 +05:30
np.sqrt(s*np.einsum('ijk->i',d**2.0))