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.
We now precise the role of the mechanical behaviours in standard displacementbased 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:
MFront
;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.
\(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 

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:
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.
If not handled directly by the calling code, appropriate pre and postprocessings allowing the use of small strain behaviours in finite strain computations can be generated. Two lagrangian finite strain strategies are currently available:
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 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
).
Apart from the domain specific languages dealing with isotropic \(J_{2}\) plastic and viscoplastic behaviours presented in the previous paragraph, MFront
also provides several generalpurpose domain specific languages:
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).RungeKutta
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.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.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 

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}\)
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 

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.
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 nonlinear) 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} \]
For state variable associated with timeindependent mechanisms, the implicit equation shall impose that the system lies on the yield surface when plastic loading occurs.
Several algorithms are available to solve the previous implicit system:
NewtonRaphson
is the standard NewtonRaphson 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 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. It can be considered as a first step toward an implicit implementation with an analytical jacobian matrix.Broyden
algorithms which do not require to 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 analitically 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 outperfom the {NewtonRaphson} algorithm.PowellDogLeg_XX}
algorithm, where XX
is one of the previous algorithm. Those trustregion algorithms implements the classical Powell dogleg method (see Chen and Stadtherr 1981) to improve the robustness of the resolution.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.
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 

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 JeanLouis Chaboche. 2001. Mécanique Non Linéaire Des Matériaux. Paris: Hermès.
Besson, Jacques, Rodolphe Le Riche, Ronald Foerch, and Georges Cailletaud. 1998. “ObjectOriented Programming Applied to the Finite Element Method Part II. Application to Material Behaviors.” Revue Européenne Des Éléments 7 (5): 567–88.
Chaboche, JeanLouis, Jean Lemaître, Ahmed Benallal, and Rodrigue Desmorat. 2009. Mécanique Des Matériaux Solides. Paris: Dunod.
Chen, HernShann, 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/00981354(81)85003X.
Doghri, Issam. 2000. Mechanics of Deformable Solids: Linear, Nonlinear, Analytical, and Computational Aspects. Berlin; New York: Springer.
EDF. 2013. “Algorithme Non Linéaire QuasiStatique STAT_NON_LINE.” Référence du Code Aster R5.03.01 révision : 10290. EDFR&D/AMA. http://www.codeaster.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/S00457825(96)011115.
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/S00457825(02)004383.
Northwest Numerics and Modeling, Inc. 2014. “ZebFront.” http://www.nwnumerics.com/Zmat/ZebFront.
Simo, J. C., and R. L. Taylor. 1985. “Consistent Tangent Operators for RateIndependent Elastoplasticity.” Computer Methods in Applied Mechanics and Engineering 48 (1): 101–18. doi:10.1016/00457825(85)900702.
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/00457825(92)90170O.
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. McGrawHill.