Equivalent stress

Hosford equivalent stress

The header TFEL/Material/Hosford1972YieldCriterion.hxx introduces three functions which are meant to compute the Hosford equivalent stress and its first and second derivatives. This header is automatically included by MFront

The Hosford equivalent stress is defined by: \[ \sigma_{\mathrm{eq}}^{H}=\sqrt[a]{\displaystyle\frac{\displaystyle 1}{\displaystyle 2}\left({\left|\sigma_{1}-\sigma_{2}\right|}^{a}+{\left|\sigma_{1}-\sigma_{3}\right|}^{a}+{\left|\sigma_{2}-\sigma_{3}\right|}^{a}\right)} \] where \(s_{1}\), \(s_{2}\) and \(s_{3}\) are the eigenvalues of the stress.

Therefore, when \(a\) goes to infinity, the Hosford stress reduces to the Tresca stress. When \(n = 2\) the Hosford stress reduces to the von Mises stress.

The following functions has been implemented:

• computeHosfordStress: return the Hosford equivalent stress
• computeHosfordStressNormal: return a tuple containing the Hosford equivalent stress and its first derivative (the normal)
• computeHosfordStressSecondDerivative: return a tuple containing the Hosford equivalent stress, its first derivative (the normal) and the second derivative.

Example

The following example computes the Hosford equivalent stress, its normal and second derivative:

stress seq;
Stensor  n;
Stensor4 dn;
std::tie(seq,n,dn) = computeHosfordStressSecondDerivative(s,a,seps);

In this example, s is the stress tensor, a is the Hosford exponent, seps is a numerical parameter used to detect when two eigenvalues are equal.

If C++-17 is available, the previous code can be made much more readable:

const auto [seq,n,dn] = computeHosfordStressSecondDerivative(s,a,seps);

Barlat equivalent stress

The header TFEL/Material/Barlat2004YieldCriterion.hxx introduces various functions which are meant to compute the Barlat equivalent stress and its first and second derivatives. This header is automatically included by MFront for orthotropic behaviours.

The Barlat equivalent stress is defined as follows (see [1]): \[ \sigma_{\mathrm{eq}}^{B}= \sqrt[a]{ \frac{1}{4}\left( \sum_{i=0}^{3} \sum_{j=0}^{3} {\left|s'_{i}-s''_{j}\right|}^{a} \right) } \]

where \(s'_{i}\) and \(s''_{i}\) are the eigenvalues of two transformed stresses \(\underline{s}'\) and \(\underline{s}''\) by two linear transformation \(\underline{\underline{\mathbf{L}}}'\) and \(\underline{\underline{\mathbf{L}}}''\): \[ \left\{ \begin{aligned} \underline{s}' &= \underline{\underline{\mathbf{L'}}} \,\colon\,\underline{\sigma}\\ \underline{s}'' &= \underline{\underline{\mathbf{L''}}}\,\colon\,\underline{\sigma}\\ \end{aligned} \right. \]

The linear transformations \(\underline{\underline{\mathbf{L}}}'\) and \(\underline{\underline{\mathbf{L}}}''\) are defined by \(9\) coefficients (each) which describe the material orthotropy. There are defined through auxiliary linear transformations \(\underline{\underline{\mathbf{C}}}'\) and \(\underline{\underline{\mathbf{C}}}''\) as follows: \[ \begin{aligned} \underline{\underline{\mathbf{L}}}' &=\underline{\underline{\mathbf{C}}}'\,\colon\,\underline{\underline{\mathbf{M}}} \\ \underline{\underline{\mathbf{L}}}''&=\underline{\underline{\mathbf{C}}}''\,\colon\,\underline{\underline{\mathbf{M}}} \end{aligned} \] where \(\underline{\underline{\mathbf{M}}}\) is the transformation of the stress to its deviator: \[ \underline{\underline{\mathbf{M}}}=\underline{\underline{\mathbf{I}}}-\displaystyle\frac{\displaystyle 1}{\displaystyle 3}\underline{I}\,\otimes\,\underline{I} \]

The linear transformations \(\underline{\underline{\mathbf{C}}}'\) and \(\underline{\underline{\mathbf{C}}}''\) of the deviator stress are defined as follows: \[ \underline{\underline{\mathbf{C}}}'= \begin{pmatrix} 0 & -c'_{12} & -c'_{13} & 0 & 0 & 0 \\ -c'_{21} & 0 & -c'_{23} & 0 & 0 & 0 \\ -c'_{31} & -c'_{32} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & c'_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & c'_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & c'_{66} \\ \end{pmatrix} \quad \text{and} \quad \underline{\underline{\mathbf{C}}}''= \begin{pmatrix} 0 & -c''_{12} & -c''_{13} & 0 & 0 & 0 \\ -c''_{21} & 0 & -c''_{23} & 0 & 0 & 0 \\ -c''_{31} & -c''_{32} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & c''_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & c''_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & c''_{66} \\ \end{pmatrix} \]

The following functions have been implemented:

• computeBarlatStress: return the Barlat equivalent stress
• computeBarlatStressNormal: return a tuple containing the Barlat equivalent stress and its first derivative (the normal)
• computeBarlatStressSecondDerivative: return a tuple containing the Barlat equivalent stress, its first derivative (the normal) and the second derivative.

Linear transformations

To define the linear transformations, the makeBarlatLinearTransformation function has been introduced. This function takes two template parameter:

• the space dimension (\(1\), \(2\), and \(3\))
• the numeric type used (automatically deduced)

This functions takes the \(9\) coefficients as arguments, as follows:

const auto l1 = makeBarlatLinearTransformation<3>(c_12,c_21,c_13,c_31,
                                                  c_23,c_32,c_44,c_55,c_66);

Note In his paper, Barlat and coworkers uses the following convention for storing symmetric tensors:

\[ \begin{pmatrix} xx & yy & zz & yz & zx & xy \end{pmatrix} \]

which is not consistent with the TFEL/Cast3M/Abaqus/Ansys conventions:

\[ \begin{pmatrix} xx & yy & zz & xy & xz & yz \end{pmatrix} \]

Therefore, if one wants to use coefficients \(c^{B}\) given by Barlat, one shall call this function as follows:

const auto l1 = makeBarlatLinearTransformation<3>(cB_12,cB_21,cB_13,cB_31,
                                                  cB_23,cB_32,cB_66,cBB_55,cBB_44);

The TFEL/Material library also provide an overload of the makeBarlatLinearTransformation which template parameters are the modelling hypothesis and the orthotropic axis conventions. The purpose of this overload is to swap appropriate coefficients to get a consistent definition of the linear transformations for all the modelling hypotheses.

General functionalities

Computation of the inverse of the Langevin function

The inverse Langevin function is used in many statistically based network behaviours describing rubber-like materials.

The Langevin function \(\mathcal{L}\) is defined as follows: \[ \mathcal{L}\left(x\right)=\displaystyle\frac{\displaystyle 1}{\displaystyle \coth\left(x\right)}-\displaystyle\frac{\displaystyle 1}{\displaystyle x} \]

The complexity of the inverse Langevin function \(\mathcal{L}^{-1}\left(x\right)\) motivated the development of various approximations [2–4].

Figure 1 compares those approximations. The approximations of Bergström and Boyce [3] and Jedynak [4] are undistinguishable. See Jedynak for a quantitative discussion of the error generated by those approximations [4]. It is worth noting that all those approximations mostly differs near the pole of inverse Langevin function \(\mathcal{L}^{-1}\left(x\right)\).

The InverseLangevinFunctionApproximations enumeration lists the approximations that have been implemented and that can be evaluated in a constexpr context:

• The value COHEN_1991 is associated with the approximation proposed by Cohen [2]: \[ \mathcal{L}^{-1}\left(x\right) \approx y\displaystyle\frac{\displaystyle 3-y^{2}}{\displaystyle 1-y^{2}}. \]
• The value JEDYNAK_2015 is associated with the approximation proposed by Jedynak [4]: \[ \mathcal{L}^{-1}\left(x\right) \approx y \, \displaystyle\frac{\displaystyle c_{0} + c_{1}\,y + c_{2}\,y^{2}}{\displaystyle 1 + d_{1}\,y + d_{2}\,y^{2}} \]
• KUHN_GRUN_1942 or MORCH_2022 are associated with a taylor expansion of \(\mathcal{L}^{-1}\left(x\right)\) at \(0\): \[ \mathcal{L}^{-1}\left(x\right) \approx y\,P\left(y^{2}\right) \quad\text{with}\quad P\left(y^{2}\right)=\sum_{i=0}^{9}c_{i}\,\left(y^{2}\right)^{i} \] where \(P\) is a \(9\)th order polynomial. Hence, the Taylor expression is of order \(19\).

The computeApproximateInverseLangevinFunction computes one approximation of the inverse Langevin function and the computeApproximateInverseLangevinFunctionAndDerivative function computes an approximation of the inverse Langevin function and its derivative. These functions have two template parameters: the approximation selected and the numeric type to be used. By default, the JEDYNAK_2015 approximation is used and the numeric type can be deduced from the type of the argument.

The approximation proposed by Bergström and Boyce [3] is given by the following function: \[ \mathcal{L}^{-1}\left(x\right) \approx \left\{ \begin{aligned} c_{1} \tan\left(c_{2} \, x\right) + c_{3} \, x &\quad\text{if}\quad \left|x\right| \leq c_{0}\\ \displaystyle\frac{\displaystyle 1}{\displaystyle \mathop{sign}\left(x\right)-x}&\quad\text{if}\quad \left|x\right| > c_{0} \end{aligned} \right. \]

The computeBergstromBoyce1998ApproximateInverseLangevinFunction function computes this approximation and the computeBergstromBoyce1998ApproximateInverseLangevinFunctionAndDerivative function computes this function and its derivative. These functions can’t be declared constexpr because of tangent function is not constexpr. These functions have one template parameter, the numeric type to be used. This template parameter can be automatically deduced from the type of the argument.

Example of usage

using ApproximationFunctions = InverseLangevinFunctionApproximations;
// compute Cohen's approximation of the inverse Langevin function
const auto v = computeApproximateInverseLangevinFunction<
            ApproximationFunctions::COHEN_1991>(y)
// compute Jedynak's approximation of the inverse Langevin function and its derivative
const auto [f, df] = computeApproximateInverseLangevinFunctionAndDerivative(y)

\(\pi\)-plane

The \(\pi\)-plane is defined in the space defined by the three eigenvalues \(S_{0}\), \(S_{1}\) and \(S_{2}\) of the stress by the following equations: \[ S_{0}+S_{1}+S_{2}=0 \]

This plane contains deviatoric stress states and is perpendicular to the hydrostatic axis. A basis of this plane is given by the following vectors: \[ \vec{n}_{0}= \frac{1}{\sqrt{2}}\, \begin{pmatrix} 1 \\ -1 \\ 0 \end{pmatrix} \quad\text{and}\quad \vec{n}_{1}= \frac{1}{\sqrt{6}}\, \begin{pmatrix} -1 \\ -1 \\ 2 \end{pmatrix} \]

This plane is used to characterize the iso-values of equivalent stresses which are not sensitive to the hydrostatic pressure.

Various functions are available:

• projectOnPiPlane: this function projects a stress state on the \(\pi\)-plane.
• buildFromPiPlane: this function builds a stress state, defined by its three eigenvalues, from its coordinate in the \(\pi\)-plane.

Orthotropic axes convention

Most finite element solver can’t have a unique definition of the orthotropic axes valid for all the modelling hypotheses.

For example, one can define a pipe using the following axes definition:

• \(\left(rr,zz,tt,...\right)\) in \(3D\), \(2D\) axysymmetric, \(1D\) axisymmetric generalised plane strain or \(1D\) axisymmetric generalised plane stress.
• \(\left(rr,tt,zz,...\right)\) in \(2D\) plane stress, \(2D\) plane strain and \(2D\) generalized plane strain.

With those conventions, named Pipe in MFront, the axial direction is either the second or the third material axis, a fact that must be taken into account when defining the stiffness tensor, the Hill tensor(s), the thermal expansion, etc.

This convention is only valid for \(3D\), \(2D\) axysymmetric, \(1D\) axisymmetric generalised plane strain or generalised plane stress. - \(\left(rr,tt,zz,...\right)\) in \(2D\) plane stress, strain, generalized plane strain.

If we were to model plates, an appropriate convention is the following:

• The first material axis is the rolling direction
• The second material axis is the in plane direction perpendicular to the rolling direction (transverse direction).
• The third material axis is the normal to the plate.

By definition, this convention, named Plate in MFront is only valid for \(3D\), \(2D\) plane stress, \(2D\) plane strain and \(2D\) generalized plane strain modelling hypotheses.

Homogeneisation

The homogenization functions are part of the namespace tfel::material::homogenization. A specialization for elasticity is defined: tfel::material::homogenization::elasticity.

Eshelby tensors

The header Eshelby.hxx introduces the function computeEshelbyTensor which computes the Eshelby tensor of an ellipsoid. If we consider a constant stress-free strain \(\underline{\varepsilon}^\mathrm{T}\) filling an ellipsoidal volume embedded in an infinite homogeneous medium whose elasticity is \(\underline{C}_0\), the strain tensor inside the ellipsoid is given by

\(\underline{\varepsilon}=\underline{S}_0:\underline{\varepsilon}^\mathrm{T}\).

where \(\underline{S}_0\) is the Eshelby tensor. Note that it is related to the Hill tensor (P_0) by

\(\underline{P}_0=\underline{S}_0:\underline{C}_0^{-1}\)

which gives the strain tensor inside the ellipsoid as a function of the polarization tensor \(\underline{\tau }= -\underline{C}_0:\underline{\varepsilon}^\mathrm{T}\) :

\(\underline{\varepsilon}=-\underline{P}_0:\underline{\tau}\).

The function computeEshelbyTensor computes the Eshelby tensor of an ellipsoid whose semi-axis lengths are a, b, c, embedded in an isotropic matrix. It returns an object of type st2tost2<3u,real>, which is the fourth-order Eshelby tensor, in a basis which is adapted to the ellipsoid.

There is also computeSphereEshelbyTensor, computeAxisymmetricalEshelbyTensor, and also computeCircularCylinderEshelbyTensor and computeEllipticCylinderEshelbyTensor for plane strain elasticity.

The expressions can be found in [5] for the axisymmetrical ellipsoid and in [6] for other cases.

When two axes are very close, the formulas for three different axes are numerically instable, hence a parameter is introduced to switch to the formulas suited for the perfect axisymmetrical case. This parameter can be modified by the user, it is called precf when using float, precd for double, and precld. for long double. In the same way, the formulas for the axisymmetrical case are instable when the aspect ratio is near one, so a parameter allows to switch to the formula for a sphere.

Strain localisation tensors

The header Eshelby.hxx also introduces three functions that compute the strain localisation tensor of an ellipsoid. If we consider an ellipsoid whose elasticity is \(\underline{C}_i\), embedded in an infinite homogeneous medium whose elasticity is \(\underline{C}_0\), submitted to a external uniform strain field at infinity \(\underline{E}\), the strain field within the ellipsoid is uniform and given by

\(\underline{\varepsilon }= \underline{A}:\underline{E}\)

where \(\underline{A}\) is the localisation tensor.

Three functions are implemented for the different possible shapes : computeEllipsoidLocalisationTensor, computeAxisymmetricalEllipsoidLocalisationTensor and computeSphereLocalisationTensor. The ellipsoid is parametrized by its semi-axis lengths \(a,b,c\) but also by its axis orientations. The functions then return the localisation tensors taking into account the orientations.

Homogenization schemes

Different schemes are implemented and return the homogenized stiffness of the material. The scheme available are Mori-Tanaka scheme and dilute scheme. The available functions are computeMoriTanakaScheme, computeDiluteScheme, computeSphereDiluteScheme, computeSphereMoriTanakaScheme.

If a distribution of ellipsoids is considered, three types of distributions are considered : isotropic, transverse isotropic and with unique orientation. The corresponding functions are computeIsotropicDiluteScheme, computeTransverseIsotropicDiluteScheme, computeOrientedDiluteScheme, computeIsotropicMoriTanakaScheme, computeTransverseIsotropicMoriTanakaScheme and computeOrientedMoriTanakaScheme.