Overview

Versatility of general purpose mechanical solver mainly relies on its ability to let the user define the material behaviour. MFront provides a high level language to write mechanical behaviours.

It can be compared to the ZebFront tool developped by the Centre des Matériaux de Mines ParisTech as part of the Zset software (see R. Foerch et al. 1997; J. Besson et al. 1998; Northwest Numerics and Modeling 2014). One major difference between ZebFront and MFront is the programming techniques used: ZebFront mostly relies on object oriented techniques where MFront relies on generic programming leading to almost orthogonal design choices.

Three kind of mechanical behaviours are currently considered with MFront: small and finite strain behaviours and cohesive zone models.

Mechanical behaviour role

We now precise the role of the mechanical behaviours in standard displacement-based finite element solver (see Zienkiewicz 1977; J. Besson, Cailletaud, and Chaboche 2001; EDF 2013). For the sake of simplicity, we only treat the case of small strain behaviours for the rest of this document.

At each time step, the following resolution procedure is used:

  1. a prediction of the displacement is made. Such a prediction may use the derivate of the stress tensor with respect of the strain tensor \({{\displaystyle \frac{\displaystyle \partial {\underline{\sigma}}}{\displaystyle \partial {\underline{\epsilon}^{\mathrm{to}}}}}}\) or an approximation of it. This prediction step will not be discussed in this article but can also be handled by behaviour implementations made with MFront;
  2. an iterative process similar to the Newton-Raphson algorithm used to find a displacement that will satisfy the mechanical equilibrium. Given an estimation of the displacement at the end the time step, one computes at each integration point an estimation of the increment of the deformation tensor \(\Delta\,{\underline{\epsilon}^{\mathrm{to}}}\). The mechanical behaviour is then called and provides an associated estimation of the stress tensor \({\underline{\sigma}}\) and the values of some internal state variables \(Y_{i}\). In the solver requires it, the mechanical behaviour may also provide an estimation of the tangent consistant operator \({{\displaystyle \frac{\displaystyle \partial \Delta\,{\underline{\sigma}}}{\displaystyle \partial \Delta\,{\underline{\epsilon}^{\mathrm{to}}}}}}\) (see J. C. Simo and Taylor (1985)) which is used to estimate a more accurate displacement field.

A mechanical behaviour can thus be viewed as functional:

\[ {\left({\left.{\underline{\sigma}}\right|_{t+\Delta\,t}},{\left.Y_{i}\right|_{t+\Delta\,t}},{{\displaystyle \frac{\displaystyle \partial \Delta\,{\underline{\sigma}}}{\displaystyle \partial \Delta\,{\underline{\epsilon}^{\mathrm{to}}}}}}\right)}= \mathcal{F}{\left({\left.{\underline{\sigma}}\right|_{t}},{\left.Y_{i}\right|_{t}},\Delta\,{\underline{\epsilon}^{\mathrm{to}}},\Delta\,t,\ldots\right)} \]

The dots \(\ldots\) means that the behaviour may also depend of external state variables evolutions with time, namely the temperature, the irradiation damage, and so on.

Isotropic \(J_{2}\) plastic/vsicplastic behaviours. Example of finite strain pre- and post-processing

\(4\) domain specific languages address the case of small strain isotropic \(J_{2}\) plastic and/or viscoplastic behaviours which are of common use and for which efficient implicit scalar radial return mapping algorithms exist (see Juan C Simo and Hughes (1998)).

The following listing shows how a simple plastic behaviour can be implemented:

1
2
3
4
5
6
7
8
9
@Parser IsotropicPlasticMisesFlow; //< domain specific language
@Behaviour Plasticity;             //< name of the behaviour
@Parameter  H  = 22e9;             //< hardening slope
@Parameter s0 = 200e6;             //< elasticity limit
@FlowRule{                         //< flow rule
  f       = seq-H*p-s0;
  df_dseq = 1;
  df_dp   = -H;
}

The plastic behaviour is governed by the following yield surface (see J. Besson, Cailletaud, and Chaboche (2001);Chaboche et al. (2009)): \[ f{\left({\sigma_{\mathrm{eq}}},p\right)}={\sigma_{\mathrm{eq}}}-H\,p-\sigma_{0} \leq 0 \] where \({\sigma_{\mathrm{eq}}}\) is the Von Mises stress, \(p\) the equivalent plastic strain, \(R\) the isotropic hardening slope and \(\sigma_{0}\) the initial elasticity limit.

The generated code represent a total amount of \(1\,512\) lines of ̀C++ code and provides:

  1. optimised implementations of the behaviour for various modelling hypotheses (axisymmetrical generalised plane strain, plane strain, plane stress, generalised plane strain, axisymmetry, tridimensional) thanks to template metaprogramming and template specialisations. An small overview of the programming techniques used can be found in~.
  2. the computation of the prediction operator;
  3. meta-data about the required material properties, the number of states variables, etc... that can be retrieved dynamically (For the sake of simplicity, no glossary name was specified in this example). This mechanism is used by applications to appropriately call the behaviour;
  4. the computation of a tangent matrix operator (various choice are possible: elastic, secant, consistent);
  5. dynamically loadable functions allowing the user to change various parameters of the behaviour (the convergence criterion of the implicit algorithm, the \(\theta\) parameter of the implicit algorithm, the hardening slope \(H\) and the initial elasticity limit \(\sigma_{0}\), etc.). Those functions by-pass the standard behaviour call and are an extremely light-weight manner to dynamically modify a behaviour (almost no runtime cost).

Local divergence

Local divergence of the implicit algorithm can be handled through an appropriate substepping procedure. This feature is not enabled by default, but appropriate keywords gives to the end user explicit control on this procedure.

Finite strain strategies

If not handled directly by the calling code, appropriate pre- and post-processings allowing the use of small strain behaviours in finite strain computations can be generated. Two lagrangian finite strain strategies are currently available:

  1. finite rotations, small strains. This method allows the re-use a behaviour whose material parameters \(H\) and \(\sigma_{0}\) were identified through small strain computations in the context of finite rotations without any re-identification. The physical meaning of the pre- and post-processing stages are discussed by Doghri (see Doghri (2000));
  2. lagrangian logarithmic strains as proposed by Miehe, Apel and Lambrecht (see Miehe, Apel, and Lambrecht (2002);edf_modeles_2013). The use of the logarithmic strains has several advantages:
  1. it preserves the small-strain classical meaning of the variables, which is truly appreciated by engineers;
  2. it may can also be used for arbitrary complex models (kinematic hardening, initial or induced anisotropy, etc.).

The following figure shows how our example can be used to model a notched specimen under a tensile test:

Von Mises stress

Von Mises stress

In this case, the material parameters \(H\) and \(\sigma_{0}\) of the behaviour must be identified on tests implying finite strains. As an additional remark, the results found using logarithmic strains were remarkably closed to those obtained by the classical formulation based on an \(F_{e}\,F_{p}\) decomposition proposed by Simo and Miehe (see J. Simo and Miehe (1992)) (that was also implemented using MFront).

Generic domain specific languages

Apart from the domain specific languages dealing with isotropic \(J_{2}\) plastic and viscoplastic behaviours presented in the previous paragraph, MFront also provides several general-purpose domain specific languages:

  1. the Default domain specific language allows the user the write its own integration algorithm. This is very useful for explicit behaviours such as the classical cohesive zone model (see Tvergaard 1990).
  2. the Runge-Kutta domain specific language allows the user to write the constitutive equations given as a system of ordinary time differential equations. Using those algorithms is generally less efficient than using implicit integration to be described. Various algorithms are however available.
  3. the Implicit domain specific language allows the user to perform the local integration using an implicit algorithm. An introduction to those algorithms is given in the next paragraph.

Example of a cohesive zone model

The implementation of the Tvergaard cohesive zone model using the Default domain specific language is given below:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
@Parser DefaultCZMParser;     // domain specific language
@Behaviour Tvergaard;         // name of the behaviour
@MaterialProperty stress kn;  // normal stiffness
@MaterialProperty stress ks;  // tangential elastic stiffness
@MaterialProperty stress smax;// maximal stress
@MaterialProperty real delta; // maximal normal opening displacement
@StateVariable real d;        // damage variable
@Integrator{
  const real C = real(27)/real(4);
  t_t = ks*(u_t+du_t);        // tangential behaviour
  if(u_n+du_n<0){             // normal behaviour in compression
    t_n = kn*(u_n+du_n);
  } else {                    // normal behaviour in traction
    const real rod = (u_n+du_n)/delta;   // reduced opening displacement
    const real d_1 = d;             // previous damage
    d   = min(max(d,rod),0.99);     // damage indicator
    const real K1 = C*smax/delta;   // initial stiffness
    const real K  = K1*(1-d)*(1-d); // secant stiffness
    t_n = K*(u_n+du_n);
  }
} // end of @Integrator

Details about the computation of the consistent tangent operator were eluded. The opening displacement \(\vec{u}\) is automatically decomposed into the normal opening displacement \(u_{n}\) and its tangential opening displacement \(\vec{u}_{t}\)

Explicit algorithm example

The implementation of a generalisation of the Norton creep law for anisotropic materials is given below:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
@Parser    RungeKutta;           // domain specific language
@Behaviour OrthotropicCreep;     // name of the behaviour
@OrthotropicBehaviour;           // treating an orthotropic behaviou
@RequireStiffnessTensor;         // requires the stiffness tensor to be computed
@StateVariable Stensor evp;      // viscoplastic strain
@StateVariable strain p;         // Equivalent viscoplastic strain
@Includes{                       // header files
#include<TFEL/Material/Hill.hxx>
}
@ComputeStress{                  /* stress computation */
  sig = D*eel;
}
@Derivative{                     /* constitutive equations */
  st2tost2<N,real> H;            // Hill Tensor
  H = hillTensor<N,real>(0.371,0.629,4.052,1.5,1.5,1.5);
  stress sigeq = sqrt(sig|H*sig);  // equivalent Hill stress
  if(sigeq>1e9){                 // automatic sub-stepping
    return false;
  }
  Stensor  n(0.);                // flow direction
  if(sigeq > 10.e-7){
    n    = H*sig/sigeq;
  }
  dp   = 8.e-67*pow(sigeq,8.2); // evolution of p
  devp = dp*n;                  // evolution of the viscoplastic strains
  deel = deto - devp;           // evolution of the elastic strains
}

Integration is performed in the material referential. The elastic strain state variable \({\underline{\epsilon}^{\mathrm{el}}}\) is automatically declared. For each state variable Y, its time derivative dY is automatically declared.

Implicit integration

If the evolution of the state variables, grouped into a single vector \(Y\) whose components \(Y_{i}\) may be scalars or symmetric tensors, is given by the following system of differential equations: \[ \dot{Y}=G{\left(Y,t\right)} \quad\quad \Leftrightarrow\quad\quad \left\{ \begin{aligned} \dot{Y}_{0} &= g_{Y_{0}}{\left(Y,t\right)}\\ &\ldots \\ \dot{Y}_{i} &= g_{Y_{i}}{\left(Y,t\right)}\\ &\ldots \\ \dot{Y}_{N} &= g_{Y_{N}}{\left(Y,t\right)}\\ \end{aligned} \right. \]

where the dependency with respect to time stands for the evolution of some external state variables and the evolution of strains (for small strain behaviours) which are supposed to evolve linearly during the time step.

The integration of this ordinary differential equation over a time step \(\Delta\,t\) using an implicit algorithm leads to the (generally non-linear) system of equations (see Besson and Desmorat 2004): \[ F{\left(\Delta\,Y\right)}=0\quad \Leftrightarrow\quad \left\{ \begin{aligned} f_{Y_{0}}&=\Delta\,Y_{0}-\Delta\,t\,g_{y_{0}}{\left({\left.Y\right|_{t+\theta\,\Delta\,t}},t\right)}&=0 \\ &\ldots \\ f_{Y_{i}}&=\Delta\,Y_{i}-\Delta\,t\,g_{y_{i}}{\left({\left.Y\right|_{t+\theta\,\Delta\,t}},t\right)}&=0 \\ &\ldots \\ f_{Y_{N}}&=\Delta\,Y_{N}-\Delta\,t\,g_{y_{N}}{\left({\left.Y\right|_{t+\theta\,\Delta\,t}},t\right)}&=0 \\ \end{aligned} \right. \]

where the unknowns are the state variables increments \(\Delta\,Y_{i}\), and we introduced the following notation: \[ {\left.Y\right|_{t+\theta\,\Delta\,t}}=Y+\theta\,\Delta\,Y \]

Algorithms used to solve this system of equations may require the jacobian matrix \(J\) of \(F\) which can be computed by blocks (see J. Besson, Cailletaud, and Chaboche 2001): \[ J={{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,Y}}}= \begin{pmatrix} {{\displaystyle \frac{\displaystyle \partial f_{y_{1}}}{\displaystyle \partial \Delta\,y_{1}}}} & \ldots & \ldots & \ldots & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & {{\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\,y_{j}}}} & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & {{\displaystyle \frac{\displaystyle \partial f_{y_{N}}}{\displaystyle \partial \Delta\,y_{N}}}} \\ \end{pmatrix} \]

Time independent mechanisms

For state variable associated with time-independent mechanisms, the implicit equation shall impose that the system lies on the yield surface when plastic loading occurs.

Available algorithms

Several algorithms are available to solve the previous implicit system:

Consistent tangent operator

For most small strain behaviours, algorithms providing the jacobian matrix \(J\) of the implicit system have a significant advantage: the consistent tangent operator \({{\displaystyle \frac{\displaystyle \partial \Delta\,\sigma}{\displaystyle \partial \Delta\,{\underline{\epsilon}^{\mathrm{to}}}}}}\) can be computed almost automatically with only a small numerical cost.

Example

The Norton creep law can be implemented as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
@Parser Implicit;
@Behaviour ImplicitNorton;
@Includes{
#include"TFEL/Material/Lame.hxx"
}

@MaterialProperty stress young; /* mandatory for castem */
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;    /* mandatory for castem */
nu.setGlossaryName("PoissonRatio");

@LocalVariable real     lambda;
@LocalVariable real     mu;

// store for the Von Mises stress 
// for the tangent operator
@LocalVariable real seq;
// store the derivative of the creep function
// for the tangent operator
@LocalVariable real df_dseq;
// store the normal tensor
// for the tangent operator
@LocalVariable Stensor n;

@StateVariable real    p;
@PhysicalBounds p in [0:*[;

/* Initialize Lame coefficients */
@InitLocalVars{
  using namespace tfel::material::lame;
  lambda = computeLambda(young,nu);
  mu = computeMu(young,nu);
} // end of @InitLocalVars

@IsTangentOperatorSymmetric true;
@TangentOperator{
  using namespace tfel::material::lame;
  if((smt==ELASTIC)||(smt==SECANTOPERATOR)||
     (smt==TANGENTOPERATOR)){
    computeAlteredElasticStiffness<hypothesis,Type>::exe(Dt,lambda,mu);
  } else if (smt==CONSISTENTTANGENTOPERATOR){
    StiffnessTensor Hooke;
    Stensor4 Je;
    computeElasticStiffness<N,Type>::exe(Hooke,lambda,mu);
    getPartialJacobianInvert(Je);
    Dt = Hooke*Je;
  } else {
    return false;
  }
}

@ComputeStress{
  sig = lambda*trace(eel)*Stensor::Id()+2*mu*eel;
} // end of @ComputeStresss

@Integrator{
  const real A = 8.e-67;
  const real E = 8.2;
  seq = sigmaeq(sig);
  const real tmp = A*pow(seq,E-1.);
  df_dseq = E*tmp;
  real inv_seq(0);
  n = Stensor(0.);
  if(seq > 1.e-8*young){
    inv_seq = 1/seq;
    n       = 1.5*deviator(sig)*inv_seq;
  }
  feel += dp*n-deto;
  fp   -= tmp*seq*dt;
  // jacobian
  dfeel_ddeel += 2.*mu*theta*dp*inv_seq*(Stensor4::M()-(n^n));
  dfeel_ddp    = n;
  dfp_ddeel    = -2*mu*theta*df_dseq*dt*n;
} // end of @Integrator

References

Besson, J., and D. Desmorat. 2004. “Numerical Implementation of Constitutive Models.” In Local Approach to Fracture, edited by J. Besson. École des Mines de Paris - les presses.

Besson, Jacques, Georges Cailletaud, and Jean-Louis Chaboche. 2001. Mécanique Non Linéaire Des Matériaux. Paris: Hermès.

Besson, Jacques, Rodolphe Le Riche, Ronald Foerch, and Georges Cailletaud. 1998. “Object-Oriented Programming Applied to the Finite Element Method Part II. Application to Material Behaviors.” Revue Européenne Des Éléments 7 (5): 567–88.

Chaboche, Jean-Louis, Jean Lemaître, Ahmed Benallal, and Rodrigue Desmorat. 2009. Mécanique Des Matériaux Solides. Paris: Dunod.

Chen, Hern-Shann, and Mark A. Stadtherr. 1981. “A Modification of Powell’s Dogleg Method for Solving Systems of Nonlinear Equations.” Computers & Chemical Engineering 5 (3): 143–50. doi:10.1016/0098-1354(81)85003-X.

Doghri, Issam. 2000. Mechanics of Deformable Solids: Linear, Nonlinear, Analytical, and Computational Aspects. Berlin; New York: Springer.

EDF. 2013. “Algorithme Non Linéaire Quasi-Statique STAT_NON_LINE.” Référence du Code Aster R5.03.01 révision : 10290. EDF-R&D/AMA. http://www.code-aster.org.

Foerch, R., J. Besson, G. Cailletaud, and P. Pilvin. 1997. “Polymorphic Constitutive Equations in Finite Element Codes.” Computer Methods in Applied Mechanics and Engineering 141 (3–4): 355–72. doi:10.1016/S0045-7825(96)01111-5.

Miehe, C., N. Apel, and M. Lambrecht. 2002. “Anisotropic Additive Plasticity in the Logarithmic Strain Space: Modular Kinematic Formulation and Implementation Based on Incremental Minimization Principles for Standard Materials.” Computer Methods in Applied Mechanics and Engineering 191 (47–48): 5383–5425. doi:10.1016/S0045-7825(02)00438-3.

Northwest Numerics and Modeling, Inc. 2014. “ZebFront.” http://www.nwnumerics.com/Z-mat/ZebFront.

Simo, J. C., and R. L. Taylor. 1985. “Consistent Tangent Operators for Rate-Independent Elastoplasticity.” Computer Methods in Applied Mechanics and Engineering 48 (1): 101–18. doi:10.1016/0045-7825(85)90070-2.

Simo, J.C., and C. Miehe. 1992. “Associative Coupled Thermoplasticity at Finite Strains: Formulation, Numerical Analysis and Implementation.” Computer Methods in Applied Mechanics and Engineering 98 (1): 41–104. doi:10.1016/0045-7825(92)90170-O.

Simo, Juan C, and Thomas J. R Hughes. 1998. Computational Inelasticity. New York: Springer.

Tvergaard, V. 1990. “Effect of Fibre Debonding in a Whisker Reinforced Metal.” Mater. Sci. Eng. A125: pp. 203–13.

Zienkiewicz, O.C. 1977. The Finite Element Method. McGraw-Hill.