Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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        2384 : ComputeLagrangianObjectiveStress::validParams()
      16             : {
      17        2384 :   InputParameters params = ComputeLagrangianStressCauchy::validParams();
      18             : 
      19        2384 :   params.addClassDescription("Stress update based on the small (engineering) stress");
      20             : 
      21        4768 :   MooseEnum objectiveRate("truesdell jaumann green_naghdi", "truesdell");
      22        4768 :   params.addParam<MooseEnum>(
      23             :       "objective_rate", objectiveRate, "Which type of objective integration to use");
      24             : 
      25        2384 :   return params;
      26        2384 : }
      27             : 
      28        1788 : ComputeLagrangianObjectiveStress::ComputeLagrangianObjectiveStress(
      29        1788 :     const InputParameters & parameters)
      30             :   : ComputeLagrangianStressCauchy(parameters),
      31        1788 :     _small_stress(declareProperty<RankTwoTensor>(_base_name + "small_stress")),
      32        3576 :     _small_stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "small_stress")),
      33        1788 :     _small_jacobian(declareProperty<RankFourTensor>(_base_name + "small_jacobian")),
      34        3576 :     _cauchy_stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "cauchy_stress")),
      35        3576 :     _mechanical_strain(getMaterialPropertyByName<RankTwoTensor>(_base_name + "mechanical_strain")),
      36        3576 :     _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
      37        1788 :     _vorticity_increment(
      38        1788 :         getMaterialPropertyByName<RankTwoTensor>(_base_name + "vorticity_increment")),
      39        3576 :     _def_grad(getMaterialPropertyByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      40        3576 :     _def_grad_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "deformation_gradient")),
      41        3576 :     _rate(getParam<MooseEnum>("objective_rate").getEnum<ObjectiveRate>()),
      42        1788 :     _polar_decomp(_rate == ObjectiveRate::GreenNaghdi),
      43        1788 :     _rotation(_polar_decomp ? &declareProperty<RankTwoTensor>(_base_name + "rotation") : nullptr),
      44        1815 :     _rotation_old(_polar_decomp ? &getMaterialPropertyOld<RankTwoTensor>(_base_name + "rotation")
      45             :                                 : nullptr),
      46        1788 :     _d_rotation_d_def_grad(
      47        3711 :         _polar_decomp ? &declareProperty<RankFourTensor>(derivativePropertyName(
      48        1842 :                             _base_name + "rotation", {_base_name + "deformation_gradient"}))
      49             :                       : nullptr),
      50        3576 :     _stretch(_polar_decomp ? &declareProperty<RankTwoTensor>(_base_name + "stretch") : nullptr)
      51             : {
      52        1815 : }
      53             : 
      54             : void
      55      446724 : ComputeLagrangianObjectiveStress::initQpStatefulProperties()
      56             : {
      57      446724 :   ComputeLagrangianStressBase::initQpStatefulProperties();
      58             : 
      59      446724 :   _small_stress[_qp].zero();
      60      446724 :   _cauchy_stress[_qp].zero();
      61             : 
      62      446724 :   if (_polar_decomp)
      63        2608 :     (*_rotation)[_qp] = RankTwoTensor::Identity();
      64      446724 : }
      65             : 
      66             : void
      67    38784580 : ComputeLagrangianObjectiveStress::computeQpCauchyStress()
      68             : {
      69    38784580 :   computeQpSmallStress();
      70             : 
      71    38784580 :   if (!_large_kinematics)
      72             :   {
      73    20870728 :     _cauchy_stress[_qp] = _small_stress[_qp];
      74    20870728 :     _cauchy_jacobian[_qp] = _small_jacobian[_qp];
      75             :   }
      76             :   else
      77             :   {
      78             :     // If large_kinematics = true, do the objective integration
      79    17913852 :     RankTwoTensor dS = _small_stress[_qp] - _small_stress_old[_qp];
      80             : 
      81    17913852 :     if (_rate == ObjectiveRate::Truesdell)
      82    16868636 :       _cauchy_stress[_qp] = objectiveUpdateTruesdell(dS);
      83     1045216 :     else if (_rate == ObjectiveRate::Jaumann)
      84      522608 :       _cauchy_stress[_qp] = objectiveUpdateJaumann(dS);
      85      522608 :     else if (_rate == ObjectiveRate::GreenNaghdi)
      86      522608 :       _cauchy_stress[_qp] = objectiveUpdateGreenNaghdi(dS);
      87             :     else
      88           0 :       mooseError("Internal error: unsupported objective rate.");
      89             :   }
      90    38784580 : }
      91             : 
      92             : RankTwoTensor
      93    16868636 : ComputeLagrangianObjectiveStress::objectiveUpdateTruesdell(const RankTwoTensor & dS)
      94             : {
      95             :   // Get the kinematic tensor
      96    16868636 :   RankTwoTensor dL = RankTwoTensor::Identity() - _inv_df[_qp];
      97             : 
      98             :   // Update the Cauchy stress
      99    16868636 :   auto [S, Jinv] = advectStress(_cauchy_stress_old[_qp] + dS, dL);
     100             : 
     101             :   // Get the appropriate tangent tensor
     102    16868636 :   RankFourTensor U = stressAdvectionDerivative(S);
     103    16868636 :   _cauchy_jacobian[_qp] = cauchyJacobian(Jinv, U);
     104             : 
     105    16868636 :   return S;
     106             : }
     107             : 
     108             : RankTwoTensor
     109      522608 : ComputeLagrangianObjectiveStress::objectiveUpdateJaumann(const RankTwoTensor & dS)
     110             : {
     111             :   usingTensorIndices(i, j, k, l);
     112             : 
     113             :   // Get the kinematic tensor
     114      522608 :   RankTwoTensor dL = RankTwoTensor::Identity() - _inv_df[_qp];
     115      522608 :   RankTwoTensor dW = 0.5 * (dL - dL.transpose());
     116             : 
     117             :   // Update the Cauchy stress
     118      522608 :   auto [S, Jinv] = advectStress(_cauchy_stress_old[_qp] + dS, dW);
     119             : 
     120             :   // Get the appropriate tangent tensor
     121             :   RankTwoTensor I = RankTwoTensor::Identity();
     122      522608 :   RankFourTensor ddW_ddL = 0.5 * (I.times<i, k, j, l>(I) - I.times<i, l, j, k>(I));
     123      522608 :   RankFourTensor U = stressAdvectionDerivative(S) * ddW_ddL;
     124      522608 :   _cauchy_jacobian[_qp] = cauchyJacobian(Jinv, U);
     125             : 
     126      522608 :   return S;
     127             : }
     128             : 
     129             : RankTwoTensor
     130      522608 : 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      522608 :   polarDecomposition();
     137             :   RankTwoTensor I = RankTwoTensor::Identity();
     138      522608 :   RankTwoTensor dR = (*_rotation)[_qp] * (*_rotation_old)[_qp].transpose() - I;
     139      522608 :   RankTwoTensor dO = dR * _inv_df[_qp];
     140             : 
     141             :   // Update the Cauchy stress
     142      522608 :   auto [S, Jinv] = advectStress(_cauchy_stress_old[_qp] + dS, dO);
     143             : 
     144             :   // Get the appropriate tangent tensor
     145      522608 :   RankFourTensor d_R_d_F = (*_d_rotation_d_def_grad)[_qp];
     146      522608 :   RankFourTensor d_F_d_dL = _inv_df[_qp].inverse().times<i, k, l, j>(_def_grad[_qp]);
     147      522608 :   RankTwoTensor T = (*_rotation_old)[_qp].transpose() * _inv_df[_qp];
     148             :   RankFourTensor d_dO_d_dL =
     149      522608 :       T.times<m, j, i, m, k, l>(d_R_d_F * d_F_d_dL) - dR.times<i, k, j, l>(I);
     150      522608 :   RankFourTensor U = stressAdvectionDerivative(S) * d_dO_d_dL;
     151      522608 :   _cauchy_jacobian[_qp] = cauchyJacobian(Jinv, U);
     152             : 
     153      522608 :   return S;
     154             : }
     155             : 
     156             : std::tuple<RankTwoTensor, RankFourTensor>
     157    17913852 : ComputeLagrangianObjectiveStress::advectStress(const RankTwoTensor & S0,
     158             :                                                const RankTwoTensor & dQ) const
     159             : {
     160    17913852 :   RankFourTensor J = updateTensor(dQ);
     161    17913852 :   RankFourTensor Jinv = J.inverse();
     162    17913852 :   RankTwoTensor S = Jinv * S0;
     163    17913852 :   return {S, Jinv};
     164             : }
     165             : 
     166             : RankFourTensor
     167    17913852 : ComputeLagrangianObjectiveStress::updateTensor(const RankTwoTensor & dQ) const
     168             : {
     169             :   auto I = RankTwoTensor::Identity();
     170             :   usingTensorIndices(i, j, k, l);
     171    35827704 :   return (1.0 + dQ.trace()) * I.times<i, k, j, l>(I) - dQ.times<i, k, j, l>(I) -
     172    35827704 :          I.times<i, k, j, l>(dQ);
     173             : }
     174             : 
     175             : RankFourTensor
     176    17913852 : ComputeLagrangianObjectiveStress::stressAdvectionDerivative(const RankTwoTensor & S) const
     177             : {
     178             :   auto I = RankTwoTensor::Identity();
     179             :   usingTensorIndices(i, j, k, l);
     180    17913852 :   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    17913852 : ComputeLagrangianObjectiveStress::cauchyJacobian(const RankFourTensor & Jinv,
     185             :                                                  const RankFourTensor & U) const
     186             : {
     187    17913852 :   return Jinv * (_small_jacobian[_qp] - U);
     188             : }
     189             : 
     190             : void
     191      522608 : ComputeLagrangianObjectiveStress::polarDecomposition()
     192             : {
     193      522608 :   FactorizedRankTwoTensor C = _def_grad[_qp].transpose() * _def_grad[_qp];
     194     1045216 :   (*_stretch)[_qp] = MathUtils::sqrt(C).get();
     195      522608 :   RankTwoTensor Uinv = MathUtils::sqrt(C).inverse().get();
     196      522608 :   (*_rotation)[_qp] = _def_grad[_qp] * Uinv;
     197             : 
     198             :   // Derivative of rotation w.r.t. the deformation gradient
     199             :   RankTwoTensor I = RankTwoTensor::Identity();
     200      522608 :   RankTwoTensor Y = (*_stretch)[_qp].trace() * I - (*_stretch)[_qp];
     201      522608 :   RankTwoTensor Z = (*_rotation)[_qp] * Y;
     202      522608 :   RankTwoTensor O = Z * (*_rotation)[_qp].transpose();
     203             :   usingTensorIndices(i, j, k, l);
     204     1045216 :   (*_d_rotation_d_def_grad)[_qp] = (O.times<i, k, l, j>(Y) - Z.times<i, l, k, j>(Z)) / Y.det();
     205      522608 : }