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 "ComputeLagrangianStressPK1.h"
      11             : 
      12             : InputParameters
      13         316 : ComputeLagrangianStressPK1::validParams()
      14             : {
      15         316 :   InputParameters params = ComputeLagrangianStressBase::validParams();
      16             : 
      17         316 :   params.addClassDescription("Stress update based on the first Piola-Kirchhoff stress");
      18             : 
      19         316 :   return params;
      20           0 : }
      21             : 
      22         237 : ComputeLagrangianStressPK1::ComputeLagrangianStressPK1(const InputParameters & parameters)
      23             :   : ComputeLagrangianStressBase(parameters),
      24         237 :     _inv_df(getMaterialPropertyByName<RankTwoTensor>(_base_name +
      25             :                                                      "inverse_incremental_deformation_gradient")),
      26         237 :     _inv_def_grad(
      27         237 :         getMaterialPropertyByName<RankTwoTensor>(_base_name + "inverse_deformation_gradient")),
      28         711 :     _F(getMaterialPropertyByName<RankTwoTensor>(_base_name + "deformation_gradient"))
      29             : {
      30         237 : }
      31             : 
      32             : void
      33     3700552 : ComputeLagrangianStressPK1::computeQpStressUpdate()
      34             : {
      35     3700552 :   computeQpPK1Stress();
      36     3700552 :   computeQpCauchyStress(); // This could be "switched"
      37     3700552 : }
      38             : 
      39             : void
      40      160600 : ComputeLagrangianStressPK1::computeQpCauchyStress()
      41             : {
      42             :   // Wrap PK1 -> sigma using the *unstabilized* F. F-bar enters only through the constitutive
      43             :   // PK1's strain dependence (via the F-bar'd `_f_inv`); the geometric wrap uses F_ust so
      44             :   // the residual matches OLD.
      45             :   //
      46             :   // All RankFourTensor algebra below builds Jacobian-only quantities
      47             :   // (`_cauchy_jacobian`, the post-multiplied `_pk1_jacobian`, and `_pk1_jacobian_bypass_fbar`).
      48             :   // Skip them on residual-only sweeps.
      49      160600 :   const bool need_jacobian = _fe_problem.currentlyComputingJacobian() ||
      50      130784 :                              _fe_problem.currentlyComputingResidualAndJacobian();
      51      160600 :   if (_large_kinematics)
      52             :   {
      53      160600 :     const RankTwoTensor F_ust_inv = _F_ust[_qp].inverse();
      54      160600 :     const RankTwoTensor F_ust_invT = F_ust_inv.transpose();
      55      160600 :     const Real J_ust = _F_ust[_qp].det();
      56      160600 :     _cauchy_stress[_qp] = _pk1_stress[_qp] * _F_ust[_qp].transpose() / J_ust;
      57             : 
      58      160600 :     if (!need_jacobian)
      59             :       return;
      60             : 
      61             :     usingTensorIndices(i_, j_, m_, n_);
      62             :     const auto I = RankTwoTensor::Identity();
      63             : 
      64             :     // dsigma/d(F_ust) via the constitutive PK1's F_stab chain plus the geometric F_ust pieces.
      65             :     //   sigma = (1/J(F_ust)) * PK1(F_stab) * F_ust^T
      66             :     //   dsigma/d(F_ust) = (1/J) * dPK1/d(F_stab) * d(F_stab)/d(F_ust) * F_ust^T
      67             :     //               + (1/J) * PK1 * d(F_ust^T)/d(F_ust)
      68             :     //               - sigma (x) F_ust^{-T}.
      69             :     // The middle factor `_d_F_stab_d_F_ust` is IdentityFour when F-bar is off, so the
      70             :     // first line collapses to the original PK1-jacobian chain.
      71             :     // Common geometric pieces (independent of F-bar):
      72       29816 :     const RankFourTensor dsigma_dF_geom = _pk1_stress[_qp].times<i_, n_, j_, m_>(I) / J_ust -
      73       59632 :                                           _cauchy_stress[_qp].times<i_, j_, n_, m_>(F_ust_inv);
      74             :     // dsigma_dF with the constitutive PK1 jacobian chained via F-bar (default path).
      75             :     RankFourTensor dsigma_dF =
      76       29816 :         (_pk1_jacobian[_qp] * _d_F_stab_d_F_ust[_qp]).singleProductJ(_F_ust[_qp]) / J_ust;
      77       29816 :     dsigma_dF += dsigma_dF_geom;
      78       29816 :     _cauchy_jacobian[_qp] = dsigma_dF * _d_deformation_gradient_increment_d_F[_qp].inverse();
      79             : 
      80             :     // Bypass-F-bar variant: sigma chain without the `_d_F_stab_d_F_ust` factor. Same
      81             :     // geometric pieces; constitutive PK1 jacobian used directly.
      82       29816 :     RankFourTensor dsigma_dF_no_fbar = _pk1_jacobian[_qp].singleProductJ(_F_ust[_qp]) / J_ust;
      83       29816 :     dsigma_dF_no_fbar += dsigma_dF_geom;
      84             : 
      85             :     // Update _pk1_jacobian and publish the bypass variant.
      86       29816 :     const RankFourTensor pk1_jac_no_fbar = _pk1_jacobian[_qp];
      87       29816 :     _pk1_jacobian[_qp] = _pk1_jacobian[_qp] * _d_F_stab_d_F_ust[_qp];
      88       29816 :     _pk1_jacobian_bypass_fbar[_qp] = pk1_jac_no_fbar;
      89             :   }
      90             :   // Small deformations: sigma = PK1 (no wrap). For TL, the kernel chain wants
      91             :   // dPK1/d(F_ust) -- pk1_jacobian (constitutive's dPK1/d(F_stab)) must be chained to
      92             :   // F_ust via `_d_F_stab_d_F_ust`. For UL, cauchy_jacobian = dsigma/d(dL) = pk1_jacobian
      93             :   // (since the constitutive's dPK1/d(F_stab) is essentially dsigma/d(dL) at small kin).
      94             :   else
      95             :   {
      96           0 :     _cauchy_stress[_qp] = _pk1_stress[_qp];
      97           0 :     if (!need_jacobian)
      98             :       return;
      99           0 :     _cauchy_jacobian[_qp] = _pk1_jacobian[_qp];
     100           0 :     _pk1_jacobian_bypass_fbar[_qp] = _pk1_jacobian[_qp];
     101           0 :     _pk1_jacobian[_qp] = _pk1_jacobian[_qp] * _d_F_stab_d_F_ust[_qp];
     102             :   }
     103             : }