The aim of this report is to describe briefly the implementation of a simple self-consistent polycrystalline model based on the use of the Berveiller-Zaoui concentration rule [1]. This model provides an explicit relationship between the local stress, inside a crystallographic phase, and the macroscopic applied stress. This model has been developed for microscopically isotropic elastic behaviour and macroscopically isotropic plastic behaviour.

After a quick recall of the equations involved, an explicit implementation of a polycrystal behaviours based on this homogeneization scheme is presented in Sections 2. This choice is justified by the huge number of state variables involved which is a priori incompatible with an implicit resolution. The main trouble with explicit scheme is the lack of consistent tangent operator which degrades considerably the performance of most mechanical solvers at the structural scale, the Cast3M solver being a noticeable exception du to its specific resolution scheme [3, 4].

The implementation described here can be downloaded here.

Constitutive equations

Notations

• \(\underline{\Sigma}\): the macroscopic stress.
• \(\Sigma_{\mathrm{eq}}\): the macroscropic von Mises stress.
• \(\underline{E}^{\mathrm{to}}\): the macroscopic strain.
• \(\underline{E}^{\mathrm{vp}}\): the macroscopic viscoplastic strain.
• \(N_{g}\): number of grains/phases.
• \({\underline{\sigma}}^{(k)}\): the microscopic stress of the k\(\mbox{}^{\text{th}}\) phase.
• \({\phi}^{(k)}\): the volume fraction of the k\(\mbox{}^{\text{th}}\) phase: \[ {\phi}^{(k)}=\dfrac{{V}^{(k)}}{\sum_{l=1}^{N_{g}}{V}^{(l)}} \]
• \(M\): the macroscropic shear modulus.
• \(\nu\): the macroscropic Poisson ratio.
• \(P\): the equivalent macroscropic viscoplastic strain.
• \({g}^{(k,i)}\): plastic slip along the i\(\mbox{}^{\text{th}}\) slip system of the k\(\mbox{}^{\text{th}}\) phase.
• \({\mu}^{(k,i)}\): orientation tensor of the i\(\mbox{}^{\text{th}}\) slip system of the k\(\mbox{}^{\text{th}}\) phase.

The Berveiller-Zaoui homogeneisation scheme

In the framework of macroscopically isotropic elasto-plastic behaviors, Berveiller and Zaoui have derived a self-consistent polycrystalline model using a secant approximation for the linearization of the non-linear behavior of the material [1]. The grains of the polycrystal have equiaxed shape and the texture of the material is isotropic. The elasticity is assumed homogeneous and istropic In the case of monotonous radial loadings Berveiller and Zaoui have shown that the local stress \({\underline{\sigma}}^{(k)}\) applied in a crystallographic phase \((i)\) (all the grains with nearly the same orientation) can be expressed explicitly as a function of the local plastic strain \({\underline{\varepsilon}^{\mathrm{vp}}}^{(k)}\), the macroscopic stress \(\underline{\Sigma}\) and the macroscopic plastic strain \(\underline{E}^{\mathrm{vp}}\).

The macrosopic viscoplastic strain is linked to the microscopic local plastic strain \({\underline{\varepsilon}^{\mathrm{vp}}}^{(k)}\) by a simple mixing rule:

\[ \underline{E}^{\mathrm{vp}}=\sum_{k=1}^{N_{g}} {\phi}^{(k)}\,{\underline{\varepsilon}^{\mathrm{vp}}}^{(k)} \qquad{(1)}\]

where \({\phi}^{(k)}\) is the volume fraction of the k\(\mbox{}^{\text{th}}\) phase.

The local stress \({\underline{\sigma}}^{(k)}\) is given by the following relationship¹:

\[ {\underline{\sigma}}^{(k)}=\underline{\Sigma}+2\,M\alpha{\left(\Sigma_{\mathrm{eq}},P\right)}\,{\left(1-\beta{\left(\nu\right)}\right)}\,{\left(\underline{E}^{\mathrm{vp}}-{\underline{\varepsilon}^{\mathrm{vp}}}^{(k)}\right)}\quad\text{with}\quad\beta{\left(\nu\right)}=\dfrac{2}{15}\,\dfrac{4-5\,\nu}{1-\nu} \qquad{(2)}\]

where:

• \(\Sigma_{\mathrm{eq}}\) is the macroscropic von Mises stress.
• \(P\) is the equivalent macroscropic viscoplastic strain defined in an incremental way as follows: \[ \dot{P}=\sqrt{\dfrac{2}{3}\,\underline{\dot{E}}^{\mathrm{vp}}\,\colon\,\underline{\dot{E}}^{\mathrm{vp}}} \]

The function \(\alpha{\left(\Sigma_{\mathrm{eq}},P\right)}\) is called the accommodation function. It represents the plastic accommodation between the Homogeneous Equivalent Medium and the considered crystallographic phase (r) in inclusion (spherical inclusion). The expression of the accommodation function is as follows:

\[ \alpha{\left(\Sigma_{\mathrm{eq}},P\right)}=\dfrac{2\,\Sigma_{\mathrm{eq}}}{2\,\Sigma_{\mathrm{eq}}+3\,M\,P} \]

Using this explicit concentration rule, it is particularly simple to compute the overall response of a polycrystalline material, provided the intra-granular constitutive behaviour is known.

Intra-granular constitutive behavior

In this report, we assume that the local constitutive equations have the following form:

\[ {\underline{\dot{\varepsilon}}^{\mathrm{vp}}}^{(k)}={v}^{(k)}{\left({\underline{\sigma}}^{(k)}\right)} \qquad{(3)}\]

In this report, the microscopic plastic strain may results from the slips along the slip systems of the crystal following:

\[ {\underline{\dot{\varepsilon}}^{\mathrm{vp}}}^{(k)}=\sum_{i=1}^{N_{s}}{\dot{g}}^{(k,i)}\,{\mu}^{(k,i)}\quad\text{with}\quad{\dot{g}}^{(k,i)}=\left\langle\dfrac{{\sigma}^{(k)}\,\colon\,{\mu}^{(k,i)}-R_{0}}{K}\right\rangle^{n} \]

This constitutive equation is very simple, but extensions of the proposed integration scheme to more complex ones , such as i.e. the Méric-Cailletaud behaviour [5], is straightforward using the proposed explicit scheme.

For the example, we will follow the work of Onimus et al. with sightly modified material coefficients [2].

Explicit integration scheme

The first explicit implementation of a polycrystal based on the Berveiller and Zaoui homogeneization scheme and the Méric-Cailletaud single crystal behaviour at the grain scale [5] is due to Jean-Michel Proix (EDF) and is available in the sources of MFront². Other implementations based on more complex local behaviours inspired by results from the Dislocation Dynamics are also available [6].

This section presents the successive steps to implement the polycrystalline model in MFront using an explicit scheme.

MFront’ Implementation

Choice of the domain specific language

Domain specific languages (DSL) are of the main concept of MFront. It specifies a kind of material knowledge and a kind of numerical scheme . It also introduces a specific syntax for the implementation to be both efficient and explicit.

For a small strain behaviour integrated using an explicit scheme, the appropriate domain specific language is called RungeKutta. By convention, the declaration of the domain specific language is placed at the beginning of the MFront file:

@DSL RungeKutta;

Name of the behaviour

The @Behaviour keyword allows defining the name of the behaviour, as follows:

@Behaviour PolycrystalHCP_BZ;

Metadata

@Author F. Onimus, C. Gicquel, T. Helfer;
@Date 05/2020;
@Description{
  Implementation of a polycrystal of behaviour  for
  zirconium alloys using the Berveiller-Zaoui
  homegeneization scheme. The presented behaviour
  is inspired from Onimus 2009 which modified
  material coefficients.

  "Onimus, Fabien and Béchade, Jean-Luc."
  "A polycrystalline modeling of the mechanical behavior of"
  "neutron irradiated zirconium alloys."
  "Journal of Nuclear Materials. 15 February 2009."
  "Vol. 384, no. 2, p. 163-174."
}

Symmetry of the material

The material considered will be isotropic on the macroscopic scale. However, due to restrictions to the keywords defining slips systems (see Section 2.1.7), the behaviour must be declared orthotropic.

@OrthotropicBehaviour;

Supported modelling hypothesis

The class defining the slips systems declares \(3D\) orientation tensors (see Section 2.1.7 and 2.1.8). To avoid conversions, we restrict the behaviour to be only valid in 3D:

@ModellingHypothesis Tridimensional;

Note that such a restriction, when applicable, notably reduces the compilation times.

Choice of the algorithm

Among the different schemes available the choice is made to use the Runge-Kutta 5/4 algorithm (which is the default):

@Algorithm rk54;

Declaration of the slip systems

The considered material as an hexagonal close packed crystallographic structure and four families of slip systems:

• the first slip system corresponds to // systèmes prismatiques A
• the second slip system corresponds to // systèmes prismatiques C+1
• the third slip system corresponds to // systèmes basaux A
• the fourth slip system corresponds to // systèmes pyramidaux A

Those slip systems are declared in a straightforward manner using the functionalities provided by the NUMODIS code:

@CrystalStructure HCP;
@SlidingSystems{<1, 1, -2, 0>{1, -1, 0, 0},
                <-2, 1, 1, 3>{1, -1, 0, 1},
                <-2, 1, 1, 0>{0, 0, 0, 1},
                <1, 1, -2, 0>{1, -1, 0, 1}};

Those keywords are fully described on this page: https://thelfer.github.io/tfel/web/singlecrystal.html.

They generate a class called PolycrystalHCP_BZSlipSystems which contains the orientation tensors of the all the slip systems.

Handling textures

Altough, the PolycrystalHCP_BZSlipSystems class contains all the information for a single crystal, a polycrystaline behaviour must handle the texture of the material, i.e. the orientation and the volume fraction of each phase. It is not relevant to hard-code such information, so we rely on an external text file. For efficiency reasons, the external file must be read and processed only once at the beginning of the computation.

To do this, we rely on an external class called ExtendedPolyCrystalsSlidingSystems which implements the singleton pattern ³ [7].

This class is declared in an external header file which is included as follows:

@Includes{
#include "TFEL/Material/ExtendedPolyCrystalsSlidingSystems.hxx"
}

The description of the implementation of this class is out of the scope of the present report, so we will only focus on the main points of interest for the rest of this report.

The ExtendedPolyCrystalsSlidingSystems class has tree template arguments:

• The number of phases.
• The class describing the sliding system of one phase. This class is obviously the PolycrystalHCP_BZSlipSystems class defined previously.
• The numeric type used.

This class provides the getPolyCrystalsSlidingSystems static method as the unique way to retrieve the unique instance of this class. This method takes the name of a external text file. Each line of this file describe a phase by three Euler angles (in degree) and a volume fraction.

The unique instance of the ExtendedPolyCrystalsSlidingSystems has two data members of interest:

• mus: the tensors of directional senses, sorted by phases. Those orientation tensors are expressed in the material frame, not in frame associated with the given phase.
• volume_fractions: the volume fractions per phases.

Number of phases

We then need to indicate the number of phases (i.e. groups of grains with identical orientation).

@IntegerConstant Np  = 240;

This number must be consistent with the number of phases read in the external file given as argument to the getPolyCrystalsSlidingSystems method.

Elastic properties

The material is assumed to be isotropic. The @ComputeStiffnessTensor is used here to compute the stiffness tensor using the macroscopic Young modulus and the macroscopic Poisson ratio:

@ComputeStiffnessTensor {
  80000., 0.4
};

This keywords automatically declares and initializes a fourth order tensor named \(\underline{\underline{\mathbf{D}}}\).

Material coefficients

Using the keyword @Parameters one can then define some useful parameters such as the critical resolved shear stresses of the gliding systems or the other coefficients involved in the intra-granular constitutive laws:

//! Norton exponents per slip system family
@Parameter real m[4] = {10., 10., 10., 10.};
//! Norton stress normalisation factor
@Parameter stress K = 5.0;
//! critical resolved shear stress per slip system family
@Parameter stress τᶜ[4] = {43.0, 135.0, 75.0, 65.0};

Declaration of internal variables

We now need to define the internal variables such as macroscopic and microscopic plastic strain and its norm. This is done using the keyword @StateVariable. Those quantities are meant to evolve during the calculation :

//! macroscopic plastic strain
@StateVariable StrainStensor    epsp;
//! microscopic plastic strain
@StateVariable StrainStensor    epsg[Np]; 
//! plastic strain’s norm
@StateVariable real    pg;

Computation of the macroscopic stress

The macroscopic stress is computed from the macroscopic elastic strain and the macroscopic stiffness tensor as follows:

@ComputeStress{
  sig = D * eel;
}

This code is used to compute the stress before each evaluation of the rate of the internal state variables and and the end of the time step to compute the final stress.

Equations of the system

Using the @Derivative code block, the evolution of the internal variables can now be implemented.

@Derivative{

As described earlier, the orientation tensors of all the phases is handled by the ExtendedPolyCrystalsSlidingSystems. This class being a template, its use is tedious. For the sake of simplicity, we declare an alias PCSlidingSystems which hides some of the details associated with the use of the ExtendedPolyCrystalsSlidingSystems class:

  using PCSlidingSystems = 
    ExtendedPolyCrystalsSlidingSystems<Np, PolycrystalHCP_BZSlipSystems<real>, real>;

The following line retrieves the unique instance of the ExtendedPolyCrystalsSlidingSystems class and give the name of the file containing the description of the phases.

  const auto& gs = PCSlidingSystems::getPolyCrystalsSlidingSystems("isotrope_240_normalise.csv");

  // macroscopic equivalent stress
  const auto σᵉ = sigmaeq(σ);

The macroscopic strain rate is then initialised to zero:

  // initialisation of the macroscopic strain rate
  ∂ₜεᵛᵖ = Stensor(strainrate(0));

Then, the loop over the phases is started, as follows:

  // loop over the phases
  for (unsigned short k = 0; k != Np; ++k) {
    // volume fraction of the phase
    const auto φ = gs.volume_fractions[k];

The local stress is computed following the Berveiller-Zaoui scheme (Equation (2)):

    // local stress following the Berveiller-Zaoui scheme
    auto σᵍ = StressStensor(stress(0));
    if (2 ⋅ σᵉ + 3 ⋅ Mu ⋅ p > 0) {
      const auto α = 2 ⋅ σᵉ / (2 ⋅ σᵉ + 3 ⋅ Mu ⋅ p);
      σᵍ = σ + 2 ⋅ (1 - β) ⋅ Mu ⋅ α ⋅ (εᵛᵖ - εᵛᵖᵍ[k]);
    }

We then initialise the microscopic strain rate:

    // initialisation of the microscopic strain rate
    ∂ₜεᵛᵖᵍ[k] = Stensor(strainrate(0));

We now start the loop over the slip system of the considered phase and evaluate the contribution to the microscopic strain rate:

    for (unsigned short i = 0; i != Nss; ++i) {
      const auto idx = [&i] {
        if (i <= 2) {
          return 0;
        } else if (i >= 3 && i <= 14) {
          return 1;
        } else if (i >= 15 && i <= 17) {
          return 2;
        }
        return 3;
      }(); // return the family number 
      // orientation tensors in the material frame
      const auto& μⁱᵏ = gs.mus[k][i];
      // resolved shear stress
      const auto τ = μⁱᵏ | σᵍ;
      const auto rτ = (abs(τ) - τᶜ[idx]) / K;
      if (rτ > 0) {
        const auto sgn = τ > 0 ? 1 : -1;
        const auto ∂ₜg = sgn ⋅ pow(rτ, m[idx]);
        dεᵛᵖᵍ[k] += ∂ₜg ⋅ μⁱᵏ;
      }
    } // end of the loop over the slip systems

The macroscopic strain rate is then updated by the contribution of the considered phase:

    ∂ₜεᵛᵖ += ∂ₜεᵛᵖᵍ[k] ⋅ φ;
  } // end of the loop over the phases

After the loop over the phases, one can update the rate of the macroscopic equivalent plastic strain and the macroscopic elastic strain rate:

  // rate of the macroscopic equivalent plastic strain
  ∂ₜp = sqrt((2 ⋅ (∂ₜεᵛᵖ | ∂ₜεᵛᵖ)) / 3);
  // elastic strain rate
  ∂ₜεᵉˡ = ∂ₜεᵗᵒ - ∂ₜεᵛᵖ;

The only tangent operator available is the elastic operator. This leads to this simple implementation of the @TangentOperator code block:

@TangentOperator {
  Dt = D;
}

Unit testing

The following MTest file imposes a tensile test with a constant strain rate of \(2\,10^{-4}\, s^{-1}\) in the tensile direction:

@AccelerationAlgorithm 'Cast3M';
@StiffnessMatrixType 'Elastic' ;
@Behaviour 'src/libBehaviour.so' 'PolycrystalHCP_BZ';
@ExternalStateVariable 'Temperature' 293.15;
@Real 'tmax' 75 ;
@ImposedStrain 'EXX' {0.: 0., 75 : 0.015 };
@Times {0.,'tmax' in 1000};

Note that since the explicit behaviour does not provide a consistent tangent operator, we rely on an elastic search operator. The equilibrium iterations are accelerated using the Cast3M acceleration scheme which mimics the algorithm used by this code [3, 4].

Specifying the temperature is required, although not used by the behaviour.

The stress-strain curve obtained with the polycrystal behaviour on a simple tensile test is reported in Figure 1.

Conclusions

Extensions to the proposed implicit scheme to other well-known homogeneizations schemes are straightforward, including the work of Taylor-Lin, Kröner, Cailletaud-Pilvin (\(\beta\)-model) and their extension by Besson et al.

References