Overview

The versatility of a 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 developed by the Centre des Matériaux de Mines ParisTech as part of the Zset software [see 1, 2, 3]. 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 specify the role of the mechanical behaviour in standard displacement-based finite element solvers [see 4, 5, 6]. 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 derivative of the stress tensor with respect to 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 is 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}\). If the solver requires it, the mechanical behaviour may also provide an estimation of the tangent consistent operator \({\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\epsilon}^{\mathrm{to}}}}\) (see [7]) 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 on external state variables evolutions with time, namely the temperature, the irradiation damage, and so on.

Isotropic \(J_{2}\) plastic/viscoplastic 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 [8]).

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

@DSL 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 [5];[9]):

\[ 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, \(H\) 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 specializations. A small overview of the programming techniques used can be found in the description of the TFEL/Math library.
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 typically used by solvers based on the MGIS libray appropriately calling 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 over this procedure.

Finite strain strategies

If not handled directly by the calling code, appropriate pre- and post-processing 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 [10]);
2. lagrangian logarithmic strains as proposed by Miehe, Apel and Lambrecht (see [11];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 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:

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 close to those obtained by the classical formulation based on an \(F_{e}\,F_{p}\) decomposition proposed by Simo and Miehe (see [12]) (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 to write its own integration algorithm. This is very useful for explicit behaviours such as the classical cohesive zone model [see 13].
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:

@DSL DefaultCZM;     // 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:

@DSL    RungeKutta;              // domain specific language
@Behaviour OrthotropicCreep;     // name of the behaviour
@OrthotropicBehaviour;           // treating an orthotropic behaviour
@RequireStiffnessTensor;         // requires the stiffness tensor to be computed
@StateVariable Stensor evp;      // viscoplastic strain
@StateVariable strain p;         // Equivalent viscoplastic strain
@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(real(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

This section provides a very coarse overview of the Implicit domain specific language. The reader can refer to this page for a comprehensive description.

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 14]: \[ 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 5]: \[ 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:

• NewtonRaphson is the standard Newton-Raphson algorithm: \[ \Delta\,Y^{\left(n+1\right)}=\Delta\,Y^{\left(n\right)}-J^{-1}\,.\,F\left(\Delta\,Y^{\left(n\right)}\right) \] The user must explicitly compute the jacobian matrix, which constitutes the main difficulty of this method. For debugging purposes, MFront may generate the comparison of each block of the jacobian matrix with a numerical approximation.
• NewtonRaphson_NumericalJacobian is a variation of the standard Newton-Raphson algorithm using a jacobian matrix computed by a second order finite difference. Writing behaviour implementations using this algorithm is as easy as using the domain specific languages based on RungeKutta algorithm. It can be considered as a first step toward an implicit implementation with an analytical jacobian matrix.
• Broyden algorithms which do not require the computation of the jacobian matrix: these algorithms update an approximation of the jacobian matrix (first Broyden algorithm) or its inverse (second Broyden algorithm) at each iteration. The first Broyden algorithm can sometimes be interesting as one may compute analytically some part of the jacobian matrix and let the algorithm compute the other parts. If the computation of those other parts takes a significant amount of CPU time, this algorithm can in same cases outperform the Newton-Raphson algorithm.
• PowellDogLeg_XX algorithm, where XX is one of the previous algorithm. Those trust-region algorithms implement the classical Powell dogleg method [see 15] to improve the robustness of the resolution.

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:

@DSL Implicit;
@Behaviour ImplicitNorton;
@Brick StandardElasticity;

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

@LocalVariable stress lambda,mu;

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

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

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

References