Compute Crystal Plasticity Thermal Eigenstrain

As is described in ComputeMultipleCrystalPlasticityStress, a thermal deformation gradient $var element = document.getElementById("moose-equation-30d00089-8666-40b1-97cd-00e339417064");katex.render("\\boldsymbol{F}^{\\theta}", element, {displayMode:false,throwOnError:false});$ is introduced to account for the thermal induced deformations in a finite strain thermo-mechanical problem with crystal plasticity. Accordingly, the total deformation gradient is multiplicatively decomposed into three components as $(1)var element = document.getElementById("moose-equation-9dadb1d4-0dcd-4e1a-8b9e-57f62ad9c52e");katex.render(" \\boldsymbol{F} = \\boldsymbol{F}^e \\boldsymbol{F}^p \\boldsymbol{F}^{\\theta},", element, {displayMode:true,throwOnError:false});$ such that $var element = document.getElementById("moose-equation-11e0ab61-4a85-4243-9c57-890cfb4bfd2b");katex.render("\\text{det}\\left( \\boldsymbol{F}^e \\right) > 0", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-a756bd59-1711-48ab-a980-427df47abdfb");katex.render("\\text{det}\\left( \\boldsymbol{F}^p \\right) = 1", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-33f84b3c-fae2-467d-a349-89463c6186c9");katex.render("\\text{det}\\left( \\boldsymbol{F}^{\\theta} \\right) > 0", element, {displayMode:false,throwOnError:false});$ (see Li et al. (2019), Ozturk et al. (2016), and Meissonnier et al. (2001)).

To account for more than one thermal eigenstrain, one can optionally decompose the total thermal deformation gradient $var element = document.getElementById("moose-equation-85a59211-294e-445f-b8bf-b7c3a7817752");katex.render("\\boldsymbol{F}^{\\theta}", element, {displayMode:false,throwOnError:false});$ into multiple components as $(2)var element = document.getElementById("moose-equation-90c8047d-3deb-40c1-a8de-2266b50017c0");katex.render(" \\boldsymbol{F}^{\\theta} = \\boldsymbol{F}^1 \\boldsymbol{F}^2 \\cdots \\boldsymbol{F}^{N},", element, {displayMode:true,throwOnError:false});$ where a total of $var element = document.getElementById("moose-equation-dec06065-5e4d-41c7-a27f-4571526f61bc");katex.render("N", element, {displayMode:false,throwOnError:false});$ different thermal eigenstrains are considered.

The thermal Lagrangian strain that is associated with $var element = document.getElementById("moose-equation-43dd76d9-d35b-48d4-8471-ae342ffa59c0");katex.render("\\boldsymbol{F}^{N}", element, {displayMode:false,throwOnError:false});$ is computed as $(3)var element = document.getElementById("moose-equation-3f0637cd-d49c-4357-b28c-cfb271b5b8ad");katex.render(" \\boldsymbol{E}^{N} = \\frac{1}{2}\\left({\\boldsymbol{F}^{N}}^{\\intercal} \\boldsymbol{F}^{N} - \\boldsymbol{I}\\right).", element, {displayMode:true,throwOnError:false});$

In this class, the evolution of a typical thermal deformation gradient component, $var element = document.getElementById("moose-equation-8c1912df-bdc7-476a-9de9-b74495255b1a");katex.render("\\boldsymbol{F}^N", element, {displayMode:false,throwOnError:false});$ (see Eq. (2)) is expressed with respect to its lattice symmetry axis as $(4)var element = document.getElementById("moose-equation-b857917a-a7ce-4ca6-b3f7-8b89b634a91f");katex.render(" \\dot{ \\boldsymbol{F}} ^{N} {(\\boldsymbol{F}^{N})}^{-1} = \\dot{\\theta} \\boldsymbol{\\beta},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-3fa61ea6-8eae-40dc-9869-36d292480a7b");katex.render("\\dot{\\theta}", element, {displayMode:false,throwOnError:false});$ is the temperature rate, and $var element = document.getElementById("moose-equation-6c661cfa-1fab-47fe-82c6-bd084dcb6a6d");katex.render("\\boldsymbol{\\beta} = \\text{diag} \\left( \\beta_1, \\beta_2, \\beta_3 \\right)", element, {displayMode:false,throwOnError:false});$ is a diagonal tensor for anisotropic thermal expansion coefficients. One can either use the above evolution equation, or create a different crystal plasticity thermal eigenstrain class with a customized constitutive equation by inheriting from the ComputeCrystalPlasticityEigenstrainBase class.

note:Base Class Requirement

Any constitutive eigenstrain model developed for use within the ComputeMultipleCrystalPlasticityStress class must inherit from the ComputeCrystalPlasticityEigenstrainBase class.

We summarize the numerical implementation details that are specialized in crystal plasticity thermal eigenstrain computation in the following.

Numerical Implementation

The calculations of the thermal deformation gradient $var element = document.getElementById("moose-equation-49ffd11e-b407-462d-ad18-33c170966119");katex.render("\\boldsymbol{F}^{\\theta}", element, {displayMode:false,throwOnError:false});$ (or its component $var element = document.getElementById("moose-equation-7b3c97df-4165-465b-abcd-9979e27515e4");katex.render("\\boldsymbol{F}^{N}", element, {displayMode:false,throwOnError:false});$ ), the Lagrangian strain $var element = document.getElementById("moose-equation-61b9f258-1fc6-4905-a44a-af3d5898b01b");katex.render("\\boldsymbol{E}^{\\theta}", element, {displayMode:false,throwOnError:false});$ (or $var element = document.getElementById("moose-equation-4923ba51-6ec7-4eb5-888b-a98717482fed");katex.render("\\boldsymbol{E}^{N}", element, {displayMode:false,throwOnError:false});$ ), and the associated derivative $var element = document.getElementById("moose-equation-d3dc0bed-37ab-4c97-950f-a998b2ef4c52");katex.render("\\partial \\boldsymbol{F}^{\\theta} / \\partial \\theta", element, {displayMode:false,throwOnError:false});$ (or $var element = document.getElementById("moose-equation-c04fca81-fe6f-4bf4-b2ec-23ae65a54ed5");katex.render("\\partial \\boldsymbol{F}^{N} / \\partial \\theta", element, {displayMode:false,throwOnError:false});$ ) are included in this ComputeCrystalPlasticityThermalEigenstrain class. Meanwhile, changes are required in the residual and Jacobian calculations in ComputeMultipleCrystalPlasticityStress in order to account for the additional configuration due to thermal expansion. The implementation details are described below.

1. Time integration for $var element = document.getElementById("moose-equation-5fb8da8d-4a70-4a50-bc6d-af86e2c9eef2");katex.render("\\boldsymbol{F}^{N}", element, {displayMode:false,throwOnError:false});$

The evolution of the thermal deformation gradient component $var element = document.getElementById("moose-equation-276a5bd7-17d7-47e0-b537-ace1840431d5");katex.render("\\boldsymbol{F} ^{N}", element, {displayMode:false,throwOnError:false});$ is defined by Eq. (4). Similar to the time integration that has been implemented for $var element = document.getElementById("moose-equation-d4799d26-a9d2-4cb5-9956-c009e6bd25f7");katex.render("\\boldsymbol{F}^p", element, {displayMode:false,throwOnError:false});$ , we start with a backward difference expression as an approximation for $var element = document.getElementById("moose-equation-d96c9424-943e-4189-8ae8-aacfb8315f73");katex.render("\\boldsymbol{F}^N", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-6704a4ca-2b6e-43aa-bccf-0103c1006693");katex.render(" {\\boldsymbol{F}^{N}_{n+1}} \\approx \\frac{ \\boldsymbol{F} ^{N}_{n+1} - \\boldsymbol{F} ^{N}_{n} }{\\Delta t},", element, {displayMode:true,throwOnError:false});$ where the subscripts $var element = document.getElementById("moose-equation-7cdbbd88-9ff9-49be-bd88-cabc6578cc13");katex.render("n", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-0a997b12-6757-4074-afe4-a0562ef56e8b");katex.render("n+1", element, {displayMode:false,throwOnError:false});$ identifies the time steps. By substituting into Eq. (4) and rearranging, we have $var element = document.getElementById("moose-equation-b4a3820b-dd0f-4f0e-97c3-aba8cc0b1642");katex.render(" {(\\boldsymbol{F}^{N}_{n+1})}^{-1} = {(\\boldsymbol{F}^{N}_{n})}^{-1} \\left( \\boldsymbol{I} - \\Delta \\theta \\boldsymbol{\\beta} \\right)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-3c3a25eb-1d4b-4e0d-995e-911feec1464b");katex.render("\\Delta \\theta = \\dot{\\theta} \\Delta t", element, {displayMode:false,throwOnError:false});$ .

2. Crystal plasticity Jacobian

The stress is updated iteratively by driving the stress residual to zero. The stress residual is computed by $var element = document.getElementById("moose-equation-179bd643-ebe5-4977-ab07-e961efdd076c");katex.render(" \\mathbb{R} = \\boldsymbol{S} - \\mathbb{C} : \\boldsymbol{E}^e,", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-e310d18e-c0ac-4018-a554-fb1ac538ed55");katex.render("\\boldsymbol{S}", element, {displayMode:false,throwOnError:false});$ is the second Piola-Kirchhoff stress, the $var element = document.getElementById("moose-equation-4ce6f4db-c188-49bb-9f5c-3a60b1e9f3d8");katex.render("\\mathbb{C}", element, {displayMode:false,throwOnError:false});$ is the elasticity tensor, and $var element = document.getElementById("moose-equation-66cabde8-330a-409e-8497-024e85a535f3");katex.render("\\boldsymbol{E}^e", element, {displayMode:false,throwOnError:false});$ is the elastic part of the Cauchy-Green deformation tensor. The crystal plasticity Jacobian is computed by taking the derivative of $var element = document.getElementById("moose-equation-71bd880e-0df0-44c3-aa9a-1a4cd1339ccb");katex.render("\\mathbb{R}", element, {displayMode:false,throwOnError:false});$ with respect to $var element = document.getElementById("moose-equation-139c6394-b439-4538-b981-6edb6398c6db");katex.render("\\boldsymbol{S}", element, {displayMode:false,throwOnError:false});$ . By using the chain rule, we have $(5)var element = document.getElementById("moose-equation-73ee4dfa-d098-4508-97bb-3af364b28485");katex.render(" \\mathbb{J} = \\mathbb{I} - \\mathbb{C} : \\left\\{\\frac{\\partial\\boldsymbol{E}^e}{\\partial\\boldsymbol{F}^e}\\frac{\\partial\\boldsymbol{F}^e }{\\partial{\\boldsymbol{F}^p}^{-1}}\\frac{\\partial{\\boldsymbol{F}^p}^{-1}}{\\partial\\boldsymbol{S}}\\right\\}.", element, {displayMode:true,throwOnError:false});$ The computation of each term in the curly bracket will be described in the follows using indicial notations. The $var element = document.getElementById("moose-equation-566d6ac1-06dc-475d-961b-a68a337efcab");katex.render("{\\partial\\boldsymbol{E}^e_{ij}}/{\\partial \\boldsymbol{F}^e_{kl}}", element, {displayMode:false,throwOnError:false});$ is computed as $var element = document.getElementById("moose-equation-bc0f9898-1c96-4e93-8920-d7cc67d978f2");katex.render("\\frac{\\partial \\boldsymbol{E}^e_{ij}}{ \\partial \\boldsymbol{F}^e_{kl}} = \\frac{1}{2} \\delta_{il} \\boldsymbol{F}^e_{kj} + \\frac{1}{2} \\delta_{jl} {\\boldsymbol{F}^e_{ik}}^{\\intercal}.", element, {displayMode:true,throwOnError:false});$

From Eq. (1), $var element = document.getElementById("moose-equation-4df724ba-302c-42d2-977e-fcdc92e75a68");katex.render("\\boldsymbol{F}^e_{ij} = \\boldsymbol{F}_{im} {\\boldsymbol{F}^{\\theta}}^{-1}_{mn} {\\boldsymbol{F}^{p}}^{-1}_{nj}", element, {displayMode:false,throwOnError:false});$ . Therefore, $var element = document.getElementById("moose-equation-dac5229e-237d-49d7-953a-3a5603e22c1f");katex.render("{\\partial \\boldsymbol{F}^e_{ij} }/{\\partial{\\boldsymbol{F}^p}^{-1}_{kj}}=\\boldsymbol{F}_{im}{\\boldsymbol{F}^{\\theta}}^{-1}_{mk}", element, {displayMode:false,throwOnError:false});$ .

we evaluate the last term in Eq. (5) via chain rule: $var element = document.getElementById("moose-equation-d3405c97-6ab5-40ba-b132-c47411588686");katex.render(" \\frac{\\partial {\\boldsymbol{F}^p}^{-1} }{\\partial \\boldsymbol{S}} = \\sum_\\alpha \\left\\{ \\frac{\\partial {\\boldsymbol{F}^p}^{-1}}{\\partial \\gamma^{\\alpha} } \\frac {\\partial \\gamma^{\\alpha} }{\\partial \\tau^\\alpha} \\frac{\\partial \\tau^\\alpha}{\\partial \\boldsymbol{S}} \\right\\},", element, {displayMode:true,throwOnError:false});$ where the evaluation of $var element = document.getElementById("moose-equation-62a7f50f-dec8-47a1-823b-9c4ab1004d5d");katex.render("{\\partial {\\boldsymbol{F}^p}^{-1}}/{\\partial \\gamma^{\\alpha} }", element, {displayMode:false,throwOnError:false});$ is based on $var element = document.getElementById("moose-equation-f0e9bea8-74ac-4b04-b105-94529be686c6");katex.render("\\boldsymbol{L^p} =\\dot{\\boldsymbol{F}}^p \\boldsymbol{F}^{p-1}", element, {displayMode:false,throwOnError:false});$ by taking chain rule. The evaluation of $var element = document.getElementById("moose-equation-14d1639e-230b-4c83-bd58-97b1ab31cb58");katex.render("{\\partial \\gamma^{\\alpha} } / {\\partial \\tau^\\alpha}", element, {displayMode:false,throwOnError:false});$ varies based on the form of the flow rule that differs for different material types. The last term is evaluated as $var element = document.getElementById("moose-equation-25ba0dbe-b7db-41ba-9746-1464c20fe24c");katex.render(" \\frac{\\partial \\tau^\\alpha}{\\partial \\boldsymbol{S}} = \\text{det}(\\boldsymbol{F}^{\\theta}) {\\boldsymbol{F}^\\theta} \\boldsymbol{S}_{o}^{\\alpha} {\\boldsymbol{F}^{\\theta}}^{-1}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-2bcffbe8-1870-4e63-b38a-04b45b5f1bee");katex.render("\\boldsymbol{S}^{\\alpha}_{o}", element, {displayMode:false,throwOnError:false});$ is the schmid tensor for slip system $var element = document.getElementById("moose-equation-980ead86-3036-43cd-b5db-a4699b009723");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ .

3. Elasto-plastic tangent moduli

The elasto-plastic tangent moduli is computed via $var element = document.getElementById("moose-equation-db13a639-91be-47c3-805d-64211db0e2fb");katex.render(" \\frac{\\partial\\boldsymbol{\\sigma}}{ \\partial \\boldsymbol{F}} = \\left( \\frac{\\partial \\boldsymbol{\\sigma}}{\\partial \\boldsymbol{S}} \\frac{\\partial \\boldsymbol{S}}{\\partial \\boldsymbol{F}^{e}} + \\frac{\\partial \\boldsymbol{\\sigma}}{\\partial \\boldsymbol{F}^e} \\right) \\frac{\\partial \\boldsymbol{F}^e}{\\partial \\boldsymbol{F}},", element, {displayMode:true,throwOnError:false});$ where the $var element = document.getElementById("moose-equation-b67ddcfe-8e4e-429c-b020-f31416be747f");katex.render("{\\partial \\boldsymbol{\\sigma}}/{\\partial \\boldsymbol{S}}", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-eddd1267-51a7-414d-8ea2-792f9ebff781");katex.render("{\\partial \\boldsymbol{S}}/{\\partial \\boldsymbol{F}^{e}}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-5194c455-259a-489c-a6a2-3c883ff84678");katex.render("{\\partial \\boldsymbol{\\sigma}}/{\\partial \\boldsymbol{F}^e}", element, {displayMode:false,throwOnError:false});$ remains the same with the inclusion of the thermal deformation gradient. The last term is evaluated as $var element = document.getElementById("moose-equation-6e741037-ba62-47c3-814f-bfacf34f689d");katex.render(" \\frac{\\partial \\boldsymbol{F}^e_{ij}}{\\partial \\boldsymbol{F}_{il}} = \\boldsymbol{F}^{\\theta -1 }_{ln} \\boldsymbol{F}^{p -1 }_{nj}.", element, {displayMode:true,throwOnError:false});$

note:Off-diagonal Jacobian

The off-diagonal Jacobian contributions due to thermal expansion is currently not accounted for in the stress divergence calculation. Improving the Jacobian entries associated with thermal expansion remains an ongoing work.

Example Input File Syntax

[stress]
  type = ComputeMultipleCrystalPlasticityStress
  crystal_plasticity_models = 'trial_xtalpl'
  eigenstrain_names = "thermal_eigenstrain_1 thermal_eigenstrain_2"
  tan_mod_type = exact
  maximum_substep_iteration = 5
[]

[thermal_eigenstrain_1]
  type = ComputeCrystalPlasticityThermalEigenstrain
  eigenstrain_name = thermal_eigenstrain_1
  deformation_gradient_name = thermal_deformation_gradient_1
  temperature = temperature
  thermal_expansion_coefficients = '1e-05 2e-05 3e-05' # thermal expansion coefficients along three directions
[]

[thermal_eigenstrain_2]
  type = ComputeCrystalPlasticityThermalEigenstrain
  eigenstrain_name = thermal_eigenstrain_2
  deformation_gradient_name = thermal_deformation_gradient_2
  temperature = temperature
  thermal_expansion_coefficients = '2e-05 3e-05 4e-05' # thermal expansion coefficients along three directions
[]

Input Parameters

deformation_gradient_nameMaterial property name for the deformation gradient tensor computed by this model.
C++ Type:std::string
Controllable:No
Description:Material property name for the deformation gradient tensor computed by this model.
eigenstrain_nameMaterial property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
C++ Type:std::string
Controllable:No
Description:Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
thermal_expansion_coefficientsVector of values defining the constant second order thermal expansion coefficients, depending on the degree of anisotropy, this should be of size 1, 3, 6 or 9
C++ Type:std::vector<double>
Controllable:No
Description:Vector of values defining the constant second order thermal expansion coefficients, depending on the degree of anisotropy, this should be of size 1, 3, 6 or 9

Required Parameters

base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
temperatureCoupled temperature variable
C++ Type:std::vector<VariableName>
Controllable:No
Description:Coupled temperature variable

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

References

Jifeng Li, Ignacio Romero, and Javier Segurado. Development of a thermo-mechanically coupled crystal plasticity modeling framework: application to polycrystalline homogenization. International Journal of Plasticity, 119:313–330, 2019.

@article{li2019development,
    author = "Li, Jifeng and Romero, Ignacio and Segurado, Javier",
    title = "Development of a thermo-mechanically coupled crystal plasticity modeling framework: application to polycrystalline homogenization",
    journal = "International Journal of Plasticity",
    volume = "119",
    pages = "313--330",
    year = "2019",
    publisher = "Elsevier"
}

FT Meissonnier, EP Busso, and NP O'Dowd. Finite element implementation of a generalised non-local rate-dependent crystallographic formulation for finite strains. International Journal of Plasticity, 17(4):601–640, 2001.

@article{meissonnier2001finite,
    author = "Meissonnier, FT and Busso, EP and O'Dowd, NP",
    title = "Finite element implementation of a generalised non-local rate-dependent crystallographic formulation for finite strains",
    journal = "International Journal of Plasticity",
    volume = "17",
    number = "4",
    pages = "601--640",
    year = "2001",
    publisher = "Elsevier"
}

Deniz Ozturk, Ahmad Shahba, and Somnath Ghosh. Crystal plasticity FE study of the effect of thermo-mechanical loading on fatigue crack nucleation in titanium alloys. Fatigue & Fracture of Engineering Materials & Structures, 39(6):752–769, 2016.

@article{ozturk2016crystal,
    author = "Ozturk, Deniz and Shahba, Ahmad and Ghosh, Somnath",
    title = "Crystal plasticity {FE} study of the effect of thermo-mechanical loading on fatigue crack nucleation in titanium alloys",
    journal = "Fatigue \\& Fracture of Engineering Materials \\& Structures",
    volume = "39",
    number = "6",
    pages = "752--769",
    year = "2016",
    publisher = "Wiley Online Library"
}

(moose/modules/tensor_mechanics/test/tests/crystal_plasticity/stress_update_material_based/multiple_eigenstrains_test.i)

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  type = GeneratedMesh
  dim = 3
  elem_type = HEX8
[]

[AuxVariables]
  [temperature]
    order = FIRST
    family = LAGRANGE
  []
  [f1_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [f1_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [f1_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [f2_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [f2_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [f2_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [feig_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [feig_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [feig_zz]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Modules/TensorMechanics/Master/all]
  strain = FINITE
  add_variables = true
  generate_output = stress_zz
[]

[AuxKernels]
  [temperature]
    type = FunctionAux
    variable = temperature
    function = '300+400*t' # temperature increases at a constant rate
    execute_on = timestep_begin
  []
  [f1_xx]
    type = RankTwoAux
    variable = f1_xx
    rank_two_tensor = thermal_deformation_gradient_1
    index_j = 0
    index_i = 0
    execute_on = timestep_end
  []
  [f1_yy]
    type = RankTwoAux
    variable = f1_yy
    rank_two_tensor = thermal_deformation_gradient_1
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [f1_zz]
    type = RankTwoAux
    variable = f1_zz
    rank_two_tensor = thermal_deformation_gradient_1
    index_j = 2
    index_i = 2
    execute_on = timestep_end
  []
  [f2_xx]
    type = RankTwoAux
    variable = f2_xx
    rank_two_tensor = thermal_deformation_gradient_2
    index_j = 0
    index_i = 0
    execute_on = timestep_end
  []
  [f2_yy]
    type = RankTwoAux
    variable = f2_yy
    rank_two_tensor = thermal_deformation_gradient_2
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [f2_zz]
    type = RankTwoAux
    variable = f2_zz
    rank_two_tensor = thermal_deformation_gradient_2
    index_j = 2
    index_i = 2
    execute_on = timestep_end
  []
  [feig_xx]
    type = RankTwoAux
    variable = feig_xx
    rank_two_tensor = eigenstrain_deformation_gradient
    index_j = 0
    index_i = 0
    execute_on = timestep_end
  []
  [feig_yy]
    type = RankTwoAux
    variable = feig_yy
    rank_two_tensor = eigenstrain_deformation_gradient
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [feig_zz]
    type = RankTwoAux
    variable = feig_zz
    rank_two_tensor = eigenstrain_deformation_gradient
    index_j = 2
    index_i = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [symmy]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  []
  [symmx]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0
  []
  [symmz]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0
  []
  [tdisp]
    type = DirichletBC
    variable = disp_z
    boundary = front
    value = 0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeElasticityTensorCP
    C_ijkl = '1.684e5 1.214e5 1.214e5 1.684e5 1.214e5 1.684e5 0.754e5 0.754e5 0.754e5'
    fill_method = symmetric9
  []
  [stress]
    type = ComputeMultipleCrystalPlasticityStress
    crystal_plasticity_models = 'trial_xtalpl'
    eigenstrain_names = "thermal_eigenstrain_1 thermal_eigenstrain_2"
    tan_mod_type = exact
    maximum_substep_iteration = 5
  []
  [trial_xtalpl]
    type = CrystalPlasticityKalidindiUpdate
    number_slip_systems = 12
    slip_sys_file_name = input_slip_sys.txt
  []
  [thermal_eigenstrain_1]
    type = ComputeCrystalPlasticityThermalEigenstrain
    eigenstrain_name = thermal_eigenstrain_1
    deformation_gradient_name = thermal_deformation_gradient_1
    temperature = temperature
    thermal_expansion_coefficients = '1e-05 2e-05 3e-05' # thermal expansion coefficients along three directions
  []
  [thermal_eigenstrain_2]
    type = ComputeCrystalPlasticityThermalEigenstrain
    eigenstrain_name = thermal_eigenstrain_2
    deformation_gradient_name = thermal_deformation_gradient_2
    temperature = temperature
    thermal_expansion_coefficients = '2e-05 3e-05 4e-05' # thermal expansion coefficients along three directions
  []
[]

[Postprocessors]
  [stress_zz]
    type = ElementAverageValue
    variable = stress_zz
  []
  [f1_xx]
    type = ElementAverageValue
    variable = f1_xx
  []
  [f1_yy]
    type = ElementAverageValue
    variable = f1_yy
  []
  [f1_zz]
    type = ElementAverageValue
    variable = f1_zz
  []
  [f2_xx]
    type = ElementAverageValue
    variable = f2_xx
  []
  [f2_yy]
    type = ElementAverageValue
    variable = f2_yy
  []
  [f2_zz]
    type = ElementAverageValue
    variable = f2_zz
  []
  [feig_xx]
    type = ElementAverageValue
    variable = feig_xx
  []
  [feig_yy]
    type = ElementAverageValue
    variable = feig_yy
  []
  [feig_zz]
    type = ElementAverageValue
    variable = feig_zz
  []
  [temperature]
    type = ElementAverageValue
    variable = temperature
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' asm      2              lu            gmres     200'
  nl_abs_tol = 1e-10
  nl_rel_tol = 1e-10
  nl_abs_step_tol = 1e-10

  dt = 0.1
  dtmin = 1e-4
  end_time = 10
[]

[Outputs]
  csv = true
  [console]
    type = Console
    max_rows = 5
  []
[]