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 "ComputeLagrangianStressPK2.h"
      11             : 
      12             : InputParameters
      13         208 : ComputeLagrangianStressPK2::validParams()
      14             : {
      15         208 :   InputParameters params = ComputeLagrangianStressPK1::validParams();
      16         208 :   return params;
      17             : }
      18             : 
      19         156 : ComputeLagrangianStressPK2::ComputeLagrangianStressPK2(const InputParameters & parameters)
      20             :   : ComputeLagrangianStressPK1(parameters),
      21         156 :     _E(declareProperty<RankTwoTensor>(_base_name + "green_lagrange_strain")),
      22         156 :     _S(declareProperty<RankTwoTensor>(_base_name + "pk2_stress")),
      23         312 :     _C(declareProperty<RankFourTensor>(_base_name + "pk2_jacobian"))
      24             : {
      25         156 : }
      26             : 
      27             : void
      28     3519536 : ComputeLagrangianStressPK2::computeQpCauchyStress()
      29             : {
      30             :   // PK2 has its own sigma-wrap structure: sigma = (1/J_ust) F_ust * S * F_ust^T. The PK1 base's
      31             :   // sigma-chain (which assumes pk1_jacobian = constitutive dPK1/dF_stab + _d_F_stab_d_F_ust)
      32             :   // doesn't fit here, so compute sigma and cauchy_jacobian directly from S and dS/dE.
      33             :   //
      34             :   // The R4 chain for `_cauchy_jacobian` is Jacobian-only; the wrap to `_cauchy_stress` is
      35             :   // not. Gate accordingly.
      36     3519536 :   const bool need_jacobian = _fe_problem.currentlyComputingJacobian() ||
      37     3390440 :                              _fe_problem.currentlyComputingResidualAndJacobian();
      38     3519536 :   if (_large_kinematics)
      39             :   {
      40     1960512 :     const Real J_ust = _F_ust[_qp].det();
      41     1960512 :     _cauchy_stress[_qp] = _F_ust[_qp] * _S[_qp] * _F_ust[_qp].transpose() / J_ust;
      42             : 
      43     1960512 :     if (!need_jacobian)
      44     1837624 :       return;
      45             : 
      46             :     // cauchy_jacobian = dsigma/d(dL). sigma depends on dL only through S(E(F_stab(dL))); the
      47             :     // F_ust factors in the wrap are CONSTANT w.r.t. dL (F_ust does not depend on dL --
      48             :     // dL is computed from F_stab via the kinematic helper, and F-bar relates F_stab to
      49             :     // F_ust). So:
      50             :     //   dsigma/d(dL) = (1/J_ust) F_ust * dS/d(dL) * F_ust^T
      51             :     //            = (1/J_ust) F_ust * _C * dE/d(F_stab) * inverse(d(dL)/d(F_stab)) * F_ust^T.
      52             :     usingTensorIndices(i_, j_, k_, l_);
      53             :     const auto I2 = RankTwoTensor::Identity();
      54      122888 :     const RankFourTensor dE_dFstab = 0.5 * (I2.template times<i_, l_, j_, k_>(_F[_qp].transpose()) +
      55      245776 :                                             _F[_qp].transpose().template times<i_, k_, j_, l_>(I2));
      56      122888 :     const RankFourTensor dE_d_dL = dE_dFstab * _d_deformation_gradient_increment_d_F[_qp].inverse();
      57      122888 :     const RankFourTensor dS_d_dL = _C[_qp] * dE_d_dL;
      58      245776 :     _cauchy_jacobian[_qp] = dS_d_dL.singleProductI(_F_ust[_qp]).singleProductJ(_F_ust[_qp]) / J_ust;
      59             :   }
      60             :   else
      61             :   {
      62     1559024 :     _cauchy_stress[_qp] = _S[_qp];
      63     1559024 :     if (need_jacobian)
      64        6208 :       _cauchy_jacobian[_qp] = _C[_qp];
      65             :   }
      66             : }
      67             : 
      68             : void
      69     3519536 : ComputeLagrangianStressPK2::computeQpPK1Stress()
      70             : {
      71             :   // Green-Lagrange strain uses the F-bar-stabilized F (`_F` from the strain calc) so
      72             :   // the constitutive law receives the stabilized strain -- F-bar's purpose is precisely
      73             :   // to feed a volumetrically-corrected strain into the constitutive update.
      74     3519536 :   _E[_qp] = 0.5 * (_F[_qp].transpose() * _F[_qp] - RankTwoTensor::Identity());
      75             : 
      76             :   // PK2 update (constitutive). This populates `_S` (always needed) and `_C` (= dPK2/dE,
      77             :   // Jacobian-only). PK2-subclass implementers that want to skip `_C` on residual sweeps
      78             :   // can check `_fe_problem.currentlyComputingJacobian()` themselves; the wrap below
      79             :   // doesn't read `_C` in the residual path.
      80     3519536 :   computeQpPK2Stress();
      81             : 
      82     3519536 :   const bool need_jacobian = _fe_problem.currentlyComputingJacobian() ||
      83     3390440 :                              _fe_problem.currentlyComputingResidualAndJacobian();
      84             : 
      85             :   // PK2 -> PK1 wrap uses the *unstabilized* F so the residual matches OLD. The constitutive
      86             :   // PK2 still carries the F-bar effect via its dependence on E (Green-Lagrange of F_stab).
      87     3519536 :   if (_large_kinematics)
      88             :   {
      89     1960512 :     _pk1_stress[_qp] = _F_ust[_qp] * _S[_qp];
      90     1960512 :     if (!need_jacobian)
      91     1837624 :       return;
      92             : 
      93             :     usingTensorIndices(i_, j_, k_, l_);
      94             :     // dE/d(F_stab) (E is computed from F_stab).
      95             :     RankFourTensor dE_dFstab =
      96      122888 :         0.5 * (RankTwoTensor::Identity().times<i_, l_, j_, k_>(_F[_qp].transpose()) +
      97      368664 :                _F[_qp].transpose().times<i_, k_, j_, l_>(RankTwoTensor::Identity()));
      98             :     // Chain dE/d(F_ust) = dE/d(F_stab) * d(F_stab)/d(F_ust). With F-bar off
      99             :     // _d_F_stab_d_F_ust = IdentityFour and this collapses to dE_dFstab.
     100      122888 :     const RankFourTensor dE_dFust = dE_dFstab * _d_F_stab_d_F_ust[_qp];
     101             : 
     102             :     // dPK1/d(F_ust) = d(F_ust)/d(F_ust) * S + F_ust * dS/d(F_ust)
     103             :     //   d(F_ust)/d(F_ust)*S gives I x S^T (per the existing template, with S^T because
     104             :     //   PK1 = F*S and we differentiate F).
     105             :     //   F_ust * dS/d(F_ust) via dS/dE * dE/d(F_ust).
     106             :     const RankFourTensor termA =
     107      122888 :         RankTwoTensor::Identity().times<i_, k_, j_, l_>(_S[_qp].transpose());
     108      122888 :     _pk1_jacobian[_qp] = termA + (_C[_qp] * dE_dFust).singleProductI(_F_ust[_qp]);
     109             :     // Bypass-F-bar variant: same Term A; sigma chain uses dE_dFstab (without
     110             :     // _d_F_stab_d_F_ust). Used by specialty kernels whose coupled variable bypasses
     111             :     // F-bar.
     112      122888 :     _pk1_jacobian_bypass_fbar[_qp] = termA + (_C[_qp] * dE_dFstab).singleProductI(_F_ust[_qp]);
     113             :   }
     114             :   // Small deformations: PK1 = PK2 = sigma; PK2 chain to PK1 still needs the F-bar local
     115             :   // contribution through dE/d(F_stab) * d(F_stab)/d(F_ust). For specialty bypass paths,
     116             :   // skip the F-bar factor.
     117             :   else
     118             :   {
     119     1559024 :     _pk1_stress[_qp] = _S[_qp];
     120     1559024 :     if (!need_jacobian)
     121             :       return;
     122        6208 :     _pk1_jacobian[_qp] = _C[_qp] * _d_F_stab_d_F_ust[_qp];
     123        6208 :     _pk1_jacobian_bypass_fbar[_qp] = _C[_qp];
     124             :   }
     125             : }