TFEL
, MFront
and MTest
The page declares the new functionalities of the 3.2 version of the
TFEL
project. This version has been released along with
version 3.1.3 and benefits from all the associated fixes. See this page for details.
The TFEL
project is a collaborative development of CEA and EDF
dedicated to material knowledge manangement with special focus on
mechanical behaviours. It provides a set of libraries (including
TFEL/Math
and TFEL/Material
) and several
executables, in particular MFront
and
MTest
.
TFEL
is available on a wide variety of operating systems
and compilers.
A new page dedicated to the python
bindings of the
TFEL
libraries is available here.
TFEL
librariesThe TFEL
project provides several libraries. This
paragraph is about updates made in those libraries.
Evaluator
classThe Evaluator
class is used to interpret textual
formula, as follows:
("sin(x)");
Evalutator ev.setVariableValue("x",12);
evconst auto s = ev.getValue();
getValue
methodIn the previous example, each variable value had to be set using the
setVariableValue
method. The new overloaded version of the
getValue
method can take a map as argument as follows:
("sin(x)");
Evalutator evconst auto s = ev.getValue({{"x",12}});
Operator()
Two overloaded versions of the Evaluator::operator()
has
been introduced as a synonyms for the getValue
method:
("sin(x)");
Evalutator evconst auto s = ev({{"x",12}});
getCxxFormula
methodThe getCxxFormula
method returns a string representing
the evaluation of the formula in standard C++
. This method
takes a map as argument which describes how certain variables shall be
represented. This method can be used, as follows:
("2*sin(x)");
Evalutator evstd::cout << ev({"x":"this->x"}}) << '\n';
The previous code displays:
sin((2)*(this->x))
This function is the basis of a new functionality of the
MFront
code generator (inline material properties), see
Section 3.4.1 for details.
The following new mathematical functions have been introduced:
exp2
: returns the base-2 exponential.expm1
: returns the base-e exponential minus one.cbrt
: returns the cubic rootlog2
: computes the base-2 logarithm of the
argment.log1p
: computes the logarithm of the argment plus
one.acosh
: computes the inverse of the hyperbolic
cosine.asinh
: computes the inverse of the hyperbolic
sine.atanh
: computes the inverse of the hyperbolic
tangent.erf
: computes the error function.erfc
: computes the complementary error function.tgamma
: computes the Gamma function.lgamma
: compute the natural logarithm of the gamma
function.hypot
: returns the hypotenuse of a right-angled
triangle whose legs are x and y, i.e. computes (.atan2
: returns the principal value of the arc tangent
of \(y/x\), expressed in radians.Let \(\underline{\sigma}\) be a stress tensor. Its deviatoric part \(\underline{s}\) is:
\[ \underline{s}=\underline{\sigma}-\displaystyle\frac{\displaystyle 1}{\displaystyle 3}\,{\mathrm{tr}\left(\underline{\sigma}\right)}\,\underline{I} =\left(\underline{\underline{\mathbf{I}}}-\displaystyle\frac{\displaystyle 1}{\displaystyle 3}\,\underline{I}\,\otimes\,\underline{I}\right)\,\colon\,\underline{\sigma} \]
The deviator of a tensor can be computed using the
deviator
function.
As it is a second order tensor, the stress deviator tensor also has a set of invariants, which can be obtained using the same procedure used to calculate the invariants of the stress tensor. It can be shown that the principal directions of the stress deviator tensor \(s_{ij}\) are the same as the principal directions of the stress tensor \(\sigma_{ij}\). Thus, the characteristic equation is
\[ \left| s_{ij}- \lambda\delta_{ij} \right| = -\lambda^3+J_1\lambda^2-J_2\lambda+J_3=0, \]
where \(J_1\), \(J_2\) and \(J_3\) are the first, second, and third deviatoric stress invariants, respectively. Their values are the same (invariant) regardless of the orientation of the coordinate system chosen. These deviatoric stress invariants can be expressed as a function of the components of \(s_{ij}\) or its principal values \(s_1\), \(s_2\), and \(s_3\), or alternatively, as a function of \(\sigma_{ij}\) or its principal values \(\sigma_1\), \(\sigma_2\), and \(\sigma_3\). Thus,
\[ \begin{aligned} J_1 &= s_{kk}=0,\, \\ J_2 &= \textstyle{\frac{1}{2}}s_{ij}s_{ji} = \displaystyle\frac{\displaystyle 1}{\displaystyle 2}{\mathrm{tr}\left(\underline{s}^2\right)}\\ &= \displaystyle\frac{\displaystyle 1}{\displaystyle 2}(s_1^2 + s_2^2 + s_3^2) \\ &= \displaystyle\frac{\displaystyle 1}{\displaystyle 6}\left[(\sigma_{11} - \sigma_{22})^2 + (\sigma_{22} - \sigma_{33})^2 + (\sigma_{33} - \sigma_{11})^2 \right ] + \sigma_{12}^2 + \sigma_{23}^2 + \sigma_{31}^2 \\ &= \displaystyle\frac{\displaystyle 1}{\displaystyle 6}\left[(\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_3 - \sigma_1)^2 \right ] \\ &= \displaystyle\frac{\displaystyle 1}{\displaystyle 3}I_1^2-I_2 = \frac{1}{2}\left[{\mathrm{tr}\left(\underline{\sigma}^2\right)} - \frac{1}{3}{\mathrm{tr}\left(\underline{\sigma}\right)}^2\right],\,\\ J_3 &= \det\left(\underline{s}\right) \\ &= \displaystyle\frac{\displaystyle 1}{\displaystyle 3}s_{ij}s_{jk}s_{ki} = \displaystyle\frac{\displaystyle 1}{\displaystyle 3} {\mathrm{tr}\left(\underline{s}^3\right)}\\ &= \displaystyle\frac{\displaystyle 1}{\displaystyle 3}(s_1^3 + s_2^3 + s_3^3) \\ &= s_1s_2s_3 \\ &= \displaystyle\frac{\displaystyle 2}{\displaystyle 27}I_1^3 - \displaystyle\frac{\displaystyle 1}{\displaystyle 3}I_1 I_2 + I_3 = \displaystyle\frac{\displaystyle 1}{\displaystyle 3}\left[{\mathrm{tr}\left(\underline{\sigma}^3\right)} - {\mathrm{tr}\left(\underline{\sigma}^2\right)}{\mathrm{tr}\left(\underline{\sigma}\right)} +\displaystyle\frac{\displaystyle 2}{\displaystyle 9}{\mathrm{tr}\left(\underline{\sigma}\right)}^3\right]. \end{aligned} \]
where \(I_{1}\), \(I_{2}\) and \(I_{3}\) are the invariants of \(\underline{\sigma}\).
\(J_{2}\) and \(J_{3}\) are building blocks for many isotropic yield critera. Classically, \(J_{2}\) is directly related to the von Mises stress \(\sigma_{\mathrm{eq}}\):
\[ \sigma_{\mathrm{eq}}=\sqrt{\displaystyle\frac{\displaystyle 3}{\displaystyle 2}\,\underline{s}\,\colon\,\underline{s}}=\sqrt{3\,J_{2}} \]
The first and second derivatives of \(J_{2}\) with respect to \(\sigma\) can be trivially implemented, as follows:
constexpr const auto id = stensor<N,real>::Id();
constexpr const auto id4 = st2tost2<N,real>::Id();
// first derivative of J2
const auto dJ2 = deviator(sig);
// second derivative of J2
const auto d2J2 = eval(id4-(id^id)/3);
In comparison, the computation of the first and second derivatives of
\(J_{3}\) with respect to \(\sigma\) are more cumbersome. In previous
versions TFEL
, one had to write:
constexpr const auto id = stensor<N,real>::Id();
constexpr const auto id4 = st2tost2<N,real>::Id();
const auto I1 = trace(sig);
const auto I2 = (I1*I1-trace(square(sig)))/2;
const auto dI2 = I1*id-sig;
const auto dI3 = computeDeterminantDerivative(sig);
const auto d2I2 = (id^id)-id4;
const auto d2I3 = computeDeterminantSecondDerivative(sig);
// first derivative of J3
const auto dJ3 = eval((2*I1*I1/9)*id-(I2*id+I1*dI2)/3+dI3);
// second derivative of J3
const auto d2J3 = eval((4*I1/9)*(id^id)-((id^dI2)+(dI2^id)+i1*d2I2)/3+d2I3);
More efficient implementations are now available using the
computeDeviatorDeterminantDerivative
and
computeDeviatorDeterminantSecondDerivative
functions:
// first derivative of J3
const auto dJ3 = computeDeviatorDeterminantDerivative(sig);
// second derivative of J3
const auto d2J3 = computeDeviatorDeterminantSecondDerivative(sig);
By definition, \(J_{2}\) and \(J_{3}\) are the second and third invariants of the deviatoric part \(\underline{s}\) of the stress tensor \(\underline{\sigma}\) (see also Section 2.1.2):
\[ \left\{ \begin{aligned} J_2 &= \displaystyle\frac{\displaystyle 1}{\displaystyle 2}{\mathrm{tr}\left(\underline{s}^2\right)}\\ J_3 &= \det(\underline{s}) \\ \end{aligned} \right. \]
The first and second derivatives of \(J_{2}\) with respect to the stress tensor \(\underline{\sigma}\) are trivially computed and implemented (see Section 2.1.2).
The first and second derivatives of \(J_{2}\) with respect to the stress tensor \(\underline{\sigma}\) can be computed respectively by:
computesJ3Derivative
function, which is a synonym
for the computeDeviatorDeterminantDerivative
function
defined in the tfel::math
namespace (see Section 2.1.2 for
details).computeJ3SecondDerivative
function, which is a
synonym for the computeDeviatorDeterminantSecondDerivative
function defined in the tfel::math
namespace (see
Section 2.1.2 for details).Within the framework of the theory of representation, generalizations to anisotropic conditions of the invariants of the deviatoric stress have been proposed by Cazacu and Barlat (see [1]):
Those invariants may be used to generalize isotropic yield criteria based on \(J_{2}\) and \(J_{3}\) invariants to orthotropy.
The following functions
\(J_{2}^{0}\), \(J_{3}^{0}\) and their first and second derivatives with respect to the stress tensor \(\underline{\sigma}\) can be computed by the following functions:
computesJ2O
, computesJ2ODerivative
and
computesJ2OSecondDerivative
.computesJ3O
, computesJ3ODerivative
and
computesJ3OSecondDerivative
.Those functions take the stress tensor as first argument and each
orthotropic coefficients. Each of those functions has an overload taking
the stress tensor as its firs arguments and a tiny vector
(tfel::math::tvector
) containing the orthotropic
coefficients.
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 pression.
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.python
bindingstfel.math
modulestensor
classThe following operations are supported:
The following functions have been introduced:
makeStensor1D
: builds a \(1D\) symmetric tensor from a tuple of three
values.makeStensor2D
: builds a \(2D\) symmetric tensor from a tuple of three
values.makeStensor3D
: builds a \(3D\) symmetric tensor from a tuple of three
values.tfel.material
moduleThe following functions are available:
buildFromPiPlane
: returns a tuple containing the three
eigenvalues of the stress corresponding to the given point in the \(\pi\)-plane.projectOnPiPlane
: projects a stress state, defined its
three eigenvalues or by a symmetric tensor, on the \(\pi\)-plane.The following script shows how to build an isosurface of the von Mises equivalent stress in the \(\pi\)-plane:
from math import pi,cos,sin
import tfel.math as tmath
import tfel.material as tmaterial
= 100
nmax for a in [pi*(-1.+(2.*i)/(nmax-1)) for i in range(0,nmax)]:
= tmath.makeStensor1D(tmaterial.buildFromPiPlane(cos(a),sin(a)))
s = tmath.sigmaeq(s);
seq *= 1/seq;
s = tmaterial.projectOnPiPlane(s);
s1,s2 print(s1,s2);
The computeHosfordStress
function, which compute the
Hosford equivalent stress, is available.
The following script shows how to print an iso-surface of the Hosford equivalent stress in the \(\pi\)-plane:
from math import pi,cos,sin
import tfel.math as tmath
import tfel.material as tmaterial
= 100
nmax for a in [pi*(-1.+(2.*i)/(nmax-1)) for i in range(0,nmax)]:
= tmath.makeStensor1D(tmaterial.buildFromPiPlane(cos(a),sin(a)))
s = tmaterial.computeHosfordStress(s,8,1.e-12);
seq /= seq;
s = tmaterial.projectOnPiPlane(s);
s1,s2 print(s1,s2);
The following functions are available:
makeBarlatLinearTransformation1D
: builds a \(1D\) linear transformation of the stress
tensor.makeBarlatLinearTransformation2D
: builds a \(2D\) linear transformation of the stress
tensor.makeBarlatLinearTransformation3D
: builds a \(3D\) linear transformation of the stress
tensor.computeBarlatStress
: computes the Barlat equivalent
Barlat stress.The following script shows how to print an iso-surface of the Barlat equivalent stress for the 2090-T3 aluminum alloy in the \(\pi\)-plane (see [2]):
from math import pi,cos,sin
import tfel.math as tmath
import tfel.material as tmaterial
= 100
nmax = tmaterial.makeBarlatLinearTransformation1D(-0.069888,0.079143,0.936408,
l1 0.524741,1.00306,1.36318,
0.954322,1.06906,1.02377);
= tmaterial.makeBarlatLinearTransformation1D(0.981171,0.575316,0.476741,
l2 1.14501,0.866827,-0.079294,
1.40462,1.1471,1.05166);
for a in [pi*(-1.+(2.*i)/(nmax-1)) for i in range(0,nmax)]:
= tmath.makeStensor1D(tmaterial.buildFromPiPlane(cos(a),sin(a)))
s = tmaterial.computeBarlatStress(s,l1,l2,8,1.e-12);
seq *= 1/seq;
s = tmaterial.projectOnPiPlane(s);
s1,s2 print(s1,s2);
MFront
StandardElastoViscoPlasticity
brickThis brick is used to describe a specific class of strain based behaviours based on an additive split of the total strain \(\underline{\epsilon}^{\mathrm{to}}\) into an elastic part \(\underline{\epsilon}^{\mathrm{el}}\) and an one or several inelastic strains describing plastic (time-independent) flows and/or viscoplastic (time-dependent) flows: \[ \underline{\epsilon}^{\mathrm{to}}=\underline{\epsilon}^{\mathrm{el}} +\sum_{i_{\mathrm{p}}=0}^{n_{\mathrm{p}}}\underline{\epsilon}^{\mathrm{p}}_{i_{\mathrm{p}}} +\sum_{i_{\mathrm{vp}}=0}^{n_{\mathrm{vp}}}\underline{\epsilon}^{\mathrm{vp}}_{i_{\mathrm{vp}}} \]
This equation defines the equation associated with the elastic strain \(\underline{\epsilon}^{\mathrm{el}}\).
The brick decomposes the behaviour into two components:
Here is a complete usage:
@Brick "StandardElastoViscoPlasticity" {
// Here the stress potential is given by the Hooke law. We define:
// - the elastic properties (Young modulus and Poisson ratio).
// Here the Young modulus is a function of the temperature.
// The Poisson ratio is constant.
// - the thermal expansion coefficient
// - the reference temperature for the thermal expansion
: "Hooke" {
stress_potential : "2.e5 - (1.e5*((T - 100.)/960.)**2)",
young_modulus : 0.3,
poisson_ratio : "1.e-5 + (1.e-5 * ((T - 100.)/960.) ** 4)",
thermal_expansion : 0
thermal_expansion_reference_temperature },
// Here we define only one viscplastic flow defined by the Norton law,
// which is based:
// - the von Mises stress criterion
// - one isotorpic hardening rule based on Voce formalism
// - one kinematic hardening rule following the Armstrong-Frederick law
: "Norton" {
inelastic_flow : "Mises",
criterion : "Voce" {R0 : 200, Rinf : 100, b : 20},
isotropic_hardening : "Armstrong-Frederick" {
kinematic_hardening : "1.e6 - 98500 * (T - 100) / 96",
C : "5000 - 5* (T - 100)"
D },
: "(4200. * (T + 20.) - 3. * (T + 20.0)**2)/4900.",
K : "7. - (T - 100.) / 160.",
n : 3
Ksf }
};
cyrano
interfaceCyrano
is the state of the art fuel performance code
developed by EDF (See [3, 4]).
Cyrano
to finite strain analysisAs Cyrano, provides a mono-dimensional description of the fuel rod, its extension to finite strain follows the same treatment used for:
MTest
to pipes.Alcyone1D
to finite strain
analysis.Miehe et al. have introduced a finite strain framework based on the
Hencky strain measure (logarithmic strain) and its dual stress (see
[5]).
This framework has been introduced in Code_Aster
in 2011
(see [6]).
A behaviour is declared to follow this framework by using the
@StrainMeasure
keyword:
@StrainMeasure Hencky;
The cyrano
interface now provides support for behaviours
based on this strain measure.
The interface handles (See [7] for details):
The generic behaviour interface has been introduced to provide an interface suitable for most needs.
Whereas other interfaces target a specific solver and thus are
restricted by choices made by this specific solver, this interface has
been created for developers of homebrew solvers who are able to modify
their code to take full advantage of MFront
behaviours.
This interface is tightly linked with the
MFrontGenericInterfaceSupport
project which is available on
github
: https://github.com/thelfer/MFrontGenericInterfaceSupport.
This project has a more liberal licence which allows it to be included
in both commercial and open-source solvers/library. This licensing
choice explains why this project is not part of the
TFEL
.
Various keywords (such as @ElasticMaterialProperties
,
@ComputeThermalExpansion
, @HillTensor
, etc.)
expects one or more material properties. In previous versions, those
material properties were constants or defined by an external
MFront
file.
This new version allows those material properties to be defined by formulae, as follows:
@Parameter E0 =2.1421e11,E1 = -3.8654e7,E2 = -3.1636e4;
@ElasticMaterialProperties {"E0+(T-273.15)*(E1+E2*(T-273.15))",0.3}
As for material properties defined in external MFront
files, the material properties evaluated by formulae will be computed
for updated values of their parameters. For example, if the previous
lines were used in the Implicit
DSL, two variables
young
and young_tdt
will be automatically made
available:
young
variable will be computed for the temperature
\(T+\theta\,\Delta\,T\).young_tdt
variable will be computed for the
temperature \(T+\Delta\,T\).--list-behaviour-bricks
argumentThe --list-behaviour-bricks
argument returns the list of
all available behaviour bricks.
$ mfront --list-behaviour-bricks
DDIF2 (undocumented)
FiniteStrainSingleCrystal (undocumented)
StandardElasticity (documented)
StandardElastoViscoPlasticity (undocumented)
--help-behaviour-brick
argumentThe --help-behaviour-brick
argument returns the help
associated with a given behaviour brick.
$ mfront --help-behaviour-brick=StandardElasticity
The `StandardElasticity` brick describes the linear elastic part of
the behaviour of an isotropic or orthotropic material.
StandardElastoViscoplasticity
brickThe StandardElastoViscoplasticity
brick is based on the
concepts of:
For each concept, one can query the list of concrete implementations of this concept and the associated help, as in the following example:
$ mfront --list-stress-potentials
- DDIF2 (documented)
- Hooke (documented)
- IsotropicDamage (undocumented)
$ mfront --help-stress-potential=Hooke
The `Hooke` brick describes the linear elastic part of
the behaviour of an isotropic or orthotropic material.
.....
The following arguments are thus avaiable:
--list-isotropic-hardening-rules
--list-kinematic-hardening-rules
--list-stress-criteria
--list-stress-potentials
--help-isotropic-hardening-rule
--help-kinematic-hardening-rule
--help-stress-criterion
--help-stress-potential
MTest
Defining an event is a way to active/desactivate a constraint, (see
@ImposedStrain
, @ImposedCohesiveForce
,
@ImposedOpeningDisplacement
, @ImposedStrain
,
@ImposedDeformationGradient
,
@NonLinearConstraint
).
The @Event
keyword is followed by the name of the event
and a time or a list of times (given as an array of values).
@Event 'Stop' 1;
At the end of the definition of a constraint, one may now optionnally
set options on the active state of this contraint at the beginning of
the computation, the activating events and desactivating events using a
JSON-like structure. This structure starts with an opening curly brace
({
) and ends with a closing curly brace (}
).
An option is given by its name, an double-dot character (:
)
and the value of the option. Consecutive options are separated by a
comma ,
(see below for an example). The following options
are available:
active
: states if the constraint is active at the
beginning of the computation. This is a boolean, so the expected value
are either true
or false
.activating_event
: gives the name of the event which
activate the constraint. An array of string is expected.activating_events
: gives the list of events which
activate the constraint. An array of string is expected.desactivating_event
: gives the name of the event which
desactivate the constraint. An array of string is expected.desactivating_events
: gives the list of events which
desactivate the constraint. An array of string is expected.@ImposedStrain<evolution> 'EXX' {0.:0.,1:1e-3}{
: 'Stop'
desactivating_event };
The @UserDefinedPostProcessing
lets the user define is
own post-processings.
This keywords is followed by:
@UserDefinedPostProcessing 'myoutput.txt' {'SXX','EquivalentPlasticStrain'};
Mfront/Python
under Windows
When compiling MFront
with python bindings support using
Visual Studio 2015 Community
, the following error message
is triggered:
C:\Codes\tfel\sources\tfel\include\TFEL/Material/PiPlane.ixx(60): error C2131: l'expression n'a pas été évaluée en co
nstante [C:\Codes\tfel\sources\build-vs-2015\bindings\python\tfel\py_tfel_material.vcxproj]
0 Avertissement(s)
1 Erreur(s)
The problem is that constexpr
variables are poorly
handled by Visual Studio 2015
. The constexpr
macro was introduced to handle this compiler (it boils down to
constexpr
for others compilers). The issue can be solved by
changing the faulty line as as follows:
constexpr const auto isqrt6 = isqrt2 * isqrt3;
For more details, see: https://sourceforge.net/p/tfel/tickets/137/
The @UserDefinedPostProcessing
lets the user define is
own post-processings.
This keywords is followed by:
For more details, see: https://sourceforge.net/p/tfel/tickets/112/
The following header files have been renamed:
Hosford.hxx
has been moved in
Hosford1972YieldCriterion.hxx
.Barlat.hxx
has been moved in
Barlat2004YieldCriterion.hxx
.