This tutorial explains how to implement Cailletaud’s \(\beta\)-rule [1] in the context of elasto-viscoplastic composites using an arbitrary local flow function. In this tutorial, a local behaviour of Meric-Cailletaud type [2] is proposed.

Here, the implementation of this homogenization scheme requires the local tangent operators, obtained via the integration of a local Behaviour.

We first recall Cailletaud’s \(\beta\)-rule and present the non-linear system derived from it. We then implement it in MFront and explain how to use the keyword @BehaviourVariable for the integration of the local behaviours.

Theoretical background

Cailletaud’s beta rule

We consider a composite made of \(N\) phases indexed by \(i\). Cailletaud’s \(\beta\)-rule was proposed in the context of elasto-viscoplastic composites for which strain partition is assumed: \[\begin{aligned} \underline{\epsilon}^{\mathrm{to}}_i=\underline{\epsilon}^{\mathrm{el}}_i+\underline{\epsilon}^{\mathrm{vp}}_i \\ \end{aligned}\] where the strain \(\underline{\epsilon}^{\mathrm{el}}_i\) is given by Hooke’s law: \[\begin{aligned} \underline{\epsilon}^{\mathrm{el}}_i=\underline{\mathbf{S}}_i\mathbin{\mathord{:}}\underline{\sigma}_i \\ \end{aligned}\] where \(\underline{\mathbf{S}}_i\) is a fourth-rank compliance tensor and \(\underline{\sigma}_i\) is the stress field; and where the viscoplastic strain \(\underline{\epsilon}^{\mathrm{vp}}_i\) is given by the following flow rule: \[\begin{aligned} \underline{\dot{\epsilon}}^{\mathrm{vp}}_i=f_i(\underline{\sigma}_i) \\ \end{aligned}\]

where \(f_i\) is a non-linear function.

Assuming that \(\phi_i\) is volume fraction of phase \(i\), the macroscopic strain field is defined as \[\begin{aligned} \underline{E}= \sum_\alpha \phi_i \underline{\epsilon}^{\mathrm{to}}_i \\ \end{aligned}\]

and is given as a parameter.

In this context, the proposition of Cailletaud is to introduce an internal variable \(\underline{\beta}_i\) on each phase, whose evolution is given by \[\begin{aligned} \underline{\dot{\beta}}_i= \underline{\dot{\epsilon}}^{\mathrm{vp}}_i-D\,||\underline{\dot{\epsilon}}^{\mathrm{vp}}_i||\,\underline{\dot{\beta}}_i \\ \end{aligned}\] where \(D\in\mathbb R\) is a parameter and \[\begin{aligned} ||\underline{\dot{\epsilon}}^{\mathrm{vp}}_i||=\sqrt{\dfrac 23 \underline{\dot{\epsilon}}^{\mathrm{vp}}_i\mathbin{\mathord{:}}\underline{\dot{\epsilon}}^{\mathrm{vp}}_i}. \\ \end{aligned}\] The localisation relation is also proposed: \[\begin{aligned} \underline{\sigma}_i= \underline{\Sigma}+ c\,(\underline{B}-\underline{\beta}_i) \\ \end{aligned}\] where \(c\) is a real parameter, \(\underline{B}\) is a macroscopic variable: \[\begin{aligned} \underline{B}= \sum_i \phi_i\underline{\beta}_i\\ \end{aligned}\] and \(\underline{\Sigma}\) is the macroscopic stress, which verifies (by averaging the localisation relation) \[\begin{aligned} \underline{\Sigma}= \sum_i \phi_i\underline{\sigma}_i\\ \end{aligned}\]

Resolution of the non-linear system

We discretize the time interval and look for the increment of the variables. We choose \(\Delta\underline{\epsilon}^{\mathrm{to}}_i\), \(\Delta{\underline{\beta}}_i\) and \(\Delta{\underline{\Sigma}}\) as unknowns of our non-linear system:

\(\left\{ \begin{aligned} &\sum_{i=1}^{N}\phi_i\,\Delta\underline{\epsilon}^{\mathrm{to}}_i - \Delta\underline{E} = 0\\ &\Delta{\underline{\sigma}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)-\Delta{\underline{\Sigma}}- c\,(\Delta\underline{B}-\Delta\underline{\beta}_i) = 0\qquad\forall i\\ &\Delta\underline{\beta}_i - \Delta\underline{\epsilon}^{\mathrm{vp}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)+ D\,||\Delta\underline{\epsilon}^{\mathrm{vp}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)||\,(\underline{\beta}_i+\theta\Delta\underline{\beta}_i)=0\qquad\forall i\\ \end{aligned}\right.\)

where the following auxiliary variable is introduced:

\(\begin{aligned} &\Delta\underline{B} = \sum_{i=1}^{N}\phi_i\,\Delta\underline{\beta}_i\\ \end{aligned}\) and where \(\Delta{\underline{\sigma}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)\) and \(\Delta\underline{\epsilon}^{\mathrm{vp}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)\) are obtained with the integration of the local behaviours. Note that the local elastic strain can also be computed with this local integration.

To solve this non-linear problem in MFront, we use the ImplicitII DSL (because we do not consider \(\underline{\epsilon}^{\mathrm{el}}\) as an integration variable). We only have to precise the residues and the Jacobian matrix. The residues are:

\(\left\{ \begin{aligned} &f_{\underline{E}}=\sum_{i=1}^{N}\phi_i\,\Delta\underline{\epsilon}^{\mathrm{to}}_i - \Delta\underline{E}\\ &f_{\underline{\epsilon}^{\mathrm{to}}_i} =\frac{1}{\sigma_0}\left[\Delta{\underline{\sigma}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)-\Delta{\underline{\Sigma}}- c\,(\Delta\underline{B}-\Delta\underline{\beta}_i)\right] \qquad \forall i\\ &f_{\underline{\beta}_i} =\Delta\underline{\beta}_i - \Delta\underline{\epsilon}^{\mathrm{vp}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)+D\,||\Delta\underline{\epsilon}^{\mathrm{vp}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)||\,(\underline{\beta}_i+\theta\Delta\underline{\beta}_i) \qquad \forall i\\ \end{aligned}\right.\)

where \(\sigma_{0}\) is a parameter that give the same dimension and magnitude to the residues. The jacobian matrix is given by

\(\left\{ \begin{aligned} &{\displaystyle \frac{\displaystyle \partial f_{\underline{E}}}{\displaystyle \partial \Delta{\underline{\Sigma}}}} = \underline{\mathbf{0}}\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{E}}}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}} = \phi_j\underline{\mathbf{I}}\qquad \forall j\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{E}}}{\displaystyle \partial \Delta{\underline{\beta}}_j}} = \underline{\mathbf{0}}\qquad \forall j\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{to}}_i}}{\displaystyle \partial \Delta{\underline{\Sigma}}}} = -\frac{1}{\sigma_0}\underline{\mathbf{I}}\qquad \forall i\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{to}}_i}}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}} = \frac{\delta_{ij}}{\sigma_0}\,{\displaystyle \frac{\displaystyle \partial \underline{\sigma}_j}{\displaystyle \partial \underline{\epsilon}^{\mathrm{to}}_j}}\left(\underline{\epsilon}^{\mathrm{to}}_j,\Delta\underline{\epsilon}^{\mathrm{to}}_j\right)\qquad\forall i,j\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{\epsilon}^{\mathrm{to}}_i}}{\displaystyle \partial \Delta{\underline{\beta}}_j}} = \frac{c}{\sigma_0}\,(\delta_{ij}-\phi_j)\underline{\mathbf{I}}\qquad \forall i,j\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{\beta}_i}}{\displaystyle \partial \Delta{\underline{\Sigma}}}} = \underline{\mathbf{0}}\qquad \forall i\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{\beta}_i}}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}} = -\delta_{ij}\,{\displaystyle \frac{\displaystyle \partial \underline{\sigma}_j}{\displaystyle \partial \underline{\epsilon}^{\mathrm{to}}_j}}\left(\underline{\epsilon}^{\mathrm{to}}_j,\Delta\underline{\epsilon}^{\mathrm{to}}_j\right)+D\,(\underline{\beta}_i+\theta\Delta\underline{\beta}_i)\otimes{\displaystyle \frac{\displaystyle \partial ||\Delta\underline{\epsilon}^{\mathrm{vp}}_i||}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}}\qquad \forall i,j\\ &{\displaystyle \frac{\displaystyle \partial f_{\underline{\beta}_i}}{\displaystyle \partial \Delta{\underline{\beta}}_j}} = \left(1+\theta\,D\,||\Delta\underline{\epsilon}^{\mathrm{vp}}_i||\right)\,\delta_{ij}\,\underline{\mathbf{I}} \qquad \forall i,j\\ \end{aligned}\right.\)

where \(||\Delta\underline{\epsilon}^{\mathrm{vp}}_i||\) is used for \(||\Delta\underline{\epsilon}^{\mathrm{vp}}_i\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)||\). Moreover,

\(\left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial ||\Delta\underline{\epsilon}^{\mathrm{vp}}_i||}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}}&=\frac 23\delta_{ij}{\displaystyle \frac{\displaystyle \Delta\underline{\epsilon}^{\mathrm{vp}}_i}{\displaystyle ||\Delta\underline{\epsilon}^{\mathrm{vp}}_i||}}\mathbin{\mathord{:}}{\displaystyle \frac{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{vp}}_i}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}}\left(\underline{\epsilon}^{\mathrm{to}}_i,\Delta\underline{\epsilon}^{\mathrm{to}}_i\right)\\ &=\frac 23\delta_{ij}{\displaystyle \frac{\displaystyle \Delta\underline{\epsilon}^{\mathrm{vp}}_i}{\displaystyle ||\Delta\underline{\epsilon}^{\mathrm{vp}}_i||}}\mathbin{\mathord{:}}\left({\displaystyle \frac{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_i}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}}-{\displaystyle \frac{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{el}}_i}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}}\right)\\ &=\frac 23\delta_{ij}{\displaystyle \frac{\displaystyle \Delta\underline{\epsilon}^{\mathrm{vp}}_i}{\displaystyle ||\Delta\underline{\epsilon}^{\mathrm{vp}}_i||}}\mathbin{\mathord{:}}\left(\delta_{ij}\underline{\mathbf{I}}-\delta_{ij}\underline{\mathbf{S}}_j\mathbin{\mathord{:}}{\displaystyle \frac{\displaystyle \partial \Delta{\underline{\sigma}}_j}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}}\right)\\ \end{aligned}\right.\)

We can see the term \({\displaystyle \frac{\displaystyle \partial \Delta{\underline{\sigma}}_j}{\displaystyle \partial \Delta\underline{\epsilon}^{\mathrm{to}}_j}}\), that will be provided by the local integration.

Consistent tangent operator

We can also compute the tangent operator with a method similar as the method explained in the Implicit DSL documentation. We define \(\Delta Y\) such as \[\begin{aligned} \Delta G_i = (\Delta\underline{\epsilon}^{\mathrm{to}}_i, \Delta{\underline{\beta}}_i)\\ \Delta Y = (\Delta{\underline{\Sigma}},\Delta G_1,...,\Delta G_N)\\ \end{aligned}\] and \(H_i=(f_{\underline{\epsilon}^{\mathrm{to}}_i},f_{\underline{\beta}_i})\) and \(F=(f_{\underline{E}},H_1,...,H_N)\) so that \[\begin{aligned} F(\Delta\underline{E},\Delta Y(\Delta\underline{E}))&=0\\ \end{aligned}\]

Hence, the following equality is still verified:

\[\begin{aligned} {\displaystyle \frac{\displaystyle \partial \Delta Y}{\displaystyle \partial \Delta\underline{E}}}&=-\mathbf{J}^{-1}\cdot{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta \underline{E}}}\\ \end{aligned}\]

where \(\mathbf{J}={\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta Y}}\) is the Jacobian of the system.

Given the form of \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta \underline{E}}}\), \({\displaystyle \frac{\displaystyle \partial \Delta{\underline{\Sigma}}}{\displaystyle \partial \Delta\underline{E}}}\) can be retrieved from the \(6\times 6\) upper-left submatrix of \(\mathbf{J}^{-1}\).

Implementation in MFront

The implementation requires the integration of the local behaviours. These are carried out via the @BehaviourVariable keyword. The implementation of the local behaviour is explained here.

All files MericCailletaudSingleCrystalViscoPlasticity.mfront, BetaRule.mfront and BetaRule.mtest can be downloaded here.

Implementation of Cailletaud’s beta rule

For the example, we assume that the composites is made of only 2 phases. We hence define the integration variables:

@StateVariable Stensor Sig;
Sig.setEntryName("MacroscopicStress");
@StateVariable Stensor eto1;
eto1.setEntryName("FirstPhaseTotalStrain");
@StateVariable Stensor eto2;
eto2.setEntryName("SecondPhaseTotalStrain");
@StateVariable Stensor beta1;
beta1.setEntryName("FirstPhaseBetaStrain");
@StateVariable Stensor beta2;
beta2.setEntryName("SecondPhaseBetaStrain");

and we define the local behaviours with @BehaviourVariable:

@BehaviourVariable b1 {
file: "MericCailletaudSingleCrystalViscoPlasticity.mfront",
variables_suffix: "1",
store_gradients: false,
external_names_prefix: "FirstPhase",
shared_external_state_variables: {".+"}
};

(the same for phase 2).

We note that the gradients are not stored because the StateVariable eto1 (resp. eto2) are defined.

The integration of the local behaviour is performed as follows:

  initialize(b1);
  b1.eto=eto1;
  b1.deto=deto1;
  constexpr auto b1_smflag = TangentOperatorTraits<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR>::STANDARDTANGENTOPERATOR;
  const auto r1 = b1.integrate(b1_smflag,CONSISTENTTANGENTOPERATOR);
  StiffnessTensor Dt1 = b1.getTangentOperator();

and the following variables are stored:

  auto dsig1=b1.sig-sig1;
  auto deel1=b1.eel-eel1;
  auto devp1=deto1-deel1;
  auto ndevp_1=sqrt(real(2)/3*devp1|devp1);
  auto beta_mts_1=beta1+theta*dbeta1;

The same goes for phase 2. Here are the residues:

  //residues
  auto dB=(1-f)*dbeta1+f*dbeta2;
  fSig = -deto+ (1-f)*deto1+f*deto2;
  feto1 = (dsig1-dSig-c*(dB-dbeta1))/sig_0;
  feto2 = (dsig2-dSig-c*(dB-dbeta2))/sig_0;
  fbeta1 = dbeta1-devp1+DD*ndevp_1*beta_mts_1;
  fbeta2 = dbeta2-devp2+DD*ndevp_2*beta_mts_2;

and the terms of the analytical jacobian:

  //jacobian matrix
  auto Null = Stensor4{real{}};
  dfSig_ddSig=Null;
  dfSig_ddeto1 = (1-f)*Id;
  dfSig_ddeto2 = f*Id;
  
  dfeto1_ddSig = -Id/sig_0;
  dfeto2_ddSig = -Id/sig_0;
  
  dfeto1_ddeto1 = Dt1/sig_0;
  dfeto2_ddeto2 = Dt2/sig_0;
  
  dfeto1_ddbeta1 = c/sig_0*f*Id;
  dfeto1_ddbeta2 = -c/sig_0*f*Id;
  dfeto2_ddbeta1 = c/sig_0*(f-1)*Id;
  dfeto2_ddbeta2 = c/sig_0*(1-f)*Id;
  
  dfbeta1_ddeto1 = -(Id-S1*Dt1) + 2*DD/3*Stensor4(beta_mts_1^devp1)*Stensor4(Id-S1*Dt1)/max(ndevp_1,eeps);
  dfbeta2_ddeto2 = -(Id-S2*Dt2) + 2*DD/3*Stensor4(beta_mts_2^devp2)*Stensor4(Id-S2*Dt2)/max(ndevp_2,eeps);
  
  dfbeta1_ddbeta1 =Id+ theta*DD*ndevp_1*Id;
  dfbeta2_ddbeta2 =Id+ theta*DD*ndevp_2*Id;

Note that S1 and S2 are the compliance tensors of each phase and were defined by the elastic coefficients of each phase. The computation of the stress and of the tangent operator is straightforward:

@ComputeFinalStress{
  sig += dSig;
}

@TangentOperator{
  Stensor4 iJs;
  getPartialJacobianInvert(iJs);
  Dt=iJs;
}

Note that getPartialJacobianInvert(iJs) permits to obtain the 6x6 left-upper part of the inverse of the Jacobian. The left-upper part is associated with the first integration variable declared (here Sig).

Results

We used MTest to carry out simulations of an uniaxial tensile test on a two-grain polycrystal. The first phase is stiffer than the second, and we also choose for the example a different Norton exponent between the phases. The following values are given by the mtest file:

@ModellingHypothesis 'Tridimensional';
@Behaviour<Generic> 'src/libBehaviour.so' 'BetaRuleBehaviour';
@MaterialProperty<constant> 'FirstPhaseFraction' 0.5;
@MaterialProperty<constant> 'E_1'               208000;
@MaterialProperty<constant> 'nu_1'               0.3;
@MaterialProperty<constant> 'mu_1'               80000;
@MaterialProperty<constant> 'E_2'               104000;
@MaterialProperty<constant> 'nu_2'               0.3;
@MaterialProperty<constant> 'mu_2'               40000;
@MaterialProperty<constant> 'tau01'               36.;
@MaterialProperty<constant> 'tau02'               25.;
@MaterialProperty<constant> 'n1'               10.;
@MaterialProperty<constant> 'n2'               5.;
@ExternalStateVariable 'Temperature' 293.15;
@ImposedStrain 'EXX' {0 : 0, 5 :1e-2};
@Times {0, 3 in 400};

Note that the elastic behaviours are isotropic for the example. (the file is available in the archive). We plotted here the macroscopic and local behaviours (\(\underline{\sigma}_i\) as a function of \(\underline{\epsilon}_i\)).

Macroscopic stress as a function of local strain, uniaxial tensile test, Cailletaud’s beta rule
1.
Cailletaud, Georges and Pilvin, Philippe. Utilisation de modèles polycristallins pour le calcul par éléments finis. Revue européenne des éléments finis. 1994. Vol. 3, p. 515–541.
2.
Méric, L. and Cailletaud, Georges. Single crystal modelling for structural calculations. Journal of Engineering Material and Technology. January 1991. Vol. 113, p. 171–182.