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