Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "ComputeLagrangianObjectiveStress.h"
      11             : 
      12             : #include "FactorizedRankTwoTensor.h"
      13             : 
      14             : InputParameters
      15        2656 : ComputeLagrangianObjectiveStress::validParams()
      16             : {
      17        2656 :   InputParameters params = ComputeLagrangianStressCauchy::validParams();
      18             : 
      19        2656 :   params.addClassDescription("Stress update based on the small (engineering) stress");
      20             : 
      21        5312 :   MooseEnum objectiveRate("truesdell jaumann green_naghdi", "truesdell");
      22        5312 :   params.addParam<MooseEnum>(
      23             :       "objective_rate", objectiveRate, "Which type of objective integration to use");
      24             : 
      25        2656 :   return params;
      26        2656 : }
      27             : 
      28        1992 : ComputeLagrangianObjectiveStress::ComputeLagrangianObjectiveStress(
      29        1992 :     const InputParameters & parameters)
      30             :   : ComputeLagrangianStressCauchy(parameters),
      31        1992 :     _small_stress(declareProperty<RankTwoTensor>(_base_name + "small_stress")),
      32        3984 :     _small_stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "small_stress")),
      33        1992 :     _small_jacobian(declareProperty<RankFourTensor>(_base_name + "small_jacobian")),
      34        3984 :     _cauchy_stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "cauchy_stress")),
      35        3984 :     _mechanical_strain(getMaterialPropertyByName<RankTwoTensor>(_base_name + "mechanical_strain")),
      36        3984 :     _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
      37        1992 :     _vorticity_increment(
      38        1992 :         getMaterialPropertyByName<RankTwoTensor>(_base_name + "vorticity_increment")),
      39        3984 :     _def_grad(getMaterialPropertyByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      40        3984 :     _def_grad_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      41        3984 :     _rate(getParam<MooseEnum>("objective_rate").getEnum<ObjectiveRate>()),
      42        1992 :     _polar_decomp(_rate == ObjectiveRate::GreenNaghdi),
      43        1992 :     _rotation(_polar_decomp ? &declareProperty<RankTwoTensor>(_base_name + "rotation") : nullptr),
      44        2025 :     _rotation_old(_polar_decomp ? &getMaterialPropertyOld<RankTwoTensor>(_base_name + "rotation")
      45             :                                 : nullptr),
      46        1992 :     _d_rotation_d_def_grad(
      47        4149 :         _polar_decomp ? &declareProperty<RankFourTensor>(derivativePropertyName(
      48        2058 :                             _base_name + "rotation", {_base_name + "deformation_gradient"}))
      49             :                       : nullptr),
      50        3984 :     _stretch(_polar_decomp ? &declareProperty<RankTwoTensor>(_base_name + "stretch") : nullptr)
      51             : {
      52        2025 : }
      53             : 
      54             : void
      55      524618 : ComputeLagrangianObjectiveStress::initQpStatefulProperties()
      56             : {
      57      524618 :   ComputeLagrangianStressBase::initQpStatefulProperties();
      58             : 
      59      524618 :   _small_stress[_qp].zero();
      60      524618 :   _cauchy_stress[_qp].zero();
      61             : 
      62      524618 :   if (_polar_decomp)
      63        2120 :     (*_rotation)[_qp] = RankTwoTensor::Identity();
      64      524618 : }
      65             : 
      66             : void
      67    30057604 : ComputeLagrangianObjectiveStress::computeQpCauchyStress()
      68             : {
      69    30057604 :   computeQpSmallStress();
      70             : 
      71    30057604 :   if (!_large_kinematics)
      72             :   {
      73    15702076 :     _cauchy_stress[_qp] = _small_stress[_qp];
      74    15702076 :     _cauchy_jacobian[_qp] = _small_jacobian[_qp];
      75             :   }
      76             :   else
      77             :   {
      78             :     // If large_kinematics = true, do the objective integration
      79    14355528 :     RankTwoTensor dS = _small_stress[_qp] - _small_stress_old[_qp];
      80             : 
      81    14355528 :     if (_rate == ObjectiveRate::Truesdell)
      82    13497228 :       _cauchy_stress[_qp] = objectiveUpdateTruesdell(dS);
      83      858300 :     else if (_rate == ObjectiveRate::Jaumann)
      84      429960 :       _cauchy_stress[_qp] = objectiveUpdateJaumann(dS);
      85      428340 :     else if (_rate == ObjectiveRate::GreenNaghdi)
      86      428340 :       _cauchy_stress[_qp] = objectiveUpdateGreenNaghdi(dS);
      87             :     else
      88           0 :       mooseError("Internal error: unsupported objective rate.");
      89             :   }
      90    30057604 : }
      91             : 
      92             : RankTwoTensor
      93    13497228 : ComputeLagrangianObjectiveStress::objectiveUpdateTruesdell(const RankTwoTensor & dS)
      94             : {
      95             :   // Get the kinematic tensor
      96    13497228 :   RankTwoTensor dL = RankTwoTensor::Identity() - _inv_df[_qp];
      97             : 
      98             :   // Update the Cauchy stress
      99    13497228 :   auto [S, Jinv] = advectStress(_cauchy_stress_old[_qp] + dS, dL);
     100             : 
     101             :   // Get the appropriate tangent tensor
     102    13497228 :   RankFourTensor U = stressAdvectionDerivative(S);
     103    13497228 :   _cauchy_jacobian[_qp] = cauchyJacobian(Jinv, U);
     104             : 
     105    13497228 :   return S;
     106             : }
     107             : 
     108             : RankTwoTensor
     109      429960 : ComputeLagrangianObjectiveStress::objectiveUpdateJaumann(const RankTwoTensor & dS)
     110             : {
     111             :   usingTensorIndices(i, j, k, l);
     112             : 
     113             :   // Get the kinematic tensor
     114      429960 :   RankTwoTensor dL = RankTwoTensor::Identity() - _inv_df[_qp];
     115      429960 :   RankTwoTensor dW = 0.5 * (dL - dL.transpose());
     116             : 
     117             :   // Update the Cauchy stress
     118      429960 :   auto [S, Jinv] = advectStress(_cauchy_stress_old[_qp] + dS, dW);
     119             : 
     120             :   // Get the appropriate tangent tensor
     121             :   RankTwoTensor I = RankTwoTensor::Identity();
     122      429960 :   RankFourTensor ddW_ddL = 0.5 * (I.times<i, k, j, l>(I) - I.times<i, l, j, k>(I));
     123      429960 :   RankFourTensor U = stressAdvectionDerivative(S) * ddW_ddL;
     124      429960 :   _cauchy_jacobian[_qp] = cauchyJacobian(Jinv, U);
     125             : 
     126      429960 :   return S;
     127             : }
     128             : 
     129             : RankTwoTensor
     130      428340 : ComputeLagrangianObjectiveStress::objectiveUpdateGreenNaghdi(const RankTwoTensor & dS)
     131             : {
     132             :   usingTensorIndices(i, j, k, l, m);
     133             : 
     134             :   // The kinematic tensor for the Green-Naghdi rate is
     135             :   // Omega = dot(R) R^T
     136      428340 :   polarDecomposition();
     137             :   RankTwoTensor I = RankTwoTensor::Identity();
     138      428340 :   RankTwoTensor dR = (*_rotation)[_qp] * (*_rotation_old)[_qp].transpose() - I;
     139      428340 :   RankTwoTensor dO = dR * _inv_df[_qp];
     140             : 
     141             :   // Update the Cauchy stress
     142      428340 :   auto [S, Jinv] = advectStress(_cauchy_stress_old[_qp] + dS, dO);
     143             : 
     144             :   // Get the appropriate tangent tensor
     145      428340 :   RankFourTensor d_R_d_F = (*_d_rotation_d_def_grad)[_qp];
     146      428340 :   RankFourTensor d_F_d_dL = _inv_df[_qp].inverse().times<i, k, l, j>(_def_grad[_qp]);
     147      428340 :   RankTwoTensor T = (*_rotation_old)[_qp].transpose() * _inv_df[_qp];
     148             :   RankFourTensor d_dO_d_dL =
     149      428340 :       T.times<m, j, i, m, k, l>(d_R_d_F * d_F_d_dL) - dR.times<i, k, j, l>(I);
     150      428340 :   RankFourTensor U = stressAdvectionDerivative(S) * d_dO_d_dL;
     151      428340 :   _cauchy_jacobian[_qp] = cauchyJacobian(Jinv, U);
     152             : 
     153      428340 :   return S;
     154             : }
     155             : 
     156             : std::tuple<RankTwoTensor, RankFourTensor>
     157    14355528 : ComputeLagrangianObjectiveStress::advectStress(const RankTwoTensor & S0,
     158             :                                                const RankTwoTensor & dQ) const
     159             : {
     160    14355528 :   RankFourTensor J = updateTensor(dQ);
     161    14355528 :   RankFourTensor Jinv = J.inverse();
     162    14355528 :   RankTwoTensor S = Jinv * S0;
     163    14355528 :   return {S, Jinv};
     164             : }
     165             : 
     166             : RankFourTensor
     167    14355528 : ComputeLagrangianObjectiveStress::updateTensor(const RankTwoTensor & dQ) const
     168             : {
     169             :   auto I = RankTwoTensor::Identity();
     170             :   usingTensorIndices(i, j, k, l);
     171    28711056 :   return (1.0 + dQ.trace()) * I.times<i, k, j, l>(I) - dQ.times<i, k, j, l>(I) -
     172    28711056 :          I.times<i, k, j, l>(dQ);
     173             : }
     174             : 
     175             : RankFourTensor
     176    14355528 : ComputeLagrangianObjectiveStress::stressAdvectionDerivative(const RankTwoTensor & S) const
     177             : {
     178             :   auto I = RankTwoTensor::Identity();
     179             :   usingTensorIndices(i, j, k, l);
     180    14355528 :   return S.times<i, j, k, l>(I) - I.times<i, k, l, j>(S) - S.times<i, l, j, k>(I);
     181             : }
     182             : 
     183             : RankFourTensor
     184    14355528 : ComputeLagrangianObjectiveStress::cauchyJacobian(const RankFourTensor & Jinv,
     185             :                                                  const RankFourTensor & U) const
     186             : {
     187    14355528 :   return Jinv * (_small_jacobian[_qp] - U);
     188             : }
     189             : 
     190             : void
     191      428340 : ComputeLagrangianObjectiveStress::polarDecomposition()
     192             : {
     193      428340 :   FactorizedRankTwoTensor C = _def_grad[_qp].transpose() * _def_grad[_qp];
     194      856680 :   (*_stretch)[_qp] = MathUtils::sqrt(C).get();
     195      428340 :   RankTwoTensor Uinv = MathUtils::sqrt(C).inverse().get();
     196      428340 :   (*_rotation)[_qp] = _def_grad[_qp] * Uinv;
     197             : 
     198             :   // Derivative of rotation w.r.t. the deformation gradient
     199             :   RankTwoTensor I = RankTwoTensor::Identity();
     200      428340 :   RankTwoTensor Y = (*_stretch)[_qp].trace() * I - (*_stretch)[_qp];
     201      428340 :   RankTwoTensor Z = (*_rotation)[_qp] * Y;
     202      428340 :   RankTwoTensor O = Z * (*_rotation)[_qp].transpose();
     203             :   usingTensorIndices(i, j, k, l);
     204      856680 :   (*_d_rotation_d_def_grad)[_qp] = (O.times<i, k, l, j>(Y) - Z.times<i, l, k, j>(Z)) / Y.det();
     205      428340 : }