Crystal Behavior (simulation.cfg)

The mechanical behavior of the crystal implies the definition of the crystal structure and the definition of the mechanical behavior per se. If the polycrystal (as defined in the mesh file) contains several phases, then the properties of each of the phases must be defined, successively. The format is as follows:

# Material Description

  number_of_phases <nphases>

  phase 1

      crystal_type <ctype>
      [ c_over_a <ratio> ]

      { parameter definition }

  ...

  phase <nphases>

      crystal_type <ctype>
      [ c_over_a <ratio> ]

      { parameter definition }

where the crystal type (<ctype>) can be fcc, bcc, hcp and bct, for face-centered cubic, body-centered cubic, hexagonal close-packed and body-centered tetragonal, respectively. For hcp and bct, the \(c/a\) ratio (c_over_a) also needs to be provided. { parameter definition } represents a block of text that defines a mechanical behavior. The mechanical behavior of a crystal includes Anisotropic Elasticity and Crystal Plasticity.

Note

All variables presented are detailed by their dimensions (if applicable) instead of any specific unit. No unit system is inherently assumed by FEPX, and the chosen unit system and value magnitudes should be consistent with the chosen length scale for the domain. For example, if it is assumed that the length scale is mm and SI units are to be used, then [force/area] will be understood to be MPa. The unit for time, however, is always assumed to be seconds (s).

Anisotropic Elasticity

The stress, \(\sigma\), is related to the elastic strain, \(\epsilon\), via Hooke’s law:

\[\sigma = \cal C \, \epsilon,\]

where \(\cal C\) is the stiffness tensor. Using Voigt notation and the strength of materials convention, \(\sigma\) and \(\epsilon\) are written as vectors of components \(\sigma_{11}\), \(\sigma_{22}\), \(\sigma_{33}\), \(\sigma_{23}\), \(\sigma_{13}\), \(\sigma_{12}\), and \(\epsilon_{11}\), \(\epsilon_{22}\), \(\epsilon_{33}\), \(2\,\epsilon_{23}\), \(2\,\epsilon_{13}\), \(2\,\epsilon_{12}\). \(\cal C\) is symmetrical and has 21 independent components; however, its structure is greatly simplified by the crystal symmetry, and the number of independent components (hence necessary parameters) is greatly reduced. The moduli are expressed in [force/area].

Cubic symmetry

For cubic symmetry (bcc and fcc crystals), Hooke’s law expands as:

\[\begin{split}\begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{Bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{12} & & & \\ C_{12} & C_{11} & C_{12} & & & \\ C_{12} & C_{12} & C_{11} & & & \\ & & & C_{44} & & \\ & & & & C_{44} & \\ & & & & & C_{44} \end{bmatrix} \begin{Bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ 2\epsilon_{23} \\ 2\epsilon_{13} \\ 2\epsilon_{12} \end{Bmatrix}.\end{split}\]

The corresponding input is:

c11 <modulus>
c12 <modulus>
c44 <modulus>

Hexagonal symmetry

For hexagonal symmetry (hcp crystals), Hooke’s law expands as:

\[\begin{split}\begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{Bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & & & \\ C_{12} & C_{11} & C_{13} & & & \\ C_{13} & C_{13} & C_{33} & & & \\ & & & C_{44} & & \\ & & & & C_{44} & \\ & & & & & \left( C_{11}-C_{12}\right)/2 \end{bmatrix} \begin{Bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ 2\epsilon_{23} \\ 2\epsilon_{13} \\ 2\epsilon_{12} \end{Bmatrix},\end{split}\]

and the following must additionally be satisfied: \(C_{33} = C_{11} + C_{12} - C_{13}\) (which is why no input is expected for \(C_{33}\)) [1].

The corresponding input is:

c11 <modulus>
c12 <modulus>
c13 <modulus>
c44 <modulus>

Tetragonal symmetry

For tetragonal symmetry (bct crystals), Hooke’s law expands as:

\[\begin{split}\begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{Bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & & & \\ C_{12} & C_{11} & C_{13} & & & \\ C_{13} & C_{13} & C_{33} & & & \\ & & & C_{44} & & \\ & & & & C_{44} & \\ & & & & & C_{66} \end{bmatrix} \begin{Bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ 2\epsilon_{23} \\ 2\epsilon_{13} \\ 2\epsilon_{12} \end{Bmatrix},\end{split}\]

and the following must additionally be satisfied: \(C_{33} = C_{11} + C_{12} - C_{13}\) (which is why no input is expected for \(C_{33}\)) [1].

The corresponding input is:

c11 <modulus>
c12 <modulus>
c13 <modulus>
c44 <modulus>
c66 <modulus>

Note

Special attention must be paid to ensure that the correct stiffness values are chosen, to align with the input convention used here. For this convention, the Zener anisotropy ratio for cubic materials (which quantifies the level of elastic anisotropy, with 1 being perfectly isotropic) would be written as:

\[A = {2 C_{44} \over C_{11} - C_{12}}.\]

For example, Tungsten (W) is a nearly perfectly elastically isotropic cubic-symmetry material, with \(C_{11} = 522.4\) GPa, \(C_{12} = 204.4\) GPa, and \(C_{44} = 160.8\) GPa, which corresponds to a Zener ratio of \(\simeq 1.01\).

Crystal Plasticity

Slip kinetics

The kinematics of slip are described by a power law:

\[\dot{\gamma}^{\alpha} = \dot{\gamma}_{0} \left( \left| {\tau}^{\alpha} \right| \over g^{\alpha} \right)^{1/m} \rm sgn({\tau}^{\alpha}),\]

where \(\dot{\gamma}_0\) (gammadot_0) is the fixed-rate strain rate scaling coefficient (typically 1 by convention, and expressed in [1/s]), and \(m\) (m) is the rate sensitivity exponent.

The corresponding input is:

m <rate_sensitivity(ies)>
gammadot_0 <slip rate(s)>

One m value may be provided for all slip systems (or families), or several m values may be provided, each applying to a specific family of slip systems (see Slip Systems).

Hardening

Hardening models of increasing complexity are available, but are all evolutions of a base model, which corresponds to a Voce-type hardening law with isotropic hardening (same behavior on all systems). Each evolution of this model enrich one aspect of the behavior and are introduced compared to this model, but the evolutions can be combined with each other (leading to more general hardening laws). An example will be provided.

Base Model

The slip system strength evolution is modeled by:

\[\dot{g^{\alpha}} = h_{0} \left (g_s - g^{\alpha} \over g_s - g_0 \right)^{n} \dot{\gamma}, \quad \hbox{with} \quad \dot{\gamma} = \sum_{\alpha} \left|\dot{\gamma}^{\alpha}\right|,\]

where \(g_0\) (g_0) is the slip system initial strength (expressed in [force/area]), \(g_s\) (g_s) is the slip system saturation strength (expressed in [force/area]), \(h_0\) (h_0) is the fixed-state hardening rate scaling coefficient, and \(n\) (n) is the non-linear Voce hardening exponent. \(g_0\) may be defined by one value for all systems, one value per slip family, or one value per slip system (see Slip Systems). Note that this law implied isotropic hardening, while Anisotropic Hardening is described in the following.

The corresponding input is:

[ hardening saturation[,isotropic] ]

g_0 <strength(s)>
g_s <strength>
h_0 <strength>
n   <hardening exponent>

Saturation Strength Evolution

Compared to the base model, the slip system saturation strength evolves as a function of the slip activity (constant \(g_s\) becomes \(g_s(\dot\gamma)\)):

\[\dot{g^{\alpha}} = h_{0} \left (g_{s}(\dot{\gamma}) - g^{\alpha} \over g_{s}(\dot{\gamma}) - g_{0} \right)^{n} \dot{\gamma}, \quad \hbox{with} \quad \dot{\gamma} = \sum_{\alpha} \left|\dot{\gamma}^{\alpha}\right|,\]

where \(g_{s}(\dot{\gamma})\) evolves via:

\[g_{s}(\dot{\gamma}) = g_{s0} \left (\dot{\gamma} \over \dot{\gamma}_{s0} \right)^{m'},\]

which introduces three new variables: \(g_s\) (g_s) is the initial slip system saturation strength (expressed in [force/area]), \(m'\) (m_prime) is the saturation strength rate scaling exponent, and \(\dot{\gamma}_{s0}\) (gammadot_s0) is the initial saturation slip system shear rate (typically 1 by convention, and expressed in [1/s]).

Compared to the base model, replace g_s by the following input:

hardening saturation_evolution

g_s0        <strength>
gammadot_s0 <slip rate(s)>
m_prime     <strength>

Cyclic Hardening

Compared to the base model, \(\dot\gamma\) becomes \(f\), defined as follows [TURKMEN04]:

\[\dot{g^{\alpha}} = h_{0} \left (g_s - g^{\alpha} \over g_s - g_{0} \right)^{n} f, \quad \hbox{with} \quad f = \sum_{\beta = 0}^{n_{a}} \left|\dot{\gamma}^{\beta}\right|.\]

So, a slip system that contributes to hardening (\(n_{a}\) systems in total) is that which has a change in shear greater than a critical value:

\[\Delta\gamma_{crit} = a \, \left( \frac{g}{g_s} \right)^c.\]

\(a\) (cyclic_a) and \(c\) (cyclic_c) are model parameters for a critical value of accumulated shear strain used to modify the form of the Voce hardening law.

Compared to the base model, modify the hardening type and provide the \(a\) and \(c\) values:

hardening cyclic

cyclic_a    <strength(s)>
cyclic_c    <strength>

Precipitation Hardening

The base slip system strength is modified to consider the effects of the presence of precipitates. This is performed by replacing \(g_{0}\) by \(g_{0} + g_{p}\), where \(g_p\) is the additional strength contribution due to a precipitate phase:

\[\dot{g^{\alpha}} = h_{0} \left (g_s - g^{\alpha} \over g_s - (g_0 + g_p) \right)^{n} \dot{\gamma}, \quad \hbox{with} \quad \dot{\gamma} = \sum_{\alpha} \left|\dot{\gamma}^{\alpha}\right|,\]

There are two options for precipitate strengthening [COURTNEY90]. Below a critical average precipitate radius (the “cutting” regime), the contribution to strength is calculated via:

\[g_p = a_p \left(\frac{f_p \, r_p}{b_p}\right)^{\frac{1}{2}},\]

where \(f_{p}\) (f_p) is the precipitate volume fraction, \(r_p\) (r_p) is the average precipitate radius, \(b_p\) b_p is the precipitate’s Burgers’ vector, and \(a_p\) (a_p) is the cutting scaling coefficient.

Above a critical average precipitate radius (the “bowing” regime), the contribution to strength is calculated via:

\[g_p = \frac{c_p \, G_m \, b_m}{L - 2 \, r_p},\]

where \(G_m\) is the shear modulus of the matrix, \(c_p\) (c_p) the bowing scaling coefficient, and \(L\) is the average center to center spacing of the precipitates, calculated as:

\[L = \frac{r_p}{\sqrt{f_p}}.\]

The increase in strength due to the presence of precipitates is applied globally equally to all elements.

Compared to the base model, provide the following additional input:

hardening precipitation

a_p <strength>
f_p <volume_fraction>
r_p <length>
b_p <length>
c_p <coefficient>

To disable precipitate cutting, omit a_p. To disable precipitate bowing, omit c_p.

Anisotropic Hardening

Contrary to the case of isotropic hardening, for while that slip on given system generates the same hardening on all systems, and so the single crystal yield surface retains the same shape and grows isotropically [2], anisotropic hardening is such that slip on a given system generates different hardenings on the systems.

Compared to the base model, the hardening coefficient (\(h_0\)) is multiplied by the slip interaction matrix, \(h_{\alpha\beta}\):

\[\dot{g^{\alpha}} = h_{0} \, h_{\alpha \beta} \, \left (g_s - g^{\alpha} \over g_s - g_{0} \right)^{n} \dot{\gamma},\]

The diagonal entries of the matrix correspond to self-hardening, while the off-diagonal entries correspond to latent hardening. The special case of isotropic hardening corresponds to equal self-hardening and latent hardening, i.e. \(h_{\alpha \beta} = 1 \, \forall \, \alpha,\beta\). The coefficients generally depend on the type of interactions between systems, but two specific cases are implemented:

• 2-parameter interaction matrix, where the parameters are the diagonal (self-hardening) value, \(d\), and the off-diagonal (latent hardening) value, \(h\).

  The corresponding input is:

  hardening anisotropic
  
  interaction_matrix_parameters <diag> <h>
  

• Self-hardening + co-planar latent hardening, for which latent hardening is limited to coplanar slip systems (other interaction parameters are 0) [CARSON17]. In this case, the slip interaction matrix is defined by the diagonal entry, \(d\), and the off-diagonal entries, \(h_{1},\dots, h_{n}\) (depending on the crystal symmetry / number of slip planes).

  The corresponding input is:

  hardening anisotropic
  
  interaction_matrix_parameters <diag> <h1> <h2> <h3> <h4>                                    (for fcc)
  interaction_matrix_parameters <diag> <h1> <h2> <h3> <h4> <h5> <h6>                          (for bcc)
  interaction_matrix_parameters <diag> <h1> <h2> <h3> <h4> <h5> <h6> <h7>                     (for hcp)
  interaction_matrix_parameters <diag> <h1> <h2> <h3> <h4> <h5> <h6> <h7> <h8> <h9> <h10>     (for bct)
  

Coupling Model Evolutions

An example of using several evolutions of the Base Model is as follows. We introduce both Saturation Strength Evolution and Anisotropic Hardening. The slip law becomes:

\[\dot{g^{\alpha}} = h_{0} h_{\alpha\beta} \left (g_{s}(\dot{\gamma}) - g^{\alpha} \over g_{s}(\dot{\gamma}) - g_{0} \right)^{n} \dot{\gamma}, \quad \hbox{with} \quad \dot{\gamma} = \sum_{\alpha} \left|\dot{\gamma}^{\alpha}\right|.\]

where \(g_{s}(\dot{\gamma})\) is the function for the saturation strength, which evolves via:

\[g_{s}(\dot{\gamma}) = g_{s0} \left (\dot{\gamma} \over \dot{\gamma}_{s0} \right)^{m'}.\]

The corresponding input is:

hardening saturation_evolution,anisotropic

g_0                          <strength(s)>
g_s0                         <strength>
h_0                          <strength>
n                            <hardening exponent>
gammadot_s0                  <slip rate(s)>
m_prime                      <strength>
interaction_matrix_parameters <diag> <h>

[CARSON17]

18. Carson, M. Obstalecki, M. Miller, and P.R. Dawson. Characterizing heterogeneous intragranular deformations in polycrystalline solids using diffraction-based and mechanics-based metrics. Modelling and Simulation in Materials Science and Engineering, 25:055008, 2017.

[TURKMEN04]

H.S. Turkmen, M.P. Miller, P.R. Dawson, and J.C. Moosbrugger. A slip-based model for strength evolution during cyclic loading. Journal of Engineering Materials and Technology, 126(4):329-338, 2004. (Note that minor differences exist between the implemented model described above and the formulation described in the paper).

[COURTNEY90]

Thomas H. Courtney. Mechanical behavior of materials. eng. 2. ed. Long Grove, Illinois: Waveland Press, 2005. isbn: 9781577664253.

[1] (1,2)

This allows for the decoupling of the hydrostatic and deviatoric portions of the elastic deformation, in the problem resolution.

[2]

Note that, for crystals with different families of slip systems (HCP and BCT materials), this implies faster (absolute) hardening on the slip systems of larger slip strengths (specifically, for an HCP material, the slip rates on the basal, prismatic, and pyramidal slip systems will harden such that the ratios of slip strengths remain constant).