The UserObject based crystal plasticity system is designed to facilitate the implementation of different constitutive laws in a modular way. Both phenomenological constitutive models and dislocation-based constitutive models can be implemented through this system. This system consists of one material class FiniteStrainUObasedCP and four userobject classes, namely CrystalPlasticitySlipRate, CrystalPlasticitySlipResistance, CrystalPlasticityStateVarRateComponent and CrystalPlasticityStateVariable.

UserObject class

In general, users need to implement their own userobject classes by deriving from the base classes. Each userobject class represents a material property and it is created by the material class. To use a material property, call getMaterialProperty<TYPE>() in userobject.

  • CrystalPlasticityStateVariable Doxygen

    • This represents a state variable, such as mobile dislocation density.

    • The initial values of state variable are either provided in the input file (readInitialValueFromInline()) or read from the file (readInitialValueFromFile()).

    • The state variable evolves as $$$y_{n+1} = y_n + \dot{y}\cdot dt$$$ where $$$\dot{y}$$$ is the state variable evolution rate.

    • The $$$\dot{y}$$$ can have multiple components, such as $$$\dot{y} = a_1 r_1 + a_2 r_2 + \cdots$$$ where "$$$a_1, a_2, \cdots$$$" are the scale_factor and $$$r_1, r_2, \cdots$$$ are the individual rate components.

  • CrystalPlasticityStateVarRateComponent Doxygen

    • This represents individual component of the state variable evolution rate, such as $$$r_1, r_2$$$.

    • The rate component $$$r = r(y,\cdots, \dot{g}, \cdots)$$$ can be a function of state variables $$$y$$$ and slip rate $$$\dot{g}$$$.

  • CrystalPlasticitySlipResistance Doxygen

    • This represents slip system resistance.

    • The slip resistance $$$s=s(y,\cdots)$$$ can be a function of state variables $$$y$$$.

  • CrystalPlasticitySlipRate Doxygen

    • This represents slip rate.

    • The slip system is read into this class.

    • The Schmid tensor is generated in this class and the flow direction is calculated by calcFlowDirection().

    • The slip rate $$$\dot{g} = g(T, s\cdots, y\cdots)$$$ can be a function of PK2 stress $$$T$$$, slip resistance $$$s$$$ and state variables $$$y$$$.

    • Overwrite the function calcSlipRate() to calculate slip rate;

    • Overwrite the function calcSlipRateDerivative() to calculate the derivative of slip increment w.r.t. resolved shear stress.

Material class

The material class is based on plastic flow on individual slip systems to obtain the inelastic deformation in materials. The present formulation considers large deformation and is based on a stress update algorithm. Some of the important functions associated with this class are outlined below.

Function - computeQpStress:


  • 2$$$^{nd}$$$ Piola-Kirchhoff (PK2) stress in the intermediate configuration at previous increment (variable _pk2_old, symbol $$$T_n$$$);
  • Plastic deformation gradient at previous increment (variable _fp_old, symbol $$$F^p_n$$$);
  • State variables at previous increment;
  • Current deformation gradient (variable _dfgrd, symbol $$$F$$$).


  • Current PK2 stress (variable _pk2, symbol $$$T$$$);
  • Current plastic deformation gradient (variable _fp, symbol $$$F^p$$$);
  • Current slip system resistances.
  • State variables are solved iteratively using a predictor corrector algorithm. $$$T$$$ is solved using Newton Raphson.


i: At iteration i=0, $$$y^\alpha_i=y^\alpha_n$$$

ii: Calculate the residual r and Jacobian J=$$$\partial r/\partial T$$$ from the function calc_resid_jacob

iii: Update T as: $$$T_{i+1}=T_{i}-J^{-1}r$$$

iv: Check $$$||r||_2 <$$$ _rtol (Optional user defined parameter)

If False then go to Step ii


v: Calculate $$$y^\alpha_{i+1}$$$ from function updateSlipSystemResistanceAndStateVariable

vi: Check $$$max|s^\alpha_{i+1}-s^\alpha_{i}|/|s^\alpha_{i}| <$$$ _stol (Optional user defined parameter)

If False then go to Step ii


vii: Obtain rotation tensor R

viii: Update rotation as $$$\tilde{R}=RQ$$$ where Q is the rotation tensor from Euler angles

Exit Function

Function - calc_resid_jacob: User can override this function.


  • PK2 stress at previous iteration (variable _pk2, symbol $$$T_i$$$, i - iteration);
  • Inverse of plastic deformation gradient at previous increment (variable _fp_old_inv, symbol $$$F^{p-1}_n$$$);
  • Current deformation gradient (variable _dfgrd, symbol $$$F$$$).


  • Residual (variable resid, symbol r);
  • Jacobian (variable jac, symbol J);
  • Updated inverse of plastic deformation gradient (variable _fp_inv, symbol $$$F^{p-1}$$$)

Residual r as implemented:

i. Get slip rates from userobject.

ii. Resultant slip increment (variable eqv_slip_incr) $$$\Delta \gamma_{res}=I-\sum\limits_m\sum \limits_\alpha \left( \Delta \gamma^\alpha_m S_m^\alpha\right)$$$ where $$$S_m^\alpha$$$ is the flow direction.

iii. Current plastic component of deformation gradient $$$F^{p-1}=F^{p-1}_n\Delta \gamma_{res}$$$.

iv. Elastic component of deformation gradient in iteration i+1, $$$F^e=FF^{p-1}$$$.

vi: PK2 stress due to $$$T_i$$$ and associated slip increment in current iteration i+1 (variable pk2_new), $$$T^\prime=C:\frac{1}{2}\left(F^{eT}F^e-I\right)$$$.

vii. Residual $$$r=T_i-T^\prime$$$.

Jacobian J formulation:

i. $$$r=T_i-T^\prime$$$ and $$$J=\frac{\partial r}{\partial T_i}=I-\frac{\partial T^\prime}{\partial T_i}$$$.

ii. $$$T^\prime=C:E^e$$$ which gives $$$\frac{\partial T^\prime}{\partial T_i}=C: \frac{\partial E^e}{\partial T_i}$$$.

iii. $$$E^e=\frac{1}{2}\left(F^{eT}F^e -I\right)$$$ which gives in indicial notation $$$\frac{\partial E^e_{ij}}{\partial F^e_{kl}}=\frac{1}{2}\left( \delta_{il}F^e_{kj}+\delta_{jl}F^e_{ki}\right)$$$ where $$$\delta_{ij}$$$ is the Kroneckar delta.

iv. $$$F^e=FF^{p-1}$$$ which gives in indicial notation $$$\frac{\partial F^e_{ij}}{\partial F^{p-1}_{kl}}=F_{ik}\delta_{lj}$$$.

v. $$$F^{p-1}=F^{p-1}_n\Delta \gamma_{res}$$$ which gives $$$\frac{\partial F^{p-1}}{\partial \Delta \gamma_{res}}=F^{p-1}_n$$$.

vi. $$$\Delta \gamma_{res}=I-\sum\limits_m\sum \limits_\alpha \left( \Delta \gamma^\alpha_m S_m^\alpha\right)$$$ which gives $$$\frac{\partial \Delta \gamma_{res}}{\partial \Delta \gamma^\alpha_m}=-\sum\limits_\alpha S_m^\alpha$$$.

vii. $$$\Delta \gamma^\alpha_m=f(\tau_m^\alpha)$$$ which gives $$$\frac{\partial \Delta \gamma^\alpha_m}{\partial \tau^\alpha_m}$$$.

viii. $$$\tau^\alpha=T_i:S_m$$$ gives $$$\frac{\partial \tau^\alpha_m}{\partial T_i}=S_m^\alpha$$$.

ix. Hence $$$\frac{\partial T^\prime}{\partial T_i}=C: \frac{\partial E^e}{\partial F^e} \frac{\partial F^e}{\partial F^{p-1}} \frac{\partial F^{p-1}}{\partial \Delta \gamma_{res}}\sum\limits_m(\frac{\partial \Delta \gamma_{res}}{\partial \Delta \gamma^\alpha_m}\frac{\partial \Delta \gamma^\alpha_m}{\partial \tau^\alpha_m})\frac{\partial \tau^\alpha}{\partial T_i}$$$.

Jacobian J as implemented:

i. Variable dtaudpk2, $$$\frac{\partial \tau^\alpha}{\partial T_i}$$$.

ii. Variable dfpinvdslip, $$$\frac{\partial F^{p-1}}{\partial \Delta \gamma^\alpha}$$$.

iii. Variable dfedfpinv, $$$\frac{\partial F^e}{\partial F^{p-1}}$$$.

iv: Variable deedfe, $$$\frac{\partial E^e}{\partial F^e}$$$.

v. Variable dfpinvdpk2, $$$\frac{\partial F^{p-1}}{\partial T_i}=\frac{\partial F^{p-1}}{\partial \Delta \gamma^\alpha} \frac{\partial \Delta \gamma^\alpha}{\partial \tau^\alpha} \frac{\partial \tau^\alpha}{\partial T_i}$$$.

vi. $$$\frac{\partial \Delta \gamma^\alpha_m}{\partial \tau^\alpha_m}$$$ obtained from slip rate userobject.

vii. Variable jac, $$$J=C:\frac{\partial E^e}{\partial F^e} \frac{\partial F^e}{\partial F^{p-1}} \frac{\partial F^{p-1}}{\partial T_i}$$$, followed by $$$J=\mathfrak{I}-J$$$ where $$$\mathfrak{I}$$$ is the fourth order identity tensor.

Update Cauchy stress (variable sig) $$$\sigma=F^e T_i F^{eT}/J^e$$$.

Function - computeQpElasticityTensor: Defines tangent moduli K and can be used as preconditioner for JFNK. User can override this function.


  • Current plastic deformation gradient (variable _fp, symbol $$$F^p$$$);
  • Current deformation gradient (variable _dfgrd, symbol $$$F$$$)


  • Tangent moduli (variable _Jacobian_mult,symbol K).


i. $$$K=\frac{\partial \sigma}{\partial F}$$$. ii. $$$\sigma=F^eTF^{eT}/J^e$$$ which gives $$$\frac{\partial \sigma}{\partial F}=\frac{1}{J^e}\left(\frac{\partial F^e}{\partial F}TF^{eT}+F^e\frac{\partial T}{\partial F}F^{eT}+F^eT\frac{\partial F^{eT}}{\partial F}\right)$$$. iii. In indicial notation $$$\frac{\partial F^e_{ij}}{\partial F_{kl}}=\delta_{ki}F^{p-1}_{lj}$$$. iv. $$$T=C:E^e$$$ which gives $$$\frac{\partial T}{\partial F}=\frac{\partial T}{\partial E^e}\frac{\partial E^e}{\partial F^e}\frac{\partial F^e}{\partial F}=C:\frac{\partial E^e}{\partial F^e}\frac{\partial F^e}{\partial F}$$$. v. In indicial notation $$$\frac{\partial E^e_{ij}}{\partial F^e_{kl}}=\frac{1}{2}\left(\delta_{il}F^e_{kj}+\delta_{jl}F^e_{ki}\right)$$$.


i. Variable dfedf calculates $$$\frac{\partial F^e}{\partial F}$$$. ii. Variable deedfe calculates $$$\frac{\partial E^e}{\partial F^e}$$$. iii. Variable dsigspk2dfe calculates $$$F^{eT} C:\frac{\partial E^e}{\partial F^e} F^e$$$.


A phenomenological constitutive model is implemented and the test files can be found under /mooose//modules/tensor_mechanics/tests/cp_user_object/.

The UserObjects block looks like this

  # uo_xxx_name is used to fetch the corresponding material properties.
    type = CrystalPlasticitySlipRateGSS
    variable_size = 12
    slip_sys_file_name = input_slip_sys.txt
    num_slip_sys_flowrate_props = 2
    flowprops = '1 4 0.001 0.1 5 8 0.001 0.1 9 12 0.001 0.1'
    uo_state_var_name = state_var_gss
    type = CrystalPlasticitySlipResistanceGSS
    variable_size = 12
    # for phenomenological constitutive model, the state variable is taken as the same as slip resistance. 
    uo_state_var_name = state_var_gss
    type = CrystalPlasticityStateVariable
    variable_size = 12
    groups = '0 4 8 12'
    group_values = '60.8 60.8 60.8'
    uo_state_var_evol_rate_comp_name = state_var_evol_rate_comp_gss
    scale_factor = 1.0
    type = CrystalPlasticityStateVarRateComponentGSS
    variable_size = 12
    hprops = '1.0 541.5 109.8 2.5'
    uo_slip_rate_name = slip_rate_gss
    uo_state_var_name = state_var_gss

The Materials block looks like this

    type = FiniteStrainUObasedCP
    block = 0
    stol = 1e-2
    tan_mod_type = exact
    # slip rate userobject
    uo_slip_rates = 'slip_rate_gss'
    # slip resistance userobject
    uo_slip_resistances = 'slip_resistance_gss'
    # state variable userobject
    uo_state_vars = 'state_var_gss'
    # state variable evolution rate component userobject
    uo_state_var_evol_rate_comps = 'state_var_evol_rate_comp_gss'
    type = ComputeFiniteStrain
    block = 0
    displacements = 'ux uy uz'
    type = ComputeElasticityTensorCP
    block = 0
    C_ijkl = '1.684e5 1.214e5 1.214e5 1.684e5 1.214e5 1.684e5 0.754e5 0.754e5 0.754e5'
    fill_method = symmetric9